From 9798fea5edddffafd53ce4e4d792b362e77ffbef Mon Sep 17 00:00:00 2001
From: Nils Kohl <nils.kohl@fau.de>
Date: Mon, 12 Aug 2024 13:32:41 +0200
Subject: [PATCH] Introducing new 'blending' map for iso-, super-,
 sub-parametric operators.

---
 generate_all_operators.py    |   4 +-
 hog/blending.py              |  41 +++++++++-
 hog/fem_helpers.py           |   8 +-
 hog/integrand.py             | 145 ++++++++++++++++++++++++++---------
 hog/quadrature/quad_loop.py  |  18 +++--
 hog/quadrature/tabulation.py |  13 +++-
 6 files changed, 182 insertions(+), 47 deletions(-)

diff --git a/generate_all_operators.py b/generate_all_operators.py
index 005b4d3..282f968 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 01d5a79..d4eea4e 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 855bdd9..4857cdf 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 b6175a3..81586f3 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 5c84449..2cf1821 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 ef409a1..61436d2 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):
-- 
GitLab