toolbox_FEniCS.py 10.9 KB
Newer Older
1 2 3 4 5 6 7 8
# coding: utf-8
"""
Created on 17/01/2019
@author: Baptiste Durand, baptiste.durand@enpc.fr

Collection of tools designed to help users working with FEniCS objects.

"""
Baptiste Durand's avatar
Baptiste Durand committed
9 10
import logging

11 12
import dolfin as fe
import matplotlib.pyplot as plt
Baptiste Durand's avatar
Baptiste Durand committed
13
import numpy as np
14
from pathlib import Path, PurePath
Baptiste Durand's avatar
Baptiste Durand committed
15 16

logger = logging.getLogger(__name__)
17

18
DOLFIN_KRYLOV_METHODS = {
Baptiste Durand's avatar
Baptiste Durand committed
19 20 21 22 23 24 25
    "bicgstab": "Biconjugate gradient stabilized method",
    "cg": "Conjugate gradient method",
    "default": "default Krylov method",
    "gmres": "Generalized minimal residual method",
    "minres": "Minimal residual method",
    "richardson": "Richardson method",
    "tfqmr": "Transpose-free quasi-minimal residual method",
26 27 28
}

DOLFIN_LU_METHODS = {
Baptiste Durand's avatar
Baptiste Durand committed
29 30 31 32 33
    "default": "default LU solver",
    "mumps": "MUMPS (MUltifrontal Massively Parallel Sparse direct Solver)",
    "petsc": "PETSc built in LU solver",
    "superlu": "SuperLU",
    "umfpack": "UMFPACK (Unsymmetric MultiFrontal sparse LU factorization)",
34 35 36
}


37 38 39
SUPPORTED_MESH_SUFFIX = (".xml", ".xdmf")


40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
def get_MeshFunction_val(msh_fctn):
    """
    Get information about the set of values that a given MeshFunction outputs.

    Parameters
    ------
    msh_fctn : instance of MeshFunction

    Returns
    -------
    Tuple (nb, vals)
    nb : int
        number of different values
    vals : list
        list of the different values that the MeshFunction outputs on its definition mesh.
    """

    val = np.unique(msh_fctn.array())
    nb = len(val)
    return (nb, val)

Baptiste Durand's avatar
Baptiste Durand committed
61 62 63 64

def facet_plot2d(
    facet_func, mesh, mesh_edges=True, markers=None, exclude_val=(0,), **kargs
):
65
    """
Baptiste Durand's avatar
Baptiste Durand committed
66
    Source : https://bitbucket.org/fenics-project/dolfin/issues/951/plotting-facetfunctions-in-matplotlib-not #noqa
67
    """
Baptiste Durand's avatar
Baptiste Durand committed
68
    x_list, y_list = [], []
Baptiste Durand's avatar
Baptiste Durand committed
69
    if markers is None:
70 71 72 73 74 75 76 77 78 79 80 81 82 83
        for facet in fe.facets(mesh):
            mp = facet.midpoint()
            x_list.append(mp.x())
            y_list.append(mp.y())
        values = facet_func.array()
    else:
        i = 0
        values = []
        for facet in fe.facets(mesh):
            if facet_func[i] in markers:
                mp = facet.midpoint()
                x_list.append(mp.x())
                y_list.append(mp.y())
                values.append(facet_func[i])
Baptiste Durand's avatar
Baptiste Durand committed
84
            i += 1
85 86 87 88 89 90 91 92 93 94 95 96 97
    if exclude_val:
        filtered_data = [], [], []
        for x, y, val in zip(x_list, y_list, values):
            if val in exclude_val:
                continue
            filtered_data[0].append(x)
            filtered_data[1].append(y)
            filtered_data[2].append(val)
        x_list, y_list, values = filtered_data

    plots = [plt.scatter(x_list, y_list, s=30, c=values, linewidths=1, **kargs)]
    if mesh_edges:
        plots.append(fe.plot(facet_func.mesh()))
98 99 100
    return plots


