diff --git a/hog/forms.py b/hog/forms.py
index 5ed4b0eeea7aeecc044ad81ba042bb9ff4ad7e0f..db70846c44ea1bce2c55c183635131d3db62f992 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -329,6 +329,7 @@ def epsilon(
variable_viscosity: bool = True,
coefficient_function_space: Optional[FunctionSpace] = None,
rotation_wrapper: bool = False,
+ rigid_body_mode: bool = False,
) -> Form:
docstring = f"""
"Epsilon" operator.
@@ -353,6 +354,7 @@ where
from hog.recipes.integrands.volume.rotation import RotationType
from hog.recipes.integrands.volume.epsilon import integrand
+ from hog.recipes.integrands.volume.epsilon_rigid import integrand as integrand_rigid
from hog.recipes.integrands.volume.epsilon_affine import (
integrand as integrand_affine,
)
@@ -363,6 +365,9 @@ where
if blending == IdentityMap():
integr = integrand_affine
+ if rigid_body_mode:
+ integr = integrand_rigid
+
return process_integrand(
integr,
trial,
@@ -640,6 +645,8 @@ def full_stokes(
component_test: int = 0,
variable_viscosity: bool = True,
coefficient_function_space: Optional[FunctionSpace] = None,
+ rotation_wrapper: bool = False,
+ rigid_body_mode: bool = False,
) -> Form:
docstring = f"""
Implements the fully coupled viscous operator of the Stokes problem.
@@ -670,11 +677,17 @@ where
ε(w) := (1/2) (∇w + (∇w)ᵀ)
"""
-
+ from hog.recipes.integrands.volume.rotation import RotationType
from hog.recipes.integrands.volume.full_stokes import integrand
+ from hog.recipes.integrands.volume.full_stokes_rigid import integrand as integrand_rigid
+
+ if rigid_body_mode:
+ intgr = integrand_rigid
+ else:
+ intgr = integrand
return process_integrand(
- integrand,
+ intgr,
trial,
test,
geometry,
@@ -683,6 +696,7 @@ 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,
)
diff --git a/hog/recipes/integrands/volume/epsilon_rigid.py b/hog/recipes/integrands/volume/epsilon_rigid.py
new file mode 100644
index 0000000000000000000000000000000000000000..15cc9564d5ff4c911f5563253c1e23ce6147b984
--- /dev/null
+++ b/hog/recipes/integrands/volume/epsilon_rigid.py
@@ -0,0 +1,71 @@
+# 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 hog.recipes.common import *
+
+
+def integrand(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ u,
+ v,
+ grad_u,
+ grad_v,
+ k,
+ x,
+ scalars,
+ tabulate,
+ volume_geometry,
+ **_,
+):
+
+ grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
+ grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)
+
+ def symm_grad(w):
+ return 0.5 * (w + w.T)
+
+ symm_grad_u = symm_grad(grad_u_chain)
+ symm_grad_v = symm_grad(grad_v_chain)
+
+ c_rot_penalty = scalars("c_rot_penalty")
+
+ if volume_geometry.dimensions == 2:
+ rigid_z = sp.Matrix([[-x[1, 0]], [x[0, 0]]])
+ rigid_body_model_val = dot(rigid_z, u) * dot(rigid_z, v)
+ elif volume_geometry.dimensions == 3:
+ rigid_z = sp.Matrix([[-x[1, 0]], [x[0, 0]], [0]])
+ rigid_z = rigid_z / rigid_z.norm()
+
+ rigid_y = sp.Matrix([[0], [-x[2, 0]], [x[1, 0]]])
+ rigid_y = rigid_y / rigid_y.norm()
+
+ rigid_x = sp.Matrix([[x[2, 0]], [0], [-x[0, 0]]])
+ rigid_x = rigid_x / rigid_x.norm()
+
+ rigid_body_model_val = dot(rigid_z, u) * dot(rigid_z, v) + dot(rigid_y, u) * dot(rigid_y, v) + dot(rigid_x, u) * dot(rigid_x, v)
+
+ return (
+ (
+ double_contraction(2 * k["mu"] * symm_grad_u, symm_grad_v) +
+ c_rot_penalty * rigid_body_model_val
+ )
+ * tabulate(jac_a_abs_det)
+ * jac_b_abs_det
+ )
diff --git a/hog/recipes/integrands/volume/full_stokes_rigid.py b/hog/recipes/integrands/volume/full_stokes_rigid.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ba8cf05716467d93b410d05e2c848783a10130c
--- /dev/null
+++ b/hog/recipes/integrands/volume/full_stokes_rigid.py
@@ -0,0 +1,74 @@
+# 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 hog.recipes.common import *
+
+
+def integrand(
+ *,
+ jac_a_inv,
+ jac_a_abs_det,
+ jac_b_inv,
+ jac_b_abs_det,
+ u,
+ v,
+ grad_u,
+ grad_v,
+ k,
+ x,
+ scalars,
+ tabulate,
+ volume_geometry,
+ **_,
+):
+
+ grad_u_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)
+ grad_v_chain = jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)
+
+ def symm_grad(w):
+ return 0.5 * (w + w.T)
+
+ symm_grad_u = symm_grad(grad_u_chain)
+ symm_grad_v = symm_grad(grad_v_chain)
+
+ div_u = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_u)).trace()
+ div_v = (jac_b_inv.T * tabulate(jac_a_inv.T * grad_v)).trace()
+
+ c_rot_penalty = scalars("c_rot_penalty")
+
+ if volume_geometry.dimensions == 2:
+ rigid_z = sp.Matrix([[-x[1, 0]], [x[0, 0]]])
+ rigid_body_model_val = dot(rigid_z, u) * dot(rigid_z, v)
+ elif volume_geometry.dimensions == 3:
+ rigid_z = sp.Matrix([[-x[1, 0]], [x[0, 0]], [0]])
+ rigid_z = rigid_z / rigid_z.norm()
+
+ rigid_y = sp.Matrix([[0], [-x[2, 0]], [x[1, 0]]])
+ rigid_y = rigid_y / rigid_y.norm()
+
+ rigid_x = sp.Matrix([[x[2, 0]], [0], [-x[0, 0]]])
+ rigid_x = rigid_x / rigid_x.norm()
+
+ rigid_body_model_val = dot(rigid_z, u) * dot(rigid_z, v) + dot(rigid_y, u) * dot(rigid_y, v) + dot(rigid_x, u) * dot(rigid_x, v)
+
+ return k["mu"] * (
+ (
+ double_contraction(2 * symm_grad_u, symm_grad_v)
+ * tabulate(jac_a_abs_det)
+ * jac_b_abs_det
+ )
+ - (2.0 / 3.0) * div_u * div_v * tabulate(jac_a_abs_det) * jac_b_abs_det
+ ) + (c_rot_penalty * rigid_body_model_val) * tabulate(jac_a_abs_det) * jac_b_abs_det