Add weak form for penalising rigid body mode - angular momentum

generate_all_operators.py

@@ -47,6 +47,7 @@ from hog.element_geometry import (
 from hog.exception import HOGException
 from hog.forms import (
     diffusion,
+    rigid_body_mode_angular,
     divergence,
     gradient,
     div_k_grad,
@@ -665,6 +666,8 @@ def all_operators(
     ops.append(OperatorInfo("P2PlusBubble", "Diffusion", TrialSpace(P2PlusBubble), TestSpace(P2PlusBubble), form=diffusion,
                             type_descriptor=type_descriptor, geometries=two_d, opts=opts, blending=blending))
 
+    ops.append(OperatorInfo("P2Vector", "RigidModeAngular", TrialSpace(P2Vector), TestSpace(P2Vector), form=rigid_body_mode_angular,
+                            type_descriptor=type_descriptor, geometries=list(geometries), opts=opts, blending=blending))
     # fmt: on
 
     p2vec_epsilon = partial(

hog/forms.py

@@ -545,6 +545,49 @@ def linear_form(
 
     return mat
 
+def rigid_body_mode_angular(
+    trial: TrialSpace,
+    test: TestSpace,
+    geometry: ElementGeometry,
+    symbolizer: Symbolizer,
+    blending: GeometryMap,
+    rotation_wrapper: bool = False,
+) -> Form:
+    docstring = f"""
+Rigid Body Mode -- Angular Momentum.
+
+Geometry map: {blending}
+
+Weak formulation
+
+    u: trial function (vectorial space: {trial})
+    v: test function  (vectorial space:    {test})
+
+    ∫ ( ud · u ) ( ud · v )
+
+    ud encompasses the unit vectors pointing in the direction
+    of the three rigid body modes
+"""
+
+    from hog.recipes.integrands.volume.rigid_mode_angular import integrand
+    from hog.recipes.integrands.volume.rotation import RotationType
+
+    if isinstance(blending, IdentityMap):
+        raise HOGException("Rigid body mode only applicable for Annulus and IcosahedralShell map")
+
+    return process_integrand(
+        integrand,
+        trial,
+        test,
+        geometry,
+        symbolizer,
+        blending=blending,
+        is_symmetric=True,
+        docstring=docstring,
+        rot_type=RotationType.POST_MULTIPLY
+        if rotation_wrapper
+        else RotationType.NO_ROTATION,
+    )
 
 def divergence(
     trial: TrialSpace,
@@ -587,6 +630,46 @@ Weak formulation
         else RotationType.NO_ROTATION,
     )
 
+def divergence(
+    trial: TrialSpace,
+    test: TestSpace,
+    geometry: ElementGeometry,
+    symbolizer: Symbolizer,
+    blending: GeometryMap = IdentityMap(),
+    component_index: int = 0,
+    rotation_wrapper: bool = False,
+) -> Form:
+    docstring = f"""
+Divergence.
+
+Component:    {component_index}
+Geometry map: {blending}
+
+Weak formulation
+
+    u: trial function (vectorial space: {trial})
+    v: test function  (scalar space:    {test})
+
+    ∫ - ( ∇ · u ) v
+"""
+
+    from hog.recipes.integrands.volume.divergence import integrand
+    from hog.recipes.integrands.volume.rotation import RotationType
+
+    return process_integrand(
+        integrand,
+        trial,
+        test,
+        geometry,
+        symbolizer,
+        blending=blending,
+        component_index=component_index,
+        is_symmetric=False,
+        docstring=docstring,
+        rot_type=RotationType.POST_MULTIPLY
+        if rotation_wrapper
+        else RotationType.NO_ROTATION,
+    )
 
 def gradient(
     trial: TrialSpace,

hog/recipes/integrands/volume/rigid_mode_angular.py

+# 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_abs_det,
+    jac_b_abs_det,
+    u,
+    v,
+    x,
+    scalars,
+    tabulate,
+    volume_geometry,
+    **_,
+):
+    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 (c_rot_penalty * rigid_body_model_val) * tabulate(jac_a_abs_det) * jac_b_abs_det