diff --git a/generate_all_operators.py b/generate_all_operators.py index 005b4d30b2a125d42503e1ad733715da8b819e87..282f968421204f8f540bd6040003dc56603159ec 100644 --- a/generate_all_operators.py +++ b/generate_all_operators.py @@ -29,7 +29,7 @@ import sympy as sp from sympy.core.cache import clear_cache from tabulate import tabulate -from hog.blending import GeometryMap, IdentityMap, AnnulusMap, IcosahedralShellMap +from hog.blending import GeometryMap, IdentityMap, AnnulusMap, IcosahedralShellMap, ParametricMap from hog.cse import CseImplementation from hog.element_geometry import ( ElementGeometry, @@ -101,6 +101,8 @@ ALL_GEOMETRY_MAPS = [ IdentityMap(), AnnulusMap(), IcosahedralShellMap(), + ParametricMap(1), + ParametricMap(2), ] diff --git a/hog/blending.py b/hog/blending.py index 01d5a79c0f0513ec4c00c3783da2afacfbd659ce..d4eea4e576be79958c6d2306e77c5f38e094229b 100644 --- a/hog/blending.py +++ b/hog/blending.py @@ -49,7 +49,7 @@ class GeometryMap: def jacobian(self, x: sp.Matrix) -> sp.Matrix: """Evaluates the Jacobian of the geometry map at the passed point.""" raise HOGException("jacobian() not implemented for this map.") - + def hessian(self, x: sp.Matrix) -> List[sp.Matrix]: """Evaluates the hessian of the geometry map at the passed point.""" raise HOGException("hessian() not implemented for this map.") @@ -100,6 +100,45 @@ class ExternalMap(GeometryMap): return [LineElement(), TriangleElement(), TetrahedronElement()] +class ParametricMap(GeometryMap): + """ + This is a special map that indicates that you want a parametric mapping. + + It uses HyTeG's MicroMesh implementation and introduces a vector finite element coefficient that represents the + coordinates at each node. Locally, the Jacobians are computed through the gradients of the shape functions. + + Depending on your choice of polynomial degree, you can use this blending map to construct sub-, super-, and + isoparametric mappings. + + The affine Jacobians in all integrands will automatically be set to the identity if you use this map. + The blending Jacobians will be set to the transpose of the gradient of the shape functions. + """ + + def __init__(self, degree: int): + + if degree not in [1, 2]: + raise HOGException( + "Only first and second order parametric maps are supported." + ) + + self.degree = degree + + def supported_geometries(self) -> List[ElementGeometry]: + return [TriangleElement(), TetrahedronElement()] + + def parameter_coupling_code(self) -> str: + return "" + + def coupling_includes(self) -> List[str]: + return [] + + def __eq__(self, other: Any) -> bool: + return isinstance(other, ParametricMap) and self.degree == other.degree + + def __str__(self): + return self.__class__.__name__ + f"P{self.degree}" + + class AnnulusMap(GeometryMap): """Projects HyTeG's approximate annulus to the actual curved geometry. diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py index 855bdd9ef3443249b0b20d6b4c3213e21b809a15..4857cdf578449ac2d8ca4e97108f4732765884d8 100644 --- a/hog/fem_helpers.py +++ b/hog/fem_helpers.py @@ -22,8 +22,7 @@ from hog.blending import ( GeometryMap, ExternalMap, IdentityMap, - AnnulusMap, - IcosahedralShellMap, + ParametricMap, ) from hog.element_geometry import ( ElementGeometry, @@ -38,11 +37,9 @@ from hog.multi_assignment import MultiAssignment from hog.symbolizer import Symbolizer from hog.external_functions import ( BlendingFTriangle, - BlendingFEmbeddedTriangle, BlendingFTetrahedron, BlendingDFTetrahedron, BlendingDFTriangle, - BlendingDFInvDFTriangle, BlendingDFEmbeddedTriangle, ScalarVariableCoefficient2D, ScalarVariableCoefficient3D, @@ -209,6 +206,9 @@ def trafo_ref_to_physical( if geometry not in blending.supported_geometries(): raise HOGException("Geometry not supported by blending map.") + if isinstance(blending, ParametricMap): + raise HOGException("Evaluation not implemented for parametric maps.") + t = trafo_ref_to_affine(geometry, symbolizer) if isinstance(blending, ExternalMap): diff --git a/hog/integrand.py b/hog/integrand.py index b6175a3a7470adcaa9cd557499dd47d6d3159c1e..81586f34781bb8bda245bfa4b72aad9d83650b05 100644 --- a/hog/integrand.py +++ b/hog/integrand.py @@ -52,11 +52,17 @@ from dataclasses import dataclass, asdict, fields, field import sympy as sp from hog.exception import HOGException -from hog.function_space import FunctionSpace, TrialSpace, TestSpace +from hog.function_space import ( + FunctionSpace, + TrialSpace, + TestSpace, + TensorialVectorFunctionSpace, + LagrangianFunctionSpace, +) from hog.element_geometry import ElementGeometry from hog.quadrature import Quadrature, Tabulation from hog.symbolizer import Symbolizer -from hog.blending import GeometryMap, IdentityMap, ExternalMap +from hog.blending import GeometryMap, IdentityMap, ExternalMap, ParametricMap from hog.fem_helpers import ( create_empty_element_matrix, element_matrix_iterator, @@ -293,6 +299,28 @@ def process_integrand( s = IntegrandSymbols() + #################### + # Element geometry # + #################### + + s.volume_geometry = volume_geometry + + if boundary_geometry is not None: + if boundary_geometry.dimensions != boundary_geometry.space_dimension - 1: + raise HOGException( + "Since you are integrating over a boundary, the boundary element's space dimension should be larger " + "than its dimension." + ) + + if boundary_geometry.space_dimension != volume_geometry.space_dimension: + raise HOGException("All geometries must be embedded in the same space.") + + s.boundary_geometry = boundary_geometry + + ############## + # Tabulation # + ############## + tabulation = Tabulation(symbolizer) def _tabulate( @@ -305,6 +333,10 @@ def process_integrand( s.tabulate = _tabulate + ################################ + # Scalar and matrix parameters # + ################################ + free_symbols = set() def _scalars(symbol_names: str | List[str]) -> sp.Symbol | List[sp.Symbol]: @@ -329,46 +361,29 @@ def process_integrand( s.scalars = _scalars s.matrix = _matrix - s.volume_geometry = volume_geometry + ################### + # FE coefficients # + ################### - s.jac_a = symbolizer.jac_ref_to_affine(volume_geometry) - s.jac_a_inv = symbolizer.jac_ref_to_affine_inv(volume_geometry) - s.jac_a_abs_det = symbolizer.abs_det_jac_ref_to_affine() + fe_coefficients_modified = {k: v for k, v in fe_coefficients.items()} - if isinstance(blending, IdentityMap): - s.jac_b = sp.eye(volume_geometry.space_dimension) - s.jac_b_inv = sp.eye(volume_geometry.space_dimension) - s.jac_b_abs_det = 1 - elif isinstance(blending, ExternalMap): - s.jac_b = jac_affine_to_physical(volume_geometry, symbolizer) - s.jac_b_inv = inv(s.jac_b) - s.jac_b_abs_det = abs(det(s.jac_b)) - else: - s.jac_b = symbolizer.jac_affine_to_blending(volume_geometry.space_dimension) - s.jac_b_inv = symbolizer.jac_affine_to_blending_inv( - volume_geometry.space_dimension - ) - s.jac_b_abs_det = symbolizer.abs_det_jac_affine_to_blending() - s.hessian_b = symbolizer.hessian_blending_map(volume_geometry.dimensions) - - if boundary_geometry is not None: - if boundary_geometry.dimensions != boundary_geometry.space_dimension - 1: + special_name_of_micromesh_coeff = "micromesh" + if isinstance(blending, ParametricMap): + # We add a vector coefficient for the parametric mapping here. + if special_name_of_micromesh_coeff in fe_coefficients: raise HOGException( - "Since you are integrating over a boundary, the boundary element's space dimension should be larger " - "than its dimension." + f"You cannot use the name {special_name_of_micromesh_coeff} for your FE coefficient." + f"It is reserved." ) - - if boundary_geometry.space_dimension != volume_geometry.space_dimension: - raise HOGException("All geometries must be embedded in the same space.") - - s.boundary_geometry = boundary_geometry - s.jac_a_boundary = symbolizer.jac_ref_to_affine(boundary_geometry) - - s.x = trafo_ref_to_physical(volume_geometry, symbolizer, blending) + fe_coefficients_modified[special_name_of_micromesh_coeff] = ( + TensorialVectorFunctionSpace( + LagrangianFunctionSpace(blending.degree, symbolizer) + ) + ) s.k = dict() s.grad_k = dict() - for name, coefficient_function_space in fe_coefficients.items(): + for name, coefficient_function_space in fe_coefficients_modified.items(): if coefficient_function_space is None: k = scalar_space_dependent_coefficient( name, volume_geometry, symbolizer, blending=blending @@ -400,6 +415,66 @@ def process_integrand( s.k[name] = k s.grad_k[name] = grad_k + ############################## + # Jacobians and determinants # + ############################## + + s.jac_a = symbolizer.jac_ref_to_affine(volume_geometry) + s.jac_a_inv = symbolizer.jac_ref_to_affine_inv(volume_geometry) + s.jac_a_abs_det = symbolizer.abs_det_jac_ref_to_affine() + + if boundary_geometry is not None: + s.jac_a_boundary = symbolizer.jac_ref_to_affine(boundary_geometry) + + if isinstance(blending, IdentityMap): + s.jac_b = sp.eye(volume_geometry.space_dimension) + s.jac_b_inv = sp.eye(volume_geometry.space_dimension) + s.jac_b_abs_det = 1 + elif isinstance(blending, ExternalMap): + s.jac_b = jac_affine_to_physical(volume_geometry, symbolizer) + s.jac_b_inv = inv(s.jac_b) + s.jac_b_abs_det = abs(det(s.jac_b)) + elif isinstance(blending, ParametricMap): + s.jac_a = sp.eye(volume_geometry.space_dimension) + s.jac_a_inv = sp.eye(volume_geometry.space_dimension) + s.jac_a_abs_det = 1 + + if boundary_geometry is not None: + raise HOGException( + "Boundary integrals not tested with parametric mappings yet. " + "We have to handle/set the affine Jacobian at the boundary appropriately.\n" + "Dev note to future me: since we assume the boundary element to have the last ref coord zero, " + "I suppose we can set this thing to:\n" + " ⎡ 1 ⎤ \n" + " ⎣ 0 ⎦ \n" + "in 2D and to\n" + " ⎡ 1 0 ⎤ \n" + " | 0 1 | \n" + " ⎣ 0 0 ⎦ \n" + "in 3D." + ) + + s.jac_b = s.grad_k[special_name_of_micromesh_coeff].T + s.jac_b_inv = inv(s.jac_b) + s.jac_b_abs_det = abs(det(s.jac_b)) + + s.x = s.k[special_name_of_micromesh_coeff] + + else: + s.jac_b = symbolizer.jac_affine_to_blending(volume_geometry.space_dimension) + s.jac_b_inv = symbolizer.jac_affine_to_blending_inv( + volume_geometry.space_dimension + ) + s.jac_b_abs_det = symbolizer.abs_det_jac_affine_to_blending() + s.hessian_b = symbolizer.hessian_blending_map(volume_geometry.dimensions) + + if not isinstance(blending, ParametricMap): + s.x = trafo_ref_to_physical(volume_geometry, symbolizer, blending) + + ####################################### + # Assembling the local element matrix # + ####################################### + mat = create_empty_element_matrix(trial, test, volume_geometry) it = element_matrix_iterator(trial, test, volume_geometry) diff --git a/hog/quadrature/quad_loop.py b/hog/quadrature/quad_loop.py index 5c8444997041774c9b7971e9ad2f5d5cd5118af4..2cf1821e45e56c7533162b515832f2144ca4e546 100644 --- a/hog/quadrature/quad_loop.py +++ b/hog/quadrature/quad_loop.py @@ -37,7 +37,7 @@ from hog.fem_helpers import ( jac_blending_inv_eval_symbols, ) from hog.element_geometry import ElementGeometry -from hog.blending import GeometryMap, IdentityMap +from hog.blending import GeometryMap, IdentityMap, ParametricMap class QuadLoop: @@ -104,7 +104,9 @@ class QuadLoop: for dim, symbol in enumerate(ref_symbols) } - if not self.blending.is_affine(): + if not self.blending.is_affine() and not isinstance( + self.blending, ParametricMap + ): jac = jac_blending_evaluate( self.symbolizer, self.quadrature.geometry, self.blending ) @@ -117,7 +119,9 @@ class QuadLoop: hess = hess_blending_evaluate( self.symbolizer, self.quadrature.geometry, self.blending ) - hess_evaluated = [fast_subs(hess[idx], coord_subs_dict) for idx in range(len(hess))] + hess_evaluated = [ + fast_subs(hess[idx], coord_subs_dict) for idx in range(len(hess)) + ] quadrature_assignments += self.blending_hessian_quad_loop_assignments( self.quadrature.geometry, self.symbolizer, hess_evaluated ) @@ -229,7 +233,7 @@ class QuadLoop: ] return quadrature_assignments - + def blending_hessian_quad_loop_assignments( self, geometry: ElementGeometry, @@ -242,7 +246,11 @@ class QuadLoop: dim = geometry.dimensions quadrature_assignments += [ - ast.SympyAssignment(hess_symbols[k][i, j], hessian_blending_evaluated[k][i, j], is_const=False) + ast.SympyAssignment( + hess_symbols[k][i, j], + hessian_blending_evaluated[k][i, j], + is_const=False, + ) for i in range(dim) for j in range(dim) for k in range(dim) diff --git a/hog/quadrature/tabulation.py b/hog/quadrature/tabulation.py index ef409a10fe1862f3783a582ce10eac047fe3d041..61436d2fbee25f316389559a0682b31abf3b781a 100644 --- a/hog/quadrature/tabulation.py +++ b/hog/quadrature/tabulation.py @@ -55,12 +55,23 @@ class Tabulation: self.symbolizer = symbolizer self.tables: Dict[str, Table] = {} - def register_factor(self, factor_name: str, factor: sp.Matrix) -> sp.Matrix: + def register_factor( + self, factor_name: str, factor: sp.Matrix | sp.Expr | int | float + ) -> sp.Matrix | int | float: """Register a factor of the weak form that can be tabulated. Returns symbols replacing the expression for the factor. The symbols are returned in the same form as the factor was given. E.g. in case of a blended full Stokes operator we might encounter J_F^-1 grad phi being a matrix.""" + if isinstance(factor, (int, float)): + return factor + + if isinstance(factor, sp.Expr): + factor = sp.Matrix([factor]) + + if all(f.is_constant() for f in factor): + return factor + replacement_symbols = sp.zeros(factor.rows, factor.cols) for r in range(factor.rows): for c in range(factor.cols):