diff --git a/generate_all_operators.py b/generate_all_operators.py
index 33edc497fbe9a2f9e94e8f370ba1215c503410da..7feef9fc015ea610e29b4c3f7636e3d08b85fb03 100644
--- a/generate_all_operators.py
+++ b/generate_all_operators.py
@@ -41,6 +41,7 @@ from hog.forms import (
epsilon,
epsilon_rotation,
divergence_rotation,
+ gradient_rotation,
pspg,
full_stokes,
nonlinear_diffusion,
@@ -641,7 +642,7 @@ def all_operators(
name=f"DivergenceRotation",
trial_space=P1,
test_space=P2Vector,
- form=divergence_rotation,
+ form=partial(divergence_rotation, normals_fspace = P2Vector),
type_descriptor=type_descriptor,
geometries=list(geometries),
opts=opts,
@@ -655,7 +656,7 @@ def all_operators(
name=f"GradientRotation",
trial_space=P2Vector,
test_space=P1,
- form=partial(divergence_rotation, transpose = True),
+ form=partial(gradient_rotation, normals_fspace = P2Vector),
type_descriptor=type_descriptor,
geometries=list(geometries),
opts=opts,
@@ -663,21 +664,19 @@ def all_operators(
)
)
- ops.append(
- OperatorInfo(
- mapping=f"P1",
- name=f"PSPG",
- trial_space=P1,
- test_space=P1,
- form=pspg,
- type_descriptor=type_descriptor,
- geometries=list(geometries),
- opts=opts,
- blending=blending,
- )
- )
-
-
+ # ops.append(
+ # OperatorInfo(
+ # mapping=f"P1",
+ # name=f"PSPG",
+ # trial_space=P1,
+ # test_space=P1,
+ # form=pspg,
+ # type_descriptor=type_descriptor,
+ # geometries=list(geometries),
+ # opts=opts,
+ # blending=blending,
+ # )
+ # )
for c in [0, 1, 2]:
# fmt: off
diff --git a/hog/forms.py b/hog/forms.py
index 142073feb65af88bc302d9e70d39b538c84e6809..4436d291fb63eb37b0e11eb4e8d2e9025d955e51 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -17,7 +17,7 @@
from dataclasses import dataclass
import logging
import sympy as sp
-from typing import List, Optional, Tuple
+from typing import List, Optional, Tuple, Union
from hog.ast import Operations, count_operations
from hog.element_geometry import ElementGeometry, TriangleElement, TetrahedronElement
@@ -35,7 +35,7 @@ from hog.fem_helpers import (
fem_function_on_element,
fem_function_gradient_on_element,
)
-from hog.function_space import FunctionSpace, EnrichedGalerkinFunctionSpace, N1E1Space
+from hog.function_space import FunctionSpace, EnrichedGalerkinFunctionSpace, N1E1Space, TensorialVectorFunctionSpace
from hog.math_helpers import dot, grad, inv, abs, det, double_contraction, e_vec
from hog.quadrature import Quadrature, Tabulation
from hog.symbolizer import Symbolizer
@@ -51,7 +51,7 @@ class RotationType(Enum):
def calculate_rotation_matrix(
- u_fspace: FunctionSpace,
+ u_fspace: TensorialVectorFunctionSpace,
geometry: ElementGeometry,
n_dof_symbols: List[List[sp.Matrix]],
) -> sp.MatrixBase:
@@ -813,26 +813,26 @@ where
def epsilon_rotation(
- trial: FunctionSpace,
- test: FunctionSpace,
+ trial: TensorialVectorFunctionSpace,
+ test: TensorialVectorFunctionSpace,
geometry: ElementGeometry,
symbolizer: Symbolizer,
+ coefficient_function_space: FunctionSpace,
blending: GeometryMap = IdentityMap(),
- component_trial: int = 0,
- component_test: int = 0,
variable_viscosity: bool = True,
- coefficient_function_space: Optional[FunctionSpace] = None,
) -> Form:
+ if not trial.is_vectorial:
+ raise HOGException(
+ "epsilon_rotation only applies for vectorial spaces"
+ )
+
if trial.is_vectorial ^ test.is_vectorial:
raise HOGException(
"Either both (trial and test) spaces or none should be vectorial."
)
-
- vectorial_spaces = trial.is_vectorial
+
docstring_components = (
""
- if vectorial_spaces
- else f"\nComponent trial: {component_trial}\nComponent test: {component_test}"
)
docstring = f"""
@@ -856,16 +856,6 @@ where
raise HOGException("Constant viscosity currently not supported.")
# TODO fix issue with undeclared p_affines
- if (
- not vectorial_spaces
- and geometry.dimensions < 3
- and (component_trial > 1 or component_test > 1)
- ):
- return create_empty_element_matrix(trial, test, geometry)
-
- # if geometry.dimensions > 2:
- # raise HOGException("Currently, only testing for 2D")
-
with TimedLogger("assembling epsilon matrix", level=logging.DEBUG):
tabulation = Tabulation(symbolizer)
@@ -1026,7 +1016,7 @@ where
epsForm = Form(
mat,
tabulation,
- symmetric=vectorial_spaces or (component_trial == component_test),
+ symmetric=True,
docstring=docstring,
rotmat=rotmat,
rot_type=RotationType.PRE_AND_POST_MULTIPLY,
@@ -1366,12 +1356,31 @@ Weak formulation
return Form(mat, Tabulation(symbolizer), symmetric=False, docstring=docstring)
+def gradient(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ component_index: int = 0,
+) -> Form:
+ """See divergence form. Just calls that with the transpose argument set to True."""
+ return divergence(
+ trial,
+ test,
+ geometry,
+ symbolizer,
+ blending=blending,
+ component_index=component_index,
+ transpose=True,
+ )
def divergence_rotation(
trial: FunctionSpace,
test: FunctionSpace,
geometry: ElementGeometry,
symbolizer: Symbolizer,
+ normals_fspace: TensorialVectorFunctionSpace,
blending: GeometryMap = IdentityMap(),
transpose: bool = False,
) -> Form:
@@ -1426,12 +1435,8 @@ Weak formulation
tabulation = Tabulation(symbolizer)
- if transpose:
- normals_function_space = test
- normals_component_function_space = test.component_function_space
- else:
- normals_function_space = trial
- normals_component_function_space = trial.component_function_space
+ normals_function_space = normals_fspace
+ normals_component_function_space = normals_fspace.component_function_space
phi_eval_symbols = tabulation.register_phi_evals(
normals_component_function_space.shape(geometry)
@@ -1499,26 +1504,25 @@ Weak formulation
)
-def gradient(
+def gradient_rotation(
trial: FunctionSpace,
test: FunctionSpace,
geometry: ElementGeometry,
symbolizer: Symbolizer,
+ normals_fspace: TensorialVectorFunctionSpace,
blending: GeometryMap = IdentityMap(),
- component_index: int = 0,
) -> Form:
"""See divergence form. Just calls that with the transpose argument set to True."""
- return divergence(
+ return divergence_rotation(
trial,
test,
geometry,
symbolizer,
+ normals_fspace,
blending=blending,
- component_index=component_index,
transpose=True,
)
-
def full_stokes(
trial: FunctionSpace,
test: FunctionSpace,