From 16c30776a739e7f6ba422e662ddcf824d105b6dd Mon Sep 17 00:00:00 2001
From: Marcus Mohr <marcus.mohr@lmu.de>
Date: Fri, 10 Jan 2025 15:30:16 +0100
Subject: [PATCH] Adds generation of Mass and Diffusion operator for
P2PlusBubble in 2D
Note that this commit includes an IMHO rather not so nice hack for
the facedof_index function. The latter assumes that index is
Tuple[int, int, int]. However, the index it receives contains only
two entries.
---
generate_all_operators.py | 18 +++++++++++++++++-
hog/operator_generation/indexing.py | 9 ++++++++-
2 files changed, 25 insertions(+), 2 deletions(-)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index 6ce4a12..25ca992 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -29,7 +29,13 @@ 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,
+)
from hog.cse import CseImplementation
from hog.element_geometry import (
ElementGeometry,
@@ -46,6 +52,7 @@ from hog.forms import (
shear_heating,
epsilon,
full_stokes,
+ mass,
nonlinear_diffusion,
nonlinear_diffusion_newton_galerkin,
supg_diffusion,
@@ -64,6 +71,7 @@ from hog.function_space import (
LagrangianFunctionSpace,
N1E1Space,
TensorialVectorFunctionSpace,
+ P2PlusBubbleSpace,
TrialSpace,
TestSpace,
)
@@ -584,6 +592,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 +602,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 +657,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/operator_generation/indexing.py b/hog/operator_generation/indexing.py
index 896bcae..1680e39 100644
--- a/hog/operator_generation/indexing.py
+++ b/hog/operator_generation/indexing.py
@@ -265,7 +265,14 @@ def facedof_index(
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
+
+ # Ugly hack; why do we receive an index with only two entries,
+ # and why do we expect one with three?
+ if len(index) == 3:
+ x, y, _ = index
+ elif len(index) == 2:
+ x, y = 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)
--
GitLab