diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index d96ca5defea2c8af6a16ac3616d77dc6eb526d6f..b37446c88d8cea049feaafdbb64ba29ce217c5ed 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -34,6 +34,7 @@ from hog.element_geometry import (
from hog.function_space import (
LagrangianFunctionSpace,
N1E1Space,
+ P2PlusBubbleSpace,
TrialSpace,
TestSpace,
)
@@ -133,6 +134,8 @@ class FormInfo:
"""A compact representation of the function spaces."""
if self.trial_family == "N1E1":
return "n1e1"
+ elif self.trial_family == "P2 enhanced with Bubble":
+ return "p2_plus_bubble"
elif self.trial_degree == self.test_degree:
return f"p{self.trial_degree}"
else:
@@ -163,6 +166,8 @@ class FormInfo:
sub_dir = "p1"
if self.trial_family == "N1E1":
sub_dir = "n1e1"
+ elif self.trial_family == "P2 enhanced with Bubble":
+ return "p2_plus_bubble"
elif self.trial_degree == self.test_degree:
sub_dir = f"p{self.trial_degree}"
else:
@@ -242,6 +247,15 @@ form_infos = [
quad_schemes={2: 3, 3: 3},
blending=ExternalMap(),
),
+ FormInfo(
+ "diffusion",
+ trial_degree=2,
+ test_degree=2,
+ trial_family="P2 enhanced with Bubble",
+ test_family="P2 enhanced with Bubble",
+ quad_schemes={2: "exact"},
+ integrate_rows=[],
+ ),
FormInfo(
"mass",
trial_degree=1,
@@ -296,6 +310,15 @@ form_infos = [
blending=ExternalMap(),
integrate_rows=[],
),
+ FormInfo(
+ "mass",
+ trial_degree=2,
+ test_degree=2,
+ trial_family="P2 enhanced with Bubble",
+ test_family="P2 enhanced with Bubble",
+ quad_schemes={2: "exact"},
+ integrate_rows=[],
+ ),
FormInfo(
"curl_curl",
trial_degree=1,
@@ -1077,6 +1100,8 @@ def main():
trial: TrialSpace
if form_info.trial_family == "N1E1":
trial = TrialSpace(N1E1Space(symbolizer))
+ elif form_info.trial_family == "P2 enhanced with Bubble":
+ trial = TrialSpace(P2PlusBubbleSpace(symbolizer))
else:
trial = TrialSpace(
LagrangianFunctionSpace(form_info.trial_degree, symbolizer)
@@ -1085,6 +1110,8 @@ def main():
test: TestSpace
if form_info.test_family == "N1E1":
test = TestSpace(N1E1Space(symbolizer))
+ elif form_info.test_family == "P2 enhanced with Bubble":
+ test = TestSpace(P2PlusBubbleSpace(symbolizer))
else:
test = TestSpace(LagrangianFunctionSpace(form_info.test_degree, symbolizer))
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 6ce4a12dcab8ddddd60fb2e05107591d8415a84a..a2cae70e1450478d1f91aa65be0e0f3c5d513ca2 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -29,7 +29,14 @@ import sympy as sp
from sympy.core.cache import clear_cache
from tabulate import tabulate
-from hog.blending import GeometryMap, IdentityMap, AnnulusMap, IcosahedralShellMap, ParametricMap
+from hog.blending import (
+ GeometryMap,
+ IdentityMap,
+ AnnulusMap,
+ IcosahedralShellMap,
+ ParametricMap,
+ AffineMap2D,
+)
from hog.cse import CseImplementation
from hog.element_geometry import (
ElementGeometry,
@@ -46,6 +53,7 @@ from hog.forms import (
shear_heating,
epsilon,
full_stokes,
+ mass,
nonlinear_diffusion,
nonlinear_diffusion_newton_galerkin,
supg_diffusion,
@@ -64,6 +72,7 @@ from hog.function_space import (
LagrangianFunctionSpace,
N1E1Space,
TensorialVectorFunctionSpace,
+ P2PlusBubbleSpace,
TrialSpace,
TestSpace,
)
@@ -105,6 +114,7 @@ ALL_GEOMETRY_MAPS = [
IcosahedralShellMap(),
ParametricMap(1),
ParametricMap(2),
+ AffineMap2D(),
]
@@ -209,7 +219,7 @@ def parse_arguments():
"--blending",
default=str(IdentityMap()),
help=f"Enter a blending type for the operator. Note some geometry maps only work in 2D or 3D respectively. "
- f"Available geometry maps: {[str(m) for m in ALL_GEOMETRY_MAPS]}",
+ f"Available geometry maps: {sorted([str(m) for m in ALL_GEOMETRY_MAPS])}",
# TODO add choices list and type=string.lower
)
@@ -584,6 +594,7 @@ def all_operators(
P1Vector = TensorialVectorFunctionSpace(P1)
P2 = LagrangianFunctionSpace(2, symbolizer)
P2Vector = TensorialVectorFunctionSpace(P2)
+ P2PlusBubble = P2PlusBubbleSpace(symbolizer)
N1E1 = N1E1Space(symbolizer)
two_d = list(geometries & {TriangleElement()})
@@ -593,6 +604,7 @@ def all_operators(
# fmt: off
# TODO switch to manual specification of opts for now/developement, later use default factory
+
ops.append(OperatorInfo("N1E1", "CurlCurl", TrialSpace(N1E1), TestSpace(N1E1), form=curl_curl,
type_descriptor=type_descriptor, geometries=three_d, opts=opts, blending=blending))
ops.append(OperatorInfo("N1E1", "Mass", TrialSpace(N1E1), TestSpace(N1E1), form=mass_n1e1,
@@ -647,6 +659,12 @@ def all_operators(
form_boundary=mass_boundary, type_descriptor=type_descriptor, geometries=list(geometries),
opts=opts, blending=blending))
+ ops.append(OperatorInfo("P2PlusBubble", "Mass", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=mass,
+ type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
+
+ ops.append(OperatorInfo("P2PlusBubble", "Diffusion", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=diffusion,
+ type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
+
# fmt: on
p2vec_epsilon = partial(
diff --git a/hog/blending.py b/hog/blending.py
index 63478a11fe9b88e99f8bb5be6d44564b0aed03fd..e908feb1eed5558fdc2dd7b092775aee4e405ade 100644
--- a/hog/blending.py
+++ b/hog/blending.py
@@ -1,5 +1,5 @@
# HyTeG Operator Generator
-# Copyright (C) 2024 HyTeG Team
+# Copyright (C) 2024-2025 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
@@ -115,7 +115,6 @@ class ParametricMap(GeometryMap):
"""
def __init__(self, degree: int):
-
if degree not in [1, 2]:
raise HOGException(
"Only first and second order parametric maps are supported."
@@ -229,7 +228,7 @@ class AnnulusMap(GeometryMap):
return jac
def hessian(self, x: sp.Matrix) -> List[sp.Matrix]:
- """Evaluates the Jacobian of the geometry map at the passed point."""
+ """Evaluates the derivatives of the inverse Jacobian matrix of the geometry map at the passed point."""
if sp.shape(x) != (2, 1):
raise HOGException("Invalid input shape for AnnulusMap.")
@@ -299,6 +298,7 @@ class AnnulusMap(GeometryMap):
f"real_t {self.rayVertex[i]} = std::dynamic_pointer_cast< AnnulusMap >( face.getGeometryMap() )->rayVertex()[{i}];",
f"real_t {self.thrVertex[i]} = std::dynamic_pointer_cast< AnnulusMap >( face.getGeometryMap() )->thrVertex()[{i}];",
]
+
return "\n".join(code)
@@ -469,7 +469,7 @@ class IcosahedralShellMap(GeometryMap):
return jac
def hessian(self, x_: sp.Matrix) -> List[sp.Matrix]:
- """Evaluates the Jacobian of the geometry map at the passed point."""
+ """Evaluates the derivatives of the inverse Jacobian matrix of the geometry map at the passed point."""
if sp.shape(x_) != (3, 1):
raise HOGException("Invalid input shape for IcosahedralShellMap.")
@@ -595,3 +595,62 @@ class IcosahedralShellMap(GeometryMap):
]
return "\n".join(code)
+
+
+class AffineMap2D(GeometryMap):
+ """
+ This blending map uses a 2x2 matrix M and a shift vector v to
+ construct an affine mapping x -> M * x + v.
+
+ This is mostly intended for testing purposes, as we normally
+ would not want to have the overhead of blending in such a situation.
+ """
+
+ def __init__(self):
+ self.mat = sp.Matrix([["bMat_00", "bMat_01"], ["bMat_10", "bMat_11"]])
+ self.vec = sp.Matrix([["bVec_0"], ["bVec_1"]])
+
+ def supported_geometries(self) -> List[ElementGeometry]:
+ return [TriangleElement()]
+
+ def is_affine(self) -> bool:
+ # is_affine is supposed to indicate that the map does not have to be re-computed
+ # for each quadrature point - however, there seems to currently be a problem with
+ # the implementation of that optimization. Thus, we return False for the moment.
+ return False
+
+ def evaluate(self, x: sp.Matrix) -> sp.Matrix:
+ return self.mat * x + self.vec
+
+ def jacobian(self, x: sp.Matrix) -> sp.Matrix:
+ jac = self.mat.copy()
+ return jac
+
+ def hessian(self, x: sp.Matrix) -> List[sp.Matrix]:
+ hess = [sp.zeros(2, 2), sp.zeros(2, 2)]
+ return hess
+
+ def coupling_includes(self) -> List[str]:
+ return ["hyteg/geometry/AffineMap2D.hpp"]
+
+ def parameter_coupling_code(self) -> str:
+ code = []
+
+ code += [
+ f"WALBERLA_CHECK_NOT_NULLPTR( std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() ),"
+ f'"This operator requires the AffineMap2D to be registered as GeometryMap on every macro-face." )'
+ ]
+
+ code += [
+ f"real_t {self.mat[0,0]} = std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() )->getMatrix()(0,0);",
+ f"real_t {self.mat[0,1]} = std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() )->getMatrix()(0,1);",
+ f"real_t {self.mat[1,0]} = std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() )->getMatrix()(1,0);",
+ f"real_t {self.mat[1,1]} = std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() )->getMatrix()(1,1);",
+ ]
+
+ for i in range(2):
+ code += [
+ f"real_t {self.vec[i]} = std::dynamic_pointer_cast< AffineMap2D >( face.getGeometryMap() )->getVector()[{i}];",
+ ]
+
+ return "\n".join(code)
diff --git a/hog/forms.py b/hog/forms.py
index f67f0781d698489310f163fcd996231aecee66e0..ec630ce41999ca80a8396f1d0de0838f231fe435 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -36,6 +36,7 @@ from hog.fem_helpers import (
from hog.function_space import (
FunctionSpace,
N1E1Space,
+ P2PlusBubbleSpace,
TrialSpace,
TestSpace,
LagrangianFunctionSpace,
@@ -1097,7 +1098,7 @@ Weak formulation
∫ ((âˆ‡Ï / Ï) · u) v
"""
-
+
from hog.recipes.integrands.volume.frozen_velocity import integrand
return process_integrand(
@@ -1114,6 +1115,7 @@ Weak formulation
},
)
+
def zero_form(
trial: TrialSpace,
test: TestSpace,
diff --git a/hog/function_space.py b/hog/function_space.py
index 2d27469ee14c707c782a33f433f4cd7798eb9ca8..6384e94ad383416da02a28235b194cc4a8de28c4 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -80,14 +80,14 @@ class FunctionSpace(ABC):
) -> List[sp.MatrixBase]:
"""Returns a list containing the gradients of the shape functions on the element.
- Particularly, for each (vector- or scalar-valued) shape function N = (N_1, ..., N_n) returns the _transposed_
+ Particularly, for each (vector- or scalar-valued) shape function N = (N_1, ..., N_n) returns the *transposed*
Jacobian
⎡∂Nâ‚/∂x₠··· ∂Nâ‚™/∂xâ‚⎤
⎢ · · ⎥
grad(N) = J(N)ᵀ = ⎢ · · ⎥
⎢ · · ⎥
- ⎣∂Nâ‚/∂xâ‚– ··· ∂Nâ‚™/∂xâ‚– ⎦
+ ⎣∂Nâ‚/∂xâ‚– ··· ∂Nâ‚™/∂xₖ⎦
i.e., the returned gradient is a column vector for scalar shape functions.
@@ -326,7 +326,7 @@ class TensorialVectorFunctionSpace(FunctionSpace):
shape functions of the given scalar function space.
For instance, a tensorial function space of dim d = 3 that is based on a P1 (Lagrangian) scalar function space with
- n shape functions N_1, ..., N_n has the shape functions
+ n shape functions N_1, ..., N_n has the shape functions::
(N_1, 0, 0),
...
@@ -338,8 +338,9 @@ class TensorialVectorFunctionSpace(FunctionSpace):
...
(0, 0, N_n).
+
This class also enables specifying only a single component to get the shape functions, e.g. only for component 1
- (starting to count components at 0)
+ (starting to count components at 0)::
(0, N_1, 0),
...
@@ -592,3 +593,106 @@ class N1E1Space(FunctionSpace):
def __repr__(self):
return str(self)
+
+
+class P2PlusBubbleSpace(FunctionSpace):
+ """
+ Space of continuous piecewise quadratic functions enhanced with local bubbles vanishing on element boundaries.
+
+ In the 2D case the triplet :math:`(K,P,\Sigma)` describing the Finite Element is given by the following:
+
+ #. K represents the standard 2D-simplex, i.e. the unit triangle.
+ #. The local function space :math:`P` on :math:`K` is given by
+ :math:`P = \mathbb{P}_2 \oplus \ell_1 \ell_2 \ell_3 \mathbb{R}`.
+ Here the :math:`\ell_i` represent barycentric coordinates. Thus, any function :math:`v\in P` can uniquely be
+ composed into :math:`v = p + b`, where p is a quadratic polynomial and b a bubble function.
+ #. The set :math:`\Sigma` of linear functionals :math:`L_k : P \rightarrow \mathbb{R}, k = 1,\ldots,7` is given
+ by
+
+ .. math::
+ L_m( v ) = v(x_m), m = 1,\ldots 6 \\
+ L_7( v ) = b(x_7)
+
+ Here :math:`x_1, \ldots, x_3` represents the vertices of :math:`K`, :math:`x_4,\ldots,x_6` the midpoints
+ of its edges and :math:`x_7` its barycenter.
+ """
+
+ def __init__(self, symbolizer: Symbolizer):
+ self._symbolizer = symbolizer
+
+ @property
+ def is_vectorial(self) -> bool:
+ return False
+
+ @property
+ def is_continuous(self) -> bool:
+ return True
+
+ @property
+ def degree(self) -> int:
+ return 2
+
+ @property
+ def symbolizer(self) -> Symbolizer:
+ return self._symbolizer
+
+ @property
+ def family(self) -> str:
+ return "P2 enhanced with Bubble"
+
+ def shape(
+ self,
+ geometry: ElementGeometry,
+ domain: str = "reference",
+ dof_map: Optional[List[int]] = None,
+ ) -> List[sp.MatrixBase]:
+ """Returns a list containing the shape functions on the element.
+
+ :param dof_map: this list can be used to specify (remap) the DoF ordering of the element
+ """
+ if domain in ["ref", "reference"]:
+ symbols = self.symbolizer.ref_coords_as_list(geometry.dimensions)
+
+ basis_functions = []
+
+ if isinstance(geometry, TriangleElement):
+ L_1 = 1 - symbols[0] - symbols[1]
+ L_2 = symbols[0]
+ L_3 = symbols[1]
+ basis_functions = [
+ L_1 * (2 * L_1 - 1),
+ L_2 * (2 * L_2 - 1),
+ L_3 * (2 * L_3 - 1),
+ 4 * L_2 * L_3,
+ 4 * L_1 * L_3,
+ 4 * L_1 * L_2,
+ 27 * L_1 * L_2 * L_3,
+ ]
+
+ elif isinstance(geometry, TetrahedronElement):
+ raise HOGException(
+ f"P2PlusBubbleSpace is currently only implemented for triangle elements"
+ )
+ else:
+ raise HOGException(
+ "Basis functions not implemented for the specified element type and geometry."
+ )
+
+ if dof_map:
+ raise HOGException("DoF reordering not implemented.")
+
+ basis_functions = [sp.Matrix([b]) for b in basis_functions]
+ return basis_functions
+
+ raise HOGException(
+ "Basis functions not implemented for the specified element type and geometry."
+ )
+
+ def __eq__(self, other: Any) -> bool:
+ return type(self) == type(other)
+
+ def __str__(self):
+ return f"P2PlusBubble"
+
+ def __repr__(self):
+ return str(self)
diff --git a/hog/hyteg_form_template.py b/hog/hyteg_form_template.py
index f0d43087871aab1c97978286f5a0a4776f6eafbb..b1c7c34124dd97cb339213ac35a8b6fc11cc2b32 100644
--- a/hog/hyteg_form_template.py
+++ b/hog/hyteg_form_template.py
@@ -21,7 +21,7 @@ from hog.ast import Assignment, CodeBlock
from hog.exception import HOGException
from hog.quadrature import Quadrature
from hog.symbolizer import Symbolizer
-from hog.function_space import N1E1Space, TrialSpace, TestSpace
+from hog.function_space import N1E1Space, P2PlusBubbleSpace, TrialSpace, TestSpace
from hog.element_geometry import ElementGeometry
from hog.code_generation import code_block_from_element_matrix
from hog.multi_assignment import Member
@@ -292,6 +292,8 @@ class HyTeGFormClass:
if isinstance(self.trial, N1E1Space):
super_class = f"n1e1::N1E1Form"
+ elif isinstance(self.trial, P2PlusBubbleSpace):
+ super_class = f"P2PlusBubbleForm"
elif self.trial.degree == self.test.degree:
super_class = f"P{self.trial.degree}FormHyTeG"
else:
@@ -358,6 +360,8 @@ class HyTeGForm:
if isinstance(self.trial, N1E1Space):
super_class = f"N1E1Form"
+ elif isinstance(self.trial, P2PlusBubbleSpace):
+ super_class = f"P2PlusBubbleForm"
elif self.trial.degree == self.test.degree:
super_class = f"form_hyteg_base/P{self.trial.degree}FormHyTeG"
else:
diff --git a/hog/operator_generation/function_space_implementations/__init__.py b/hog/operator_generation/function_space_implementations/__init__.py
index b133d9c3e2a9fb91d0a0adb2c637c1069dc457c8..f18bade22e37884ba7fa6b49a68da21bb257b993 100644
--- a/hog/operator_generation/function_space_implementations/__init__.py
+++ b/hog/operator_generation/function_space_implementations/__init__.py
@@ -20,5 +20,6 @@ __all__ = [
"p0_space_impl",
"p1_space_impl",
"p2_space_impl",
+ "p2_plus_bubble_space_impl",
"n1e1_space_impl",
]
diff --git a/hog/operator_generation/function_space_implementations/function_space_impl_factory.py b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py
index c3762ee3c596c3279b34c93a2366528934cb81d5..4a6eb73c35d71080cefd1b06b0392e7d588abe36 100644
--- a/hog/operator_generation/function_space_implementations/function_space_impl_factory.py
+++ b/hog/operator_generation/function_space_implementations/function_space_impl_factory.py
@@ -19,6 +19,7 @@ from hog.exception import HOGException
from hog.function_space import (
FunctionSpace,
LagrangianFunctionSpace,
+ P2PlusBubbleSpace,
TensorialVectorFunctionSpace,
N1E1Space,
)
@@ -38,6 +39,9 @@ from hog.operator_generation.function_space_implementations.p2_space_impl import
P2FunctionSpaceImpl,
P2VectorFunctionSpaceImpl,
)
+from hog.operator_generation.function_space_implementations.p2_plus_bubble_space_impl import (
+ P2PlusBubbleFunctionSpaceImpl,
+)
from hog.operator_generation.function_space_implementations.n1e1_space_impl import (
N1E1FunctionSpaceImpl,
)
@@ -92,6 +96,8 @@ def create_impl(
)
elif isinstance(func_space, N1E1Space):
impl_class = N1E1FunctionSpaceImpl
+ elif isinstance(func_space, P2PlusBubbleSpace):
+ impl_class = P2PlusBubbleFunctionSpaceImpl
else:
raise HOGException("Unexpected function space")
diff --git a/hog/operator_generation/function_space_implementations/p2_plus_bubble_space_impl.py b/hog/operator_generation/function_space_implementations/p2_plus_bubble_space_impl.py
new file mode 100644
index 0000000000000000000000000000000000000000..9f7e33049f068cee8d63073a142485e09086f866
--- /dev/null
+++ b/hog/operator_generation/function_space_implementations/p2_plus_bubble_space_impl.py
@@ -0,0 +1,136 @@
+# 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 itertools
+import numpy as np
+import sympy as sp
+from typing import List, Tuple, Set, Union, Dict
+
+import pystencils as ps
+from pystencils import Field, FieldType
+from pystencils.backends.cbackend import CustomCodeNode
+
+from hog.element_geometry import ElementGeometry
+from hog.exception import HOGException
+from hog.function_space import (
+ FunctionSpace,
+ LagrangianFunctionSpace,
+ TensorialVectorFunctionSpace,
+ N1E1Space,
+)
+from hog.operator_generation.indexing import (
+ CellType,
+ FaceType,
+ VolumeDoFMemoryLayout,
+ micro_element_to_vertex_indices,
+ micro_vertex_to_edge_indices,
+ micro_element_to_volume_indices,
+ IndexingInfo,
+)
+from hog.operator_generation.types import HOGType
+from hog.operator_generation.function_space_implementations.function_space_impl_base import (
+ FunctionSpaceImpl,
+)
+
+from hog.operator_generation.function_space_implementations.p2_space_impl import (
+ P2FunctionSpaceImpl,
+)
+
+
+class P2PlusBubbleFunctionSpaceImpl(FunctionSpaceImpl):
+ def __init__(
+ self,
+ fe_space: FunctionSpace,
+ name: str,
+ type_descriptor: HOGType,
+ is_pointer: bool = False,
+ ):
+ super().__init__(fe_space, name, type_descriptor, is_pointer)
+ self.field = self._create_field(name)
+ self.p2Impl = P2FunctionSpaceImpl(fe_space, name, type_descriptor, is_pointer)
+
+ def pre_communication(self, dim: int) -> str:
+ return self.p2Impl.pre_communication(dim)
+
+ def zero_halos(self, dim: int) -> str:
+ return self.p2Impl.zero_halos(dim)
+
+ def post_communication(
+ self, dim: int, params: str, transform_basis: bool = True
+ ) -> str:
+ return self.p2Impl.post_communication(dim, params, transform_basis)
+
+ def pointer_retrieval(self, dim: int) -> str:
+ """C++ code for retrieving pointers to the numerical data stored in the macro primitives `face` (2d) or `cell` (3d)."""
+ Macro = {2: "Face", 3: "Cell"}[dim]
+ macro = {2: "face", 3: "cell"}[dim]
+
+ result_str = self.p2Impl.pointer_retrieval(dim)
+ result_str += f"{self.type_descriptor.pystencils_type}* _data_{self.name} = {self._deref()}.getVolumeDoFFunction().dofMemory( it.first, level );"
+
+ return result_str
+
+ def invert_elementwise(self, dim: int) -> str:
+ return f"{self._deref()}.invertElementwise( level );"
+
+ def local_dofs(
+ self,
+ geometry: ElementGeometry,
+ element_index: Tuple[int, int, int],
+ element_type: Union[FaceType, CellType],
+ indexing_info: IndexingInfo,
+ element_vertex_ordering: List[int],
+ ) -> List[Field.Access]:
+ complete_list_of_indices = self.p2Impl.local_dofs(
+ geometry,
+ element_index,
+ element_type,
+ indexing_info,
+ element_vertex_ordering,
+ )
+
+ volume_dof_indices = micro_element_to_volume_indices(
+ element_type, element_index, 1, VolumeDoFMemoryLayout.SoA
+ )
+
+ volume_array_indices = [
+ dof_idx.array_index(geometry, indexing_info)
+ for dof_idx in volume_dof_indices
+ ]
+
+ complete_list_of_indices += [
+ self.field.absolute_access((idx,), (0,)) for idx in volume_array_indices
+ ]
+
+ return complete_list_of_indices
+
+ def func_type_string(self) -> str:
+ return f"P2PlusBubbleFunction< {self.type_descriptor.pystencils_type} >"
+
+ def includes(self) -> Set[str]:
+ return {f"hyteg/experimental/P2PlusBubbleFunction.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("", [], []),
+ sp.Identity(self.fe_space.num_dofs(geometry)),
+ )
diff --git a/hog/operator_generation/indexing.py b/hog/operator_generation/indexing.py
index 896bcae8e6c4ff5ea2d04c69af5b447a2fd9b70d..218a8fef66eb9cf5dcec11e9494d2bcd92e130ba 100644
--- a/hog/operator_generation/indexing.py
+++ b/hog/operator_generation/indexing.py
@@ -27,7 +27,7 @@ from hog.symbolizer import Symbolizer
from pystencils.integer_functions import int_div
from pystencils import TypedSymbol
from enum import Enum
-from typing import Tuple, List, Union
+from typing import Tuple, List, Union, cast
import operator
import sympy as sp
@@ -259,13 +259,21 @@ def linear_macro_cell_index(width: int, x: int, y: int, z: int) -> int:
def facedof_index(
level: int,
- index: Tuple[int, int, int],
+ # Hack: see issue #46 for details
+ index: Union[Tuple[int, int], Tuple[int, int, int]],
faceType: Union[None, EdgeType, FaceType, CellType],
num_microfaces_per_face: sp.Symbol,
num_microedges_per_edge: sp.Symbol,
) -> int:
"""Indexes triangles/faces. Used to compute offsets in volume dof indexing in 2D and AoS layout."""
- x, y, _ = index
+
+ # second part of hack for #46
+ # also see https://github.com/python/mypy/issues/1178
+ if len(index) == 3:
+ x, y, _ = cast(Tuple[int, int, int], index)
+ else:
+ x, y = cast(Tuple[int, int], index)
+
# width = num_faces_per_row_by_type(level, faceType)
if faceType == FaceType.GRAY:
return linear_macro_face_index(num_microedges_per_edge, x, y)