toolbox_FEniCS.py 10.1 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 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
    return plots


def function_from_xdmf(function_space, function_name, xdmf_path):
    """Read a finite element function from a xdmf file with checkpoint format.

    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.
    function_name : str
        The name of the saved function.
    xdmf_path : pathlib.Path
        Path of the xdmf file. It can be a Path object or a string.
        Extension must be '.xmdf'

    Returns
    -------
    dolfin.Function
        New function that is identical to the previously saved function.
    """
    f = fe.Function(function_space)
    file_in = fe.XDMFFile(str(xdmf_path))
    file_in.read_checkpoint(f, function_name)
    file_in.close()
    return f


Baptiste Durand's avatar
Baptiste Durand committed
127
def function_errornorm(u, v, norm_type="L2", enable_diff_fspace=False):
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
    """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)
155
        difference.vector().axpy(-1.0, v.vector())
156
    elif enable_diff_fspace:
Baptiste Durand's avatar
Baptiste Durand committed
157 158
        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
159
        difference = fe.project(u - v, u.function_space())
160
    else:
Baptiste Durand's avatar
Baptiste Durand committed
161
        raise RuntimeError("Cannot compute error norm, Function spaces do not match.")
162 163
    return fe.norm(difference, norm_type)

164

165
def local_project(v, fspace, solver_method: str = "", **kwargs):
166
    """
Baptiste Durand's avatar
Baptiste Durand committed
167
    Info : https://comet-fenics.readthedocs.io/en/latest/tips_and_tricks.html#efficient-projection-on-dg-or-quadrature-spaces #noqa
168 169 170 171

    Parameters
    ----------
    v : [type]
172
        input field
173 174
    fspace : [type]
        [description]
Baptiste Durand's avatar
Baptiste Durand committed
175
    solver_method : str, optional
176 177 178 179
        "LU" or "Cholesky" factorization. LU method used by default.

    keyword arguments
    ----------
180 181 182 183 184
    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
185 186 187 188
        ** 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.
189 190 191 192 193

    Returns
    -------
    Function
    """
194 195 196
    metadata = kwargs.get("metadata", {})
    dX = kwargs.get("dx", fe.dx(metadata=metadata))

197 198
    dv = fe.TrialFunction(fspace)
    v_ = fe.TestFunction(fspace)
199 200 201
    a_proj = fe.inner(dv, v_) * dX
    b_proj = fe.inner(v, v_) * dX

Baptiste Durand's avatar
Baptiste Durand committed
202 203 204 205 206 207 208 209 210 211
    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)
212 213 214 215
    solver.factorize()
    u = fe.Function(fspace)
    solver.solve_local_rhs(u)
    return u
Baptiste Durand's avatar
Baptiste Durand committed
216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238


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
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 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


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")
309
        facet_regions = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc)
310 311 312 313 314 315 316 317 318
    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
319