Source code for dolfinx.fem.forms
# Copyright (C) 2017-2021 Chris N. Richardson, Garth N. Wells and Michal Habera
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
from __future__ import annotations
import collections
import collections.abc
import typing
from dolfinx.fem.function import FunctionSpace
if typing.TYPE_CHECKING:
from dolfinx.fem import function
from dolfinx.mesh import Mesh
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import jit
from petsc4py import PETSc
[docs]class FormMetaClass:
def __init__(self, form, V: list[_cpp.fem.FunctionSpace], coeffs, constants,
subdomains: dict[_cpp.mesh.MeshTags_int32, typing.Union[None, typing.Any]], mesh: _cpp.mesh.Mesh,
ffi, code):
"""A finite element form
Notes:
Forms should normally be constructed using
:func:`forms.form` and not using this class initialiser.
This class is combined with different base classes that
depend on the scalar type used in the Form.
Args:
form: Compiled UFC form
V: The argument function spaces
coeffs: Finite element coefficients that appear in the form
constants: Constants appearing in the form
subdomains: Subdomains for integrals
mesh: The mesh that the form is defined on
"""
self._code = code
self._ufcx_form = form
super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
V, coeffs, constants, subdomains, mesh) # type: ignore
@property
def ufcx_form(self):
"""The compiled ufcx_form object"""
return self._ufcx_form
@property
def code(self) -> str:
"""C code strings"""
return self._code
@property
def function_spaces(self) -> typing.List[FunctionSpace]:
"""Function spaces on which this form is defined"""
return super().function_spaces # type: ignore
@property
def dtype(self) -> np.dtype:
"""dtype of this form"""
return super().dtype # type: ignore
@property
def mesh(self) -> Mesh:
"""Mesh on which this form is defined"""
return super().mesh # type: ignore
@property
def integral_types(self):
"""Integral types in the form"""
return super().integral_types # type: ignore
form_types = typing.Union[FormMetaClass, _cpp.fem.Form_float32, _cpp.fem.Form_float64,
_cpp.fem.Form_complex64, _cpp.fem.Form_complex128]
[docs]def form(form: typing.Union[ufl.Form, typing.Iterable[ufl.Form]], dtype: np.dtype = PETSc.ScalarType,
form_compiler_options: dict = {}, jit_options: dict = {}):
"""Create a DOLFINx Form or an array of Forms
Args:
form: A UFL form or list(s) of UFL forms
dtype: Scalar type to use for the compiled form
form_compiler_options: See :func:`ffcx_jit <dolfinx.jit.ffcx_jit>`
jit_options:See :func:`ffcx_jit <dolfinx.jit.ffcx_jit>`
Returns:
Compiled finite element Form
Notes:
This function is responsible for the compilation of a UFL form
(using FFCx) and attaching coefficients and domains specific
data to the underlying C++ form. It dynamically create a
:class:`Form` instance with an appropriate base class for the
scalar type, e.g. `_cpp.fem.Form_float64`.
"""
if dtype == np.float32:
ftype = _cpp.fem.Form_float32
form_compiler_options["scalar_type"] = "float"
elif dtype == np.float64:
ftype = _cpp.fem.Form_float64
form_compiler_options["scalar_type"] = "double"
elif dtype == np.complex64:
ftype = _cpp.fem.Form_complex64
form_compiler_options["scalar_type"] = "float _Complex"
elif dtype == np.complex128:
ftype = _cpp.fem.Form_complex128
form_compiler_options["scalar_type"] = "double _Complex"
else:
raise NotImplementedError(f"Type {dtype} not supported.")
formcls = type("Form", (FormMetaClass, ftype), {})
def _form(form):
""""Compile a single UFL form"""
# Extract subdomain data from UFL form
sd = form.subdomain_data()
domain, = list(sd.keys()) # Assuming single domain
# Get subdomain data for each integral type
subdomains = {}
for integral_type, data in sd.get(domain).items():
# Check that the subdomain data for each integral of this type is
# the same
assert all([id(d) == id(data[0]) for d in data])
subdomains[integral_type] = data[0]
mesh = domain.ufl_cargo()
if mesh is None:
raise RuntimeError("Expecting to find a Mesh in the form.")
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
form_compiler_options=form_compiler_options,
jit_options=jit_options)
# For each argument in form extract its function space
V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
# Prepare coefficients data. For every coefficient in form take its
# C++ object.
original_coefficients = form.coefficients()
coeffs = [original_coefficients[ufcx_form.original_coefficient_position[i]
]._cpp_object for i in range(ufcx_form.num_coefficients)]
constants = [c._cpp_object for c in form.constants()]
# Subdomain markers (possibly None for some dimensions)
subdomains = {_cpp.fem.IntegralType.cell: subdomains.get("cell"),
_cpp.fem.IntegralType.exterior_facet: subdomains.get("exterior_facet"),
_cpp.fem.IntegralType.interior_facet: subdomains.get("interior_facet"),
_cpp.fem.IntegralType.vertex: subdomains.get("vertex")}
return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)
def _create_form(form):
"""Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
return form argument"""
if isinstance(form, ufl.Form):
return _form(form)
elif isinstance(form, collections.abc.Iterable):
return list(map(lambda sub_form: _create_form(sub_form), form))
return form
return _create_form(form)
[docs]def extract_function_spaces(forms: typing.Union[typing.Iterable[FormMetaClass], # type: ignore [return]
typing.Iterable[typing.Iterable[FormMetaClass]]],
index: int = 0) -> typing.Iterable[typing.Union[None, function.FunctionSpace]]:
"""Extract common function spaces from an array of forms. If `forms`
is a list of linear form, this function returns of list of the
corresponding test functions. If `forms` is a 2D array of bilinear
forms, for index=0 the list common test function space for each row
is returned, and if index=1 the common trial function spaces for
each column are returned."""
_forms = np.array(forms)
if _forms.ndim == 0:
raise RuntimeError("Expected an array for forms, not a single form")
elif _forms.ndim == 1:
assert index == 0
for form in _forms:
if form is not None:
assert form.rank == 1, "Expected linear form"
return [form.function_spaces[0] if form is not None else None for form in forms] # type: ignore[union-attr]
elif _forms.ndim == 2:
assert index == 0 or index == 1
extract_spaces = np.vectorize(lambda form: form.function_spaces[index] if form is not None else None)
V = extract_spaces(_forms)
def unique_spaces(V):
# Pick spaces from first column
V0 = V[:, 0]
# Iterate over each column
for col in range(1, V.shape[1]):
# Iterate over entry in column, updating if current
# space is None, or where both spaces are not None check
# that they are the same
for row in range(V.shape[0]):
if V0[row] is None and V[row, col] is not None:
V0[row] = V[row, col]
elif V0[row] is not None and V[row, col] is not None:
assert V0[row] is V[row, col], "Cannot extract unique function spaces"
return V0
if index == 0:
return list(unique_spaces(V))
elif index == 1:
return list(unique_spaces(V.transpose()))
else:
raise RuntimeError("Unsupported array of forms")