diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index e2731976571de4544c7147174f04d6116ad3ba56..e21806a73402828daceb6cb44137fefc7160833b 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -24,7 +24,7 @@ import sympy as sp
from sympy.core.cache import clear_cache
import re
-from hog.blending import GeometryMap, IdentityMap, ExternalMap
+from hog.blending import GeometryMap, IdentityMap, ExternalMap, AnnulusMap
from hog.element_geometry import TriangleElement, TetrahedronElement, ElementGeometry
from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
from hog.forms import (
@@ -38,9 +38,10 @@ from hog.forms import (
divergence,
full_stokes,
divdiv,
+ supg_diffusion,
)
from hog.forms_vectorial import mass_n1e1, curl_curl
-from hog.quadrature import Quadrature
+from hog.quadrature import Quadrature, select_quadrule
from hog.exception import HOGException
from hog.logger import TimedLogger
from hog.symbolizer import Symbolizer
@@ -419,6 +420,20 @@ form_infos = [
description="Implements a linear form of type: (k(x), psi) where psi a test function and k = k(x) a vectorial, external function.",
integrate_rows=[],
),
+ FormInfo(
+ "supg_diffusion",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={2: 4, 3: 4},
+ blending=ExternalMap(),
+ ),
+ FormInfo(
+ "supg_diffusion",
+ trial_degree=2,
+ test_degree=2,
+ quad_schemes={2: 4, 3: 4},
+ blending=AnnulusMap(),
+ )
]
for d in [1, 2]:
@@ -658,6 +673,14 @@ def form_func(
transpose=False,
blending=blending,
).integrate(quad, symbolizer)
+ elif name.startswith("supg_d"):
+ return supg_diffusion(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ ).integrate(quad, symbolizer)
else:
raise HOGException(f"Cannot call form function with name {name}.")
@@ -821,7 +844,8 @@ def main():
f"- Generating code for class {form_info.class_name(row,col)}, {geometry.dimensions}D"
):
quad = Quadrature(
- form_info.quad_schemes[geometry.dimensions],
+ select_quadrule(form_info.quad_schemes[geometry.dimensions], geometry),
+ # form_info.quad_schemes[geometry.dimensions],
geometry,
inline_values=form_info.inline_quad,
)
diff --git a/generate_all_operators.py b/generate_all_operators.py
index a9c58a69a59bba041aa3b4c95e04e84220a26cab..8806bc6ff8746aa8e4c35b944d352248cc9e71e8 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -40,6 +40,7 @@ from hog.forms import (
full_stokes,
nonlinear_diffusion,
nonlinear_diffusion_newton_galerkin,
+ supg_diffusion,
)
from hog.forms_vectorial import curl_curl, curl_curl_plus_mass, mass_n1e1
from hog.function_space import FunctionSpace, LagrangianFunctionSpace, N1E1Space
@@ -552,6 +553,10 @@ def all_operators(
coefficient_function_space=P1, onlyNewtonGalerkinPartOfForm=False),
type_descriptor=type_descriptor, opts=opts, blending=blending))
+ ops.append(OperatorInfo(mapping="P2", name="SUPGDiffusion", trial_space=P2, test_space=P2,
+ form=partial(supg_diffusion, velocity_function_space=P2, diffusivityXdelta_function_space=P2),
+ type_descriptor=type_descriptor, opts=opts, blending=blending))
+
# fmt: on
for c in [0, 1, 2]:
diff --git a/hog/forms.py b/hog/forms.py
index 4c7c6b217e935c2f0fb95554f1c19df80cc35f10..568e0ef3f2e8d583764e68eb72986d1d4bef90ec 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -20,12 +20,15 @@ import sympy as sp
from typing import List, Optional, Tuple
from hog.ast import Operations, count_operations
-from hog.element_geometry import ElementGeometry
+from hog.element_geometry import ElementGeometry, TriangleElement, TetrahedronElement
from hog.exception import HOGException
from hog.fem_helpers import (
trafo_ref_to_affine,
jac_ref_to_affine,
jac_affine_to_physical,
+ hessian_ref_to_affine,
+ hessian_affine_to_blending,
+ jacinvjac_affine_to_physical,
create_empty_element_matrix,
element_matrix_iterator,
scalar_space_dependent_coefficient,
@@ -1375,6 +1378,180 @@ Weak formulation
docstring=docstring,
)
+def supg_diffusion(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ velocity_function_space: FunctionSpace = None,
+ diffusivityXdelta_function_space: FunctionSpace = None
+) -> Form:
+ docstring = f"""
+Second derivative operator for testing.
+
+Geometry map: {blending}
+
+Weak formulation
+
+ T: trial function (space: {trial})
+ s: test function (space: {test})
+ u: velocity function (space: {velocity_function_space})
+ k𝛿: FE function representing k·𝛿 (space: {diffusivityXdelta_function_space})
+
+ For OpGen,
+
+ ∫ k(ΔT) · 𝛿(u · ∇s)
+
+ -------------------
+
+ For ExternalMap (only for testing),
+
+ ∫ (ΔT) s
+"""
+
+ if trial != test:
+ raise HOGException(
+ "Trial space must be equal to test space to assemble second derivative matrix."
+ )
+
+ with TimedLogger("assembling second derivative matrix", level=logging.DEBUG):
+ tabulation = Tabulation(symbolizer)
+
+ jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
+ jac_affine_inv = symbolizer.jac_ref_to_affine_inv(geometry.dimensions)
+ jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
+
+ if isinstance(blending, ExternalMap):
+ TimedLogger("ExternalMap is not tested well").log()
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ jacinvjac_blending = jacinvjac_affine_to_physical(geometry, symbolizer)
+ else:
+ affine_coords = trafo_ref_to_affine(geometry, symbolizer)
+ jac_blending = blending.jacobian(affine_coords)
+ if not isinstance(blending, IdentityMap):
+ hessian_blending_map = blending.hessian(affine_coords)
+
+ jac_blending_det = abs(det(jac_blending))
+ with TimedLogger("inverting blending Jacobian", level=logging.DEBUG):
+ jac_blending_inv = inv(jac_blending)
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+ it = element_matrix_iterator(trial, test, geometry)
+
+ if velocity_function_space != None and diffusivityXdelta_function_space != None:
+ u_eval_symbols = tabulation.register_phi_evals(
+ velocity_function_space.shape(geometry)
+ )
+
+ ux, _ = fem_function_on_element(
+ velocity_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="ux",
+ basis_eval=u_eval_symbols,
+ )
+
+ uy, _ = fem_function_on_element(
+ velocity_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uy",
+ basis_eval=u_eval_symbols,
+ )
+
+ if isinstance(geometry, TetrahedronElement):
+ uz, _ = fem_function_on_element(
+ velocity_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uz",
+ basis_eval=u_eval_symbols,
+ )
+ u = sp.Matrix([[ux], [uy], [uz]])
+ else:
+ u = sp.Matrix([[ux], [uy]])
+
+ kdelta_eval_symbols = tabulation.register_phi_evals(
+ diffusivityXdelta_function_space.shape(geometry)
+ )
+
+ kdelta, _ = fem_function_on_element(
+ diffusivityXdelta_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="kdelta",
+ basis_eval=kdelta_eval_symbols,
+ )
+
+ for data in it:
+ psi = data.test_shape
+ grad_phi = data.trial_shape_grad
+ grad_psi = data.test_shape_grad
+ hessian_phi = data.trial_shape_hessian
+ hessian_affine = hessian_ref_to_affine(geometry, hessian_phi, jac_affine_inv)
+
+ hessian_affine_symbols = tabulation.register_factor(
+ "hessian_affine",
+ hessian_affine,
+ )
+
+ jac_affine_inv_T_grad_phi_symbols = tabulation.register_factor(
+ "jac_affine_inv_T_grad_phi",
+ jac_affine_inv.T * grad_phi,
+ )
+
+ jac_affine_inv_T_grad_psi_symbols = tabulation.register_factor(
+ "jac_affine_inv_T_grad_psi",
+ jac_affine_inv.T * grad_psi,
+ )
+
+ jac_blending_inv_T_jac_affine_inv_T_grad_psi_symbols = tabulation.register_factor(
+ "jac_affine_inv_T_grad_psi",
+ jac_blending_inv.T * jac_affine_inv_T_grad_psi_symbols,
+ )
+
+ if isinstance(blending, ExternalMap):
+ hessian_blending = hessian_affine_to_blending(geometry, hessian_affine, None, jac_blending_inv.T, jac_affine_inv_T_grad_phi_symbols, jacinvjac_blending)
+
+ laplacian = sum([hessian_blending[i, i] for i in range(geometry.dimensions)])
+
+ form = (
+ laplacian * psi
+ # * dot(u, jac_blending_inv_T_jac_affine_inv_T_grad_psi_symbols)
+ # * kdelta
+ * jac_affine_det
+ * jac_blending_det
+ )
+ elif isinstance(blending, IdentityMap):
+ laplacian = sum([hessian_affine_symbols[i, i] for i in range(geometry.dimensions)])
+ form = (
+ laplacian
+ * dot(u, jac_affine_inv_T_grad_psi_symbols)
+ * kdelta
+ * jac_affine_det
+ )
+ else:
+ hessian_blending = hessian_affine_to_blending(geometry, hessian_affine, hessian_blending_map, jac_blending_inv.T, jac_affine_inv_T_grad_phi_symbols)
+
+ laplacian = sum([hessian_blending[i, i, 0] for i in range(geometry.dimensions)])
+
+ form = (
+ laplacian
+ * dot(u, jac_blending_inv_T_jac_affine_inv_T_grad_psi_symbols)
+ * kdelta
+ * jac_affine_det
+ * jac_blending_det
+ )
+ # HOGException("Only for testing with Blending map")
+
+ mat[data.row, data.col] = form
+
+ return Form(mat, tabulation, symmetric=False, docstring=docstring)
def zero_form(
trial: FunctionSpace, test: FunctionSpace, geometry: ElementGeometry