diff --git a/hog/forms.py b/hog/forms.py
index 8ae046393a7a8ad2b0330b12c960bb6eccfa7c1b..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
+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
@@ -321,6 +327,7 @@ def epsilon(
component_test: int = 0,
variable_viscosity: bool = True,
coefficient_function_space: Optional[FunctionSpace] = None,
+ rotation_wrapper: bool = False,
) -> Form:
docstring = f"""
"Epsilon" operator.
@@ -339,9 +346,11 @@ 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,
@@ -363,6 +372,9 @@ where
is_symmetric=trial == test,
docstring=docstring,
fe_coefficients={"mu": coefficient_function_space},
+ rot_type=RotationType.PRE_AND_POST_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -540,6 +552,7 @@ def divergence(
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
component_index: int = 0,
+ rotation_wrapper: bool = False,
) -> Form:
docstring = f"""
Divergence.
@@ -556,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,
@@ -566,6 +580,9 @@ Weak formulation
blending=blending,
is_symmetric=False,
docstring=docstring,
+ rot_type=RotationType.POST_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -576,6 +593,7 @@ def gradient(
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
component_index: int = 0,
+ rotation_wrapper: bool = False,
) -> Form:
docstring = f"""
Gradient.
@@ -592,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,
@@ -602,6 +621,9 @@ def gradient(
blending=blending,
is_symmetric=False,
docstring=docstring,
+ rot_type=RotationType.PRE_MULTIPLY
+ if rotation_wrapper
+ else RotationType.NO_ROTATION,
)
@@ -767,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,
@@ -844,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,
@@ -888,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 81586f34781bb8bda245bfa4b72aad9d83650b05..07313705041256cc2c0ffb8f2f35d62ee448dbbc 100644
--- a/hog/integrand.py
+++ b/hog/integrand.py
@@ -73,6 +73,7 @@ from hog.fem_helpers import (
trafo_ref_to_physical,
)
from hog.math_helpers import inv, det
+from hog.recipes.integrands.volume.rotation import RotationType
@dataclass
@@ -230,6 +231,7 @@ def process_integrand(
boundary_geometry: ElementGeometry | None = None,
fe_coefficients: Dict[str, Union[FunctionSpace, None]] | None = None,
is_symmetric: bool = False,
+ rot_type: RotationType = RotationType.NO_ROTATION,
docstring: str = "",
) -> Form:
"""
@@ -291,6 +293,8 @@ def process_integrand(
finite-element function coefficients, they will be available to the callable as `k`
supply None as the FunctionSpace for a std::function-type coeff (only works for old forms)
:param is_symmetric: whether the bilinear form is symmetric - this is exploited by the generator
+ :param rot_type: whether the operator has to be wrapped with rotation matrix and the type of rotation that needs
+ to be applied, only applicable for Vectorial spaces
:param docstring: documentation of the integrand/bilinear form - will end up in the docstring of the generated code
"""
@@ -375,10 +379,10 @@ def process_integrand(
f"You cannot use the name {special_name_of_micromesh_coeff} for your FE coefficient."
f"It is reserved."
)
- fe_coefficients_modified[special_name_of_micromesh_coeff] = (
- TensorialVectorFunctionSpace(
- LagrangianFunctionSpace(blending.degree, symbolizer)
- )
+ fe_coefficients_modified[
+ special_name_of_micromesh_coeff
+ ] = TensorialVectorFunctionSpace(
+ LagrangianFunctionSpace(blending.degree, symbolizer)
)
s.k = dict()
@@ -491,6 +495,108 @@ def process_integrand(
free_symbols_sorted = sorted(list(free_symbols), key=lambda x: str(x))
+ 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 = rot_space.component_function_space
+
+ normals = (
+ ["nx_rotation", "ny_rotation"]
+ if volume_geometry.dimensions == 2
+ else ["nx_rotation", "ny_rotation", "nz_rotation"]
+ )
+
+ n_dof_symbols = []
+
+ for normal in normals:
+ if normal in fe_coefficients:
+ raise HOGException(
+ f"You cannot use the name {normal} for your FE coefficient."
+ f"It is reserved."
+ )
+ if normal_fspace is None:
+ raise HOGException("Invalid normal function space")
+ else:
+ _, nc_dof_symbols = fem_function_on_element(
+ normal_fspace,
+ volume_geometry,
+ symbolizer,
+ domain="reference",
+ function_id=normal,
+ )
+
+ n_dof_symbols.append(nc_dof_symbols)
+
+ rotmat = rotation_matrix(
+ rot_space.num_dofs(volume_geometry),
+ int(rot_space.num_dofs(volume_geometry) / volume_geometry.dimensions),
+ n_dof_symbols,
+ volume_geometry,
+ )
+
+ rot_doc_string = ""
+
+ if rot_type == RotationType.PRE_AND_POST_MULTIPLY:
+ mat = rotmat * mat * rotmat.T
+ rot_doc_string = "RKRᵀ uᵣ"
+ elif rot_type == RotationType.PRE_MULTIPLY:
+ mat = rotmat * mat
+ rot_doc_string = "RK uᵣ"
+ elif rot_type == RotationType.POST_MULTIPLY:
+ mat = mat * rotmat.T
+ rot_doc_string = "KRᵀ uᵣ"
+
+ docstring += f"""
+
+And the assembled FE matrix (K) is wrapped with a Rotation matrix (R) locally as below,
+
+ {rot_doc_string}
+
+where
+ R : Rotation matrix calculated with the normal vector (n̂) at the DoF
+ uᵣ: FE function but the components rotated at the boundaries according to the normal FE function passed
+
+ n̂ : normals (vectorial space: {normal_fspace})
+ * The passed normal vector must be normalized
+ * The radial component of the rotated vector will be pointing in the given normal direction
+ * If the normals are zero at a DoF, the rotation matrix is identity matrix
+ """
+
return Form(
mat,
tabulation,
diff --git a/hog/recipes/integrands/volume/rotation.py b/hog/recipes/integrands/volume/rotation.py
new file mode 100644
index 0000000000000000000000000000000000000000..f648988540bbfb867752a41513ebb4b58b98ece2
--- /dev/null
+++ b/hog/recipes/integrands/volume/rotation.py
@@ -0,0 +1,130 @@
+# 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 enum import Enum
+from typing import List
+from hog.recipes.common import *
+from hog.exception import HOGException
+from hog.element_geometry import ElementGeometry
+
+
+class RotationType(Enum):
+ """
+ This decides how the operator must be wrapped
+
+ The global stiffness matrix (K) assembled from an operator would result in,
+
+ NO_ROTATION - K
+ PRE_MULTIPLY - R K
+ POST_MULTIPLY - K R^T
+ PRE_AND_POST_MULTIPLY - R K R^T
+ """
+
+ NO_ROTATION = 0
+ PRE_MULTIPLY = 1
+ POST_MULTIPLY = 2
+ PRE_AND_POST_MULTIPLY = 3
+
+
+def rotation_matrix(
+ mat_size: int,
+ mat_comp_size: int,
+ n_dof_symbols: List[List[sp.Symbol]],
+ geometry: ElementGeometry,
+) -> sp.Matrix:
+ """
+ The Rotation matrix is used to wrap the operators, especially to apply freeslip boundary conditions.
+ The matrix itself is calculated from the normal vectors according to the method from [Engelmann 1982]
+
+ :param mat_size: Number of local DoF
+ :param mat_comp_size: Number of local DoF of each vector component
+ :param n_dof_symbols: Normal vector symbols packed as a list 2D/3D for each DoF
+ :param geometry: the geometry of the volume element
+ """
+
+ rotmat = sp.zeros(mat_size, mat_size)
+
+ for row in range(mat_comp_size):
+ for col in range(mat_comp_size):
+ if row != col:
+ continue
+
+ nx_ = n_dof_symbols[0][row]
+ ny_ = n_dof_symbols[1][row]
+
+ if geometry.dimensions == 2:
+ n_ = sp.Matrix([[nx_], [ny_]])
+ nnorm = sp.sqrt(nx_ * nx_ + ny_ * ny_)
+ rotmat_ = sp.Matrix([[-ny_, nx_], [nx_, ny_]])
+ elif geometry.dimensions == 3:
+ nz_ = n_dof_symbols[2][row]
+ n_ = sp.Matrix([[nx_], [ny_], [nz_]])
+ nnorm = sp.sqrt(nx_ * nx_ + ny_ * ny_ + nz_ * nz_)
+
+ ex = sp.Matrix([[1.0], [0.0], [0.0]])
+ ey = sp.Matrix([[0.0], [1.0], [0.0]])
+ ez = sp.Matrix([[0.0], [0.0], [1.0]])
+
+ ncross = [n_.cross(ex), n_.cross(ey), n_.cross(ez)]
+ ncrossnorm = [nc_.norm() for nc_ in ncross]
+
+ t1all = [
+ n_.cross(ex).normalized(),
+ n_.cross(ey).normalized(),
+ n_.cross(ez).normalized(),
+ ]
+
+ getrotmat = lambda t1: (t1.row_join(n_.cross(t1)).row_join(n_)).T
+
+ machine_eps = 1e-10
+ rotmat_ = sp.eye(geometry.dimensions)
+ for iRot in range(geometry.dimensions):
+ for jRot in range(geometry.dimensions):
+ rotmat_[iRot, jRot] = sp.Piecewise(
+ (
+ getrotmat(t1all[0])[iRot, jRot],
+ sp.And(
+ ncrossnorm[0] + machine_eps > ncrossnorm[1],
+ ncrossnorm[0] + machine_eps > ncrossnorm[2],
+ ),
+ ),
+ (
+ getrotmat(t1all[1])[iRot, jRot],
+ sp.And(
+ ncrossnorm[1] + machine_eps > ncrossnorm[0],
+ ncrossnorm[1] + machine_eps > ncrossnorm[2],
+ ),
+ ),
+ (getrotmat(t1all[2])[iRot, jRot], True),
+ )
+
+ else:
+ raise HOGException(
+ "We have not reached your dimensions, supreme being or little one"
+ )
+
+ rotmatid_ = sp.eye(geometry.dimensions)
+
+ ndofs = mat_comp_size
+
+ for iRot in range(geometry.dimensions):
+ for jRot in range(geometry.dimensions):
+ rotmat[row + iRot * ndofs, col + jRot * ndofs] = sp.Piecewise(
+ (rotmat_[iRot, jRot], nnorm > 0.5),
+ (rotmatid_[iRot, jRot], True),
+ )
+
+ return rotmat