diff --git a/hog/forms.py b/hog/forms.py
index 081b05854ac98850d1a6621498542ee5e0cfedaf..4a9439ad920f1796126b7f3085311ebe056adf36 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -33,7 +33,13 @@ from hog.fem_helpers import (
fem_function_on_element,
fem_function_gradient_on_element,
)
-from hog.function_space import FunctionSpace, N1E1Space, TrialSpace, TestSpace, LagrangianFunctionSpace
+from hog.function_space import (
+ FunctionSpace,
+ N1E1Space,
+ TrialSpace,
+ TestSpace,
+ LagrangianFunctionSpace,
+)
from hog.math_helpers import dot, inv, abs, det, double_contraction
from hog.quadrature import Quadrature, Tabulation
from hog.symbolizer import Symbolizer
@@ -322,7 +328,6 @@ def epsilon(
variable_viscosity: bool = True,
coefficient_function_space: Optional[FunctionSpace] = None,
rotation_wrapper: bool = False,
- normal_fspace: Optional[FunctionSpace] = None,
) -> Form:
docstring = f"""
"Epsilon" operator.
@@ -342,10 +347,10 @@ where
ε(w) := (1/2) (∇w + (∇w)ᵀ)
"""
-
if not variable_viscosity:
raise HOGException("Constant viscosity currently not supported.")
+ from hog.recipes.integrands.volume.rotation import RotationType
from hog.recipes.integrands.volume.epsilon import integrand
from hog.recipes.integrands.volume.epsilon_affine import (
integrand as integrand_affine,
@@ -357,8 +362,6 @@ where
if blending == IdentityMap():
integr = integrand_affine
-
-
return process_integrand(
integr,
trial,
@@ -369,7 +372,9 @@ where
is_symmetric=trial == test,
docstring=docstring,
fe_coefficients={"mu": coefficient_function_space},
- rotation_wrapper=rotation_wrapper,
+ rot_type=RotationType.PRE_AND_POST_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -547,6 +552,7 @@ def divergence(
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
component_index: int = 0,
+ rotation_wrapper: bool = False,
) -> Form:
docstring = f"""
Divergence.
@@ -563,6 +569,7 @@ Weak formulation
"""
from hog.recipes.integrands.volume.divergence import integrand
+ from hog.recipes.integrands.volume.rotation import RotationType
return process_integrand(
integrand,
@@ -573,6 +580,9 @@ Weak formulation
blending=blending,
is_symmetric=False,
docstring=docstring,
+ rot_type=RotationType.POST_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -583,6 +593,7 @@ def gradient(
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
component_index: int = 0,
+ rotation_wrapper: bool = False,
) -> Form:
docstring = f"""
Gradient.
@@ -599,6 +610,7 @@ def gradient(
"""
from hog.recipes.integrands.volume.gradient import integrand
+ from hog.recipes.integrands.volume.rotation import RotationType
return process_integrand(
integrand,
@@ -609,6 +621,9 @@ def gradient(
blending=blending,
is_symmetric=False,
docstring=docstring,
+ rot_type=RotationType.PRE_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -774,6 +789,7 @@ The resulting matrix must be multiplied with a vector of ones to be used as the
docstring=docstring,
)
+
def divdiv(
trial: TrialSpace,
test: TestSpace,
@@ -851,11 +867,10 @@ Weak formulation
blending=blending,
is_symmetric=trial == test,
docstring=docstring,
- fe_coefficients={
- "k": coefficient_function_space
- }
+ fe_coefficients={"k": coefficient_function_space},
)
+
def advection(
trial: TrialSpace,
test: TestSpace,
@@ -895,7 +910,9 @@ Weak formulation
"ux": velocity_function_space,
"uy": velocity_function_space,
"uz": velocity_function_space,
- } if constant_cp else {
+ }
+ if constant_cp
+ else {
"ux": velocity_function_space,
"uy": velocity_function_space,
"uz": velocity_function_space,
diff --git a/hog/integrand.py b/hog/integrand.py
index 7f82fed228c8de2b2e27c10d762b74f647ec3890..b43c3721b4d8675859838c62be39214636610053 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -88,8 +88,6 @@ class Form:
symmetric: bool
free_symbols: List[sp.Symbol] = field(default_factory=lambda: list())
docstring: str = ""
- rotmat: sp.MatrixBase = sp.Matrix([0])
- rot_type: RotationType = RotationType.NO_ROTATION
def integrate(self, quad: Quadrature, symbolizer: Symbolizer) -> sp.Matrix:
"""Integrates the form using the passed quadrature directly, i.e. without tabulations or loops."""
@@ -233,7 +231,7 @@ def process_integrand(
boundary_geometry: ElementGeometry | None = None,
fe_coefficients: Dict[str, Union[FunctionSpace, None]] | None = None,
is_symmetric: bool = False,
- rotation_wrapper: bool = False,
+ rot_type: RotationType = RotationType.NO_ROTATION,
docstring: str = "",
) -> Form:
"""
@@ -496,13 +494,45 @@ def process_integrand(
free_symbols_sorted = sorted(list(free_symbols), key=lambda x: str(x))
- if rotation_wrapper:
- if not trial.is_vectorial:
- raise HOGException("Rotation wrapper can only work with vectorial spaces")
+ if not rot_type == RotationType.NO_ROTATION:
+ if rot_type == RotationType.PRE_AND_POST_MULTIPLY:
+ if not trial == test:
+ HOGException(
+ "Trial and Test spaces must be the same for RotationType.PRE_AND_POST_MULTIPLY"
+ )
+
+ if not trial.is_vectorial:
+ raise HOGException(
+ "Rotation wrapper can only work with vectorial spaces"
+ )
+
+ rot_space = TensorialVectorFunctionSpace(
+ LagrangianFunctionSpace(trial.degree, symbolizer)
+ )
+
+ elif rot_type == RotationType.PRE_MULTIPLY:
+ if not test.is_vectorial:
+ raise HOGException(
+ "Rotation wrapper can only work with vectorial spaces"
+ )
+
+ rot_space = TensorialVectorFunctionSpace(
+ LagrangianFunctionSpace(test.degree, symbolizer)
+ )
+
+ elif rot_type == RotationType.POST_MULTIPLY:
+ if not trial.is_vectorial:
+ raise HOGException(
+ "Rotation wrapper can only work with vectorial spaces"
+ )
+
+ rot_space = TensorialVectorFunctionSpace(
+ LagrangianFunctionSpace(trial.degree, symbolizer)
+ )
from hog.recipes.integrands.volume.rotation import rotation_matrix
- normal_fspace = LagrangianFunctionSpace(trial.degree, symbolizer)
+ normal_fspace = rot_space.component_function_space
normals = (
["nx_rotation", "ny_rotation"]
@@ -532,12 +562,19 @@ def process_integrand(
n_dof_symbols.append(nc_dof_symbols)
rotmat = rotation_matrix(
- trial.num_dofs(volume_geometry),
- int(trial.num_dofs(volume_geometry) / volume_geometry.dimensions),
+ rot_space.num_dofs(volume_geometry),
+ int(rot_space.num_dofs(volume_geometry) / volume_geometry.dimensions),
n_dof_symbols,
volume_geometry,
)
+ if rot_type == RotationType.PRE_AND_POST_MULTIPLY:
+ mat = rotmat * mat * rotmat.T
+ elif rot_type == RotationType.PRE_MULTIPLY:
+ mat = rotmat * mat
+ elif rot_type == RotationType.POST_MULTIPLY:
+ mat = mat * rotmat.T
+
docstring += f"""
And the assembled FE matrix (K) is wrapped with a Rotation matrix (R) locally as below,
@@ -560,6 +597,4 @@ where
symmetric=is_symmetric,
free_symbols=free_symbols_sorted,
docstring=docstring,
- rotmat=rotmat if rotation_wrapper else sp.Matrix([0]),
- rot_type=RotationType.PRE_AND_POST_MULTIPLY if rotation_wrapper else RotationType.NO_ROTATION,
)
diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py
index 3bc0c03b617bc31185326040fe3fe1ea876038c0..78181e54045231302d4eded19c20da780c36704c 100644
--- a/hog/operator_generation/operators.py
+++ b/hog/operator_generation/operators.py
@@ -398,17 +398,6 @@ class HyTeGElementwiseOperator:
mat[row, col], self.symbolizer, blending
)
- if not form.rot_type == RotationType.NO_ROTATION:
- if form.rot_type == RotationType.PRE_AND_POST_MULTIPLY:
- mat = form.rotmat * mat * form.rotmat.T
- elif form.rot_type == RotationType.PRE_MULTIPLY:
- mat = form.rotmat * mat
- elif form.rot_type == RotationType.POST_MULTIPLY:
- mat = mat * form.rotmat.T
- else:
- raise HOGException("Not implemented")
-
-
if volume_geometry.space_dimension not in self.integration_infos:
self.integration_infos[volume_geometry.space_dimension] = []