101
def function_from_xdmf(function_space, function_name="", xdmf_path=""):
102
    """Read a finite element function from a xdmf file with checkpoint format.
103
        Multiple functions can also be read at once.
104 105 106 107 108 109

    Parameters
    ----------
    function_space : dolfin.FunctionSpace
        Function space appropriate for the previously saved function.
        The mesh must be identical to the one used for the saved function.
110 111 112
    function_name : str or list of str
        The name(s) of the saved function.
    xdmf_path : pathlib.Path or str
113 114 115
        Path of the xdmf file. It can be a Path object or a string.
        Extension must be '.xmdf'

116 117 118 119 120
    If multiple functions are requested: #TODO : test
    the first argument must be a list of tuples (FunctionSpace, function name)
    the second argument must be empty.
    Ex: function_from_xdmf([(V,"u"),(W,"eps")], xdmf_path=checkpoint_path)

121 122
    Returns
    -------
123 124
    dolfin.Function, list of dolfin.Function
        New function(s) that is(are) identical to the previously saved function(s).
125
    """
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
    if isinstance(function_space, list):
        all_f = list()
        with fe.XDMFFile(str(xdmf_path)) as f_in:
            for fspace, fname in function_space:
                f = fe.Function(fspace)
                f_in.read_checkpoint(f, fname)
                f.rename(fname, "")
                all_f.append(f)
        return all_f
    else:
        f = fe.Function(function_space)
        f_in = fe.XDMFFile(str(xdmf_path))
        f_in.read_checkpoint(f, function_name)
        f_in.close()
        f.rename(function_name, "")
        return f
142 143


Baptiste Durand's avatar
Baptiste Durand committed
144
def function_errornorm(u, v, norm_type="L2", enable_diff_fspace=False):
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
    """Compute the difference between two functions
    defined on the same functionspace with the given norm.

    If the two objects are not functions defined in the same functionspace
    the FEniCS function errornorm should be used.
    Alternatively the constraint of function space equality can be relaxed
    with caution with the enable_diff_fspace flag.

    Based of the FEniCS functions errornorm and norm.

    Parameters
    ----------
    u, v : dolfin.functions.function.Function
    norm_type : string
        Type of norm. The :math:`L^2` -norm is default.
        For other norms, see :py:func:`norm <dolfin.fem.norms.norm>`.
    enable_diff_fspace: bool
        Relax the constraint of function space equality

    Returns
    -------
    float
        Norm of the difference
    """
    if u.function_space() == v.function_space():
        difference = fe.Function(u.function_space())
        difference.assign(u)
172
        difference.vector().axpy(-1.0, v.vector())
173
    elif enable_diff_fspace:
Baptiste Durand's avatar
Baptiste Durand committed
174 175
        logger.warning("Function spaces not equals.")
        logger.warning(f"Projection to compute the difference between {u} and {v}")
Baptiste Durand's avatar
Baptiste Durand committed
176
        difference = fe.project(u - v, u.function_space())
177
    else:
Baptiste Durand's avatar
Baptiste Durand committed
178
        raise RuntimeError("Cannot compute error norm, Function spaces do not match.")
179 180
    return fe.norm(difference, norm_type)

181

182
def local_project(v, fspace, solver_method: str = "", **kwargs):
183
    """
Baptiste Durand's avatar
Baptiste Durand committed
184
    Info : https://comet-fenics.readthedocs.io/en/latest/tips_and_tricks.html#efficient-projection-on-dg-or-quadrature-spaces #noqa
185 186 187 188

    Parameters
    ----------
    v : [type]
189
        input field
190 191
    fspace : [type]
        [description]
Baptiste Durand's avatar
Baptiste Durand committed
192
    solver_method : str, optional
193 194 195 196
        "LU" or "Cholesky" factorization. LU method used by default.

    keyword arguments
    ----------
197 198 199 200 201
    metadata : dict, optional
        This can be used to define different quadrature degrees for different
        terms in a form, and to override other form compiler specific options
        separately for different terms. By default {}
        See UFL user manual for more information
202 203 204 205
        ** WARNING** : May be over ignored if dx argument is also used. (Non étudié)
    dx : Measure
        Impose the dx measure for the projection.
        This is for example useful to do some kind of interpolation on DG functionspaces.
206 207 208 209 210

    Returns
    -------
    Function
    """
211 212 213
    metadata = kwargs.get("metadata", {})
    dX = kwargs.get("dx", fe.dx(metadata=metadata))

214 215
    dv = fe.TrialFunction(fspace)
    v_ = fe.TestFunction(fspace)
216 217 218
    a_proj = fe.inner(dv, v_) * dX
    b_proj = fe.inner(v, v_) * dX

