|
1 | | -# Copyright (C) 2024 Garth N. Wells |
| 1 | +# Copyright (C) 2024 Garth N. Wells and Paul T. Kühner |
2 | 2 | # |
3 | 3 | # This file is part of DOLFINx (https://www.fenicsproject.org) |
4 | 4 | # |
|
12 | 12 | import numpy.typing as npt |
13 | 13 |
|
14 | 14 | import basix |
| 15 | +import ufl |
| 16 | +import ufl.finiteelement |
15 | 17 | from dolfinx import cpp as _cpp |
16 | 18 |
|
17 | 19 |
|
@@ -93,7 +95,7 @@ def pull_back( |
93 | 95 | ``shape=(num_points, geometrical_dimension)``. |
94 | 96 | cell_geometry: Physical coordinates describing the cell, |
95 | 97 | shape ``(num_of_geometry_basis_functions, geometrical_dimension)`` |
96 | | - They can be created by accessing `geometry.x[geometry.dofmap.cell_dofs(i)]`, |
| 98 | + They can be created by accessing ``geometry.x[geometry.dofmap.cell_dofs(i)]``, |
97 | 99 |
|
98 | 100 | Returns: |
99 | 101 | Reference coordinates of the physical points ``x``. |
@@ -160,3 +162,190 @@ def _(e: basix.finite_element.FiniteElement): |
160 | 162 | return CoordinateElement(_cpp.fem.CoordinateElement_float32(e._e)) |
161 | 163 | except TypeError: |
162 | 164 | return CoordinateElement(_cpp.fem.CoordinateElement_float64(e._e)) |
| 165 | + |
| 166 | + |
| 167 | +class FiniteElement: |
| 168 | + _cpp_object: typing.Union[_cpp.fem.FiniteElement_float32, _cpp.fem.FiniteElement_float64] |
| 169 | + |
| 170 | + def __init__( |
| 171 | + self, |
| 172 | + cpp_object: typing.Union[_cpp.fem.FiniteElement_float32, _cpp.fem.FiniteElement_float64], |
| 173 | + ): |
| 174 | + """Creates a Python wrapper for the exported finite element class. |
| 175 | +
|
| 176 | + Note: |
| 177 | + Do not use this constructor directly. Instead use :func:``finiteelement``. |
| 178 | +
|
| 179 | + Args: |
| 180 | + The underlying cpp instance that this object will wrap. |
| 181 | + """ |
| 182 | + self._cpp_object = cpp_object |
| 183 | + |
| 184 | + def __eq__(self, other): |
| 185 | + return self._cpp_object == other._cpp_object |
| 186 | + |
| 187 | + @property |
| 188 | + def dtype(self) -> np.dtype: |
| 189 | + """Geometry type of the Mesh that the FunctionSpace is defined on.""" |
| 190 | + return self._cpp_object.dtype |
| 191 | + |
| 192 | + @property |
| 193 | + def basix_element(self) -> basix.finite_element.FiniteElement: |
| 194 | + """Return underlying Basix C++ element (if it exists). |
| 195 | +
|
| 196 | + Raises: |
| 197 | + Runtime error if Basix element does not exist. |
| 198 | + """ |
| 199 | + return self._cpp_object.basix_element |
| 200 | + |
| 201 | + @property |
| 202 | + def num_sub_elements(self) -> int: |
| 203 | + """Number of sub elements (for a mixed or blocked element).""" |
| 204 | + return self._cpp_object.num_sub_elements |
| 205 | + |
| 206 | + @property |
| 207 | + def value_shape(self) -> npt.NDArray[np.integer]: |
| 208 | + """Value shape of the finite element field. |
| 209 | +
|
| 210 | + The value shape describes the shape of the finite element field, e.g. ``{}`` for a scalar, |
| 211 | + ``{2}`` for a vector in 2D, ``{3, 3}`` for a rank-2 tensor in 3D, etc. |
| 212 | + """ |
| 213 | + return self._cpp_object.value_shape |
| 214 | + |
| 215 | + @property |
| 216 | + def interpolation_points(self) -> npt.NDArray[np.floating]: |
| 217 | + """Points on the reference cell at which an expression needs to be evaluated in order to |
| 218 | + interpolate the expression in the finite element space. |
| 219 | +
|
| 220 | + Interpolation point coordinates on the reference cell, returning the coordinates data |
| 221 | + (row-major) storage with shape ``(num_points, tdim)``. |
| 222 | +
|
| 223 | + Note: |
| 224 | + For Lagrange elements the points will just be the nodal positions. For other elements |
| 225 | + the points will typically be the quadrature points used to evaluate moment degrees of |
| 226 | + freedom. |
| 227 | + """ |
| 228 | + return self._cpp_object.interpolation_points |
| 229 | + |
| 230 | + @property |
| 231 | + def interpolation_ident(self) -> bool: |
| 232 | + """Check if interpolation into the finite element space is an identity operation given the |
| 233 | + evaluation on an expression at specific points, i.e. the degree-of-freedom are equal to |
| 234 | + point evaluations. The function will return `true` for Lagrange elements.""" |
| 235 | + return self._cpp_object.interpolation_ident |
| 236 | + |
| 237 | + @property |
| 238 | + def space_dimension(self) -> int: |
| 239 | + """Dimension of the finite element function space (the number of degrees-of-freedom for the |
| 240 | + element). |
| 241 | +
|
| 242 | + For 'blocked' elements, this function returns the dimension of the full element rather than |
| 243 | + the dimension of the base element. |
| 244 | + """ |
| 245 | + return self._cpp_object.space_dimension |
| 246 | + |
| 247 | + @property |
| 248 | + def needs_dof_transformations(self) -> bool: |
| 249 | + """Check if DOF transformations are needed for this element. |
| 250 | +
|
| 251 | + DOF transformations will be needed for elements which might not be continuous when two |
| 252 | + neighbouring cells disagree on the orientation of a shared sub-entity, and when this cannot |
| 253 | + be corrected for by permuting the DOF numbering in the dofmap. |
| 254 | +
|
| 255 | + For example, Raviart-Thomas elements will need DOF transformations, as the neighbouring |
| 256 | + cells may disagree on the orientation of a basis function, and this orientation cannot be |
| 257 | + corrected for by permuting the DOF numbers on each cell. |
| 258 | + """ |
| 259 | + return self._cpp_object.needs_dof_transformations |
| 260 | + |
| 261 | + @property |
| 262 | + def signature(self) -> str: |
| 263 | + """String identifying the finite element.""" |
| 264 | + return self._cpp_object.signature |
| 265 | + |
| 266 | + def T_apply(self, x: npt.NDArray[np.floating], cell_permutations: np.int32, dim: int) -> None: |
| 267 | + """Transform basis functions from the reference element ordering and orientation to the |
| 268 | + globally consistent physical element ordering and orientation. |
| 269 | +
|
| 270 | + Args: |
| 271 | + x: Data to transform (in place). The shape is ``(m, n)``, where `m` is the number of |
| 272 | + dgerees-of-freedom and the storage is row-major. |
| 273 | + cell_permutations: Permutation data for the cell. |
| 274 | + dim: Number of columns in ``data``. |
| 275 | +
|
| 276 | + Note: |
| 277 | + Exposed for testing. Function is not vectorised across multiple cells. Please see |
| 278 | + `basix.numba_helpers` for performant versions. |
| 279 | + """ |
| 280 | + self._cpp_object.T_apply(x, cell_permutations, dim) |
| 281 | + |
| 282 | + def Tt_apply(self, x: npt.NDArray[np.floating], cell_permutations: np.int32, dim: int) -> None: |
| 283 | + """Apply the transpose of the operator applied by T_apply(). |
| 284 | +
|
| 285 | + Args: |
| 286 | + x: Data to transform (in place). The shape is ``(m, n)``, where `m` is the number of |
| 287 | + dgerees-of-freedom and the storage is row-major. |
| 288 | + cell_permutations: Permutation data for the cell. |
| 289 | + dim: Number of columns in `data`. |
| 290 | +
|
| 291 | + Note: |
| 292 | + Exposed for testing. Function is not vectorised across multiple cells. Please see |
| 293 | + `basix.numba_helpers` for performant versions. |
| 294 | + """ |
| 295 | + self._cpp_object.Tt_apply(x, cell_permutations, dim) |
| 296 | + |
| 297 | + def Tt_inv_apply( |
| 298 | + self, x: npt.NDArray[np.floating], cell_permutations: np.int32, dim: int |
| 299 | + ) -> None: |
| 300 | + """Apply the inverse transpose of the operator applied by T_apply(). |
| 301 | +
|
| 302 | + Args: |
| 303 | + x: Data to transform (in place). The shape is ``(m, n)``, where ``m`` is the number of |
| 304 | + dgerees-of-freedom and the storage is row-major. |
| 305 | + cell_permutations: Permutation data for the cell. |
| 306 | + dim: Number of columns in `data`. |
| 307 | +
|
| 308 | + Note: |
| 309 | + Exposed for testing. Function is not vectorised across multiple cells. Please see |
| 310 | + ``basix.numba_helpers`` for performant versions. |
| 311 | + """ |
| 312 | + self._cpp_object.Tt_apply(x, cell_permutations, dim) |
| 313 | + |
| 314 | + |
| 315 | +def finiteelement( |
| 316 | + cell_type: _cpp.mesh.CellType, |
| 317 | + ufl_e: ufl.finiteelement, |
| 318 | + FiniteElement_dtype: np.dtype, |
| 319 | +) -> FiniteElement: |
| 320 | + """Create a DOLFINx element from a basix.ufl element. |
| 321 | +
|
| 322 | + Args: |
| 323 | + cell_type: Element cell type, see ``mesh.CellType`` |
| 324 | + ufl_e: UFL element, holding quadrature rule and other properties of the selected element. |
| 325 | + FiniteElement_dtype: Geometry type of the element. |
| 326 | + """ |
| 327 | + if np.issubdtype(FiniteElement_dtype, np.float32): |
| 328 | + CppElement = _cpp.fem.FiniteElement_float32 |
| 329 | + elif np.issubdtype(FiniteElement_dtype, np.float64): |
| 330 | + CppElement = _cpp.fem.FiniteElement_float64 |
| 331 | + else: |
| 332 | + raise ValueError(f"Unsupported dtype: {FiniteElement_dtype}") |
| 333 | + |
| 334 | + if ufl_e.is_mixed: |
| 335 | + elements = [ |
| 336 | + finiteelement(cell_type, e, FiniteElement_dtype)._cpp_object for e in ufl_e.sub_elements |
| 337 | + ] |
| 338 | + return FiniteElement(CppElement(elements)) |
| 339 | + elif ufl_e.is_quadrature: |
| 340 | + return FiniteElement( |
| 341 | + CppElement( |
| 342 | + cell_type, |
| 343 | + ufl_e.custom_quadrature()[0], |
| 344 | + ufl_e.reference_value_shape, |
| 345 | + ufl_e.is_symmetric, |
| 346 | + ) |
| 347 | + ) |
| 348 | + else: |
| 349 | + basix_e = ufl_e.basix_element._e |
| 350 | + value_shape = ufl_e.reference_value_shape if ufl_e.block_size > 1 else None |
| 351 | + return FiniteElement(CppElement(basix_e, value_shape, ufl_e.is_symmetric)) |
0 commit comments