diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index e59ef72ae8e198132ab6fd19282af32bd3a1ea7b..7f9341b5b1079a54b0ffabe628bf9917ad3da426 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -29,7 +29,6 @@ from hog.blending import GeometryMap, IdentityMap, ExternalMap, AnnulusMap
from hog.element_geometry import (
TriangleElement,
TetrahedronElement,
- EmbeddedTriangle,
ElementGeometry,
)
from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
@@ -463,12 +462,12 @@ form_infos = [
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={2: 1}),
FormInfo(
"manifold_mass",
trial_degree=1,
test_degree=1,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
blending=ExternalMap(),
),
]
@@ -609,7 +608,7 @@ for blending in [IdentityMap(), ExternalMap()]:
"manifold_vector_mass",
trial_degree=2,
test_degree=2,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=3,
is_implemented=is_implemented_for_vector_to_vector,
@@ -622,7 +621,7 @@ for blending in [IdentityMap(), ExternalMap()]:
"manifold_normal_penalty",
trial_degree=2,
test_degree=2,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=3,
is_implemented=is_implemented_for_vector_to_vector,
@@ -635,7 +634,7 @@ for blending in [IdentityMap(), ExternalMap()]:
"manifold_epsilon",
trial_degree=2,
test_degree=2,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=3,
is_implemented=is_implemented_for_vector_to_vector,
@@ -648,7 +647,7 @@ for blending in [IdentityMap(), ExternalMap()]:
"manifold_epsilon",
trial_degree=2,
test_degree=2,
- quad_schemes={3: 6},
+ quad_schemes={2: 6},
row_dim=3,
col_dim=3,
is_implemented=is_implemented_for_vector_to_vector,
@@ -661,7 +660,7 @@ for blending in [IdentityMap(), ExternalMap()]:
"vector_laplace_beltrami",
trial_degree=2,
test_degree=2,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=3,
is_implemented=is_implemented_for_vector_to_vector,
@@ -677,7 +676,7 @@ for trial_deg, test_deg, transpose in [(1, 2, True), (2, 1, False)]:
"manifold_div",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=1,
col_dim=3,
is_implemented=is_implemented_for_vector_to_scalar,
@@ -690,7 +689,7 @@ for trial_deg, test_deg, transpose in [(1, 2, True), (2, 1, False)]:
"manifold_divt",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=1,
is_implemented=is_implemented_for_scalar_to_vector,
@@ -711,7 +710,7 @@ for trial_deg, test_deg, transpose in [
"manifold_vector_div",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=1,
col_dim=3,
is_implemented=is_implemented_for_vector_to_scalar,
@@ -724,7 +723,7 @@ for trial_deg, test_deg, transpose in [
"manifold_vector_divt",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={3: 3},
+ quad_schemes={2: 3},
row_dim=3,
col_dim=1,
is_implemented=is_implemented_for_scalar_to_vector,
@@ -1056,7 +1055,7 @@ def main():
geometries = [TetrahedronElement()]
elif args.geometry == "embedded_triangle":
logger.info(f"- selected geometry: embedded triangle")
- geometries = [EmbeddedTriangle()]
+ geometries = [TriangleElement(space_dimension=3)]
else:
logger.info(f"- selected geometries: triangle, tetrahedron")
geometries = [TriangleElement(), TetrahedronElement()]
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 2f81aeefce7ebb8c46278109dcb12d4b29bb5e3b..dd8efae68b8016c338bdac7f1114981062ae2a97 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -704,7 +704,6 @@ def generate_elementwise_op(
operator = HyTeGElementwiseOperator(
name,
symbolizer,
- opts=optimizations,
kernel_wrapper_types=op_info.kernel_types,
type_descriptor=type_descriptor,
)
@@ -727,15 +726,14 @@ def generate_elementwise_op(
blending=blending, # type: ignore[call-arg] # kw-args are not supported by Callable
)
- operator.add_integral(
+ operator.add_volume_integral(
name="".join(name.split()),
- dim=geometry.dimensions,
- geometry=geometry,
- integration_domain=MacroIntegrationDomain.VOLUME,
+ volume_geometry=geometry,
quad=quad,
blending=blending,
form=form,
loop_strategy=loop_strategy,
+ optimizations=optimizations,
)
dir_path = os.path.join(args.output, op_info.name.split("_")[0])
diff --git a/hog/blending.py b/hog/blending.py
index 7a8f884e0924f6882fdf65593e416c38fe0c1bf8..01d5a79c0f0513ec4c00c3783da2afacfbd659ce 100644
--- a/hog/blending.py
+++ b/hog/blending.py
@@ -148,7 +148,7 @@ class AnnulusMap(GeometryMap):
"""Evaluates the Jacobian of the geometry map at the passed point."""
if sp.shape(x) != (2, 1):
- raise HOGException("Invalid input shape for AnnulusMap.")
+ raise HOGException(f"Invalid input shape {sp.shape(x)} for AnnulusMap.")
radRefVertex = self.radRefVertex
radRayVertex = self.radRayVertex
diff --git a/hog/code_generation.py b/hog/code_generation.py
index e20e3197b26b988b002ed119d87ff6e4c9a0fd09..afa5aaea38a5a0a1a11b1ec72239c4cb747b9234 100644
--- a/hog/code_generation.py
+++ b/hog/code_generation.py
@@ -160,9 +160,9 @@ def replace_multi_assignments(
# Actually filling the dict.
for ma in multi_assignments:
replacement_symbol = replacement_symbols[ma.output_arg()]
- multi_assignments_replacement_symbols[
- ma.unique_identifier
- ] = replacement_symbol
+ multi_assignments_replacement_symbols[ma.unique_identifier] = (
+ replacement_symbol
+ )
if multi_assignments_replacement_symbols:
with TimedLogger(
@@ -238,9 +238,15 @@ def jacobi_matrix_assignments(
assignments = []
- jac_aff_symbol = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_aff_inv_symbol = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
- jac_aff_det_symbol = symbolizer.abs_det_jac_ref_to_affine()
+ # If the reference space dimension (geometry.dimension) is not equal to the space dimension
+ # the Jacobian of the affine map is not square. Processing its inverse and determinant is
+ # thus not meaningful.
+ is_affine_jac_square = geometry.dimensions == geometry.space_dimension
+
+ jac_aff_symbol = symbolizer.jac_ref_to_affine(geometry)
+ if is_affine_jac_square:
+ jac_aff_inv_symbol = symbolizer.jac_ref_to_affine_inv(geometry)
+ jac_aff_det_symbol = symbolizer.abs_det_jac_ref_to_affine()
free_symbols = element_matrix.free_symbols | {
free_symbol
@@ -250,27 +256,28 @@ def jacobi_matrix_assignments(
}
# Steps 1 and 2.
- jac_affine_inv_in_expr = set(jac_aff_inv_symbol).intersection(free_symbols)
- abs_det_jac_affine_in_expr = jac_aff_det_symbol in free_symbols
-
- # Just an early exit. Not strictly required, but might accelerate this process in some cases.
- if jac_affine_inv_in_expr:
- jac_aff_inv_expr = jac_aff_symbol.inv()
- for s_ij, e_ij in zip(jac_aff_inv_symbol, jac_aff_inv_expr):
- if s_ij in jac_affine_inv_in_expr:
- assignments.append(SympyAssignment(s_ij, e_ij))
+ if is_affine_jac_square:
+ jac_affine_inv_in_expr = set(jac_aff_inv_symbol).intersection(free_symbols)
+ abs_det_jac_affine_in_expr = jac_aff_det_symbol in free_symbols
+
+ if jac_affine_inv_in_expr:
+ jac_aff_inv_expr = jac_aff_symbol.inv()
+ for s_ij, e_ij in zip(jac_aff_inv_symbol, jac_aff_inv_expr):
+ if s_ij in jac_affine_inv_in_expr:
+ assignments.append(SympyAssignment(s_ij, e_ij))
- if abs_det_jac_affine_in_expr:
- assignments.append(
- SympyAssignment(jac_aff_det_symbol, sp.Abs(jac_aff_symbol.det()))
- )
+ if abs_det_jac_affine_in_expr:
+ assignments.append(
+ SympyAssignment(jac_aff_det_symbol, sp.Abs(jac_aff_symbol.det()))
+ )
- # Collecting all expressions to parse for step 3.
- free_symbols |= {free_symbol for a in assignments for free_symbol in a.rhs.atoms()}
+ # Collecting all expressions to parse for step 3.
+ free_symbols |= {
+ free_symbol for a in assignments for free_symbol in a.rhs.atoms()
+ }
jac_affine_in_expr = set(jac_aff_symbol).intersection(free_symbols)
- # Just an early exit. Not strictly required, but might accelerate this process in some cases.
if jac_affine_in_expr:
jac_aff_expr = jac_ref_to_affine(geometry, symbolizer, affine_points)
for s_ij, e_ij in zip(jac_aff_symbol, jac_aff_expr):
@@ -335,7 +342,6 @@ def blending_jacobi_matrix_assignments(
jac_blending_in_expr = set(jac_blend_symbol).intersection(free_symbols)
- # Just an early exit. Not strictly required, but might accelerate this process in some cases.
if jac_blending_in_expr:
if isinstance(blending, ExternalMap):
HOGException("Not implemented or cannot be?")
@@ -344,7 +350,7 @@ def blending_jacobi_matrix_assignments(
# Collecting all expressions to parse for step 3.
spat_coord_subs = {}
for idx, symbol in enumerate(
- symbolizer.ref_coords_as_list(geometry.dimensions)
+ symbolizer.ref_coords_as_list(quad_info.geometry.dimensions)
):
spat_coord_subs[symbol] = point[idx]
jac_blend_expr_sub = jac_blend_expr.subs(spat_coord_subs)
diff --git a/hog/element_geometry.py b/hog/element_geometry.py
index b69c3e95945c5cdf59ee4e05e46f7d807b9067cd..9fe91fc1b7a7e0f886003a20beec765e4f0862c7 100644
--- a/hog/element_geometry.py
+++ b/hog/element_geometry.py
@@ -14,69 +14,68 @@
# 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 hog.exception import HOGException
+
class ElementGeometry:
- def __init__(self, dimensions: int, num_vertices: int):
+ def __init__(self, dimensions: int, num_vertices: int, space_dimension: int):
+
+ if space_dimension < dimensions:
+ raise HOGException(
+ "The space dimension should be larger or equal to the dimension of the geometry."
+ )
+
self.dimensions = dimensions
self.num_vertices = num_vertices
+ self.space_dimension = space_dimension
def __str__(self):
- return f"ElementGeometry(dim: {self.dimensions}, vertices: {self.num_vertices})"
+ return f"ElementGeometry(dim: {self.dimensions}, vertices: {self.num_vertices}), space dim: {self.space_dimension}"
def __repr__(self):
return str(self)
def __hash__(self):
- return hash((self.dimensions, self.num_vertices))
+ return hash((self.dimensions, self.num_vertices, self.space_dimension))
def __eq__(self, other):
if isinstance(other, ElementGeometry):
return (
self.dimensions == other.dimensions
and self.num_vertices == other.num_vertices
+ and self.space_dimension == other.space_dimension
)
return False
class LineElement(ElementGeometry):
- def __init__(self):
- super().__init__(1, 2)
+ def __init__(self, space_dimension: int = 1):
+ super().__init__(1, 2, space_dimension=space_dimension)
def __str__(self):
- return f"line, dim: 1, vertices: 2"
+ return f"line, dim: 1, vertices: 2, spacedim: {self.space_dimension}"
def __repr__(self):
return str(self)
class TriangleElement(ElementGeometry):
- def __init__(self):
- super().__init__(2, 3)
-
- def __str__(self):
- return f"triangle, dim: 2, vertices: 3"
-
- def __repr__(self):
- return str(self)
-
-
-class EmbeddedTriangle(ElementGeometry):
- def __init__(self):
- super().__init__(3, 3)
+ def __init__(self, space_dimension: int = 2):
+ super().__init__(2, 3, space_dimension=space_dimension)
def __str__(self):
- return f"embedded triangle, dim: 3, vertices: 3"
+ return f"triangle, dim: 2, vertices: 3, spacedim: {self.space_dimension}"
def __repr__(self):
return str(self)
class TetrahedronElement(ElementGeometry):
- def __init__(self):
- super().__init__(3, 4)
+ def __init__(self, space_dimension: int = 3):
+ super().__init__(3, 4, space_dimension=space_dimension)
def __str__(self):
- return f"tetrahedron, dim: 3, vertices: 4"
+ return f"tetrahedron, dim: 3, vertices: 4, spacedim: {self.space_dimension}"
def __repr__(self):
return str(self)
diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py
index a6a525fbb850ea8cdc317ec794c428d6c40c8bbe..3cbb2be5dd2efcc32a1f8e47d0759bdad8f53c09 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -28,7 +28,6 @@ from hog.blending import (
from hog.element_geometry import (
ElementGeometry,
TriangleElement,
- EmbeddedTriangle,
TetrahedronElement,
LineElement,
)
@@ -118,12 +117,10 @@ 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
+ geometry.space_dimension, geometry.num_vertices
)
else:
if len(affine_points) != geometry.num_vertices:
@@ -135,12 +132,12 @@ def trafo_ref_to_affine(
affine_points[p][d] - affine_points[0][d]
for p in range(1, geometry.num_vertices)
]
- for d in range(geometry.dimensions)
+ for d in range(geometry.space_dimension)
]
)
trafo = trafo * ref_symbols_vector
trafo = trafo + sp.Matrix(
- [[affine_points[0][d]] for d in range(geometry.dimensions)]
+ [[affine_points[0][d]] for d in range(geometry.space_dimension)]
)
return trafo
@@ -199,8 +196,6 @@ 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)
@@ -220,8 +215,8 @@ def trafo_ref_to_physical(
blending_class: type[MultiAssignment]
if isinstance(geometry, TriangleElement):
blending_class = BlendingFTriangle
- elif isinstance(geometry, EmbeddedTriangle):
- blending_class = BlendingDFEmbeddedTriangle
+ if geometry.space_dimension == 3:
+ blending_class = BlendingDFEmbeddedTriangle
elif isinstance(geometry, TetrahedronElement):
blending_class = BlendingFTetrahedron
else:
@@ -246,15 +241,15 @@ def jac_affine_to_physical(
blending_class: type[MultiAssignment]
if isinstance(geometry, TriangleElement):
blending_class = BlendingDFTriangle
- elif isinstance(geometry, EmbeddedTriangle):
- blending_class = BlendingDFEmbeddedTriangle
+ if geometry.space_dimension == 3:
+ blending_class = BlendingDFEmbeddedTriangle
elif isinstance(geometry, TetrahedronElement):
blending_class = BlendingDFTetrahedron
else:
raise HOGException("Blending not implemented for the passed element geometry.")
t = trafo_ref_to_affine(geometry, symbolizer)
- jac = sp.zeros(geometry.dimensions, geometry.dimensions)
+ jac = sp.zeros(geometry.space_dimension, geometry.space_dimension)
rows, cols = jac.shape
for row in range(rows):
for col in range(cols):
diff --git a/hog/forms.py b/hog/forms.py
index 8f316d6f6dfa1b085a763c5c7a805713775d8082..8f7ed9f7be3e15c7ecd6cfb1ca1dbbf3ba89ab6d 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -24,6 +24,7 @@ from hog.element_geometry import ElementGeometry, TriangleElement, TetrahedronEl
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,
hessian_ref_to_affine,
@@ -96,7 +97,7 @@ Weak formulation
with TimedLogger("assembling diffusion matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -249,8 +250,8 @@ Weak formulation
with TimedLogger("assembling div-k-grad matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -360,8 +361,8 @@ Note: :math:`a(c) = 1/8 + u^2` is currently hard-coded and the form is intended
with TimedLogger("assembling non-linear diffusion matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -444,8 +445,8 @@ Note: :math:`a(c) = 1/8 + u^2` is currently hard-coded and the form is intended
):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -571,8 +572,8 @@ where
with TimedLogger("assembling epsilon matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -1011,8 +1012,8 @@ Weak formulation
f"assembling divergence {'transpose' if transpose else ''} matrix",
level=logging.DEBUG,
):
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -1122,8 +1123,8 @@ where
with TimedLogger("assembling full stokes matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -1353,8 +1354,8 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
with TimedLogger("assembling shear heating matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -1388,7 +1389,9 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
basis_eval=phi_eval_symbols,
)
else:
- raise HOGException("scalar_space_dependent_coefficient currently not supported in opgen.")
+ raise HOGException(
+ "scalar_space_dependent_coefficient currently not supported in opgen."
+ )
# mu = scalar_space_dependent_coefficient(
# "mu", geometry, symbolizer, blending=blending
# )
@@ -1433,7 +1436,6 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
dof_symbols=dof_symbols_uy,
)
-
# if geometry.dimensions > 2:
uz, dof_symbols_uz = fem_function_on_element(
velocity_function_space,
@@ -1488,21 +1490,17 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
for data in it:
phi = data.trial_shape
psi = data.test_shape
-
+
if blending != IdentityMap():
affine_factor = (
tabulation.register_factor(
"affine_factor_symbol",
sp.Matrix([phi * psi * jac_affine_det]),
)
- )[
- 0
- ]
+ )[0]
form = (
mu[0]
- * (
- double_contraction(tau, grad_u)[0]
- )
+ * (double_contraction(tau, grad_u)[0])
* jac_blending_det
* affine_factor
)
@@ -1510,19 +1508,10 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
shear_heating_det_symbol = (
tabulation.register_factor(
"shear_heating_det_symbol",
- (
- double_contraction(tau, grad_u)
- )
- * phi
- * psi
- * jac_affine_det,
+ (double_contraction(tau, grad_u)) * phi * psi * jac_affine_det,
)
- )[
- 0
- ]
- form = (
- mu[0] * shear_heating_det_symbol
- )
+ )[0]
+ form = mu[0] * shear_heating_det_symbol
mat[data.row, data.col] = form
@@ -1534,7 +1523,6 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
)
-
def divdiv(
trial: FunctionSpace,
test: FunctionSpace,
@@ -1680,8 +1668,8 @@ Weak formulation
with TimedLogger("assembling second derivative matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -1695,7 +1683,9 @@ Weak formulation
jac_blending_det = symbolizer.abs_det_jac_affine_to_blending()
if not isinstance(blending, IdentityMap):
# hessian_blending_map = blending.hessian(affine_coords)
- hessian_blending_map = symbolizer.hessian_blending_map(geometry.dimensions)
+ hessian_blending_map = symbolizer.hessian_blending_map(
+ geometry.dimensions
+ )
# jac_blending_det = abs(det(jac_blending))
# with TimedLogger("inverting blending Jacobian", level=logging.DEBUG):
diff --git a/hog/forms_boundary.py b/hog/forms_boundary.py
new file mode 100644
index 0000000000000000000000000000000000000000..edff32d2c0a481aa46148dd4d50a6cb2e39b6993
--- /dev/null
+++ b/hog/forms_boundary.py
@@ -0,0 +1,104 @@
+# 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 dataclasses import dataclass
+import logging
+import sympy as sp
+from typing import List, Optional, Tuple
+
+from hog.ast import Operations, count_operations
+from hog.element_geometry import ElementGeometry, TriangleElement, TetrahedronElement
+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,
+ hessian_ref_to_affine,
+ hessian_affine_to_blending,
+ create_empty_element_matrix,
+ element_matrix_iterator,
+ scalar_space_dependent_coefficient,
+ vector_space_dependent_coefficient,
+ fem_function_on_element,
+ fem_function_gradient_on_element,
+)
+from hog.function_space import FunctionSpace, EnrichedGalerkinFunctionSpace, N1E1Space
+from hog.math_helpers import dot, grad, inv, abs, det, double_contraction, e_vec
+from hog.quadrature import Quadrature, Tabulation
+from hog.symbolizer import Symbolizer
+from hog.logger import TimedLogger, get_logger
+from hog.blending import GeometryMap, ExternalMap, IdentityMap
+from hog.forms import Form
+
+
+def mass_boundary(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ volume_geometry: ElementGeometry,
+ boundary_geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+) -> Form:
+ docstring = f"""
+Mass operator.
+
+Geometry map: {blending}
+
+Weak formulation
+
+ u: trial function (space: {trial})
+ v: test function (space: {test})
+
+ ∫ uv ds
+"""
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble mass matrix."
+ )
+
+ 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."
+ )
+
+ with TimedLogger("assembling mass matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ jac_affine = symbolizer.jac_ref_to_affine(boundary_geometry)
+ jac_blending = symbolizer.jac_affine_to_blending(volume_geometry.dimensions)
+
+ fundamental_form_det = abs(
+ det(jac_affine.T * jac_blending.T * jac_blending * jac_affine)
+ )
+
+ mat = create_empty_element_matrix(trial, test, volume_geometry)
+ it = element_matrix_iterator(trial, test, volume_geometry)
+
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for data in it:
+ phi = data.trial_shape
+ psi = data.test_shape
+
+ form = phi * psi * fundamental_form_det**0.5
+
+ mat[data.row, data.col] = form
+
+ return Form(mat, tabulation, symmetric=True, docstring=docstring)
diff --git a/hog/forms_vectorial.py b/hog/forms_vectorial.py
index a5619145cbe36452bd936a8525be1cf4e5da61b4..237b3619db6ce8e13bba2d0084e8ca581953391c 100644
--- a/hog/forms_vectorial.py
+++ b/hog/forms_vectorial.py
@@ -159,8 +159,8 @@ def mass_n1e1(
raise HOGException("mass_n1e1 form is only implemented for N1E1.")
with TimedLogger("assembling mass matrix", level=logging.DEBUG):
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -287,7 +287,7 @@ def curl_curl(
raise HOGException("curl-curl form is only implemented for N1E1.")
with TimedLogger("assembling curl-curl matrix", level=logging.DEBUG):
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
@@ -353,8 +353,8 @@ Strong formulation
with TimedLogger("assembling curl_curl_plus_mass matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
- jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
- jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine = symbolizer.jac_ref_to_affine(geometry)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry)
jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
if isinstance(blending, ExternalMap):
diff --git a/hog/function_space.py b/hog/function_space.py
index 8ee718ed5551fb85631fc5a49066fe16fa974889..3e1b048daa84cb68385b61d24c05eaab31f46144 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -21,7 +21,6 @@ import sympy as sp
from hog.element_geometry import (
ElementGeometry,
TriangleElement,
- EmbeddedTriangle,
TetrahedronElement,
)
from hog.exception import HOGException
@@ -89,8 +88,6 @@ 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)
@@ -201,6 +198,7 @@ class LagrangianFunctionSpace(FunctionSpace):
elif (
isinstance(geometry, TriangleElement)
+ and geometry.dimensions == geometry.space_dimension
and self.family in ["Lagrange"]
and self._degree == 2
):
@@ -216,14 +214,16 @@ class LagrangianFunctionSpace(FunctionSpace):
]
elif (
- isinstance(geometry, EmbeddedTriangle)
+ isinstance(geometry, TriangleElement)
+ and geometry.dimensions == geometry.space_dimension - 1
and self.family in ["Lagrange"]
and self._degree == 0
):
basis_functions = [sp.sympify(1)]
elif (
- isinstance(geometry, EmbeddedTriangle)
+ isinstance(geometry, TriangleElement)
+ and geometry.dimensions == geometry.space_dimension - 1
and self.family in ["Lagrange"]
and self._degree == 1
):
@@ -234,7 +234,8 @@ class LagrangianFunctionSpace(FunctionSpace):
]
elif (
- isinstance(geometry, EmbeddedTriangle)
+ isinstance(geometry, TriangleElement)
+ and geometry.dimensions == geometry.space_dimension - 1
and self.family in ["Lagrange"]
and self._degree == 2
):
diff --git a/hog/manifold_forms.py b/hog/manifold_forms.py
index f53fcfed70272186bfd12aa2cc7afb418bc7bb9c..81f92c5b6941735829caab8151b67bbf980ad2c8 100644
--- a/hog/manifold_forms.py
+++ b/hog/manifold_forms.py
@@ -17,8 +17,7 @@
import logging
import sympy as sp
-from hog.element_geometry import ElementGeometry
-from hog.element_geometry import EmbeddedTriangle
+from hog.element_geometry import ElementGeometry, TriangleElement
from hog.exception import HOGException
from hog.fem_helpers import (
trafo_ref_to_affine,
@@ -65,8 +64,10 @@ Weak formulation
"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")
+ if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
with TimedLogger("assembling laplace beltrami matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
@@ -139,8 +140,10 @@ Weak formulation
"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.")
+ if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
with TimedLogger("assembling laplace beltrami matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
@@ -210,8 +213,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
projection = face_projection(geometry, symbolizer, blending=blending)
@@ -281,8 +288,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
normal = embedded_normal(geometry, symbolizer, blending)
@@ -351,8 +362,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
jac_affine = jac_ref_to_affine(geometry, symbolizer)
@@ -425,8 +440,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
projection_mat = face_projection(geometry, symbolizer, blending=blending)
@@ -509,8 +528,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
projection_mat = face_projection(geometry, symbolizer, blending=blending)
@@ -601,8 +624,12 @@ Weak formulation
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.")
+ if not (
+ isinstance(geometry, TriangleElement) and geometry.space_dimension == 3
+ ):
+ raise HOGException(
+ "Laplace Beltrami only works for triangles embedded in 3D space."
+ )
projection_mat = face_projection(geometry, symbolizer, blending=blending)
diff --git a/hog/manifold_helpers.py b/hog/manifold_helpers.py
index 4d2bf663ab74c27ceb6f1ed96d9b5914e4a3a855..125e36a957d1b08d2bdf17696da538354ca99084 100644
--- a/hog/manifold_helpers.py
+++ b/hog/manifold_helpers.py
@@ -19,7 +19,7 @@ import sympy as sp
from hog.exception import HOGException
-from hog.element_geometry import ElementGeometry, EmbeddedTriangle
+from hog.element_geometry import ElementGeometry, TriangleElement
from hog.symbolizer import Symbolizer
from hog.blending import GeometryMap, IdentityMap
from hog.fem_helpers import jac_affine_to_physical
@@ -33,13 +33,13 @@ def embedded_normal(
) -> sp.Matrix:
"""Returns an unoriented unit normal vector for an embedded triangle."""
- if not isinstance(geometry, EmbeddedTriangle):
+ if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
raise HOGException(
- "Embedded normal vectors are only defined for embedded triangles."
+ "Embedded normal vectors are only defined for triangles embedded in 3D space."
)
vert_points = symbolizer.affine_vertices_as_vectors(
- geometry.dimensions, geometry.num_vertices
+ geometry.space_dimension, geometry.num_vertices
)
span0 = vert_points[1] - vert_points[0]
span1 = vert_points[2] - vert_points[0]
@@ -67,9 +67,9 @@ def face_projection(
) -> sp.Matrix:
"""Returns a projection matrix for an embedded triangle."""
- if not isinstance(geometry, EmbeddedTriangle):
+ if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
raise HOGException(
- "Projection matrices are only defined for embedded triangles."
+ "Projection matrices are only defined for triangles embedded in 3D space."
)
normal = embedded_normal(geometry, symbolizer, blending=blending)
diff --git a/hog/operator_generation/function_space_impls.py b/hog/operator_generation/function_space_impls.py
index 130539f6040a59defd1ed11b668adcffa775ed5f..8421f9ff013e66a58bdd21d67cd943cd9449bc06 100644
--- a/hog/operator_generation/function_space_impls.py
+++ b/hog/operator_generation/function_space_impls.py
@@ -166,8 +166,15 @@ class FunctionSpaceImpl(ABC):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
- """Returns a list of local dof values on the current element."""
+ """
+ Returns a list of local dof values on the current element.
+
+ The element_vertex_ordering is a list that specifies the ordering of the reference vertices.
+ The ordering in which the DoFs are returned depends on this list. The "default" ordering is
+ [0, 1, ..., num_vertices - 1].
+ """
...
@abstractmethod
@@ -186,6 +193,7 @@ class FunctionSpaceImpl(ABC):
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
"""Returns HyTeG code that computes the basis/DoF transformation and a symbolic expression of the result.
@@ -214,6 +222,10 @@ class FunctionSpaceImpl(ABC):
assembling operators into matrices, these transformations must be
"baked into" the matrix because vectors are assembled locally and our
communication routine is not performed during the operator application.
+
+ The element_vertex_ordering is a list that specifies the ordering of the reference vertices.
+ The ordering in which the DoFs are returned depends on this list. The "default" ordering is
+ [0, 1, ..., num_vertices - 1].
"""
...
@@ -296,10 +308,14 @@ class P1FunctionSpaceImpl(FunctionSpaceImpl):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
vertex_dof_indices = micro_element_to_vertex_indices(
geometry, element_type, element_index
)
+
+ vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering]
+
vertex_array_indices = [
dof_idx.array_index(geometry, indexing_info)
for dof_idx in vertex_dof_indices
@@ -320,6 +336,7 @@ class P1FunctionSpaceImpl(FunctionSpaceImpl):
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
return (
CustomCodeNode("", [], []),
@@ -422,10 +439,14 @@ class P2FunctionSpaceImpl(FunctionSpaceImpl):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
vertex_dof_indices = micro_element_to_vertex_indices(
geometry, element_type, element_index
)
+
+ vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering]
+
edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices)
vrtx_array_idcs = [
@@ -454,6 +475,7 @@ class P2FunctionSpaceImpl(FunctionSpaceImpl):
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
return (
CustomCodeNode("", [], []),
@@ -557,6 +579,7 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
"""
Returns the element-local DoFs (field accesses) in a list (i.e., linearized).
@@ -606,6 +629,9 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
vertex_dof_indices = micro_element_to_vertex_indices(
geometry, element_type, element_index
)
+
+ vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering]
+
vertex_array_indices = [
dof_idx.array_index(geometry, indexing_info)
for dof_idx in vertex_dof_indices
@@ -628,6 +654,7 @@ class P1VectorFunctionSpaceImpl(FunctionSpaceImpl):
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
return (
CustomCodeNode("", [], []),
@@ -756,6 +783,7 @@ class P2VectorFunctionSpaceImpl(FunctionSpaceImpl):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
"""
Returns the element-local DoFs (field accesses) in a list (i.e., linearized).
@@ -772,6 +800,8 @@ class P2VectorFunctionSpaceImpl(FunctionSpaceImpl):
geometry, element_type, element_index
)
+ vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering]
+
edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices)
vertex_array_indices = [
@@ -801,14 +831,17 @@ class P2VectorFunctionSpaceImpl(FunctionSpaceImpl):
return f"P2VectorFunction< {self.type_descriptor.pystencils_type} >"
def includes(self) -> Set[str]:
- return {f"hyteg/p2functionspace/P2VectorFunction.hpp",
- f"hyteg/p2functionspace/P2Function.hpp"}
+ return {
+ f"hyteg/p2functionspace/P2VectorFunction.hpp",
+ f"hyteg/p2functionspace/P2Function.hpp",
+ }
def dof_transformation(
self,
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
return (
CustomCodeNode("", [], []),
@@ -872,10 +905,14 @@ class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
) -> List[Field.Access]:
vertex_dof_indices = micro_element_to_vertex_indices(
geometry, element_type, element_index
)
+
+ vertex_dof_indices = [vertex_dof_indices[i] for i in element_vertex_ordering]
+
edge_dof_indices = micro_vertex_to_edge_indices(geometry, vertex_dof_indices)
edge_array_indices = [
dof_idx.array_index(geometry, indexing_info) for dof_idx in edge_dof_indices
@@ -897,7 +934,14 @@ class N1E1FunctionSpaceImpl(FunctionSpaceImpl):
geometry: ElementGeometry,
element_index: Tuple[int, int, int],
element_type: Union[FaceType, CellType],
+ element_vertex_ordering: List[int],
) -> Tuple[CustomCodeNode, sp.MatrixBase]:
+
+ if element_vertex_ordering != [0, 1, 2, 3]:
+ raise HOGException(
+ "Element vertex re-ordering not supported for Nédélec elements (yet)."
+ )
+
Macro = {2: "Face", 3: "Cell"}[geometry.dimensions]
macro = {2: "face", 3: "cell"}[geometry.dimensions]
name = "basisTransformation"
diff --git a/hog/operator_generation/indexing.py b/hog/operator_generation/indexing.py
index dd445364aec0157b14a58ede2d6f61af509fd749..b29e566c5e6459c9f4170b1c11edc2ce4cf51d65 100644
--- a/hog/operator_generation/indexing.py
+++ b/hog/operator_generation/indexing.py
@@ -334,7 +334,7 @@ class IndexingInfo:
)
-def all_element_types(dimensions: int) -> Union[List[FaceType], List[CellType]]:
+def all_element_types(dimensions: int) -> List[Union[FaceType, CellType]]:
if dimensions == 2:
return [FaceType.GRAY, FaceType.BLUE]
if dimensions == 3:
diff --git a/hog/operator_generation/kernel_types.py b/hog/operator_generation/kernel_types.py
index 9050623c9d65c5a87c88ce944bf9e01fbc615d87..474a9b45d55a797fcd62f2469755a02e92dcee13 100644
--- a/hog/operator_generation/kernel_types.py
+++ b/hog/operator_generation/kernel_types.py
@@ -114,6 +114,7 @@ class KernelType(ABC):
element_type: Union[FaceType, CellType],
src_vecs_accesses: List[List[Field.Access]],
dst_vecs_accesses: List[List[Field.Access]],
+ element_vertex_ordering: List[int],
) -> List[ps.astnodes.Node]:
"""
Operations to be executed after the access resolution and common subexpression elimination on the symbols
@@ -158,6 +159,7 @@ class Apply(KernelType):
element_type: Union[FaceType, CellType],
src_vecs_accesses: List[List[Field.Access]],
dst_vecs_accesses: List[List[Field.Access]],
+ element_vertex_ordering: List[int],
) -> List[ps.astnodes.Node]:
tmp_symbols = sp.numbered_symbols(self.result_prefix)
@@ -197,6 +199,7 @@ class AssembleDiagonal(KernelType):
element_type: Union[FaceType, CellType],
src_vecs_accesses: List[List[Field.Access]],
dst_vecs_accesses: List[List[Field.Access]],
+ element_vertex_ordering: List[int],
) -> List[ps.astnodes.Node]:
tmp_symbols = sp.numbered_symbols(self.result_prefix)
@@ -243,6 +246,7 @@ class Assemble(KernelType):
element_type: Union[FaceType, CellType],
src_vecs_accesses: List[List[Field.Access]],
dst_vecs_accesses: List[List[Field.Access]],
+ element_vertex_ordering: List[int],
) -> List[ps.astnodes.Node]:
src, dst = dst_vecs_accesses
el_mat = sp.Matrix(
@@ -254,10 +258,10 @@ class Assemble(KernelType):
# apply basis/dof transformations
transform_src_code, transform_src_mat = self.src.dof_transformation(
- geometry, element_index, element_type
+ geometry, element_index, element_type, element_vertex_ordering
)
transform_dst_code, transform_dst_mat = self.dst.dof_transformation(
- geometry, element_index, element_type
+ geometry, element_index, element_type, element_vertex_ordering
)
transformed_el_mat = transform_dst_mat.T * el_mat * transform_src_mat
diff --git a/hog/operator_generation/loop_strategies.py b/hog/operator_generation/loop_strategies.py
index 9eb28a48cd4277af001d5de9052753a80c3ddcf8..0c370e42b7901a15b0fea556bf75365f608078cf 100644
--- a/hog/operator_generation/loop_strategies.py
+++ b/hog/operator_generation/loop_strategies.py
@@ -16,7 +16,7 @@
from abc import ABC, abstractmethod
import re
-from typing import Type, Union
+from typing import Type, Union, Dict
from pystencils import TypedSymbol
from pystencils.astnodes import (
@@ -24,6 +24,7 @@ from pystencils.astnodes import (
Conditional,
ResolvedFieldAccess,
SourceCodeComment,
+ LoopOverCoordinate,
)
from pystencils.sympyextensions import fast_subs
from pystencils.typing import FieldPointerSymbol
@@ -33,6 +34,7 @@ import sympy as sp
from hog.exception import HOGException
from hog.operator_generation.pystencils_extensions import (
loop_over_simplex,
+ loop_over_simplex_facet,
get_innermost_loop,
create_field_access,
create_micro_element_loops,
@@ -96,12 +98,14 @@ class LoopStrategy(ABC):
class CUBES(LoopStrategy):
- """For the "cubes" loop strategy, we want to update all micro-elements in the current "cube", i.e. loop over all
- # element types for the current micro-element index before incrementing.
- # This way we only have one loop, but need to insert conditionals that take care of those element types that
- # are not existing at the loop boundary.
- # The hope is that this strategy induces better cache locality, and that the conditionals can be automatically
- # transformed by the loop cutting features of pystencils."""
+ """
+ For the "cubes" loop strategy, we want to update all micro-elements in the current "cube", i.e. loop over all
+ element types for the current micro-element index before incrementing.
+ This way we only have one loop, but need to insert conditionals that take care of those element types that
+ are not existing at the loop boundary.
+ The hope is that this strategy induces better cache locality, and that the conditionals can be automatically
+ transformed by the loop cutting features of pystencils.
+ """
def __init__(self):
super(CUBES, self).__init__()
@@ -310,3 +314,120 @@ class FUSEDROWS(LoopStrategy):
def __str__(self):
return "FUSEDROWS"
+
+
+class BOUNDARY(LoopStrategy):
+ """
+ Special loop strategy that only loops over elements of a specified boundary.
+
+ Concretely, this means that it iterates over all elements with facets (edges in 2D, faces in 3D) that fully overlap
+ with the specified boundary (one of 3 macro-edges in 2D, one of 4 macro-faces in 3D).
+
+ To loop over multiple boundaries, just construct multiple loops.
+
+ The loop is constructed using conditionals, see the CUBES loop strategy.
+ """
+
+ def __init__(self, facet_id: int):
+ """
+ Constructs and initializes the BOUNDARY loop strategy.
+
+ :param facet_id: in [0, 2] for 2D, in [0, 3] for 3D
+ """
+ super(BOUNDARY, self).__init__()
+ self.facet_id = facet_id
+ self.element_loops: Dict[Union[FaceType, CellType], LoopOverCoordinate] = dict()
+
+ def create_loop(self, dim, element_index, micro_edges_per_macro_edge):
+
+ if dim == 2:
+ if self.facet_id not in [0, 1, 2]:
+ raise HOGException("Invalid facet ID for BOUNDARY loop strategy in 2D.")
+
+ self.element_loops = {
+ FaceType.GRAY: loop_over_simplex_facet(
+ dim, micro_edges_per_macro_edge, self.facet_id
+ ),
+ }
+
+ elif dim == 3:
+ if self.facet_id not in [0, 1, 2, 3]:
+ raise HOGException("Invalid facet ID for BOUNDARY loop strategy in 3D.")
+
+ second_cell_type = {
+ 0: CellType.BLUE_UP,
+ 1: CellType.GREEN_UP,
+ 2: CellType.BLUE_DOWN,
+ 3: CellType.GREEN_DOWN,
+ }
+
+ self.element_loops = {
+ CellType.WHITE_UP: loop_over_simplex_facet(
+ dim, micro_edges_per_macro_edge, self.facet_id
+ ),
+ second_cell_type[self.facet_id]: loop_over_simplex_facet(
+ dim, micro_edges_per_macro_edge - 1, self.facet_id
+ ),
+ }
+
+ return Block(
+ [
+ Block([SourceCodeComment(str(element_type)), loop])
+ for element_type, loop in self.element_loops.items()
+ ]
+ )
+
+ def add_body_to_loop(self, loop, body, element_type):
+ """Register a given loop body to innermost loop of the outer loop corresponding to element_type"""
+ if element_type not in self.element_loops:
+ return
+ innermost_loop = get_innermost_loop(self.element_loops[element_type])
+ body = Block(body)
+ innermost_loop[0].body = body
+ body.parent = innermost_loop[0]
+
+ def get_inner_bodies(self, loop_body):
+ return [
+ inner_loop.body
+ for inner_loop in get_innermost_loop(loop_body, return_all_inner=True)
+ ]
+
+ def add_preloop_for_loop(self, loops, preloop_stmts, element_type):
+ """add given list of statements directly in front of the loop corresponding to element_type."""
+
+ if element_type not in self.element_loops:
+ return loops
+
+ preloop_stmts_lhs_subs = {
+ stmt.lhs: get_element_replacement(stmt.lhs, element_type)
+ for stmt in preloop_stmts
+ }
+
+ if not isinstance(loops, Block):
+ loops = Block(loops)
+
+ blocks = loops.take_child_nodes()
+ new_blocks = []
+ for block in blocks:
+ idx_to_slice_at = -1
+ for idx, stmt in enumerate(block.args):
+ if stmt == self.element_loops[element_type]:
+ idx_to_slice_at = idx
+ break
+
+ if idx_to_slice_at == -1:
+ new_blocks.append(block)
+ else:
+ new_stmts = (
+ block.args[0:idx_to_slice_at]
+ + preloop_stmts
+ + block.args[idx_to_slice_at:]
+ )
+ new_block = Block(new_stmts)
+ new_block.fast_subs(preloop_stmts_lhs_subs)
+ new_blocks.append(new_block)
+
+ return new_blocks
+
+ def __str__(self):
+ return "BOUNDARY"
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index 2eb16fc268b28e80f33601fe340d2b5fc87eba76..d8fd1167c70af2feb3ab02e54db1104b5eeaf20c 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -14,10 +14,10 @@
# 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 dataclasses import dataclass
+from dataclasses import dataclass, field
from enum import auto, Enum
import logging
-from typing import Dict, List, Optional, Set, Tuple
+from typing import Dict, List, Optional, Set, Tuple, Union
import os
from textwrap import indent
@@ -76,15 +76,27 @@ from hog.blending import GeometryMap
import hog.code_generation
import hog.cse
from hog.dof_symbol import DoFSymbol
-from hog.element_geometry import ElementGeometry
+from hog.element_geometry import (
+ ElementGeometry,
+ TriangleElement,
+ TetrahedronElement,
+)
from hog.exception import HOGException
from hog.operator_generation.indexing import (
all_element_types,
element_vertex_coordinates,
IndexingInfo,
+ FaceType,
+ CellType,
+)
+from hog.operator_generation.kernel_types import KernelWrapperType
+from hog.operator_generation.kernel_types import Assemble, KernelType
+from hog.operator_generation.loop_strategies import (
+ LoopStrategy,
+ SAWTOOTH,
+ BOUNDARY,
+ CUBES,
)
-from hog.operator_generation.kernel_types import KernelWrapperType, KernelType, Assemble
-from hog.operator_generation.loop_strategies import LoopStrategy, SAWTOOTH, CUBES
from hog.operator_generation.optimizer import Optimizer, Opts
from hog.quadrature import QuadLoop, Quadrature
from hog.symbolizer import Symbolizer
@@ -94,9 +106,43 @@ from hog.operator_generation.types import HOGType
class MacroIntegrationDomain(Enum):
"""Enum type to specify where to integrate."""
- # Integration over the volume element
+ # Integration over the volume element.
VOLUME = "Volume"
+ # Integration over the boundary of the domain (for forms like ∫ ... d(∂Ω)).
+ #
+ # Note: Having one flag for all boundaries of one "type" is not very flexible.
+ # At some point there should be additional logic that uses HyTeG's BoundaryUIDs.
+ # To avoid ambiguity, the operator should (at least if there are boundary integrals) have a
+ # hyteg::BoundaryCondition constructor parameter, and for each boundary integral one hyteg::BoundaryUID
+ # parameter. Then the BC is tested by comparing those for each facet.
+ #
+ # Like so (application pseudocode in HyTeG):
+ #
+ # // generating an operator with one volume and one free-slip boundary integral
+ # // ∫ F dx + ∫ G ds
+ # // where ds corresponds to integrating over parts of the boundary that are marked with a specific
+ # // BoundaryUID
+ #
+ # // see HyTeG's documentation and/or BC tutorial
+ # BoundaryCondition someBC( ... );
+ #
+ # // setting up the BCs
+ # BoundaryUID freeslipBC = someBC.createFreeslipBC( ... );
+ #
+ # // the generated operator
+ # MyFancyFreeSlipOperator op( storage,
+ # ...,
+ # someBC, // this is what will be tested against - no ambiguity because src
+ # // and dst func might have different BCs!
+ # freeslipBC // this is the BCUID that will be used for the boundary integral
+ # // if there are more boundary integrals there are more UIDs in the
+ # // constructor (this BCUID is linked to 'ds' above)
+ # );
+ #
+
+ DOMAIN_BOUNDARY = "Domain boundary"
+
@dataclass
class IntegrationInfo:
@@ -106,6 +152,7 @@ class IntegrationInfo:
integration_domain: (
MacroIntegrationDomain # entity of geometry to integrate over, e.g. facet
)
+
quad: Quadrature # quadrature over integration domain of geometry, e.g. triangle
blending: GeometryMap
@@ -115,8 +162,11 @@ class IntegrationInfo:
loop_strategy: LoopStrategy
+ optimizations: Set[Opts]
+
name: str = "unknown_integral"
docstring: str = ""
+ boundary_uid_name: str = ""
def _str_(self):
return f"Integration Info: {self.name}, {self.geometry}, {self.integration_domain}, mat shape {self.mat.shape}, quad degree {self.quad.degree}, blending {self.blending}"
@@ -144,6 +194,66 @@ class CppClassFiles(Enum):
HEADER_IMPL_AND_VARIANTS = auto()
+def micro_vertex_permutation_for_facet(
+ volume_geometry: ElementGeometry,
+ element_type: Union[FaceType, CellType],
+ facet_id: int,
+) -> List[int]:
+ """
+ Provides a re-ordering of the micro-vertices such that the facet of the micro-element that coincides with the
+ macro-facet (which is given by the facet_id parameter) is spanned by the first three returned vertex positions.
+
+ The reordering can then for instance be executed by
+
+ ```python
+ element_vertex_order = shuffle_order_for_element_micro_vertices( ... )
+
+ el_vertex_coordinates = [
+ el_vertex_coordinates[i] for i in element_vertex_order
+ ]
+ ```
+
+ """
+
+ if volume_geometry == TriangleElement():
+
+ if element_type == FaceType.BLUE:
+ return [0, 1, 2]
+
+ shuffle_order_gray = {
+ 0: [0, 1, 2],
+ 1: [0, 2, 1],
+ 2: [1, 2, 0],
+ }
+
+ return shuffle_order_gray[facet_id]
+
+ elif volume_geometry == TetrahedronElement():
+
+ if element_type == CellType.WHITE_DOWN:
+ return [0, 1, 2, 3]
+
+ # All element types but WHITE_UP only overlap with a single macro-facet.
+ # It's a different element type for each facet. WHITE_DOWN is never at the boundary.
+ shuffle_order: Dict[Union[FaceType, CellType], Dict[int, List[int]]] = {
+ CellType.WHITE_UP: {
+ 0: [0, 1, 2, 3],
+ 1: [0, 1, 3, 2],
+ 2: [0, 2, 3, 1],
+ 3: [1, 2, 3, 0],
+ },
+ CellType.BLUE_UP: {0: [0, 1, 2, 3]},
+ CellType.GREEN_UP: {1: [0, 2, 3, 1]},
+ CellType.BLUE_DOWN: {2: [0, 1, 3, 2]},
+ CellType.GREEN_DOWN: {3: [1, 2, 3, 0]},
+ }
+
+ return shuffle_order[element_type][facet_id]
+
+ else:
+ raise HOGException("Not implemented.")
+
+
class HyTeGElementwiseOperator:
"""
This class handles the code generation of HyTeG-type 'elementwise' operators.
@@ -157,6 +267,8 @@ class HyTeGElementwiseOperator:
"hyteg/edgedofspace/EdgeDoFMacroCell.hpp",
"hyteg/primitivestorage/PrimitiveStorage.hpp",
"hyteg/LikwidWrapper.hpp",
+ "hyteg/boundary/BoundaryConditions.hpp",
+ "hyteg/types/types.hpp",
}
VAR_NAME_MICRO_EDGES_PER_MACRO_EDGE = "micro_edges_per_macro_edge"
@@ -171,7 +283,6 @@ class HyTeGElementwiseOperator:
name: str,
symbolizer: Symbolizer,
kernel_wrapper_types: List[KernelWrapperType],
- opts: Set[Opts],
type_descriptor: HOGType,
):
self.name = name
@@ -181,9 +292,6 @@ class HyTeGElementwiseOperator:
# Each IntegrationInfo object represents one integral of the weak formulation.
self.integration_infos: Dict[int, List[IntegrationInfo]] = {}
- self._optimizer = Optimizer(opts)
- self._optimizer_no_vec = Optimizer(opts - {Opts.VECTORIZE, Opts.VECTORIZE512})
-
# set the precision in which the operations are to be performed
self._type_descriptor = type_descriptor
@@ -192,31 +300,33 @@ class HyTeGElementwiseOperator:
# implementations for each kernel, generated at a later stage
self.operator_methods: List[OperatorMethod] = []
- def add_integral(
+ def _add_integral(
self,
name: str,
- dim: int,
- geometry: ElementGeometry,
+ volume_geometry: ElementGeometry,
integration_domain: MacroIntegrationDomain,
quad: Quadrature,
blending: GeometryMap,
form: Form,
loop_strategy: LoopStrategy,
+ boundary_uid_name: str,
+ optimizations: Set[Opts],
) -> None:
"""
- Use this method to add integrals to the operator.
+ Use this method to add integrals to the operator if you know what you are doing. There are helper methods for
+ adding integrals that are a little simpler to use.
:param name: some name for this integral (no spaces please)
- :param dim: the dimensionality of the domain, i.e. the volume - it may be that this does not match the geometry
- of the element, e.g., when facet integrals are required, note that only one routine per dim is
- created
- :param geometry: geometry that shall be integrated over
+ :param volume_geometry: the volume element (even for boundary integrals pass the element with dim == space_dim)
:param integration_domain: where to integrate - see MacroIntegrationDomain
- :param quad: the employed quadrature scheme - should match what has been used to integrate the weak form
+ :param quad: the employed quadrature scheme
:param blending: the same geometry map that has been passed to the form
:param form: the integrand
:param loop_strategy: loop pattern over the refined macro-volume - must somehow be compatible with the
integration domain
+ :param boundary_uid_name: string that defines the name of the boundary UID if this is a boundary integral
+ the parameter is ignored for volume integrals
+ :param optimizations: optimizations that shall be applied to this integral
"""
if "".join(name.split()) != name:
@@ -224,13 +334,15 @@ class HyTeGElementwiseOperator:
"Please give the integral an identifier without white space."
)
- if dim not in [2, 3]:
- raise HOGException("Only supporting 2D and 3D. Dim should be in [2, 3]")
+ if volume_geometry.space_dimension in self.integration_infos:
+ if name in [
+ ii.name
+ for ii in self.integration_infos[volume_geometry.space_dimension]
+ ]:
+ raise HOGException(f"Integral with name {name} already added!")
- if integration_domain != MacroIntegrationDomain.VOLUME:
- raise HOGException("Only volume integrals supported as of now.")
- if dim != geometry.dimensions:
- raise HOGException("Only volume integrals supported as of now.")
+ if volume_geometry.space_dimension not in [2, 3]:
+ raise HOGException("Only supporting 2D and 3D. Dim should be in [2, 3]")
if integration_domain == MacroIntegrationDomain.VOLUME and not (
isinstance(loop_strategy, SAWTOOTH) or isinstance(loop_strategy, CUBES)
@@ -241,14 +353,14 @@ class HyTeGElementwiseOperator:
quad_loop = None
mat = form.integrand
- if self._optimizer[Opts.TABULATE]:
+ if Opts.TABULATE in optimizations:
mat = form.tabulation.resolve_table_accesses(mat, self._type_descriptor)
with TimedLogger(f"constructing tables", level=logging.DEBUG):
tables = form.tabulation.construct_tables(quad, self._type_descriptor)
else:
mat = form.tabulation.inline_tables(mat)
- if self._optimizer[Opts.QUADLOOPS]:
+ if Opts.QUADLOOPS in optimizations:
quad_loop = QuadLoop(
self.symbolizer,
quad,
@@ -272,13 +384,13 @@ class HyTeGElementwiseOperator:
mat[row, col], self.symbolizer, blending
)
- if dim not in self.integration_infos:
- self.integration_infos[dim] = []
+ if volume_geometry.space_dimension not in self.integration_infos:
+ self.integration_infos[volume_geometry.space_dimension] = []
- self.integration_infos[dim].append(
+ self.integration_infos[volume_geometry.space_dimension].append(
IntegrationInfo(
name=name,
- geometry=geometry,
+ geometry=volume_geometry,
integration_domain=integration_domain,
quad=quad,
blending=blending,
@@ -287,17 +399,125 @@ class HyTeGElementwiseOperator:
mat=mat,
docstring=form.docstring,
loop_strategy=loop_strategy,
+ boundary_uid_name=boundary_uid_name,
+ optimizations=optimizations,
)
)
+ def add_volume_integral(
+ self,
+ name: str,
+ volume_geometry: ElementGeometry,
+ quad: Quadrature,
+ blending: GeometryMap,
+ form: Form,
+ loop_strategy: LoopStrategy,
+ optimizations: Union[None, Set[Opts]] = None,
+ ) -> None:
+ """
+ Adds a volume integral to the operator. Wrapper around _add_integral() for volume integrals.
+
+ :param name: some name for this integral (no spaces please)
+ :param volume_geometry: the volume element
+ :param quad: the employed quadrature scheme
+ :param blending: the same geometry map that has been passed to the form
+ :param form: the integrand
+ :param loop_strategy: loop pattern over the refined macro-volume - must somehow be compatible with the
+ integration domain
+ :param optimizations: optimization applied to this integral
+ """
+ if optimizations is None:
+ optimizations = set()
+
+ if volume_geometry.dimensions != quad.geometry.dimensions:
+ raise HOGException(
+ "The quadrature geometry does not match the volume geometry."
+ )
+
+ self._add_integral(
+ name,
+ volume_geometry,
+ MacroIntegrationDomain.VOLUME,
+ quad,
+ blending,
+ form,
+ loop_strategy,
+ "",
+ optimizations,
+ )
+
+ def add_boundary_integral(
+ self,
+ name: str,
+ volume_geometry: ElementGeometry,
+ quad: Quadrature,
+ blending: GeometryMap,
+ form: Form,
+ optimizations: Union[None, Set[Opts]] = None,
+ ) -> None:
+ """
+ Adds a boundary integral to the operator. Wrapper around _add_integral() for boundary integrals.
+
+ :param name: some name for this integral (no spaces please)
+ :param volume_geometry: the volume element (not the geometry of the boundary, also no embedded elements, just
+ the volume element geometry)
+ :param quad: the employed quadrature scheme - this must use the embedded geometry (e.g., for boundary integrals
+ in 2D, this should use a LineElement(space_dimension=2))
+ :param blending: the same map that has been passed to the form
+ :param form: the integrand
+ :param optimizations: optimization applied to this integral
+ """
+
+ if optimizations is None:
+ optimizations = set()
+
+ allowed_boundary_optimizations = {Opts.MOVECONSTANTS}
+ if optimizations - allowed_boundary_optimizations != set():
+ raise HOGException(
+ f"Only allowed (aka tested and working) optimizations for boundary integrals are "
+ f"{allowed_boundary_optimizations}."
+ )
+
+ if volume_geometry not in [TriangleElement(), TetrahedronElement()]:
+ raise HOGException(
+ "Boundary integrals only implemented for triangle and tetrahedral elements."
+ )
+
+ if volume_geometry.dimensions - 1 != quad.geometry.dimensions:
+ raise HOGException(
+ "The quadrature geometry does not match the boundary geometry."
+ )
+
+ # Since we will only integrate over the reference facet that lies on the x-line (2D) or xy-plane (3D) we need to
+ # set the last reference coordinate to zero since it will otherwise appear as a free, uninitialized variable.
+ #
+ # This has to be repeated later before the qudrature is applied in case we are working with symbols.
+ #
+ form.integrand = form.integrand.subs(
+ self.symbolizer.ref_coords_as_list(volume_geometry.dimensions)[-1], 0
+ )
+
+ for facet_id in range(volume_geometry.num_vertices):
+ self._add_integral(
+ name + f"_facet_id_{facet_id}",
+ volume_geometry,
+ MacroIntegrationDomain.DOMAIN_BOUNDARY,
+ quad,
+ blending,
+ form,
+ BOUNDARY(facet_id=facet_id),
+ name + "_boundary_uid",
+ optimizations,
+ )
+
def coefficients(self) -> List[FunctionSpaceImpl]:
"""Returns all coefficients sorted by name.
During generation coefficents are detected in the element matrix and
stored in the `coeffs` field. Being a Python dictionary, iterating over
it yields the coefficients in an arbitrary order. Whenever generating
- code for all coefficents it is a good idea to access them in a well
- defined order. Most importantly, the order of constructor arguments must
+ code for all coefficents it is a good idea to access them in a well-defined
+ order. Most importantly, the order of constructor arguments must
be deterministic.
"""
return sorted(self.coeffs.values(), key=lambda c: c.name)
@@ -311,17 +531,18 @@ class HyTeGElementwiseOperator:
"""
Invokes the code generation process, writing the full operator C++ code to file.
- :param dir_path: directory where to write the files - the file names are built automatically
- :param loop_strategy: iteration pattern
- :param header_only: if True, the entire class (incl. implementation) is written into a single file
- :clang_format_binary: path and/or name of binary for clang-format, defaults to None, which turns
- off formatting
+ :param dir_path: directory where to write the files - the file names are built automatically
+ :param class_files: determines whether header and or impl files are generated
+ :param clang_format_binary: path and/or name of binary for clang-format, defaults to None, which turns
+ off formatting
"""
- # Asking the optimizer if optimizations are valid.
- self._optimizer.check_opts_validity()
- # Generate each kernel type (apply, gemv, ...).
- self.generate_kernels()
+ with TimedLogger(
+ f"Generating kernels for operator {self.name}", level=logging.INFO
+ ):
+
+ # Generate each kernel type (apply, gemv, ...).
+ self.generate_kernels()
# Setting up the final C++ class.
operator_cpp_class = CppClass(
@@ -427,6 +648,30 @@ class HyTeGElementwiseOperator:
)
)
+ # Let's now check whether we need ctor arguments and member variables for boundary integrals.
+ boundary_condition_vars = []
+ for integration_infos in self.integration_infos.values():
+ if not all(
+ ii.integration_domain == MacroIntegrationDomain.VOLUME
+ for ii in integration_infos
+ ):
+ bc_var = CppVariable(name="boundaryCondition", type="BoundaryCondition")
+ if bc_var not in boundary_condition_vars:
+ boundary_condition_vars.append(bc_var)
+
+ for ii in integration_infos:
+ if ii.integration_domain == MacroIntegrationDomain.DOMAIN_BOUNDARY:
+ bcuid_var = CppVariable(
+ name=ii.boundary_uid_name, type="BoundaryUID"
+ )
+ if bcuid_var not in boundary_condition_vars:
+ boundary_condition_vars.append(bcuid_var)
+
+ boundary_condition_vars_members = [
+ CppVariable(name=bcv.name + "_", type=bcv.type)
+ for bcv in boundary_condition_vars
+ ]
+
# Finally we know what fields we need and can build the constructors, member variables, and includes.
# Constructors ...
@@ -450,9 +695,16 @@ class HyTeGElementwiseOperator:
is_reference=True,
)
for coeff in self.coefficients()
- ],
+ ]
+ + boundary_condition_vars,
initializer_list=["Operator( storage, minLevel, maxLevel )"]
- + [f"{coeff.name}( _{coeff.name} )" for coeff in self.coefficients()],
+ + [f"{coeff.name}( _{coeff.name} )" for coeff in self.coefficients()]
+ + [
+ f"{bcv[0].name}( {bcv[1].name} )"
+ for bcv in zip(
+ boundary_condition_vars_members, boundary_condition_vars
+ )
+ ],
)
)
@@ -468,7 +720,11 @@ class HyTeGElementwiseOperator:
)
)
- os.makedirs(dir_path, exist_ok=True) # Create path if it doesn't exist
+ for bcv in boundary_condition_vars_members:
+ operator_cpp_class.add(CppMemberVariable(bcv, visibility="private"))
+
+ # Create path if it doesn't exist
+ os.makedirs(dir_path, exist_ok=True)
output_path_header = os.path.join(dir_path, f"{self.name}.hpp")
output_path_impl = os.path.join(dir_path, f"{self.name}.cpp")
@@ -625,7 +881,10 @@ class HyTeGElementwiseOperator:
# This list is only filled if we want to vectorize.
loop_counter_custom_code_nodes = []
- if not (self._optimizer[Opts.VECTORIZE] or self._optimizer[Opts.VECTORIZE512]):
+ if (
+ Opts.VECTORIZE not in integration_info.optimizations
+ and Opts.VECTORIZE512 not in integration_info.optimizations
+ ):
# The Jacobians are loop-counter dependent, and we do not care about vectorization.
# So we just use the indices. pystencils will handle casting them to float.
el_matrix_element_index = element_index.copy()
@@ -689,7 +948,9 @@ class HyTeGElementwiseOperator:
]
# Let's fill the array.
- float_ctr_array_size = 8 if self._optimizer[Opts.VECTORIZE512] else 4
+ float_ctr_array_size = (
+ 8 if Opts.VECTORIZE512 in integration_info.optimizations else 4
+ )
custom_code = ""
custom_code += f"const int64_t phantom_ctr_0 = ctr_0;\n"
@@ -756,6 +1017,9 @@ class HyTeGElementwiseOperator:
rows, cols = mat.shape
+ optimizer = Optimizer(integration_info.optimizations)
+ optimizer.check_opts_validity()
+
kernel_config = CreateKernelConfig(
default_number_float=self._type_descriptor.pystencils_type,
data_type=self._type_descriptor.pystencils_type,
@@ -794,7 +1058,7 @@ class HyTeGElementwiseOperator:
# Common subexpression elimination.
with TimedLogger("cse on kernel operation", logging.DEBUG):
- cse_impl = self._optimizer.cse_impl()
+ cse_impl = optimizer.cse_impl()
kernel_op_assignments = hog.cse.cse(
kernel_op_assignments,
cse_impl,
@@ -810,7 +1074,7 @@ class HyTeGElementwiseOperator:
with TimedLogger("constructing quadrature loops"):
quad_loop = integration_info.quad_loop.construct_quad_loop(
- accessed_mat_entries, self._optimizer.cse_impl()
+ accessed_mat_entries, optimizer.cse_impl()
)
else:
quad_loop = []
@@ -861,7 +1125,52 @@ class HyTeGElementwiseOperator:
# element coordinates, jacobi matrix, tabulations, etc.
preloop_stmts = {}
- for element_type in all_element_types(geometry.dimensions):
+ # Deciding on which element types we want to iterate over.
+ # We skip certain element types for macro-volume boundary integrals.
+ element_types: List[Union[FaceType, CellType]] = all_element_types(
+ geometry.dimensions
+ )
+ if isinstance(integration_info.loop_strategy, BOUNDARY):
+ element_types = list(integration_info.loop_strategy.element_loops.keys())
+
+ for element_type in element_types:
+
+ # Re-ordering micro-element vertices for the handling of domain boundary integrals.
+ #
+ # Boundary integrals are handled by looping over all (volume-)elements that have a facet at one of the
+ # macro-element boundaries. There are three such boundaries in 2D and four in 3D. For each case, a separate
+ # kernel has to be generated. The logic that decides which of these kernels are called on which
+ # macro-volumes is handled by the HyTeG operator that is generated around these kernels.
+ #
+ # Instead of integrating over the micro-volume, we need to integrate over one of the micro-element facets.
+ # For that, the micro-element vertices are re-ordered such that the transformation from the affine element
+ # to the reference element results in the integration (boundary-)domain being mapped to the reference
+ # domain.
+ #
+ # E.g., in 2D, if the integration over the xy-boundary of the macro-volume (== macro-face) shall be
+ # generated this is signalled here by loop_strategy.facet_id == 2.
+ # In 2D all boundaries are only touched by the GRAY micro-elements (this is a little more complicated in 3D
+ # where at each macro-volume boundary, two types of elements overlap).
+ # Thus, we
+ # a) Shuffle the affine element vertices such that the xy-boundary of the GRAY elements is mapped to the
+ # (0, 1) line, which is the reference line for integration.
+ # b) Only iterate over the GRAY elements (handled later below).
+
+ # Default order.
+ element_vertex_order = list(range(geometry.num_vertices))
+
+ # Shuffling vertices if a boundary integral is requested.
+ if (
+ integration_info.integration_domain
+ == MacroIntegrationDomain.DOMAIN_BOUNDARY
+ and isinstance(integration_info.loop_strategy, BOUNDARY)
+ ):
+ element_vertex_order = micro_vertex_permutation_for_facet(
+ volume_geometry=geometry,
+ element_type=element_type,
+ facet_id=integration_info.loop_strategy.facet_id,
+ )
+
# Create array accesses to the source and destination vector(s) for
# the kernel.
src_vecs_accesses = [
@@ -870,6 +1179,7 @@ class HyTeGElementwiseOperator:
element_index, # type: ignore[arg-type] # list of sympy expressions also works
element_type,
indexing_info,
+ element_vertex_order,
)
for src_field in src_fields
]
@@ -879,6 +1189,7 @@ class HyTeGElementwiseOperator:
element_index, # type: ignore[arg-type] # list of sympy expressions also works
element_type,
indexing_info,
+ element_vertex_order,
)
for dst_field in dst_fields
]
@@ -896,6 +1207,7 @@ class HyTeGElementwiseOperator:
element_type,
src_vecs_accesses,
dst_vecs_accesses,
+ element_vertex_order,
)
# Load DoFs of coefficients. Those appear whenever a form is
@@ -906,6 +1218,16 @@ class HyTeGElementwiseOperator:
for ass in kernel_op_assignments + quad_loop
for a in ass.atoms(DoFSymbol)
}
+
+ if (
+ integration_info.integration_domain != MacroIntegrationDomain.VOLUME
+ and len(dof_symbols_set) > 0
+ ):
+ raise HOGException(
+ "Boundary integrals and FE coefficients cannot be combined at the moment."
+ "Dev note: the coefficients need to be sorted."
+ )
+
dof_symbols = sorted(dof_symbols_set, key=lambda ds: ds.name)
coeffs = dict(
(
@@ -925,7 +1247,11 @@ class HyTeGElementwiseOperator:
SympyAssignment(
sp.Symbol(dof_symbol.name),
coeffs[dof_symbol.function_id].local_dofs(
- geometry, element_index, element_type, indexing_info
+ geometry,
+ element_index,
+ element_type,
+ indexing_info,
+ element_vertex_order,
)[dof_symbol.dof_id],
)
)
@@ -953,6 +1279,13 @@ class HyTeGElementwiseOperator:
self.symbolizer,
)
+ # Re-ordering the element vertex coordinates
+ # (for all computations but affine transformation Jacobians - those are re-ordered later).
+ # See comment on boundary integrals above.
+ el_vertex_coordinates = [
+ el_vertex_coordinates[i] for i in element_vertex_order
+ ]
+
coords_assignments = [
SympyAssignment(
element_vertex_coordinates_symbols[vertex][component],
@@ -962,7 +1295,7 @@ class HyTeGElementwiseOperator:
for component in range(geometry.dimensions)
]
- if integration_info.blending.is_affine() or self._optimizer[Opts.QUADLOOPS]:
+ if integration_info.blending.is_affine() or optimizer[Opts.QUADLOOPS]:
blending_assignments = []
else:
blending_assignments = (
@@ -988,7 +1321,7 @@ class HyTeGElementwiseOperator:
)
with TimedLogger("cse on blending operation", logging.DEBUG):
- cse_impl = self._optimizer.cse_impl()
+ cse_impl = optimizer.cse_impl()
blending_assignments = hog.cse.cse(
blending_assignments,
cse_impl,
@@ -1006,7 +1339,25 @@ class HyTeGElementwiseOperator:
+ kernel_op_post_assignments
)
- if not self._optimizer[Opts.QUADLOOPS]:
+ if (
+ integration_info.integration_domain
+ == MacroIntegrationDomain.DOMAIN_BOUNDARY
+ and isinstance(integration_info.loop_strategy, BOUNDARY)
+ ):
+ with TimedLogger(
+ "boundary integrals: setting unused reference coordinate to 0",
+ logging.DEBUG,
+ ):
+ for node in body:
+ node.subs(
+ {
+ self.symbolizer.ref_coords_as_list(geometry.dimensions)[
+ -1
+ ]: 0
+ }
+ )
+
+ if not optimizer[Opts.QUADLOOPS]:
# Only now we replace the quadrature points and weights - if there are any.
# We also setup sympy assignments in body
with TimedLogger(
@@ -1060,6 +1411,12 @@ class HyTeGElementwiseOperator:
self.symbolizer,
)
+ # Re-ordering the element vertex coordinates for the Jacobians.
+ # See comment on boundary integrals above.
+ el_vertex_coordinates = [
+ el_vertex_coordinates[i] for i in element_vertex_order
+ ]
+
coords_assignments = [
SympyAssignment(
coord_symbols_for_jac_affine[vertex][component],
@@ -1072,7 +1429,7 @@ class HyTeGElementwiseOperator:
elem_dependent_stmts = (
hog.cse.cse(
coords_assignments + jacobi_assignments,
- self._optimizer.cse_impl(),
+ optimizer.cse_impl(),
"tmp_coords_jac",
return_type=SympyAssignment,
)
@@ -1141,7 +1498,7 @@ class HyTeGElementwiseOperator:
# generate AST of kernel loop
with TimedLogger(
- f"Generating kernel {integration_info.name} (wrapper: {kernel_wrapper_type.name} in {dim}D",
+ f"Generating kernel {integration_info.name} ({kernel_wrapper_type.name}, {dim}D)",
logging.INFO,
):
@@ -1174,11 +1531,12 @@ class HyTeGElementwiseOperator:
f"Optimizing kernel: {kernel_function.function_name} in {dim}D",
logging.INFO,
):
- optimizer = (
- self._optimizer
- if not isinstance(kernel_wrapper_type.kernel_type, Assemble)
- else self._optimizer_no_vec
- )
+ optimizer = Optimizer(integration_info.optimizations)
+ optimizer.check_opts_validity()
+
+ if isinstance(kernel_wrapper_type.kernel_type, Assemble):
+ optimizer = optimizer.copy_without_vectorization()
+
platform_dep_kernels.append(
optimizer.apply_to_kernel(
kernel_function, dim, integration_info.loop_strategy
@@ -1242,7 +1600,37 @@ class HyTeGElementwiseOperator:
# Kernel function call(s).
kernel_function_call_strings = []
- for kernel_function in kernel_functions:
+ for kernel_function, integration_info in zip(
+ kernel_functions, integration_infos
+ ):
+
+ pre_call_code = ""
+ post_call_code = ""
+
+ if (
+ integration_info.integration_domain
+ == MacroIntegrationDomain.DOMAIN_BOUNDARY
+ ):
+
+ if not isinstance(integration_info.loop_strategy, BOUNDARY):
+ raise HOGException(
+ "The loop strategy should be BOUNDARY for boundary integrals."
+ )
+
+ facet_type = "Edge" if dim == 2 else "Face"
+
+ neighbor_facet = (
+ f"getStorage()->get{facet_type}( "
+ f"{macro_type[dim]}.getLowerDimNeighbors()"
+ f"[{integration_info.loop_strategy.facet_id}] )"
+ )
+
+ pre_call_code = (
+ f"if ( boundaryCondition_.getBoundaryUIDFromMeshFlag( "
+ f"{neighbor_facet}->getMeshBoundaryFlag() ) == {integration_info.boundary_uid_name}_ ) {{"
+ )
+ post_call_code = "}"
+
kernel_parameters = kernel_function.get_parameters()
kernel_function_call_parameter_string = ",\n".join(
@@ -1251,8 +1639,10 @@ class HyTeGElementwiseOperator:
kernel_function_call_strings.append(
f"""
+ {pre_call_code}\n
{kernel_function.function_name}(\n
- {kernel_function_call_parameter_string});"""
+ {kernel_function_call_parameter_string});\n
+ {post_call_code}\n"""
)
kernel_wrapper_type.substitute(
diff --git a/hog/operator_generation/optimizer.py b/hog/operator_generation/optimizer.py
index 6d433c4b027553fd681954c9411f47fdc450c75d..e9922a24ea34fd3a628c5e37dd6f9a81c8345a93 100644
--- a/hog/operator_generation/optimizer.py
+++ b/hog/operator_generation/optimizer.py
@@ -18,7 +18,7 @@ from copy import deepcopy
import enum
import logging
import sympy as sp
-from typing import Dict, Iterable, List, Set
+from typing import Dict, Iterable, List, Set, Any
from pystencils import TypedSymbol
from pystencils.astnodes import (
@@ -108,6 +108,11 @@ class Optimizer:
def __getitem__(self, opt):
return opt in self._opts
+ def copy_without_vectorization(
+ self,
+ ) -> Any: # Optimizer and Self are not accepted :/
+ return Optimizer(self._opts - {Opts.VECTORIZE, Opts.VECTORIZE512})
+
def check_opts_validity(self) -> None:
"""Checks if the desired optimizations are valid."""
diff --git a/hog/operator_generation/pystencils_extensions.py b/hog/operator_generation/pystencils_extensions.py
index 9aa6dbe4a1d71fb89bd9a1d692cf49c39a90a8b2..8ee01d446f1d28cd2e8460d914d62bdb1f8b4512 100644
--- a/hog/operator_generation/pystencils_extensions.py
+++ b/hog/operator_generation/pystencils_extensions.py
@@ -18,7 +18,7 @@ from typing import Dict, List, Tuple, Union
import sympy as sp
from pystencils import FieldType, Field
-from pystencils.astnodes import Block, LoopOverCoordinate, Node
+from pystencils.astnodes import Block, LoopOverCoordinate, Node, get_next_parent_of_type
import pystencils as ps
from pystencils.astnodes import (
Block,
@@ -88,6 +88,60 @@ def loop_over_simplex(
return loops[dim - 1]
+def loop_over_simplex_facet(dim: int, width: int, facet_id: int) -> LoopOverCoordinate:
+ """
+ Loops over one boundary facet of a simplex (e.g., one edge in 2D, one face in 3D).
+
+ The facet is specified by an integer. It is required that
+
+ 0 ≤ facet_id ≤ dim
+
+ Let [x_0, x_1, ..., x_(dim-1)] be the coordinate of one element that is looped over.
+
+ For facet_id < dim, we have that
+
+ x_(dim - 1 - facet_id) = 0
+
+ and the remaining boundary is selected with facet_id == dim.
+
+ So in 2D for example, we get
+
+ facet_id = 0: (x_0, 0 )
+ facet_id = 1: (0, x_1)
+ facet_id = 2: the xy- (or "diagonal") boundary
+ """
+ loop = loop_over_simplex(dim, width)
+ innermost_loops = get_innermost_loop(loop)
+ if len(innermost_loops) != 1:
+ raise HOGException("There should be only one innermost loop.")
+ innermost_loop: LoopOverCoordinate = innermost_loops[0]
+
+ if facet_id not in range(0, dim + 1):
+ raise HOGException(f"Bad facet_id ({facet_id}) for dim {dim}.")
+
+ if facet_id == dim:
+ # For the "diagonal" loop we can just iterate as usual but skip all but the last element of the innermost loop.
+ # I hope.
+ innermost_loop.start = innermost_loop.stop - 1
+ return loop
+
+ # For the other facet_ids we need to find the corresponding loop and set that counter to 0.
+ # I am doing this here by just traversing the loops from the inside out, assuming that no loop cutting etc.
+ # occurred.
+ loop_with_counter_to_be_set_to_zero = innermost_loop
+ for d in range(dim - 1 - facet_id):
+ loop_with_counter_to_be_set_to_zero = get_next_parent_of_type(
+ loop_with_counter_to_be_set_to_zero, LoopOverCoordinate
+ )
+ if loop_with_counter_to_be_set_to_zero is None:
+ raise HOGException("There was no parent loop. This should not happen. I think.")
+
+ loop_with_counter_to_be_set_to_zero.start = 0
+ loop_with_counter_to_be_set_to_zero.stop = 1
+
+ return loop
+
+
def create_micro_element_loops(
dim: int, micro_edges_per_macro_edge: int
) -> Dict[Union[FaceType, CellType], LoopOverCoordinate]:
diff --git a/hog/quadrature/quadrature.py b/hog/quadrature/quadrature.py
index b51dd0bfce6014ac3e10c0ce05c89694295d5c20..dc0d39b4e7d49fe816ae05f331bac1a8554915fb 100644
--- a/hog/quadrature/quadrature.py
+++ b/hog/quadrature/quadrature.py
@@ -26,7 +26,6 @@ import sympy as sp
from hog.element_geometry import (
ElementGeometry,
TriangleElement,
- EmbeddedTriangle,
TetrahedronElement,
LineElement,
)
@@ -51,7 +50,6 @@ def select_quadrule(
"""Checks for availability of a specified quadrature rule and chooses a rule with minimal points
if only a degree is given."""
- # TODO for now, leave out line elements as we have no use for them without DG
logger = get_logger()
# quadrule given by name, just check if it exists and return it
@@ -96,24 +94,43 @@ def select_quadrule(
warnings.simplefilter("ignore")
all_schemes = []
- if isinstance(geometry, TriangleElement) or isinstance(
- geometry, EmbeddedTriangle
- ):
+ if isinstance(geometry, LineElement):
+ # Since line integrals can be approximated with arbitrary order, registering schemes is not
+ # meaningful. I doubt that we will need quadrature with order > 20, so we simply limit it here.
+ # If we require higher orders, simply bump this number.
+ line_quad_degree_limit = 20
+ # Also, for now let's restrict ourselves to Gauss-Legendre polynomials.
+ schemes = {
+ f"gauss_legendre_{d}": quadpy.c1.gauss_legendre(d)
+ for d in range(1, line_quad_degree_limit + 1)
+ }
+ elif isinstance(geometry, TriangleElement):
schemes = quadpy.t2.schemes
elif isinstance(geometry, TetrahedronElement):
schemes = quadpy.t3.schemes
for key, s in schemes.items():
try:
- scheme = s()
+ if isinstance(geometry, LineElement):
+ # For some reason the quadpy implementation is a little different for line segments.
+ scheme = s
+ else:
+ scheme = s()
+
except TypeError as e:
pass
all_schemes.append(scheme)
def select_degree(x: TnScheme) -> int:
- if x.degree >= degree:
- return x.points.shape[1]
+ if isinstance(geometry, LineElement):
+ if x.degree >= degree:
+ return x.degree
+ else:
+ return 10**10000 # just a large number
else:
- return 10**10000 # just a large number
+ if x.degree >= degree:
+ return x.points.shape[1]
+ else:
+ return 10**10000 # just a large number
return min(all_schemes, key=lambda x: select_degree(x))
@@ -121,9 +138,11 @@ def select_quadrule(
else:
raise HOGException(f"Unexpected {scheme_info}")
- logger.info(
- f"Integrating over {geometry} with rule: {scheme.name} (degree: {scheme.degree}, #points: {scheme.points.shape[1]})."
- )
+ if isinstance(geometry, LineElement):
+ num_points = len(scheme.points)
+ else:
+ num_points = scheme.points.shape[1]
+
return scheme
@@ -226,8 +245,6 @@ 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)
@@ -245,10 +262,12 @@ class Quadrature:
for i, (point, weight) in enumerate(zip(inline_points, inline_weights)):
spat_coord_subs = {}
for idx, symbol in enumerate(ref_symbols):
+ if not isinstance(point, sp.Matrix):
+ point = sp.Matrix([point])
spat_coord_subs[symbol] = point[idx]
if not blending.is_affine():
for symbol in symbolizer.quadpoint_dependent_free_symbols(
- self._geometry.dimensions
+ self._geometry.space_dimension
):
spat_coord_subs[symbol] = sp.Symbol(symbol.name + f"_q_{i}")
f_sub = fast_subs(f, spat_coord_subs)
@@ -291,9 +310,6 @@ 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.")
diff --git a/hog/symbolizer.py b/hog/symbolizer.py
index 721e82000c6707a5fe5ff7069c891fad25b427b2..d31c7524cbeb8feaefc7b47cf61b770d1dfc7ce7 100644
--- a/hog/symbolizer.py
+++ b/hog/symbolizer.py
@@ -18,6 +18,7 @@ from typing import List
import sympy as sp
import string
from hog.exception import HOGException
+from hog.element_geometry import ElementGeometry
class Symbolizer:
@@ -167,25 +168,27 @@ class Symbolizer:
def tmp_prefix(self) -> str:
return self._tmp_prefix
- def jac_ref_to_affine(self, dimensions: int) -> sp.Matrix:
+ def jac_ref_to_affine(self, geometry: ElementGeometry) -> sp.Matrix:
return sp.Matrix(
[
[
sp.Symbol(f"{self._symbol_jac_affine}_{i}_{j}")
- for j in range(dimensions)
+ for j in range(geometry.dimensions)
]
- for i in range(dimensions)
+ for i in range(geometry.space_dimension)
]
)
- def jac_ref_to_affine_inv(self, dimensions: int) -> sp.Matrix:
+ def jac_ref_to_affine_inv(self, geometry: ElementGeometry) -> sp.Matrix:
+ if geometry.dimensions != geometry.space_dimension:
+ raise HOGException("Cannot invert Jacobian for embedded elements.")
return sp.Matrix(
[
[
sp.Symbol(f"{self._symbol_jac_affine_inv}_{i}_{j}")
- for j in range(dimensions)
+ for j in range(geometry.dimensions)
]
- for i in range(dimensions)
+ for i in range(geometry.dimensions)
]
)
@@ -216,7 +219,7 @@ class Symbolizer:
def abs_det_jac_affine_to_blending(self, q_pt: str = "") -> sp.Symbol:
return sp.Symbol(f"{self._symbol_abs_det_jac_blending}{q_pt}")
-
+
def hessian_blending_map(self, dimensions: int, q_pt: str = "") -> List[sp.Matrix]:
return [
sp.Matrix(
@@ -227,7 +230,8 @@ class Symbolizer:
]
for i in range(dimensions)
]
- ) for k in range(dimensions)
+ )
+ for k in range(dimensions)
]
def blending_parameter_symbols(self, num_symbols: int) -> List[sp.Symbol]:
diff --git a/hog/sympy_extensions.py b/hog/sympy_extensions.py
index 3f92ee99ba1dcd4802d435817a208d850cd1043e..5d6467975da0349fafb80d8d84f46857f3e12924 100644
--- a/hog/sympy_extensions.py
+++ b/hog/sympy_extensions.py
@@ -22,7 +22,9 @@ T = TypeVar("T")
def fast_subs(
- expression: T, substitutions: Dict[sp.Expr, sp.Expr], skip: Optional[Callable[[sp.Expr], bool]] = None
+ expression: T,
+ substitutions: Dict[sp.Expr, sp.Expr],
+ skip: Optional[Callable[[sp.Expr], bool]] = None,
) -> T:
"""Similar to sympy subs function.
diff --git a/hog_tests/operator_generation/test_boundary_loop.py b/hog_tests/operator_generation/test_boundary_loop.py
new file mode 100644
index 0000000000000000000000000000000000000000..34a504e5eef0f153c382e00964fc1af26b812d74
--- /dev/null
+++ b/hog_tests/operator_generation/test_boundary_loop.py
@@ -0,0 +1,107 @@
+# 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
+
+from sympy.core.cache import clear_cache
+
+from hog.blending import AnnulusMap, IcosahedralShellMap
+from hog.element_geometry import LineElement, TriangleElement, TetrahedronElement
+from hog.function_space import LagrangianFunctionSpace
+from hog.operator_generation.operators import (
+ HyTeGElementwiseOperator,
+)
+from hog.symbolizer import Symbolizer
+from hog.quadrature import Quadrature, select_quadrule
+from hog.operator_generation.kernel_types import ApplyWrapper
+from hog.operator_generation.types import hyteg_type
+from hog.forms_boundary import mass_boundary
+from hog.logger import TimedLogger
+
+
+def test_boundary_loop():
+
+ # TimedLogger.set_log_level(logging.DEBUG)
+
+ dims = [2, 3]
+
+ clear_cache()
+
+ symbolizer = Symbolizer()
+
+ name = f"P2MassBoundary"
+
+ trial = LagrangianFunctionSpace(2, symbolizer)
+ test = LagrangianFunctionSpace(2, symbolizer)
+
+ type_descriptor = hyteg_type()
+
+ kernel_types = [
+ ApplyWrapper(
+ test,
+ trial,
+ type_descriptor=type_descriptor,
+ dims=dims,
+ )
+ ]
+
+ operator = HyTeGElementwiseOperator(
+ name,
+ symbolizer=symbolizer,
+ kernel_wrapper_types=kernel_types,
+ type_descriptor=type_descriptor,
+ )
+
+ for dim in dims:
+
+ if dim == 2:
+
+ volume_geometry = TriangleElement()
+ boundary_geometry = LineElement(space_dimension=2)
+ blending = AnnulusMap()
+
+ else:
+
+ volume_geometry = TetrahedronElement()
+ boundary_geometry = TriangleElement(space_dimension=3)
+ blending = IcosahedralShellMap()
+
+ quad = Quadrature(select_quadrule(5, boundary_geometry), boundary_geometry)
+
+ form = mass_boundary(
+ trial,
+ test,
+ volume_geometry,
+ boundary_geometry,
+ symbolizer,
+ blending=blending,
+ )
+
+ operator.add_boundary_integral(
+ name=f"boundary_mass",
+ volume_geometry=volume_geometry,
+ quad=quad,
+ blending=blending,
+ form=form,
+ )
+
+ operator.generate_class_code(
+ ".",
+ clang_format_binary="clang-format",
+ )
+
+
+if __name__ == "__main__":
+ test_boundary_loop()
diff --git a/hog_tests/operator_generation/test_opgen_smoke.py b/hog_tests/operator_generation/test_opgen_smoke.py
index b328de9311ecd64a9ad32fb435e2b0d599212314..747fd53d2b2498766db6c9cd54afc3e89bd6d17c 100644
--- a/hog_tests/operator_generation/test_opgen_smoke.py
+++ b/hog_tests/operator_generation/test_opgen_smoke.py
@@ -16,21 +16,18 @@
from sympy.core.cache import clear_cache
-from hog.blending import IdentityMap
from hog.operator_generation.loop_strategies import CUBES
from hog.operator_generation.optimizer import Opts
-from hog.element_geometry import TriangleElement
+from hog.element_geometry import LineElement, TriangleElement, TetrahedronElement
from hog.function_space import LagrangianFunctionSpace
-from hog.operator_generation.operators import (
- HyTeGElementwiseOperator,
- MacroIntegrationDomain,
-)
+from hog.operator_generation.operators import HyTeGElementwiseOperator
from hog.symbolizer import Symbolizer
from hog.quadrature import Quadrature, select_quadrule
-from hog.forms import div_k_grad, mass
-from hog.operator_generation.kernel_types import ApplyWrapper
+from hog.forms import div_k_grad
+from hog.forms_boundary import mass_boundary
+from hog.operator_generation.kernel_types import ApplyWrapper, AssembleWrapper
from hog.operator_generation.types import hyteg_type
-from hog.blending import AnnulusMap
+from hog.blending import AnnulusMap, IcosahedralShellMap
def test_opgen_smoke():
@@ -41,26 +38,24 @@ def test_opgen_smoke():
We are generating a matvec method here for
- ∫ k ∇u · ∇v dx + ∫ uv dx
+ ∫ k ∇u · ∇v dx + ∫ uv ds
with either integral being evaluated in their own kernel.
That may not be reasonable but tests some features.
"""
+
clear_cache()
symbolizer = Symbolizer()
- volume_geometry = TriangleElement()
- name = f"P2DivKGradBlendingPlusMass"
+ name = f"P2DivKGradBlendingPlusBoundaryMass"
+
+ dims = [2, 3]
trial = LagrangianFunctionSpace(2, symbolizer)
test = LagrangianFunctionSpace(2, symbolizer)
coeff = LagrangianFunctionSpace(2, symbolizer)
- quad = Quadrature(select_quadrule(2, volume_geometry), volume_geometry)
-
- divkgrad = div_k_grad(trial, test, volume_geometry, symbolizer, AnnulusMap(), coeff)
- m = mass(trial, test, volume_geometry, symbolizer, AnnulusMap())
type_descriptor = hyteg_type()
@@ -69,41 +64,72 @@ def test_opgen_smoke():
test,
trial,
type_descriptor=type_descriptor,
- dims=[2],
- )
+ dims=dims,
+ ),
+ AssembleWrapper(
+ test,
+ trial,
+ type_descriptor=type_descriptor,
+ dims=dims,
+ ),
]
- opts = {Opts.MOVECONSTANTS, Opts.VECTORIZE, Opts.TABULATE, Opts.QUADLOOPS}
-
operator = HyTeGElementwiseOperator(
name,
symbolizer=symbolizer,
kernel_wrapper_types=kernel_types,
- opts=opts,
type_descriptor=type_descriptor,
)
- operator.add_integral(
- name="div_k_grad",
- dim=volume_geometry.dimensions,
- geometry=volume_geometry,
- integration_domain=MacroIntegrationDomain.VOLUME,
- quad=quad,
- blending=AnnulusMap(),
- form=divkgrad,
- loop_strategy=CUBES(),
- )
+ opts_volume = {Opts.MOVECONSTANTS, Opts.VECTORIZE, Opts.TABULATE, Opts.QUADLOOPS}
+ opts_boundary = {Opts.MOVECONSTANTS}
- operator.add_integral(
- name="mass",
- dim=volume_geometry.dimensions,
- geometry=volume_geometry,
- integration_domain=MacroIntegrationDomain.VOLUME,
- quad=quad,
- blending=AnnulusMap(),
- form=m,
- loop_strategy=CUBES(),
- )
+ for d in dims:
+
+ if d == 2:
+ volume_geometry = TriangleElement()
+ boundary_geometry = LineElement(space_dimension=2)
+ blending_map = AnnulusMap()
+ else:
+ volume_geometry = TetrahedronElement()
+ boundary_geometry = TriangleElement(space_dimension=3)
+ blending_map = IcosahedralShellMap()
+
+ quad_volume = Quadrature(select_quadrule(2, volume_geometry), volume_geometry)
+ quad_boundary = Quadrature(
+ select_quadrule(5, boundary_geometry), boundary_geometry
+ )
+
+ divkgrad = div_k_grad(
+ trial, test, volume_geometry, symbolizer, blending_map, coeff
+ )
+ mass_b = mass_boundary(
+ trial,
+ test,
+ volume_geometry,
+ boundary_geometry,
+ symbolizer,
+ blending_map,
+ )
+
+ operator.add_volume_integral(
+ name="div_k_grad",
+ volume_geometry=volume_geometry,
+ quad=quad_volume,
+ blending=blending_map,
+ form=divkgrad,
+ loop_strategy=CUBES(),
+ optimizations=opts_volume,
+ )
+
+ operator.add_boundary_integral(
+ name="mass_boundary",
+ volume_geometry=volume_geometry,
+ quad=quad_boundary,
+ blending=blending_map,
+ form=mass_b,
+ optimizations=opts_boundary,
+ )
operator.generate_class_code(
".",