diff --git a/generate_all_hyteg_forms.py b/generate_all_hyteg_forms.py
index 867d0727a17725f7a74d62af935325e8673fa047..e33e68759cb9302fa4daf6e11ffa8869d52174b5 100644
--- a/generate_all_hyteg_forms.py
+++ b/generate_all_hyteg_forms.py
@@ -57,7 +57,9 @@ from hog.manifold_forms import (
vector_laplace_beltrami
)
from hog.manifold_facet_forms import (
- manifold_epsilon_sip_facet
+ manifold_epsilon_sip_facet,
+ manifold_epsilon_nip_facet,
+ manifold_divergence_facet
)
from hog.forms_vectorial import mass_n1e1, curl_curl
from hog.quadrature import Quadrature, select_quadrule
@@ -136,7 +138,7 @@ class FormInfo:
if self.trial_family == "EG" and self.test_family == "Lagrange":
return f"eg_to_p{self.test_degree}"
if self.trial_family == "Lagrange" and self.test_family == "EG":
- return f"p{self.test_degree}_to_eg"
+ return f"p{self.trial_degree}_to_eg"
elif self.trial_degree == self.test_degree:
return f"p{self.trial_degree}"
else:
@@ -697,6 +699,20 @@ for blending in [IdentityMap(), ExternalMap()]:
)
)
+form_infos.append(
+ FormInfo(
+ "manifold_normal_penalty",
+ trial_degree=1,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 3},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
form_infos.append(
FormInfo(
"manifold_vector_mass",
@@ -787,6 +803,8 @@ for trial_deg, test_deg, transpose in [
(2, 1, False),
(0, 2, True),
(2, 0, False),
+ (0, 1, True),
+ (1, 0, False)
]:
for blending in [IdentityMap(), ExternalMap()]:
if not transpose:
@@ -795,7 +813,7 @@ for trial_deg, test_deg, transpose in [
"manifold_vector_div",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={2: 3},
+ quad_schemes={2: 6},
row_dim=1,
col_dim=3,
is_implemented=is_implemented_for_vector_to_scalar,
@@ -808,7 +826,7 @@ for trial_deg, test_deg, transpose in [
"manifold_vector_divt",
trial_degree=trial_deg,
test_degree=test_deg,
- quad_schemes={2: 3},
+ quad_schemes={2: 6},
row_dim=3,
col_dim=1,
is_implemented=is_implemented_for_scalar_to_vector,
@@ -816,6 +834,34 @@ for trial_deg, test_deg, transpose in [
)
)
+form_infos.append(
+ FormInfo(
+ "manifold_vector_div",
+ trial_degree=1,
+ test_degree=0,
+ trial_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=ExternalMap(),
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_vector_divt",
+ trial_degree=0,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=ExternalMap(),
+ )
+)
+
form_infos.append(
FormInfo(
"manifold_epsilon",
@@ -872,6 +918,8 @@ form_infos.append(
)
)
+# --- Manifold facet epsilon forms ---
+
form_infos.append(
FormInfo(
"manifold_epsilon_sip_inner_facet",
@@ -958,6 +1006,202 @@ form_infos.append(
)
)
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_inner_facet",
+ trial_degree=1,
+ test_degree=1,
+ trial_family="EG",
+ quad_schemes={2: 3},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_inner_facet",
+ trial_degree=1,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 3},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_inner_facet",
+ trial_degree=1,
+ test_degree=1,
+ trial_family="EG",
+ test_family="EG",
+ quad_schemes={2: 3},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_outer_facet",
+ trial_degree=1,
+ test_degree=1,
+ trial_family="EG",
+ quad_schemes={2: 3},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_outer_facet",
+ trial_degree=1,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 3},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_epsilon_nip_outer_facet",
+ trial_degree=1,
+ test_degree=1,
+ trial_family="EG",
+ test_family="EG",
+ quad_schemes={2: 3},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_vector,
+ blending=blending,
+ )
+)
+
+# --- Manifold facet divergence forms ---
+
+form_infos.append(
+ FormInfo(
+ "manifold_div_inner_facet",
+ trial_degree=1,
+ test_degree=0,
+ quad_schemes={2: 6},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_div_outer_facet",
+ trial_degree=1,
+ test_degree=0,
+ quad_schemes={2: 6},
+ row_dim=3,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_divt_inner_facet",
+ trial_degree=0,
+ test_degree=1,
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_divt_outer_facet",
+ trial_degree=0,
+ test_degree=1,
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=3,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_div_inner_facet",
+ trial_degree=1,
+ test_degree=0,
+ trial_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_div_outer_facet",
+ trial_degree=1,
+ test_degree=0,
+ trial_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_vector_to_scalar,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_divt_inner_facet",
+ trial_degree=0,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=ExternalMap()
+ )
+)
+
+form_infos.append(
+ FormInfo(
+ "manifold_divt_outer_facet",
+ trial_degree=0,
+ test_degree=1,
+ test_family="EG",
+ quad_schemes={2: 6},
+ row_dim=1,
+ col_dim=1,
+ is_implemented=is_implemented_for_scalar_to_vector,
+ blending=ExternalMap()
+ )
+)
+
def form_func(
name: str,
row: int,
@@ -1094,6 +1338,54 @@ def form_func(
component_trial=col,
component_test=row,
).integrate(quad, symbolizer)
+ elif name.startswith("manifold_divt_inner_facet"):
+ return manifold_divergence_facet(
+ "inner",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ transpose=True,
+ blending=blending,
+ component_index=col
+ )
+ elif name.startswith("manifold_divt_outer_facet"):
+ return manifold_divergence_facet(
+ "outer",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ transpose=True,
+ blending=blending,
+ component_index=col
+ )
+ elif name.startswith("manifold_div_inner_facet"):
+ return manifold_divergence_facet(
+ "inner",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ transpose=False,
+ blending=blending,
+ component_index=row
+ )
+ elif name.startswith("manifold_div_outer_facet"):
+ return manifold_divergence_facet(
+ "outer",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ transpose=False,
+ blending=blending,
+ component_index=row
+ )
elif name.startswith("manifold_divt"):
return manifold_divergence(
trial,
@@ -1158,6 +1450,30 @@ def form_func(
trial_component=col,
test_component=row,
)
+ elif name.startswith("manifold_epsilon_nip_inner_facet"):
+ return manifold_epsilon_nip_facet(
+ "inner",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ blending,
+ trial_component=col,
+ test_component=row,
+ )
+ elif name.startswith("manifold_epsilon_nip_outer_facet"):
+ return manifold_epsilon_nip_facet(
+ "outer",
+ test,
+ trial,
+ geometry,
+ quad,
+ symbolizer,
+ blending,
+ trial_component=col,
+ test_component=row,
+ )
elif name.startswith("manifold_epsilon"):
return manifold_epsilon(
trial,
diff --git a/hog/__init__.py b/hog/__init__.py
index 3ab7a424466b0205f0a3eb903ca2bba6beea202a..5de3b234acab2511ecabb820269b3ab52e132386 100644
--- a/hog/__init__.py
+++ b/hog/__init__.py
@@ -32,6 +32,7 @@ __all__ = [
"hyteg_code_generation",
"hyteg_form_template",
"logger",
+ "manifold_facet_forms",
"math_helpers",
"multi_assignment",
"operator_generation",
diff --git a/hog/fem_helpers.py b/hog/fem_helpers.py
index 94b039ba1cdf64f076558c4dc74d850404f00d94..d476e3e2410443bbb3d01e90fe2dd63936a1579a 100644
--- a/hog/fem_helpers.py
+++ b/hog/fem_helpers.py
@@ -126,7 +126,7 @@ def trafo_shifted_ref_to_affine(
if affine_points is None:
affine_points = symbolizer.affine_vertices_as_vectors(
- geometry.dimensions, geometry.num_vertices
+ geometry.space_dimension, geometry.num_vertices
)
else:
if len(affine_points) != geometry.num_vertices:
@@ -284,6 +284,60 @@ def trafo_ref_to_physical(
return phy
+def trafo_shifted_ref_to_physical( #TODO: Remove?
+ geometry: ElementGeometry, symbolizer: Symbolizer, shift: sp.Matrix, blending: GeometryMap, affine_points: Optional[List[sp.Matrix]] = None
+) -> sp.Matrix:
+ """Returns the transformation of a point in the reference space to the physical element."""
+
+ if geometry not in blending.supported_geometries():
+ raise HOGException("Geometry not supported by blending map.")
+
+ t = trafo_shifted_ref_to_affine(geometry, symbolizer, shift, affine_points=affine_points)
+
+ if isinstance(blending, ExternalMap):
+ blending_class: type[MultiAssignment]
+ if isinstance(geometry, LineElement) and geometry.space_dimension == 3:
+ blending_class = BlendingFEmbeddedLine
+ if isinstance(geometry, TriangleElement):
+ blending_class = BlendingFTriangle
+ if geometry.space_dimension == 3:
+ blending_class = BlendingDFEmbeddedTriangle
+ elif isinstance(geometry, TetrahedronElement):
+ blending_class = BlendingFTetrahedron
+ else:
+ raise HOGException(
+ "Blending not implemented for the passed element geometry."
+ )
+
+ phy = sp.zeros(geometry.dimensions, 1)
+ for coord in range(geometry.dimensions):
+ phy[coord] = blending_class(sp.Symbol("blend"), coord, *t)
+
+ else:
+ phy = blending.evaluate(t)
+
+ return phy
+
+def jac_shifted_ref_to_physical(
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ shift: sp.Matrix,
+ affine_points: Optional[List[sp.Matrix]] = None,
+) -> sp.Matrix:
+ """Returns the Jacobian of the transformation from the shifted reference to the affine (computational) element.
+
+ :param geometry: the element geometry
+ :param symbolizer: a symbolizer object
+ :param shift: a shift on the reference space
+ :param affine_points: the symbols of the vertices of the affine element - if omitted, they default to the first d+1
+ vector symbols, useful for example if the trafo of two or more different elements is required
+ """
+
+ jac_affine = jac_ref_to_affine(geometry, symbolizer, affine_points=affine_points)
+ jac_blending = jac_shifted_affine_to_physical(geometry, symbolizer, shift, affine_points=affine_points)
+
+ return (jac_blending * jac_affine)
+
def jac_affine_to_physical(
geometry: ElementGeometry, symbolizer: Symbolizer, affine_points: Optional[List[sp.Matrix]] = None
@@ -319,7 +373,7 @@ def jac_affine_to_physical(
return jac
def jac_shifted_affine_to_physical(
- geometry: ElementGeometry, symbolizer: Symbolizer, shift: sp.Matrix
+ geometry: ElementGeometry, symbolizer: Symbolizer, shift: sp.Matrix, affine_points: Optional[List[sp.Matrix]] = None
) -> sp.Matrix:
"""Returns the Jacobian of the transformation from the affine (computational) to the physical element."""
blending_class: type[MultiAssignment]
@@ -332,7 +386,7 @@ def jac_shifted_affine_to_physical(
else:
raise HOGException("Blending not implemented for the passed element geometry.")
- t = trafo_ref_to_affine(geometry, symbolizer)
+ t = trafo_shifted_ref_to_affine(geometry, symbolizer, shift, affine_points=affine_points)
jac = sp.zeros(geometry.space_dimension, geometry.space_dimension)
rows, cols = jac.shape
for row in range(rows):
diff --git a/hog/function_space.py b/hog/function_space.py
index 3e1b048daa84cb68385b61d24c05eaab31f46144..0f3a76d596c9f3941d36fd4f5d52c682f8bbcfff 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -493,7 +493,8 @@ class EnrichedGalerkinFunctionSpace(FunctionSpace):
dof_map: Optional[List[int]] = None,
) -> List[sp.MatrixBase]:
# return [self._affine]
- return sp.Matrix([[1, 0], [0, 1]])
+ result = sp.Matrix([[1, 0], [0, 1]])
+ return [result]
def num_dofs(self, geometry: ElementGeometry) -> int:
"""Returns the number of DoFs per element."""
diff --git a/hog/manifold_facet_forms.py b/hog/manifold_facet_forms.py
index 5217af977b37ea98203839566e03a3114644888c..99edb4c7c4af932d24694da49f9f66a1fc318d5c 100644
--- a/hog/manifold_facet_forms.py
+++ b/hog/manifold_facet_forms.py
@@ -30,14 +30,17 @@ from hog.fem_helpers import (
element_matrix_iterator,
)
from hog.function_space import FunctionSpace, EnrichedGalerkinFunctionSpace, LagrangianFunctionSpace
-from hog.math_helpers import dot, inv, abs, det, double_contraction, e_vec
+from hog.math_helpers import dot, inv, abs, det, double_contraction, e_vec, vol
from hog.quadrature import Quadrature
from hog.symbolizer import Symbolizer
from hog.logger import TimedLogger
from hog.blending import GeometryMap, ExternalMap, IdentityMap
from hog.forms import Form
+from hog.external_functions import (
+ ScalarVariableCoefficient2D,
+)
-from hog.manifold_helpers import face_projection, embedded_normal, embedded_normal_grad
+from hog.manifold_helpers import face_projection, embedded_normal, embedded_normal_grad, grad_jac_ref_to_physical
from hog.forms_facets import _affine_element_vertices, trafo_ref_interface_to_ref_element, vol
@@ -144,15 +147,26 @@ def manifold_epsilon_sip_facet(
jac_total_I = jac_blending_I * jac_affine_I
fundamental_I = jac_total_I.T * jac_total_I
volume_interface = abs(det(fundamental_I))**0.5
- discontinuity_penalty = 3.0
+
+ #TODO: Make this from affine to physical
+ trafo_ref_interface_to_affine_interface = trafo_ref_to_affine(
+ volume_element_geometry,
+ symbolizer,
+ affine_points=affine_points_I + [affine_point_E_opposite],
+ )
+ penalty_coeff_class = ScalarVariableCoefficient2D
+ discontinuity_penalty = penalty_coeff_class(
+ sp.Symbol("discontinuity_penalty"), 0, *trafo_ref_interface_to_affine_interface
+ )
interface_affine_tangential = affine_points_I[1] - affine_points_I[0]
interface_tangential = jac_blending_I * interface_affine_tangential
- manifold_normal = embedded_normal(volume_element_geometry, symbolizer, blending)
+ manifold_normal_E1 = embedded_normal(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E1)
+ manifold_normal_E2 = embedded_normal(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E2)
outward_normal = sp.Matrix([
- interface_tangential[1] * manifold_normal[2] - interface_tangential[2] * manifold_normal[1],
- interface_tangential[2] * manifold_normal[0] - interface_tangential[0] * manifold_normal[2],
- interface_tangential[0] * manifold_normal[1] - interface_tangential[1] * manifold_normal[0],
+ interface_tangential[1] * manifold_normal_E1[2] - interface_tangential[2] * manifold_normal_E1[1],
+ interface_tangential[2] * manifold_normal_E1[0] - interface_tangential[0] * manifold_normal_E1[2],
+ interface_tangential[0] * manifold_normal_E1[1] - interface_tangential[1] * manifold_normal_E1[0],
])
outward_normal_orientation = -sp.sign(dot(external_outward_normal, affine_points_E1[0] + affine_points_E1[1] + affine_points_E1[2] - 1.5 * affine_points_I[0] - 1.5 * affine_points_I[1])[0, 0])
@@ -172,7 +186,13 @@ def manifold_epsilon_sip_facet(
#TODO: Implement variable viscousity
- projection_mat = face_projection(volume_element_geometry, symbolizer, blending)
+ projection_mat_E1 = face_projection(volume_element_geometry, symbolizer, blending, affine_points_E1)
+ projection_mat_E2 = face_projection(volume_element_geometry, symbolizer, blending, affine_points_E2)
+ # normal = embedded_normal(volume_element_geometry, symbolizer, blending)
+ normal_grad_E1 = embedded_normal_grad(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E1)
+ normal_grad_E2 = embedded_normal_grad(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E2)
+ grad_jac_total_E1 = grad_jac_ref_to_physical(volume_element_geometry, symbolizer, blending, affine_points_E1)
+ grad_jac_total_E2 = grad_jac_ref_to_physical(volume_element_geometry, symbolizer, blending, affine_points_E2)
with TimedLogger(
f"integrating {mat.shape[0] * mat.shape[1]} expressions",
@@ -230,33 +250,82 @@ def manifold_epsilon_sip_facet(
if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace):
if interface_type == "outer":
+ #grad_phi = jac_total_E2 * grad_phi
+ grad_phi_x = projection_mat_E2 * (jac_total_E2 * grad_phi.col(0) + grad_jac_total_E2[0] * phi)
+ grad_phi_y = projection_mat_E2 * (jac_total_E2 * grad_phi.col(1) + grad_jac_total_E2[1] * phi)
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
phi = jac_total_E2 * phi
- grad_phi = jac_total_E2 * grad_phi
else:
+ #grad_phi = jac_total_E1 * grad_phi
+ grad_phi_x = projection_mat_E1 * (jac_total_E1 * grad_phi.col(0) + grad_jac_total_E1[0] * phi)
+ grad_phi_y = projection_mat_E1 * (jac_total_E1 * grad_phi.col(1) + grad_jac_total_E1[1] * phi)
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
phi = jac_total_E1 * phi
- grad_phi = jac_total_E1 * grad_phi
elif isinstance(trial_element_2, LagrangianFunctionSpace):
- phi = projection_mat * phi
- grad_phi = projection_mat * grad_phi
+ #grad_phi = projection_mat * grad_phi
+ if interface_type == "outer":
+ grad_phi_x = projection_mat_E2 * (
+ grad_phi.col(0) - normal_grad_E2[0] * manifold_normal_E2.T * phi
+ )
+ grad_phi_y = projection_mat_E2 * (
+ grad_phi.col(1) - normal_grad_E2[1] * manifold_normal_E2.T * phi
+ )
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
+ phi = projection_mat_E2 * phi
+ else:
+ grad_phi_x = projection_mat_E1 * (
+ grad_phi.col(0) - normal_grad_E1[0] * manifold_normal_E1.T * phi
+ )
+ grad_phi_y = projection_mat_E1 * (
+ grad_phi.col(1) - normal_grad_E1[1] * manifold_normal_E1.T * phi
+ )
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
+ phi = projection_mat_E1 * phi
if isinstance(test_element_1, EnrichedGalerkinFunctionSpace):
+ #grad_psi = jac_total_E1 * grad_psi
+ grad_psi_x = projection_mat_E1 * (jac_total_E1 * grad_psi.col(0) + grad_jac_total_E1[0] * psi)
+ grad_psi_y = projection_mat_E1 * (jac_total_E1 * grad_psi.col(1) + grad_jac_total_E1[1] * psi)
+ grad_psi = (
+ grad_psi_x.row_join(grad_psi_y)
+ ).T
psi = jac_total_E1 * psi
- grad_psi = jac_total_E1 * grad_psi
elif isinstance(test_element_1, LagrangianFunctionSpace):
- psi = projection_mat * psi
- grad_psi = projection_mat * grad_psi
+ #grad_psi = projection_mat * grad_psi
+ grad_psi_x = projection_mat_E1 * (
+ grad_psi.col(0) - normal_grad_E1[0] * manifold_normal_E1.T * psi
+ )
+ grad_psi_y = projection_mat_E1 * (
+ grad_psi.col(1) - normal_grad_E1[1] * manifold_normal_E1.T * psi
+ )
+ grad_psi = (
+ grad_psi_x.row_join(grad_psi_y)
+ ).T
+ psi = projection_mat_E1 * psi
if interface_type == "inner":
+ scaled_grad_psi = jac_total_E1 * fundamental_inv_E1 * grad_psi
+ scaled_grad_phi = jac_total_E1 * fundamental_inv_E1 * grad_phi
+ epsilon_psi = 0.5 * (scaled_grad_psi + scaled_grad_psi.T)
+ epsilon_phi = 0.5 * (scaled_grad_phi + scaled_grad_phi.T)
form = (
(
- -1.0
+ -0.5
* dot(
- manifold_symm_grad(grad_psi, jac_total_E1, fundamental_inv_E1) * outward_normal,
+ epsilon_psi * outward_normal,
phi,
)[0, 0]
- - 1.0
+ - 0.5
* dot(
- manifold_symm_grad(grad_phi, jac_total_E1, fundamental_inv_E1) * outward_normal,
+ epsilon_phi * outward_normal,
psi,
)[0, 0]
+ (discontinuity_penalty / volume_interface)
@@ -265,19 +334,23 @@ def manifold_epsilon_sip_facet(
* volume_interface
)
elif interface_type == "outer":
+ scaled_grad_psi = jac_total_E1 * fundamental_inv_E1 * grad_psi
+ scaled_grad_phi = jac_total_E2 * fundamental_inv_E2 * grad_phi
+ epsilon_psi = 0.5 * (scaled_grad_psi + scaled_grad_psi.T)
+ epsilon_phi = 0.5 * (scaled_grad_phi + scaled_grad_phi.T)
form = (
(
- 1.0
+ 0.5
* dot(
- manifold_symm_grad(grad_psi, jac_total_E1, fundamental_inv_E1) * outward_normal,
+ epsilon_psi * outward_normal,
phi,
)[0, 0]
- - 1.0
+ - 0.5
* dot(
- manifold_symm_grad(grad_phi, jac_total_E2, fundamental_inv_E2) * outward_normal,
+ epsilon_phi * outward_normal,
psi,
)[0, 0]
- + (discontinuity_penalty / volume_interface)
+ - (discontinuity_penalty / volume_interface)
* dot(phi, psi)[0, 0]
)
* volume_interface
@@ -288,4 +361,510 @@ def manifold_epsilon_sip_facet(
mat[data.row, data.col] = facet_quad.integrate(form, symbolizer).subs(reference_symbols[volume_element_geometry.dimensions - 1], 0)
- return mat
\ No newline at end of file
+ return mat
+
+def manifold_epsilon_nip_facet(
+ interface_type: str,
+ test_element_1: FunctionSpace,
+ trial_element_2: FunctionSpace,
+ volume_element_geometry: ElementGeometry,
+ facet_quad: Quadrature,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ trial_component: Optional[int] = None,
+ test_component: Optional[int] = None
+) -> sp.Matrix:
+ """
+ See; diffusion_sip_facet_vectorial
+ """
+
+ if not (interface_type == "inner" or interface_type == "outer"):
+ raise HOGException("The manifold facet form is only implemented for interfaces of type inner and outer.")
+
+ if not isinstance(volume_element_geometry, TriangleElement) or not volume_element_geometry.space_dimension == 3:
+ raise HOGException("The manifold facet form is only implemented for embedded triangles.")
+
+ if not (isinstance(test_element_1, LagrangianFunctionSpace) or isinstance(test_element_1, EnrichedGalerkinFunctionSpace)):
+ raise HOGException("The manifold facet form is only implemented for Lagrangian function spaces or enriched Galerkin function spaces.")
+
+ if not (isinstance(trial_element_2, LagrangianFunctionSpace) or isinstance(trial_element_2, EnrichedGalerkinFunctionSpace)):
+ raise HOGException("The manifold facet form is only implemented for Lagrangian function spaces or enriched Galerkin function spaces.")
+
+ if isinstance(test_element_1, LagrangianFunctionSpace) and test_component is None:
+ raise HOGException("The manifold facet form needs a component index when working with Lagrangian basis functions.")
+
+ if isinstance(trial_element_2, LagrangianFunctionSpace) and trial_component is None:
+ raise HOGException("The manifold facet form needs a component index when working with Lagrangian basis functions.")
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace) and test_component is not None:
+ logging.info("INFO: Manifold_epsilon_facet form received EG basis functions and a component index. The index will be ignored.")
+
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace) and trial_component is not None:
+ logging.info("INFO: Manifold_epsilon_facet form received EG basis functions and a component index. The index will be ignored.")
+
+ logging.info("INFO: Manifold_epsilon_facet_form form ignores the gradient of the projection matrix for simplicity.")
+
+ with TimedLogger(
+ "assembling interface diffusion matrix with SIP", level=logging.DEBUG
+ ):
+ # Grabbing the symbols for the vertices of both elements, the interface, and the outward normal.
+ (
+ affine_points_E1,
+ affine_points_E2,
+ affine_points_I,
+ affine_point_E_opposite,
+ _,
+ external_outward_normal,
+ ) = _affine_element_vertices(volume_element_geometry, symbolizer)
+
+ # These Jacobians are for the trafo from the element-local reference space to the affine element space.
+ # Both are required for application of the chain rule.
+
+ jac_affine_E1 = jac_ref_to_affine(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E1
+ )
+ jac_affine_E2 = jac_ref_to_affine(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E2
+ )
+
+ jac_blending_E1 = jac_affine_to_physical(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E1
+ )
+ jac_blending_E2 = jac_affine_to_physical(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E2
+ )
+
+ jac_total_E1 = jac_blending_E1 * jac_affine_E1
+ jac_total_E2 = jac_blending_E2 * jac_affine_E2
+
+ # These trafos are required to evaluate correct points in the reference space of the elements.
+ # The reference coordinates of the interface (where the quadrature happens) are mapped to element ref space.
+
+ trafo_ref_interface_to_ref_element_E1 = trafo_ref_interface_to_ref_element(
+ volume_element_geometry, symbolizer, "inner"
+ )
+
+ trafo_ref_interface_to_ref_element_E2 = trafo_ref_interface_to_ref_element(
+ volume_element_geometry, symbolizer, "outer"
+ )
+
+ jac_affine_I = jac_ref_to_affine(
+ LineElement(space_dimension=3), symbolizer, affine_points=affine_points_I
+ )
+ jac_blending_I = jac_affine_to_physical(
+ LineElement(space_dimension=3), symbolizer, affine_points=affine_points_I
+ )
+ jac_total_I = jac_blending_I * jac_affine_I
+ fundamental_I = jac_total_I.T * jac_total_I
+ volume_interface = abs(det(fundamental_I))**0.5
+ #volume_interface = abs(vol(affine_points_I)) #TODO: Add blending
+ discontinuity_penalty = 100.0
+
+ #TODO: Make this from affine to physical
+ trafo_ref_interface_to_affine_interface = trafo_ref_to_affine(
+ volume_element_geometry,
+ symbolizer,
+ affine_points=affine_points_I + [affine_point_E_opposite],
+ )
+ penalty_coeff_class = ScalarVariableCoefficient2D
+ discontinuity_penalty = penalty_coeff_class(
+ sp.Symbol("discontinuity_penalty"), 0, *trafo_ref_interface_to_affine_interface
+ )
+
+ interface_affine_tangential = affine_points_I[1] - affine_points_I[0]
+ interface_tangential = jac_blending_I * interface_affine_tangential
+ manifold_normal = embedded_normal(volume_element_geometry, symbolizer, blending)
+ outward_normal = sp.Matrix([
+ interface_tangential[1] * manifold_normal[2] - interface_tangential[2] * manifold_normal[1],
+ interface_tangential[2] * manifold_normal[0] - interface_tangential[0] * manifold_normal[2],
+ interface_tangential[0] * manifold_normal[1] - interface_tangential[1] * manifold_normal[0],
+ ])
+
+ outward_normal_orientation = -sp.sign(dot(external_outward_normal, affine_points_E1[0] + affine_points_E1[1] + affine_points_E1[2] - 1.5 * affine_points_I[0] - 1.5 * affine_points_I[1])[0, 0])
+ outward_normal_norm = outward_normal_orientation * (outward_normal[0]**2 + outward_normal[1]**2 + outward_normal[2]**2)**0.5
+ outward_normal = outward_normal / outward_normal_norm
+
+ mat = create_empty_element_matrix(
+ trial_element_2, test_element_1, volume_element_geometry
+ )
+ it = element_matrix_iterator(
+ trial_element_2, test_element_1, volume_element_geometry
+ )
+
+ reference_symbols = symbolizer.ref_coords_as_list(
+ volume_element_geometry.dimensions
+ )
+
+ #TODO: Implement variable viscousity
+
+ projection_mat = face_projection(volume_element_geometry, symbolizer, blending)
+ normal = embedded_normal(volume_element_geometry, symbolizer, blending)
+ normal_grad = embedded_normal_grad(volume_element_geometry, symbolizer, blending)
+ grad_jac_total_E1 = grad_jac_ref_to_physical(volume_element_geometry, symbolizer, blending, affine_points_E1)
+ grad_jac_total_E2 = grad_jac_ref_to_physical(volume_element_geometry, symbolizer, blending, affine_points_E2)
+
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for data in it:
+ # TODO: fix this by introducing extra symbols for the shape functions
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace):
+ phi = data.trial_shape
+ grad_phi = data.trial_shape_grad
+ elif isinstance(trial_element_2, LagrangianFunctionSpace):
+ phi_val = data.trial_shape
+ phi = e_vec(volume_element_geometry.space_dimension, trial_component) * phi_val
+ grad_phi = (e_vec(volume_element_geometry.space_dimension, trial_component) * phi_val).jacobian(reference_symbols)
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace):
+ psi = data.test_shape
+ grad_psi = data.test_shape_grad
+ elif isinstance(test_element_1, LagrangianFunctionSpace):
+ psi_val = data.test_shape
+ psi = e_vec(volume_element_geometry.space_dimension, test_component) * psi_val
+ grad_psi = (e_vec(volume_element_geometry.space_dimension, test_component) * psi_val).jacobian(reference_symbols)
+
+ shape_symbols = ["xi_shape_0", "xi_shape_1", "xi_shape_2"][
+ : volume_element_geometry.dimensions
+ ]
+ phi = phi.subs(zip(reference_symbols, shape_symbols))
+ psi = psi.subs(zip(reference_symbols, shape_symbols))
+ grad_phi = grad_phi.subs(zip(reference_symbols, shape_symbols))
+ grad_psi = grad_psi.subs(zip(reference_symbols, shape_symbols))
+
+ # The evaluation is performed on the reference space of the _elements_ and _not_ in the reference space
+ # of the interface. This way we do not need to sort DoFs. So apply the trafo to the trial functions.
+ if interface_type == "outer":
+ phi = phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E2)
+ )
+ grad_phi = grad_phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E2)
+ )
+ else:
+ phi = phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+ grad_phi = grad_phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+
+ psi = psi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+ grad_psi = grad_psi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace):
+ if interface_type == "outer":
+ #grad_phi = jac_total_E2 * grad_phi
+ grad_phi_x = projection_mat * (jac_total_E2 * grad_phi.col(0) + grad_jac_total_E2[0] * phi)
+ grad_phi_y = projection_mat * (jac_total_E2 * grad_phi.col(1) + grad_jac_total_E2[1] * phi)
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
+ grad_phi = (jac_total_E2 * grad_phi).T #DEBUG
+ phi = jac_total_E2 * phi
+ else:
+ #grad_phi = jac_total_E1 * grad_phi
+ grad_phi_x = projection_mat * (jac_total_E1 * grad_phi.col(0) + grad_jac_total_E1[0] * phi)
+ grad_phi_y = projection_mat * (jac_total_E1 * grad_phi.col(1) + grad_jac_total_E1[1] * phi)
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
+ grad_phi = (jac_total_E1 * grad_phi).T #DEBUG
+ phi = jac_total_E1 * phi
+ elif isinstance(trial_element_2, LagrangianFunctionSpace):
+ #grad_phi = projection_mat * grad_phi
+ grad_phi_x = projection_mat * (
+ grad_phi.col(0) - normal_grad[0] * normal.T * phi
+ )
+ grad_phi_y = projection_mat * (
+ grad_phi.col(1) - normal_grad[1] * normal.T * phi
+ )
+ grad_phi = (
+ grad_phi_x.row_join(grad_phi_y)
+ ).T
+ phi = projection_mat * phi
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace):
+ #grad_psi = jac_total_E1 * grad_psi
+ grad_psi_x = projection_mat * (jac_total_E1 * grad_psi.col(0) + grad_jac_total_E1[0] * psi)
+ grad_psi_y = projection_mat * (jac_total_E1 * grad_psi.col(1) + grad_jac_total_E1[1] * psi)
+ grad_psi = (
+ grad_psi_x.row_join(grad_psi_y)
+ ).T
+ grad_psi = (jac_total_E1 * grad_psi).T #DEBUG
+ psi = jac_total_E1 * psi
+ elif isinstance(test_element_1, LagrangianFunctionSpace):
+ #grad_psi = projection_mat * grad_psi
+ grad_psi_x = projection_mat * (
+ grad_psi.col(0) - normal_grad[0] * normal.T * psi
+ )
+ grad_psi_y = projection_mat * (
+ grad_psi.col(1) - normal_grad[1] * normal.T * psi
+ )
+ grad_psi = (
+ grad_psi_x.row_join(grad_psi_y)
+ ).T
+ psi = projection_mat * psi
+
+ if interface_type == "inner":
+ form = (
+ (
+ + (discontinuity_penalty / volume_interface**0.5)
+ * dot(phi, psi)[0, 0]
+ )
+ * volume_interface
+ )
+ elif interface_type == "outer":
+ form = (
+ (
+ - (discontinuity_penalty / volume_interface**0.5)
+ * dot(phi, psi)[0, 0]
+ )
+ * volume_interface
+ )
+
+ elif interface_type == "dirichlet":
+ raise HOGException("Not implemented")
+
+ mat[data.row, data.col] = facet_quad.integrate(form, symbolizer).subs(reference_symbols[volume_element_geometry.dimensions - 1], 0)
+ return mat
+
+def manifold_divergence_facet(
+ interface_type: str,
+ test_element_1: FunctionSpace,
+ trial_element_2: FunctionSpace,
+ volume_element_geometry: ElementGeometry,
+ facet_quad: Quadrature,
+ symbolizer: Symbolizer,
+ transpose: bool,
+ blending: GeometryMap = IdentityMap(),
+ component_index: Optional[int] = None,
+) -> sp.Matrix:
+ """
+ See; diffusion_sip_facet_vectorial
+ """
+
+ #TODO: Make sure all these checks are really necessary and correct.
+
+ if not (interface_type == "inner" or interface_type == "outer"):
+ raise HOGException("The manifold facet form is only implemented for interfaces of type inner and outer.")
+
+ if not isinstance(volume_element_geometry, TriangleElement) or not volume_element_geometry.space_dimension == 3:
+ raise HOGException("The manifold facet form is only implemented for embedded triangles.")
+
+ if not (isinstance(test_element_1, LagrangianFunctionSpace) or isinstance(test_element_1, EnrichedGalerkinFunctionSpace)):
+ raise HOGException("The manifold facet form is only implemented for Lagrangian function spaces or enriched Galerkin function spaces.")
+
+ if not (isinstance(trial_element_2, LagrangianFunctionSpace) or isinstance(trial_element_2, EnrichedGalerkinFunctionSpace)):
+ raise HOGException("The manifold facet form is only implemented for Lagrangian function spaces or enriched Galerkin function spaces.")
+
+ if isinstance(test_element_1, LagrangianFunctionSpace) and component_index is None:
+ raise HOGException("The manifold facet form needs a component index when working with Lagrangian basis functions.")
+
+ if isinstance(trial_element_2, LagrangianFunctionSpace) and component_index is None:
+ raise HOGException("The manifold facet form needs a component index when working with Lagrangian basis functions.")
+
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace) and transpose:
+ raise HOGException("The transposed form can only use scalar valued functions for trial space.")
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace) and not transpose:
+ raise HOGException("The non-transposed form can only use sclar valued functions for test space.")
+
+ logging.info("INFO: Manifold_epsilon_facet_form form ignores the gradient of the projection matrix for simplicity.")
+
+ with TimedLogger(
+ "assembling interface diffusion matrix with SIP", level=logging.DEBUG
+ ):
+ # Grabbing the symbols for the vertices of both elements, the interface, and the outward normal.
+ (
+ affine_points_E1,
+ affine_points_E2,
+ affine_points_I,
+ affine_point_E_opposite,
+ _,
+ external_outward_normal,
+ ) = _affine_element_vertices(volume_element_geometry, symbolizer)
+
+ # These Jacobians are for the trafo from the element-local reference space to the affine element space.
+ # Both are required for application of the chain rule.
+
+ jac_affine_E1 = jac_ref_to_affine(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E1
+ )
+ jac_affine_E2 = jac_ref_to_affine(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E2
+ )
+
+ jac_blending_E1 = jac_affine_to_physical(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E1
+ )
+ jac_blending_E2 = jac_affine_to_physical(
+ volume_element_geometry, symbolizer, affine_points=affine_points_E2
+ )
+
+ jac_total_E1 = jac_blending_E1 * jac_affine_E1
+ jac_total_E2 = jac_blending_E2 * jac_affine_E2
+
+ # These trafos are required to evaluate correct points in the reference space of the elements.
+ # The reference coordinates of the interface (where the quadrature happens) are mapped to element ref space.
+
+ trafo_ref_interface_to_ref_element_E1 = trafo_ref_interface_to_ref_element(
+ volume_element_geometry, symbolizer, "inner"
+ )
+
+ trafo_ref_interface_to_ref_element_E2 = trafo_ref_interface_to_ref_element(
+ volume_element_geometry, symbolizer, "outer"
+ )
+
+ jac_affine_I = jac_ref_to_affine(
+ LineElement(space_dimension=3), symbolizer, affine_points=affine_points_I
+ )
+ jac_blending_I = jac_affine_to_physical(
+ LineElement(space_dimension=3), symbolizer, affine_points=affine_points_I
+ )
+ jac_total_I = jac_blending_I * jac_affine_I
+ fundamental_I = jac_total_I.T * jac_total_I
+ volume_interface = abs(det(fundamental_I))**0.5
+ #volume_interface = abs(vol(affine_points_I)) #TODO: Add blending
+
+ interface_affine_tangential = affine_points_I[1] - affine_points_I[0]
+ interface_tangential = jac_blending_I * interface_affine_tangential
+ manifold_normal = embedded_normal(volume_element_geometry, symbolizer, blending)
+ outward_normal = sp.Matrix([
+ interface_tangential[1] * manifold_normal[2] - interface_tangential[2] * manifold_normal[1],
+ interface_tangential[2] * manifold_normal[0] - interface_tangential[0] * manifold_normal[2],
+ interface_tangential[0] * manifold_normal[1] - interface_tangential[1] * manifold_normal[0],
+ ])
+
+ outward_normal_orientation = -sp.sign(dot(external_outward_normal, affine_points_E1[0] + affine_points_E1[1] + affine_points_E1[2] - 1.5 * affine_points_I[0] - 1.5 * affine_points_I[1])[0, 0])
+ outward_normal_norm = outward_normal_orientation * (outward_normal[0]**2 + outward_normal[1]**2 + outward_normal[2]**2)**0.5
+ outward_normal = outward_normal / outward_normal_norm
+
+ mat = create_empty_element_matrix(
+ trial_element_2, test_element_1, volume_element_geometry
+ )
+ it = element_matrix_iterator(
+ trial_element_2, test_element_1, volume_element_geometry
+ )
+
+ reference_symbols = symbolizer.ref_coords_as_list(
+ volume_element_geometry.dimensions
+ )
+
+ #TODO: Implement variable viscousity
+
+ projection_mat_E1 = face_projection(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E1)
+ projection_mat_E2 = face_projection(volume_element_geometry, symbolizer, blending, affine_points=affine_points_E2)
+
+ with TimedLogger(
+ f"integrating {mat.shape[0] * mat.shape[1]} expressions",
+ level=logging.DEBUG,
+ ):
+ for data in it:
+ # TODO: fix this by introducing extra symbols for the shape functions
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace):
+ phi = data.trial_shape
+ elif isinstance(trial_element_2, LagrangianFunctionSpace):
+ if not transpose:
+ phi = e_vec(volume_element_geometry.space_dimension, component_index) * data.trial_shape
+ else:
+ phi = data.trial_shape
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace):
+ psi = data.test_shape
+ elif isinstance(test_element_1, LagrangianFunctionSpace):
+ if transpose:
+ psi = e_vec(volume_element_geometry.space_dimension, component_index) * data.test_shape
+ else:
+ psi = data.test_shape
+
+
+ shape_symbols = ["xi_shape_0", "xi_shape_1", "xi_shape_2"][
+ : volume_element_geometry.dimensions
+ ]
+ phi = phi.subs(zip(reference_symbols, shape_symbols))
+ psi = psi.subs(zip(reference_symbols, shape_symbols))
+
+ # The evaluation is performed on the reference space of the _elements_ and _not_ in the reference space
+ # of the interface. This way we do not need to sort DoFs. So apply the trafo to the trial functions.
+ if interface_type == "outer":
+ phi = phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E2)
+ )
+ else:
+ phi = phi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+
+ psi = psi.subs(
+ zip(shape_symbols, trafo_ref_interface_to_ref_element_E1)
+ )
+
+ if isinstance(trial_element_2, EnrichedGalerkinFunctionSpace):
+ if interface_type == "outer":
+ phi = jac_total_E2 * phi
+ else:
+ phi = jac_total_E1 * phi
+ elif isinstance(trial_element_2, LagrangianFunctionSpace) and not transpose:
+ if interface_type == "outer":
+ phi = projection_mat_E2 * phi
+ else:
+ phi = projection_mat_E1 * phi
+
+ if isinstance(test_element_1, EnrichedGalerkinFunctionSpace):
+ psi = jac_total_E1 * psi
+ elif isinstance(test_element_1, LagrangianFunctionSpace) and transpose:
+ psi = projection_mat_E1 * psi
+
+ #TODO: Signs should be flipped, but manifold_vector_divergence also has "wrong" signs so this is necessary to make everything compatible.
+ if interface_type == "inner":
+ if not transpose:
+ form = (
+ -0.5 * dot(
+ phi,
+ outward_normal
+ )
+ * psi
+ * volume_interface
+ )
+ else:
+ form = (
+ -0.5 * dot(
+ psi,
+ outward_normal
+ )
+ * phi
+ * volume_interface
+ )
+ elif interface_type == "outer":
+ if not transpose:
+ form = (
+ 0.5 * dot(
+ phi,
+ outward_normal
+ )
+ * psi
+ * volume_interface
+ )
+ else:
+ form = (
+ -0.5 * dot(
+ psi,
+ outward_normal
+ )
+ * phi
+ * volume_interface
+ )
+
+ elif interface_type == "dirichlet":
+ raise HOGException("Not implemented")
+
+ mat[data.row, data.col] = facet_quad.integrate(form, symbolizer, blending).subs(reference_symbols[volume_element_geometry.dimensions - 1], 0)
+
+ return mat
+
diff --git a/hog/manifold_forms.py b/hog/manifold_forms.py
index aba0aa9bbd41a67c3d705429e5eeebde6d7809d2..6575af9c347fdf71ebc12a86b0973d1039da849e 100644
--- a/hog/manifold_forms.py
+++ b/hog/manifold_forms.py
@@ -35,7 +35,7 @@ from hog.logger import TimedLogger
from hog.blending import GeometryMap, ExternalMap, IdentityMap
from hog.forms import Form
-from hog.manifold_helpers import face_projection, embedded_normal, embedded_normal_grad
+from hog.manifold_helpers import face_projection, embedded_normal, embedded_normal_grad, grad_jac_ref_to_physical
def laplace_beltrami(
@@ -251,10 +251,7 @@ Weak formulation
projected_psi = projection * psi_vec
form = sp.Matrix(
- tabulation.register_factor(
- "manifold_vector_mass_symbol",
- dot(projected_phi, projected_psi) * fundamental_form_det**0.5,
- )
+ dot(projected_phi, projected_psi) * fundamental_form_det**0.5,
)[0]
mat[data.row, data.col] = form
@@ -325,10 +322,7 @@ Weak formulation
psi_normal = dot(psi_vec, normal)
form = sp.Matrix(
- tabulation.register_factor(
- "manifold_normal_penalty_symbol",
- phi_normal * psi_normal * fundamental_form_det**0.5,
- )
+ phi_normal * psi_normal * fundamental_form_det**0.5,
)[0]
mat[data.row, data.col] = form
@@ -452,6 +446,7 @@ Weak formulation
projection_mat = face_projection(geometry, symbolizer, blending)
normal = embedded_normal(geometry, symbolizer, blending)
normal_grad = embedded_normal_grad(geometry, symbolizer, blending)
+ grad_jac_total = grad_jac_ref_to_physical(geometry, symbolizer, blending)
jac_affine = jac_ref_to_affine(geometry, symbolizer)
@@ -469,40 +464,87 @@ Weak formulation
it = element_matrix_iterator(trial, test, geometry)
for data in it:
- ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions - 1)
+ ref_symbols_list = symbolizer.ref_coords_as_list(geometry.dimensions)
if not transpose:
- phi_vec_jac = (
- e_vec(geometry.dimensions, component_index) * data.trial_shape
- ).jacobian(ref_symbols_list)
- phi_vec_raw = (
- e_vec(geometry.dimensions, component_index) * data.trial_shape
- )
- phi = data.test_shape
+ if isinstance(trial, EnrichedGalerkinFunctionSpace):
+ raw_phi_grad = data.trial_shape_grad
+ raw_phi = data.trial_shape
+ phi_grad_x = projection_mat * (
+ jac_total * raw_phi_grad.col(0) + grad_jac_total[0] * raw_phi
+ )
+ phi_grad_y = projection_mat * (
+ jac_total * raw_phi_grad.col(1) + grad_jac_total[1] * raw_phi
+ )
+ phi = (
+ phi_grad_x.row_join(phi_grad_y)
+ ).T
+
+ phi = (jac_total * raw_phi_grad).T #DEBUG
+ # phi = (jac_total * raw_phi_grad).T
+ else:
+ raw_phi_grad = (
+ e_vec(geometry.space_dimension, component_index) * data.trial_shape
+ ).jacobian(ref_symbols_list)
+ raw_phi = (
+ e_vec(geometry.space_dimension, component_index) * data.trial_shape
+ )
+ phi_grad_x = projection_mat * (
+ raw_phi_grad.col(0) - normal_grad[0] * normal.T * raw_phi #- normal_grad[0] * normal.T * raw_phi
+ )
+ phi_grad_y = projection_mat * (
+ raw_phi_grad.col(1) - normal_grad[1] * normal.T * raw_phi #- normal_grad[1] * normal.T * raw_phi
+ )
+ phi = (
+ phi_grad_x.row_join(phi_grad_y)
+ ).T
+ psi = data.test_shape
else:
phi = data.trial_shape
- phi_vec_jac = (
- e_vec(geometry.dimensions, component_index) * data.test_shape
- ).jacobian(ref_symbols_list)
- phi_vec_raw = (
- e_vec(geometry.dimensions, component_index) * data.test_shape
- )
-
- phi_vec_x = projection_mat * (
- phi_vec_jac.col(0) - normal_grad[0] * normal.T * phi_vec_raw
- )
- phi_vec_y = projection_mat * (
- phi_vec_jac.col(1) - normal_grad[1] * normal.T * phi_vec_raw
- )
-
- phi_vec = (
- phi_vec_x.row_join(phi_vec_y)
- ).T
+ if isinstance(test, EnrichedGalerkinFunctionSpace):
+ raw_psi_grad = data.test_shape_grad
+ raw_psi = data.test_shape
+ psi_grad_x = projection_mat * (
+ jac_total * raw_psi_grad.col(0) + grad_jac_total[0] * raw_psi
+ )
+ psi_grad_y = projection_mat * (
+ jac_total * raw_psi_grad.col(1) + grad_jac_total[1] * raw_psi
+ )
+ psi = (
+ psi_grad_x.row_join(psi_grad_y)
+ ).T
+
+ psi = (jac_total * raw_psi_grad).T # DEBUG
+ #psi = (jac_total * raw_psi_grad).T
+ else:
+ raw_psi_grad = (
+ e_vec(geometry.space_dimension, component_index) * data.test_shape
+ ).jacobian(ref_symbols_list)
+ raw_psi = (
+ e_vec(geometry.space_dimension, component_index) * data.test_shape
+ )
+ psi_grad_x = projection_mat * (
+ raw_psi_grad.col(0) - normal_grad[0] * normal.T * raw_psi
+ )
+ psi_grad_y = projection_mat * (
+ raw_psi_grad.col(1) - normal_grad[1] * normal.T * raw_psi
+ )
+ psi = (
+ psi_grad_x.row_join(psi_grad_y)
+ ).T
- form = (
- (jac_total * fundamental_form_inv * phi_vec).trace()
- * phi
- * fundamental_form_det**0.5
- )
+ if not transpose:
+ form = (
+ (jac_total * fundamental_form_inv * phi).trace()
+ * psi
+ * fundamental_form_det**0.5
+ )
+ else:
+ form = (
+ (jac_total * fundamental_form_inv * psi).trace()
+ * phi
+ * fundamental_form_det**0.5
+ )
+ #form = normal_grad[0][0] + normal_grad[0][1] + normal_grad[0][2]
mat[data.row, data.col] = form
@@ -606,8 +648,7 @@ Weak formulation
phi_epsilon = 0.5 * (phi_projected_grad + phi_projected_grad.T)
psi_epsilon = 0.5 * (psi_projected_grad + psi_projected_grad.T)
- form = tabulation.register_factor(
- "epsilon_epsilon_prod",
+ form = (
double_contraction(phi_epsilon, psi_epsilon)
* fundamental_form_det**0.5,
)[0]
diff --git a/hog/manifold_helpers.py b/hog/manifold_helpers.py
index 1e09a8f02bce431aeb7ab675fc7462de64aea5ac..6193ae5fbf103611dc8fd8b2946e943fa64ae727 100644
--- a/hog/manifold_helpers.py
+++ b/hog/manifold_helpers.py
@@ -18,11 +18,12 @@ from typing import List
import sympy as sp
from hog.exception import HOGException
+from typing import Optional
from hog.element_geometry import ElementGeometry, TriangleElement
from hog.symbolizer import Symbolizer
from hog.blending import GeometryMap, IdentityMap
-from hog.fem_helpers import jac_affine_to_physical, jac_shifted_affine_to_physical
+from hog.fem_helpers import jac_affine_to_physical, jac_shifted_affine_to_physical, jac_shifted_ref_to_physical
from hog.external_functions import BlendingFEmbeddedTriangle
@@ -30,6 +31,7 @@ def embedded_normal(
geometry: ElementGeometry,
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
+ affine_points: Optional[List[sp.Matrix]] = None
) -> sp.Matrix:
"""Returns an unoriented unit normal vector for an embedded triangle."""
@@ -38,16 +40,19 @@ def embedded_normal(
"Embedded normal vectors are only defined for triangles embedded in 3D space."
)
- vert_points = symbolizer.affine_vertices_as_vectors(
- geometry.space_dimension, geometry.num_vertices
- )
+ if affine_points is None:
+ vert_points = symbolizer.affine_vertices_as_vectors(
+ geometry.space_dimension, geometry.num_vertices
+ )
+ else:
+ vert_points = affine_points
span0 = vert_points[1] - vert_points[0]
span1 = vert_points[2] - vert_points[0]
if not isinstance(blending, IdentityMap):
- blending_jac = jac_affine_to_physical(geometry, symbolizer)
- span0 = blending_jac.T * span0
- span1 = blending_jac.T * span1
+ blending_jac = jac_affine_to_physical(geometry, symbolizer, affine_points=vert_points)
+ span0 = blending_jac * span0
+ span1 = blending_jac * span1
normal = sp.Matrix(
[
@@ -59,7 +64,11 @@ def embedded_normal(
normal = normal / (normal[0] ** 2 + normal[1] ** 2 + normal[2] ** 2) ** 0.5
return normal
-def embedded_normal_grad(geometry: ElementGeometry, symbolizer: Symbolizer, blending: GeometryMap = IdentityMap()) -> List[sp.Matrix]:
+def embedded_normal_grad(
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ affine_points: Optional[List[sp.Matrix]] = None) -> List[sp.Matrix]:
"""Returns an approximation for partial_x and partial_y of the unoriented unit normal vector for an embedded triangle."""
if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
@@ -72,20 +81,23 @@ def embedded_normal_grad(geometry: ElementGeometry, symbolizer: Symbolizer, blen
shiftX = sp.Matrix([[shiftSize],[0.0]])
shiftY = sp.Matrix([[0.0],[shiftSize]])
- vert_points = symbolizer.affine_vertices_as_vectors(geometry.dimensions, geometry.num_vertices)
+ if affine_points is None:
+ vert_points = symbolizer.affine_vertices_as_vectors(geometry.space_dimension, geometry.num_vertices)
+ else:
+ vert_points = affine_points
span0 = vert_points[1] - vert_points[0]
span1 = vert_points[2] - vert_points[0]
- blending_jac_here = jac_affine_to_physical(geometry, symbolizer)
- blending_jac_x = jac_shifted_affine_to_physical(geometry, symbolizer, shiftX)
- blending_jac_y = jac_shifted_affine_to_physical(geometry, symbolizer, shiftY)
+ blending_jac_here = jac_affine_to_physical(geometry, symbolizer, affine_points=vert_points)
+ blending_jac_x = jac_shifted_affine_to_physical(geometry, symbolizer, shiftX, affine_points=vert_points)
+ blending_jac_y = jac_shifted_affine_to_physical(geometry, symbolizer, shiftY, affine_points=vert_points)
- span0_here = blending_jac_here.T * span0
- span1_here = blending_jac_here.T * span1
- span0_x = blending_jac_x.T * span0
- span1_x = blending_jac_x.T * span1
- span0_y = blending_jac_y.T * span0
- span1_y = blending_jac_y.T * span1
+ span0_here = blending_jac_here * span0
+ span1_here = blending_jac_here * span1
+ span0_x = blending_jac_x * span0
+ span1_x = blending_jac_x * span1
+ span0_y = blending_jac_y * span0
+ span1_y = blending_jac_y * span1
normal_here = sp.Matrix([
span0_here[1] * span1_here[2] - span0_here[2] * span1_here[1],
@@ -117,6 +129,7 @@ def face_projection(
geometry: ElementGeometry,
symbolizer: Symbolizer,
blending: GeometryMap = IdentityMap(),
+ affine_points: Optional[List[sp.Matrix]] = None
) -> sp.Matrix:
"""Returns a projection matrix for an embedded triangle."""
@@ -125,7 +138,7 @@ def face_projection(
"Projection matrices are only defined for triangles embedded in 3D space."
)
- normal = embedded_normal(geometry, symbolizer, blending=blending)
+ normal = embedded_normal(geometry, symbolizer, blending=blending, affine_points=affine_points)
projection = sp.Matrix(
[
[1.0 - normal[0] ** 2, -normal[0] * normal[1], -normal[0] * normal[2]],
@@ -134,3 +147,29 @@ def face_projection(
]
)
return projection
+
+def grad_jac_ref_to_physical(
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap,
+ affine_points: Optional[List[sp.Matrix]] = None
+ ) -> List[sp.Matrix]:
+ """Returns an approximation for partial_x and partial_y of the Jacobi Matrix for the trafo from reference space to physical space."""
+
+ if not (isinstance(geometry, TriangleElement) and geometry.space_dimension == 3):
+ pass
+ #TODO: Check other cases or output warning.
+
+ shiftSize = 0.0001
+ null_shift = sp.Matrix([[0.0],[0.0]])
+ x_shift = sp.Matrix([[shiftSize],[0.0]])
+ y_shift = sp.Matrix([[0.0],[shiftSize]])
+
+ jac_here = jac_shifted_ref_to_physical(geometry, symbolizer, null_shift, affine_points=affine_points)
+ jac_x = jac_shifted_ref_to_physical(geometry, symbolizer, x_shift, affine_points=affine_points)
+ jac_y = jac_shifted_ref_to_physical(geometry, symbolizer, y_shift, affine_points=affine_points)
+
+ grad_x = (jac_x - jac_here) / shiftSize
+ grad_y = (jac_y - jac_here) / shiftSize
+
+ return [grad_x, grad_y]