Baptiste Durand's avatar
Baptiste Durand committed
219 220 221 222 223 224 225 226 227 228
    if solver_method == "LU":
        solver = fe.LocalSolver(
            a_proj, b_proj, solver_type=fe.cpp.fem.LocalSolver.SolverType.LU
        )
    elif solver_method == "Cholesky":
        solver = fe.LocalSolver(
            a_proj, b_proj, solver_type=fe.cpp.fem.LocalSolver.SolverType.Cholesky
        )
    else:
        solver = fe.LocalSolver(a_proj, b_proj)
229 230 231 232
    solver.factorize()
    u = fe.Function(fspace)
    solver.solve_local_rhs(u)
    return u
Baptiste Durand's avatar
Baptiste Durand committed
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255


def _wrap_in_list(obj, name, types=type):
    """
    Transform single object or a collection of objects into a list.

    Source
    ------
    python/dolfin/fem/assembling.py, commit 4c72333
    """

    if obj is None:
        lst = []
    elif hasattr(obj, "__iter__"):
        lst = list(obj)
    else:
        lst = [obj]
    for obj in lst:
        if not isinstance(obj, types):
            raise TypeError(
                "expected a (list of) %s as '%s' argument" % (str(types), name)
            )
    return lst
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325


def xdmf_mesh(mesh_file, import_subdomains=False, facet_file="", physical_file=""):
    """Create a FeniCS mesh from a mesh file with xdmf format.

    Parameters
    ----------
    mesh_file : str or Path
        Path of the file mesh (xdmf file)
    import_subdomains : bool, optional
        True if information about subdomains have to be imported.
        The paths of the auxiliary files that contains information about subdomains
        can be indicated with facet_file and physical_file.
        The paths used by default are :
            - "<mesh path>_facet_region.xdmf" and
            - "<mesh path>_physical_region.xdmf" (for subdomains)
    facet_file : str or Path, optional
        Path to the mesh auxiliary file that contains subdomains data.
        Defaults to "" i.e. the default path will be used.
    physical_file : str or Path, optional
        Path to the mesh auxiliary file that contains facet regions data.
        Defaults to "" i.e. the default path will be used.

    Returns
    -------
    Mesh / tuple
        If subdomains are not requested :
            - The Mesh instance
        If subdomains are requested :
            - The Mesh instance,
            - The MeshFunction for subdomains if it exists else None,
            - The MeshFunction for facets if it exists else None,

    Source
    ------
    Gist meshtagging_mvc.py, June 2018, Michal Habera
    https://gist.github.com/michalhabera/bbe8a17f788192e53fd758a67cbf3bed
    """
    if not isinstance(mesh_file, PurePath):
        mesh_file = Path(mesh_file)
    if not mesh_file.suffix == ".xdmf":
        raise TypeError("Wrong suffix for the path to the mesh.")
    mesh = fe.Mesh()
    with fe.XDMFFile(str(mesh_file)) as f_in:
        f_in.read(mesh)

    if not import_subdomains:
        return mesh

    dim = mesh.geometric_dimension()

    if not facet_file:
        facet_file = mesh_file.with_name(f"{mesh_file.stem}_facet_region.xdmf")
    else:
        facet_file = Path(facet_file)
        if not facet_file.suffix == ".xdmf":
            raise TypeError("Wrong suffix for the path to facet regions.")
    if not physical_file:
        physical_file = mesh_file.with_name(f"{mesh_file.stem}_physical_region.xdmf")
    else:
        physical_file = Path(physical_file)
        if not physical_file.suffix == ".xdmf":
            raise TypeError("Wrong suffix for the path to subdomains.")

    if not facet_file.exists():
        facet_regions = None
    else:
        facet_vc = fe.MeshValueCollection("size_t", mesh, dim - 1)
        with fe.XDMFFile(str(facet_file)) as f_in:
            f_in.read(facet_vc, "facet_data")
326
        facet_regions = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc)
327 328 329 330 331 332 333 334 335
    if not physical_file.exists():
        subdomains = None
    else:
        cell_vc = fe.MeshValueCollection("size_t", mesh, dim - 1)
        with fe.XDMFFile(str(physical_file)) as f_in:
            f_in.read(cell_vc, "cell_data")
        subdomains = fe.cpp.mesh.MeshFunctionSizet(mesh, cell_vc)

    return mesh, subdomains, facet_regions