diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index 313c5a185c90b50c7f510317793eb1b4e4c0eba8..e59ef72ae8e198132ab6fd19282af32bd3a1ea7b 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -23,9 +23,15 @@ import itertools
import sympy as sp
from sympy.core.cache import clear_cache
import re
+import quadpy
from hog.blending import GeometryMap, IdentityMap, ExternalMap, AnnulusMap
-from hog.element_geometry import TriangleElement, TetrahedronElement, ElementGeometry
+from hog.element_geometry import (
+ TriangleElement,
+ TetrahedronElement,
+ EmbeddedTriangle,
+ ElementGeometry,
+)
from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
from hog.forms import (
mass,
@@ -40,6 +46,16 @@ from hog.forms import (
divdiv,
supg_diffusion,
)
+from hog.manifold_forms import (
+ laplace_beltrami,
+ manifold_mass,
+ manifold_vector_mass,
+ manifold_normal_penalty,
+ manifold_divergence,
+ manifold_vector_divergence,
+ manifold_epsilon,
+ vector_laplace_beltrami,
+)
from hog.forms_vectorial import mass_n1e1, curl_curl
from hog.quadrature import Quadrature, select_quadrule
from hog.exception import HOGException
@@ -433,7 +449,28 @@ form_infos = [
test_degree=2,
quad_schemes={2: 4, 3: 4},
blending=AnnulusMap(),
- )
+ ),
+ FormInfo(
+ "laplace_beltrami",
+ trial_degree=1,
+ test_degree=1,
+ quad_schemes={3: 1},
+ ),
+ FormInfo(
+ "laplace_beltrami",
+ trial_degree=1,
+ test_degree=1,
+ quad_schemes={3: 3},
+ blending=ExternalMap(),
+ ),
+ FormInfo("manifold_mass", trial_degree=1, test_degree=1, quad_schemes={3: 1}),
+ FormInfo(
+ "manifold_mass",
+ trial_degree=1,
+ test_degree=1,
+ quad_schemes={3: 3},
+ blending=ExternalMap(),
+ ),
]
for d in [1, 2]:
@@ -566,6 +603,135 @@ for trial_deg, test_deg, transpose in [
)
)
+for blending in [IdentityMap(), ExternalMap()]:
+ form_infos.append(
+ FormInfo(
+ "manifold_vector_mass",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+ )
+
+ form_infos.append(
+ FormInfo(
+ "manifold_normal_penalty",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+ )
+
+ form_infos.append(
+ FormInfo(
+ "manifold_epsilon",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+ )
+
+ form_infos.append(
+ FormInfo(
+ "manifold_epsilon",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={3: 6},
+ row_dim=3,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+ )
+
+ form_infos.append(
+ FormInfo(
+ "vector_laplace_beltrami",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+ )
+
+for trial_deg, test_deg, transpose in [(1, 2, True), (2, 1, False)]:
+ for blending in [IdentityMap(), ExternalMap()]:
+ if not transpose:
+ form_infos.append(
+ FormInfo(
+ "manifold_div",
+ trial_degree=trial_deg,
+ test_degree=test_deg,
+ quad_schemes={3: 3},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=blending,
+ )
+ )
+ else:
+ form_infos.append(
+ FormInfo(
+ "manifold_divt",
+ trial_degree=trial_deg,
+ test_degree=test_deg,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=blending,
+ )
+ )
+
+for trial_deg, test_deg, transpose in [
+ (1, 2, True),
+ (2, 1, False),
+ (0, 2, True),
+ (2, 0, False),
+]:
+ for blending in [IdentityMap(), ExternalMap()]:
+ if not transpose:
+ form_infos.append(
+ FormInfo(
+ "manifold_vector_div",
+ trial_degree=trial_deg,
+ test_degree=test_deg,
+ quad_schemes={3: 3},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=blending,
+ )
+ )
+ else:
+ form_infos.append(
+ FormInfo(
+ "manifold_vector_divt",
+ trial_degree=trial_deg,
+ test_degree=test_deg,
+ quad_schemes={3: 3},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=blending,
+ )
+ )
+
def form_func(
name: str,
@@ -675,6 +841,94 @@ def form_func(
).integrate(quad, symbolizer)
elif name.startswith("supg_d"):
raise HOGException(f"SUPG Diffusion is not supported for form generation")
+ elif name.startswith("laplace_beltrami"):
+ return laplace_beltrami(
+ trial, test, geometry, symbolizer, blending=blending
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_mass"):
+ return manifold_mass(
+ trial, test, geometry, symbolizer, blending=blending
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_vector_mass"):
+ return manifold_vector_mass(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_trial=col,
+ component_test=row,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_normal_penalty"):
+ return manifold_normal_penalty(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_trial=col,
+ component_test=row,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_divt"):
+ return manifold_divergence(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_index=row,
+ transpose=True,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_div"):
+ return manifold_divergence(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_index=col,
+ transpose=False,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_vector_divt"):
+ return manifold_vector_divergence(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_index=row,
+ transpose=True,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_vector_div"):
+ return manifold_vector_divergence(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_index=col,
+ transpose=False,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("manifold_epsilon"):
+ return manifold_epsilon(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_trial=col,
+ component_test=row,
+ ).integrate(quad, symbolizer)
+ elif name.startswith("vector_laplace_beltrami"):
+ return vector_laplace_beltrami(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_trial=col,
+ component_test=row,
+ ).integrate(quad, symbolizer)
else:
raise HOGException(f"Cannot call form function with name {name}.")
@@ -702,7 +956,7 @@ def parse_arguments():
type=str,
help="build form(s) only for triangle (2D) or tetrahedron (3D) elements; if not speficied we do both",
nargs="?",
- choices=["triangle", "tetrahedron", "both"],
+ choices=["triangle", "tetrahedron", "embedded_triangle", "both"],
default="both",
)
parser.add_argument(
@@ -800,6 +1054,9 @@ def main():
elif args.geometry == "tetrahedron":
logger.info(f"- selected geometry: tetrahedron")
geometries = [TetrahedronElement()]
+ elif args.geometry == "embedded_triangle":
+ logger.info(f"- selected geometry: embedded triangle")
+ geometries = [EmbeddedTriangle()]
else:
logger.info(f"- selected geometries: triangle, tetrahedron")
geometries = [TriangleElement(), TetrahedronElement()]
@@ -838,11 +1095,14 @@ def main():
f"- Generating code for class {form_info.class_name(row,col)}, {geometry.dimensions}D"
):
quad = Quadrature(
- select_quadrule(form_info.quad_schemes[geometry.dimensions], geometry),
- # form_info.quad_schemes[geometry.dimensions],
+ select_quadrule(
+ form_info.quad_schemes[geometry.dimensions],
+ geometry,
+ ),
geometry,
inline_values=form_info.inline_quad,
)
+
mat = form_func(
form_info.form_name,
row,
diff --git a/hog/element_geometry.py b/hog/element_geometry.py
index 0dc03539006370c33d12d1b12ebf6a8e84afb4c4..b69c3e95945c5cdf59ee4e05e46f7d807b9067cd 100644
--- a/hog/element_geometry.py
+++ b/hog/element_geometry.py
@@ -14,6 +14,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
class ElementGeometry:
def __init__(self, dimensions: int, num_vertices: int):
self.dimensions = dimensions
@@ -59,6 +60,17 @@ class TriangleElement(ElementGeometry):
return str(self)
+class EmbeddedTriangle(ElementGeometry):
+ def __init__(self):
+ super().__init__(3, 3)
+
+ def __str__(self):
+ return f"embedded triangle, dim: 3, vertices: 3"
+
+ def __repr__(self):
+ return str(self)
+
+
class TetrahedronElement(ElementGeometry):
def __init__(self):
super().__init__(3, 4)
@@ -67,4 +79,4 @@ class TetrahedronElement(ElementGeometry):
return f"tetrahedron, dim: 3, vertices: 4"
def __repr__(self):
- return str(self)
\ No newline at end of file
+ return str(self)
diff --git a/hog/external_functions.py b/hog/external_functions.py
index c22943568e72627f142cdb35cc116b95c8c5c8c6..0180c9cc6ef912980de4ded5eebcc15989212aab 100644
--- a/hog/external_functions.py
+++ b/hog/external_functions.py
@@ -35,6 +35,25 @@ geometryMap_->evalF( in, out );
*out_1 = out[1];"""
+class BlendingFEmbeddedTriangle(MultiAssignment):
+ def function_name(self):
+ return "Blending_F_EmbeddedTriangle"
+
+ def num_input_args(self):
+ return 3
+
+ def num_output_args(self):
+ return 3
+
+ def implementation(self):
+ return """Point3D in( {in_0, in_1, in_2} );
+Point3D out;
+geometryMap_->evalF( in, out );
+*out_0 = out[0];
+*out_1 = out[1];
+*out_2 = out[2];"""
+
+
class BlendingFTetrahedron(MultiAssignment):
def function_name(self):
return "Blending_F_Tetrahedron"
@@ -76,6 +95,33 @@ geometryMap_->evalDF( mappedPt, DPsi );
*out_3 = DPsi( 1, 1 );"""
+class BlendingDFEmbeddedTriangle(MultiAssignment):
+ def function_name(self):
+ return "Blending_DF_EmbeddedTriangle"
+
+ @classmethod
+ def num_input_args(cls):
+ return 3
+
+ @classmethod
+ def num_output_args(cls):
+ return 3 * 3
+
+ def implementation(self):
+ return """Point3D mappedPt( {in_0, in_1, in_2} );
+Matrix3r DPsi;
+geometryMap_->evalDF( mappedPt, DPsi );
+*out_0 = DPsi( 0, 0 );
+*out_1 = DPsi( 0, 1 );
+*out_2 = DPsi( 0, 2 );
+*out_3 = DPsi( 1, 0 );
+*out_4 = DPsi( 1, 1 );
+*out_5 = DPsi( 1, 2 );
+*out_6 = DPsi( 2, 0 );
+*out_7 = DPsi( 2, 1 );
+*out_8 = DPsi( 2, 2 );"""
+
+
class BlendingDFInvDFTriangle(MultiAssignment):
def function_name(self):
return "Blending_DFInvDF_Triangle"
diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py
index 13c5dff89aad186e7034e229812d13a5678e1f5b..a6a525fbb850ea8cdc317ec794c428d6c40c8bbe 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -28,6 +28,7 @@ from hog.blending import (
from hog.element_geometry import (
ElementGeometry,
TriangleElement,
+ EmbeddedTriangle,
TetrahedronElement,
LineElement,
)
@@ -38,10 +39,12 @@ 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,
VectorVariableCoefficient3D,
@@ -115,6 +118,9 @@ def trafo_ref_to_affine(
vector symbols, useful for example if the trafo of two or more different elements is required
"""
ref_symbols_vector = symbolizer.ref_coords_as_vector(geometry.dimensions)
+ if isinstance(geometry, EmbeddedTriangle):
+ ref_symbols_vector = symbolizer.ref_coords_as_vector(geometry.dimensions - 1)
+
if affine_points is None:
affine_points = symbolizer.affine_vertices_as_vectors(
geometry.dimensions, geometry.num_vertices
@@ -193,6 +199,9 @@ def jac_ref_to_affine(
vector symbols, useful for example if the trafo of two or more different elements is required
"""
ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions)
+ if isinstance(geometry, EmbeddedTriangle):
+ ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions - 1)
+
trafo = trafo_ref_to_affine(geometry, symbolizer, affine_points=affine_points)
return trafo.jacobian(ref_symbols_list)
@@ -211,6 +220,8 @@ def trafo_ref_to_physical(
blending_class: type[MultiAssignment]
if isinstance(geometry, TriangleElement):
blending_class = BlendingFTriangle
+ elif isinstance(geometry, EmbeddedTriangle):
+ blending_class = BlendingDFEmbeddedTriangle
elif isinstance(geometry, TetrahedronElement):
blending_class = BlendingFTetrahedron
else:
@@ -235,6 +246,8 @@ def jac_affine_to_physical(
blending_class: type[MultiAssignment]
if isinstance(geometry, TriangleElement):
blending_class = BlendingDFTriangle
+ elif isinstance(geometry, EmbeddedTriangle):
+ blending_class = BlendingDFEmbeddedTriangle
elif isinstance(geometry, TetrahedronElement):
blending_class = BlendingDFTetrahedron
else:
@@ -267,6 +280,7 @@ def jac_blending_evaluate(
jac = blending.jacobian(trafo_ref_to_affine(geometry, symbolizer, affine_points))
return jac
+
def hess_blending_evaluate(
symbolizer: Symbolizer, geometry: ElementGeometry, blending: GeometryMap
) -> List[sp.Matrix]:
@@ -276,6 +290,7 @@ def hess_blending_evaluate(
hess = blending.hessian(trafo_ref_to_affine(geometry, symbolizer, affine_points))
return hess
+
def abs_det_jac_blending_eval_symbols(
geometry: ElementGeometry, symbolizer: Symbolizer, q_pt: str = ""
) -> sp.Expr:
diff --git a/hog/function_space.py b/hog/function_space.py
index 7d0823a103fb946f7fda56c8dd33657a405e3cb3..8ee718ed5551fb85631fc5a49066fe16fa974889 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -18,7 +18,12 @@ from pyclbr import Function
from typing import Any, List, Optional, Protocol
import sympy as sp
-from hog.element_geometry import ElementGeometry, TriangleElement, TetrahedronElement
+from hog.element_geometry import (
+ ElementGeometry,
+ TriangleElement,
+ EmbeddedTriangle,
+ TetrahedronElement,
+)
from hog.exception import HOGException
from hog.math_helpers import grad, hessian
from hog.symbolizer import Symbolizer
@@ -84,6 +89,8 @@ class FunctionSpace(Protocol):
"""
if domain in ["ref", "reference"]:
symbols = self.symbolizer.ref_coords_as_list(geometry.dimensions)
+ if isinstance(geometry, EmbeddedTriangle):
+ symbols = self.symbolizer.ref_coords_as_list(geometry.dimensions - 1)
basis_functions_gradients = [
grad(f, symbols)
for f in self.shape(geometry, domain=domain, dof_map=dof_map)
@@ -93,7 +100,7 @@ class FunctionSpace(Protocol):
raise HOGException(
f"Gradient of shape function not available for domain type {domain}"
)
-
+
def hessian_shape(
self,
geometry: ElementGeometry,
@@ -208,6 +215,40 @@ class LagrangianFunctionSpace(FunctionSpace):
-4 * x**2 - 4 * x * y + 4 * x,
]
+ elif (
+ isinstance(geometry, EmbeddedTriangle)
+ and self.family in ["Lagrange"]
+ and self._degree == 0
+ ):
+ basis_functions = [sp.sympify(1)]
+
+ elif (
+ isinstance(geometry, EmbeddedTriangle)
+ and self.family in ["Lagrange"]
+ and self._degree == 1
+ ):
+ basis_functions = [
+ 1 - symbols[0] - symbols[1],
+ symbols[0],
+ symbols[1],
+ ]
+
+ elif (
+ isinstance(geometry, EmbeddedTriangle)
+ and self.family in ["Lagrange"]
+ and self._degree == 2
+ ):
+ x = symbols[0]
+ y = symbols[1]
+ basis_functions = [
+ 2 * x**2 + 4 * x * y - 3 * x + 2 * y**2 - 3 * y + 1,
+ 2 * x**2 - x,
+ 2 * y**2 - y,
+ 4 * x * y,
+ -4 * x * y - 4 * y**2 + 4 * y,
+ -4 * x**2 - 4 * x * y + 4 * x,
+ ]
+
elif (
isinstance(geometry, TetrahedronElement)
and self.family in ["Lagrange"]
diff --git a/hog/manifold_forms.py b/hog/manifold_forms.py
new file mode 100644
index 0000000000000000000000000000000000000000..f53fcfed70272186bfd12aa2cc7afb418bc7bb9c
--- /dev/null
+++ b/hog/manifold_forms.py
@@ -0,0 +1,661 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024 HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+import logging
+import sympy as sp
+
+from hog.element_geometry import ElementGeometry
+from hog.element_geometry import EmbeddedTriangle
+from hog.exception import HOGException
+from hog.fem_helpers import (
+ trafo_ref_to_affine,
+ trafo_ref_to_physical,
+ jac_ref_to_affine,
+ jac_affine_to_physical,
+ create_empty_element_matrix,
+ element_matrix_iterator,
+)
+from hog.function_space import FunctionSpace
+from hog.math_helpers import dot, inv, abs, det, double_contraction, e_vec
+from hog.quadrature import Tabulation
+from hog.symbolizer import Symbolizer
+from hog.logger import TimedLogger
+from hog.blending import GeometryMap, ExternalMap, IdentityMap
+from hog.forms import Form
+
+from hog.manifold_helpers import face_projection, embedded_normal
+
+
+def laplace_beltrami(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+Laplace beltrami operator.
+
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+ G: first fundamental form of manifold
+
+ ∫ ∇u · G^(-1) · ∇v · (det(G))^0.5
+"""
+
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble laplace beltrami matrix."
+ )
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Laplace Beltrami only works for embedded triangles")
+
+ with TimedLogger("assembling laplace beltrami matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = blending.jacobian(trafo_ref_to_affine(geometry, symbolizer))
+
+ fundamental_form = jac_affine.T * jac_blending.T * jac_blending * jac_affine
+ with TimedLogger("inverting first fundamental form", level=logging.DEBUG):
+ fundamental_form_inv = inv(fundamental_form)
+ fundamental_form_det = abs(det(fundamental_form))
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for data in it:
+ with TimedLogger(
+ f"Integrating row = {data.row} , col = {data.col}",
+ level=logging.DEBUG,
+ ):
+ grad_phi = data.trial_shape_grad
+ grad_psi = data.test_shape_grad
+ laplace_beltrami_symbol = sp.Matrix(
+ tabulation.register_factor(
+ "laplace_beltrami_symbol",
+ dot(
+ grad_phi,
+ fundamental_form_inv * grad_psi,
+ )
+ * (fundamental_form_det**0.5),
+ )
+ )
+ form = laplace_beltrami_symbol[0]
+ mat[data.row, data.col] = form
+
+ return Form(mat, tabulation, symmetric=True, docstring=docstring)
+
+
+def manifold_mass(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+Manifold mass operator.
+
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+ G: first fundamental form of manifold
+
+ ∫ uv · (det(G))^0.5
+"""
+
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble laplace beltrami matrix."
+ )
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ with TimedLogger("assembling laplace beltrami matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = blending.jacobian(trafo_ref_to_affine(geometry, symbolizer))
+
+ fundamental_form_det = abs(
+ det(jac_affine.T * jac_blending.T * jac_blending * jac_affine)
+ )
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for data in it:
+ with TimedLogger(
+ f"Integrating row = {data.row} , col = {data.col}",
+ level=logging.DEBUG,
+ ):
+ phi = data.trial_shape
+ psi = data.test_shape
+ manifold_mass_symbol = sp.Matrix(
+ tabulation.register_factor(
+ "manifold_mass_symbol",
+ sp.Matrix([phi * psi * fundamental_form_det**0.5]),
+ )
+ )
+ form = manifold_mass_symbol[0]
+ mat[data.row, data.col] = form
+
+ return Form(mat, tabulation, symmetric=True, docstring=docstring)
+
+
+def manifold_vector_mass(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+) -> Form:
+ docstring = f"""
+Manifold vector mass operator operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ P: projection matrix onto face
+
+ ∫ Pu · Pv · (det(G))^0.5
+
+"""
+ with TimedLogger("assembling manifold vector mass matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ projection = face_projection(geometry, symbolizer, blending=blending)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ fundamental_form_det = abs(
+ det(jac_affine.T * jac_blending.T * jac_blending * jac_affine)
+ )
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ phi = data.trial_shape
+ psi = data.test_shape
+ phi_vec = e_vec(geometry.dimensions, component_trial) * phi
+ psi_vec = e_vec(geometry.dimensions, component_test) * psi
+ projected_phi = projection * phi_vec
+ projected_psi = projection * psi_vec
+
+ form = sp.Matrix(
+ tabulation.register_factor(
+ "manifold_vector_mass_symbol",
+ dot(projected_phi, projected_psi) * fundamental_form_det**0.5,
+ )
+ )[0]
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=component_trial == component_test,
+ docstring=docstring,
+ )
+
+
+def manifold_normal_penalty(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+) -> Form:
+ docstring = f"""
+Manifold normal penaly operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ n: normal vector of the face
+
+ ∫ (n · u)(n · v) · (det(G))^0.5
+
+"""
+ with TimedLogger("assembling manifold normal penalty matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ normal = embedded_normal(geometry, symbolizer, blending)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ fundamental_form_det = abs(
+ det(jac_affine.T * jac_blending.T * jac_blending * jac_affine)
+ )
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ phi = data.trial_shape
+ psi = data.test_shape
+ phi_vec = e_vec(geometry.dimensions, component_trial) * phi
+ psi_vec = e_vec(geometry.dimensions, component_test) * psi
+ phi_normal = dot(phi_vec, normal)
+ psi_normal = dot(psi_vec, normal)
+
+ form = sp.Matrix(
+ tabulation.register_factor(
+ "manifold_normal_penalty_symbol",
+ phi_normal * psi_normal * fundamental_form_det**0.5,
+ )
+ )[0]
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=component_trial == component_test,
+ docstring=docstring,
+ )
+
+
+def manifold_divergence(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_index: int = 0,
+ transpose: bool = False,
+) -> Form:
+ docstring = f"""
+Manifold divergence operator.
+
+Component index: {component_index}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (space: {test})
+
+ ∫ u · grad_Gamma(v) · (det(G))^0.5
+
+"""
+
+ with TimedLogger("assembling manifold div matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ jac_total = jac_blending * jac_affine
+ fundamental_form = jac_total.T * jac_total
+ fundamental_form_inv = inv(fundamental_form)
+ fundamental_form_det = abs(det(fundamental_form))
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ if not transpose:
+ phi = data.trial_shape
+ phi_grad = data.test_shape_grad
+ else:
+ phi = data.test_shape
+ phi_grad = data.trial_shape_grad
+
+ form = (
+ (jac_total * fundamental_form_inv * phi_grad)[component_index]
+ * phi
+ * fundamental_form_det**0.5
+ )
+
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=False,
+ docstring=docstring,
+ )
+
+
+def manifold_vector_divergence(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_index: int = 0,
+ transpose: bool = False,
+) -> Form:
+ docstring = f"""
+Manifold vector divergence operator.
+
+Component index: {component_index}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (space: {test})
+ P: projection matrix
+
+ ∫ div_Gamma(Pu) · v · (det(G))^0.5
+
+"""
+
+ with TimedLogger("assembling manifold div matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ logging.info(
+ f"WARNING: Manifold vector divergence does NOT compute derivative of matrix P yet. Generated form might not work as intended."
+ )
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ projection_mat = face_projection(geometry, symbolizer, blending=blending)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ jac_total = jac_blending * jac_affine
+ fundamental_form = jac_total.T * jac_total
+ fundamental_form_inv = inv(fundamental_form)
+ fundamental_form_det = abs(det(fundamental_form))
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions - 1)
+ if not transpose:
+ phi_vec = (
+ projection_mat
+ * (
+ e_vec(geometry.dimensions, component_index) * data.trial_shape
+ ).jacobian(ref_symbols_list)
+ ).T
+ phi = data.test_shape
+ else:
+ phi = data.trial_shape
+ phi_vec = (
+ projection_mat
+ * (
+ e_vec(geometry.dimensions, component_index) * data.test_shape
+ ).jacobian(ref_symbols_list)
+ ).T
+
+ form = (
+ (jac_total * fundamental_form_inv * phi_vec).trace()
+ * phi
+ * fundamental_form_det**0.5
+ )
+
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=False,
+ docstring=docstring,
+ )
+
+
+def manifold_epsilon(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+) -> Form:
+ docstring = f"""
+Manifold epsilon operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ G: 2D manifold embedded in 3D space
+
+ ∫ epsilon_Gamma(Pu) : epsilon_Gamma(Pv) · (det(G))^0.5
+
+"""
+
+ with TimedLogger("assembling manifold epsilon matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ projection_mat = face_projection(geometry, symbolizer, blending=blending)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ jac_total = jac_blending * jac_affine
+ fundamental_form = jac_total.T * jac_total
+ fundamental_form_inv = inv(fundamental_form)
+ fundamental_form_det = abs(det(fundamental_form))
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions - 1)
+ phi = data.trial_shape
+ psi = data.test_shape
+ unscaled_phi_projected_grad = (
+ projection_mat
+ * (e_vec(geometry.dimensions, component_trial) * phi).jacobian(
+ ref_symbols_list
+ )
+ ).T
+ unscaled_psi_projected_grad = (
+ projection_mat
+ * (e_vec(geometry.dimensions, component_test) * psi).jacobian(
+ ref_symbols_list
+ )
+ ).T
+
+ phi_projected_grad = (
+ jac_total * fundamental_form_inv * unscaled_phi_projected_grad
+ )
+ psi_projected_grad = (
+ jac_total * fundamental_form_inv * unscaled_psi_projected_grad
+ )
+
+ phi_epsilon = 0.5 * (phi_projected_grad + phi_projected_grad.T)
+ psi_epsilon = 0.5 * (psi_projected_grad + psi_projected_grad.T)
+
+ form = tabulation.register_factor(
+ "epsilon_epsilon_prod",
+ double_contraction(phi_epsilon, psi_epsilon)
+ * fundamental_form_det**0.5,
+ )[0]
+
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=component_trial == component_test,
+ docstring=docstring,
+ )
+
+
+def vector_laplace_beltrami(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_trial: int = 0,
+ component_test: int = 0,
+) -> Form:
+ docstring = f"""
+Manifold vector laplace beltrami operator.
+
+Component trial: {component_trial}
+Component test: {component_test}
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (vectorial space: {trial})
+ v: test function (vectorial space: {test})
+ G: 2D manifold embedded in 3D space
+
+ ∫ grad_Gamma(u) : grad_Gamma(v) · (det(G))^0.5
+
+"""
+
+ with TimedLogger("assembling vector laplace beltrami matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException("Manifold forms only work for embedded triangles.")
+
+ projection_mat = face_projection(geometry, symbolizer, blending=blending)
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer)
+
+ if isinstance(blending, ExternalMap):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ jac_blending = sp.eye(3)
+
+ jac_total = jac_blending * jac_affine
+ fundamental_form = jac_total.T * jac_total
+ fundamental_form_inv = inv(fundamental_form)
+ fundamental_form_det = abs(det(fundamental_form))
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions - 1)
+ phi = data.trial_shape
+ psi = data.test_shape
+ unscaled_phi_projected_grad = (
+ projection_mat
+ * (e_vec(geometry.dimensions, component_trial) * phi).jacobian(
+ ref_symbols_list
+ )
+ ).T
+ unscaled_psi_projected_grad = (
+ projection_mat
+ * (e_vec(geometry.dimensions, component_test) * psi).jacobian(
+ ref_symbols_list
+ )
+ ).T
+
+ phi_projected_grad = (
+ jac_total * fundamental_form_inv * unscaled_phi_projected_grad
+ )
+ psi_projected_grad = (
+ jac_total * fundamental_form_inv * unscaled_psi_projected_grad
+ )
+
+ form = tabulation.register_factor(
+ "phi_psi_prod",
+ double_contraction(phi_projected_grad, psi_projected_grad)
+ * fundamental_form_det**0.5,
+ )[0]
+
+ mat[data.row, data.col] = form
+
+ return Form(
+ mat,
+ tabulation,
+ symmetric=component_trial == component_test,
+ docstring=docstring,
+ )
diff --git a/hog/manifold_helpers.py b/hog/manifold_helpers.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d2bf663ab74c27ceb6f1ed96d9b5914e4a3a855
--- /dev/null
+++ b/hog/manifold_helpers.py
@@ -0,0 +1,83 @@
+# HyTeG Operator Generator
+# Copyright (C) 2024 HyTeG Team
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+from typing import List
+import sympy as sp
+
+from hog.exception import HOGException
+
+from hog.element_geometry import ElementGeometry, EmbeddedTriangle
+from hog.symbolizer import Symbolizer
+from hog.blending import GeometryMap, IdentityMap
+from hog.fem_helpers import jac_affine_to_physical
+from hog.external_functions import BlendingFEmbeddedTriangle
+
+
+def embedded_normal(
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+) -> sp.Matrix:
+ """Returns an unoriented unit normal vector for an embedded triangle."""
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException(
+ "Embedded normal vectors are only defined for embedded triangles."
+ )
+
+ vert_points = symbolizer.affine_vertices_as_vectors(
+ geometry.dimensions, geometry.num_vertices
+ )
+ span0 = vert_points[1] - vert_points[0]
+ span1 = vert_points[2] - vert_points[0]
+
+ if not isinstance(blending, IdentityMap):
+ blending_jac = jac_affine_to_physical(geometry, symbolizer)
+ span0 = blending_jac.T * span0
+ span1 = blending_jac.T * span1
+
+ normal = sp.Matrix(
+ [
+ span0[1] * span1[2] - span0[2] * span1[1],
+ span0[2] * span1[0] - span0[0] * span1[2],
+ span0[0] * span1[1] - span0[1] * span1[0],
+ ]
+ )
+ normal = normal / (normal[0] ** 2 + normal[1] ** 2 + normal[2] ** 2) ** 0.5
+ return normal
+
+
+def face_projection(
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+) -> sp.Matrix:
+ """Returns a projection matrix for an embedded triangle."""
+
+ if not isinstance(geometry, EmbeddedTriangle):
+ raise HOGException(
+ "Projection matrices are only defined for embedded triangles."
+ )
+
+ normal = embedded_normal(geometry, symbolizer, blending=blending)
+ projection = sp.Matrix(
+ [
+ [1.0 - normal[0] ** 2, -normal[0] * normal[1], -normal[0] * normal[2]],
+ [-normal[0] * normal[1], 1.0 - normal[1] ** 2, -normal[1] * normal[2]],
+ [-normal[0] * normal[2], -normal[1] * normal[2], 1.0 - normal[2] ** 2],
+ ]
+ )
+ return projection
diff --git a/hog/quadrature/quadrature.py b/hog/quadrature/quadrature.py
index 5526cc700162ecade98047372725c16500cb55e5..b51dd0bfce6014ac3e10c0ce05c89694295d5c20 100644
--- a/hog/quadrature/quadrature.py
+++ b/hog/quadrature/quadrature.py
@@ -26,6 +26,7 @@ import sympy as sp
from hog.element_geometry import (
ElementGeometry,
TriangleElement,
+ EmbeddedTriangle,
TetrahedronElement,
LineElement,
)
@@ -95,7 +96,9 @@ def select_quadrule(
warnings.simplefilter("ignore")
all_schemes = []
- if isinstance(geometry, TriangleElement):
+ if isinstance(geometry, TriangleElement) or isinstance(
+ geometry, EmbeddedTriangle
+ ):
schemes = quadpy.t2.schemes
elif isinstance(geometry, TetrahedronElement):
schemes = quadpy.t3.schemes
@@ -223,6 +226,8 @@ class Quadrature:
f"Cannot apply quadrature rule to matrix of shape {f.shape}: {f}."
)
ref_symbols = symbolizer.ref_coords_as_list(self._geometry.dimensions)
+ if isinstance(self._geometry, EmbeddedTriangle):
+ ref_symbols = symbolizer.ref_coords_as_list(self._geometry.dimensions - 1)
if self._scheme_name == "exact":
mat_entry = integrate_exact_over_reference(f, self._geometry, symbolizer)
@@ -286,6 +291,9 @@ class Quadrature:
elif isinstance(geometry, LineElement):
vertices = np.asarray([[0.0], [1.0]])
degree = scheme.degree
+ elif isinstance(geometry, EmbeddedTriangle):
+ vertices = np.asarray([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
+ degree = scheme.degree
else:
raise HOGException("Invalid geometry for quadrature.")