diff --git a/hog/forms.py b/hog/forms.py
index c4aa334890e1cf8237154ad63db53761e72c1f66..14f9c06178b27332f70c1973431a1c9d2b7ec5c9 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -902,7 +902,7 @@ for details.
# we would want tau = -pow(volume, 2.0 / 3.0) / 12.0
# however pow(volume, 2.0 / 3.0) is not vectorized
# do three iterations of newton's method with starting value x0 = volume instead
- tau = -sp.sqrt(volume)/sp.S(162) - volume/sp.S(324) - volume/(sp.S(18)*sp.sqrt(sp.S(2)*sp.sqrt(volume)/sp.S(9) + volume/sp.S(9) + sp.S(2)*volume/(sp.S(3)*sp.sqrt(sp.S(2)*sp.sqrt(volume)/sp.S(3) + volume/sp.S(3))))) - volume/(sp.S(54)*sp.sqrt(sp.S(2)*sp.sqrt(volume)/sp.S(3) + volume/sp.S(3)))
+ tau = -sp.sqrt(volume)*sp.Rational(1,162) - volume*sp.Rational(1,324) - volume/(18*sp.sqrt(2*sp.sqrt(volume)*sp.Rational(1,9) + volume*sp.Rational(1,9) + 2*volume/(3*sp.sqrt(2*sp.sqrt(volume)*sp.Rational(1,3) + volume*sp.Rational(1,3))))) - volume/(54*sp.sqrt(2*sp.sqrt(volume)*sp.Rational(1,3) + volume*sp.Rational(1,3)))
mat = create_empty_element_matrix(trial, test, geometry)
@@ -1520,6 +1520,15 @@ def grad_rho_rho(
domain="reference",
function_id="rho",
)
+
+ # invRho (density)
+ invRho,invRho_values,invRho_basis,invRho_len,invRho_derivatives,invRho_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="invRho",
+ )
else:
raise HOGException(
"grad_rho_rho only defined for FEM function coefficients"
@@ -1541,7 +1550,7 @@ def grad_rho_rho(
grad_rho_component_trafo = (jac_blending_inv.T * jac_affine_inv.T * rho_grad)[component_index]
form = (
- grad_rho_component_trafo/rho[0] * affine_factor * jac_blending_det
+ invRho[0] * grad_rho_component_trafo * affine_factor * jac_blending_det
)
mat[data.row, data.col] = form
@@ -1877,13 +1886,13 @@ def shear_heating(
vRef = symbolizer.ref_coords_as_list(dimensions=geometry.dimensions)
dim = geometry.dimensions
- # rho (density)
- rho, _ = fem_function_on_element(
+ # invRho (density)
+ invRho, _ = fem_function_on_element(
rho_function_space,
geometry,
symbolizer,
domain="reference",
- function_id="rho",
+ function_id="invRho",
)
# eta (viscosity)
@@ -1960,13 +1969,13 @@ def shear_heating(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta[0] * epsilon_u
+ tau_u = 2 * eta[0] * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
# scaling
- scaling = sp.S(1) / rho[0]
+ scaling = invRho[0]
affine_factor = tabulation.register_factor("affine_factor_symbol", sp.Matrix([phi * psi * jac_affine_det]))[0]
@@ -2060,13 +2069,13 @@ def shear_heating_supg(
vRef = symbolizer.ref_coords_as_list(dimensions=geometry.dimensions)
dim = geometry.dimensions
- # rho (density)
- rho, _ = fem_function_on_element(
+ # invRho (density)
+ invRho, _ = fem_function_on_element(
rho_function_space,
geometry,
symbolizer,
domain="reference",
- function_id="rho",
+ function_id="invRho",
)
# eta (viscosity)
@@ -2169,13 +2178,13 @@ def shear_heating_supg(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta[0] * epsilon_u
+ tau_u = 2 * eta[0] * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
# scaling
- scaling = sp.S(1) / rho[0]
+ scaling = invRho[0]
jac_affine_inv_T_grad_psi = tabulation.register_factor("jac_affine_inv_T_grad_psi", sp.Matrix([jac_affine_inv.T * grad_psi]))
affine_factor = tabulation.register_factor("affine_factor_symbol", sp.Matrix([phi * jac_affine_det]))[0]
@@ -2252,6 +2261,15 @@ def diffusion_inv_rho(
domain="reference",
function_id="rho",
)
+
+ # invRho (density)
+ invRho,invRho_values,invRho_basis,invRho_len,invRho_derivatives,invRho_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="invRho",
+ )
else:
raise HOGException(
"diffusion_inv_rho only defined for FEM function coefficients"
@@ -2273,7 +2291,7 @@ def diffusion_inv_rho(
grad_rho_trafo = (jac_blending_inv.T * jac_affine_inv.T * rho_grad)
grad_phi_trafo = (jac_blending_inv.T * jac_affine_inv_T_grad_phi)
- scaling = psi / (rho[0]*rho[0])
+ scaling = psi * invRho[0] * invRho[0]
form = (
-scaling * (grad_phi_trafo.T * grad_rho_trafo)[0] * jac_blending_det
@@ -2903,13 +2921,13 @@ def diffusion_supg(
vRef_tuple = tuple([vRef[i] for i in range(geometry.dimensions)])
dim = geometry.dimensions
- # rho (density)
- rho, _ = fem_function_on_element(
+ # invRho (density)
+ invRho, _ = fem_function_on_element(
rho_function_space,
geometry,
symbolizer,
domain="reference",
- function_id="rho",
+ function_id="invRho",
)
# uX (velocity x direction)
@@ -3005,7 +3023,7 @@ def diffusion_supg(
supg = (u_vec.T * grad_psi_trafo)[0]
# scaling
- scaling = sp.S(1) / rho[0]
+ scaling = invRho[0]
form = (
-delta * scaling * laplace_phi_trafo * supg * jac_affine_det * jac_blending_det
@@ -3080,10 +3098,10 @@ def centroid_test_mass(
if dim > 2:
# expr = x*z + x - 2y + 3
- k = vCentroid[0] * vCentroid[2] + vCentroid[0] - 2*vCentroid[1] + sp.S(3)
+ k = vCentroid[0] * vCentroid[2] + vCentroid[0] - 2*vCentroid[1] + 3
else:
# expr = x*y + x - 2y + 3
- k = vCentroid[0] * vCentroid[1] + vCentroid[0] - 2*vCentroid[1] + sp.S(3)
+ k = vCentroid[0] * vCentroid[1] + vCentroid[0] - 2*vCentroid[1] + 3
mat = create_empty_element_matrix(trial, test, geometry)
@@ -4605,7 +4623,7 @@ def shear_heating_simple_visc(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta * epsilon_u
+ tau_u = 2 * eta * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
@@ -4811,7 +4829,7 @@ def shear_heating_supg_simple_visc(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta * epsilon_u
+ tau_u = 2 * eta * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
@@ -5276,7 +5294,7 @@ def shear_heating_frank_kamenetskii_simple_visc_base(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta * epsilon_u
+ tau_u = 2 * eta * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
@@ -5512,7 +5530,7 @@ def shear_heating_supg_frank_kamenetskii_simple_visc_base(
sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
epsilon_u = sigma_u + xi_u
- tau_u = sp.S(2) * eta * epsilon_u
+ tau_u = 2 * eta * epsilon_u
# double contraction
double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
@@ -5532,4 +5550,421 @@ def shear_heating_supg_frank_kamenetskii_simple_visc_base(
mat[data.row, data.col] = form
+ return Form(mat, tabulation, symmetric=False, docstring=docstring)
+
+
+def shear_heating_no_surface(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ coefficient_function_space: Optional[FunctionSpace] = None,
+ eta_function_space: Optional[FunctionSpace] = None,
+ rho_function_space: Optional[FunctionSpace] = None
+) -> Form:
+ docstring = f"""
+ Shear heating operator for the TALA
+
+ Intended for RHS use.
+
+ For the two dimensional case this is the pseudo-3D form of the operator.
+ Otherwise eps(u) would need to be := 1/2 ∇u + 1/2 (∇u)^T - 1/2 (∇ · u) I.
+
+ Typical usage sets T = 1, i.e. applying the operator to a function containing only ones.
+ TODO: Linear forms in hyteg
+
+ Suggested quad degrees:
+ Maximum: {(((1 if coefficient_function_space is None else coefficient_function_space.degree))-1)*2 + trial.degree + test.degree + 1 + (0 if isinstance(blending, IdentityMap) else 1)}
+ Reduced: {trial.degree + test.degree + 1 + (0 if isinstance(blending, IdentityMap) else 1)}
+ Minimum: {trial.degree + test.degree + 1}
+
+ Geometry map: {blending}
+
+ Weak formulation
+
+ T: trial function (scalar space: {trial})
+ w: test function (scalar space: {test})
+ u: coefficient (vector space: {coefficient_function_space})
+ rho: coefficient (scalar space: {rho_function_space})
+ eta: coefficient (scalar space: {eta_function_space})
+
+ ∫ 1/rho ( tau(u,eta) : eps(u) ) T * w
+
+ with
+
+ tau(u,eta) = 2 eta eps(u)
+ eps(u) := 1/2 ∇u + 1/2 (∇u)^T - 1/dim (∇ · u) I
+ I := Identity Matrix
+ """
+
+ with TimedLogger(
+ f"assembling shear_heating 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):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ affine_coords = trafo_ref_to_affine(geometry, symbolizer)
+ jac_blending = blending.jacobian(affine_coords)
+
+ jac_blending_det = abs(det(jac_blending))
+ with TimedLogger("inverting blending Jacobian", level=logging.DEBUG):
+ jac_blending_inv = inv(jac_blending)
+
+ if coefficient_function_space and eta_function_space and rho_function_space:
+ vRef = symbolizer.ref_coords_as_list(dimensions=geometry.dimensions)
+ dim = geometry.dimensions
+
+ # invRho (density)
+ invRho, _ = fem_function_on_element(
+ rho_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="invRho",
+ )
+
+ # eta (viscosity)
+ eta, _ = fem_function_on_element(
+ eta_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="eta",
+ )
+
+ # uX (velocity x direction)
+ uX,uX_values,uX_basis,uX_len,uX_derivatives,uX_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uX",
+ )
+
+ # uY (velocity y direction)
+ uY,uY_values,uY_basis,uY_len,uY_derivatives,uY_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uY",
+ )
+
+ # lists
+ u = [uX, uY]
+ u_values = [uX_values,uY_values]
+ u_derivatives = [uX_derivatives,uY_derivatives]
+ u_len = [uX_len, uY_len]
+
+ if dim > 2:
+ # uZ (velocity z direction)
+ uZ,uZ_values,uZ_basis,uZ_len,uZ_derivatives,uZ_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uZ",
+ )
+
+ u.append(uZ)
+ u_values.append(uZ_values)
+ u_derivatives.append(uZ_derivatives)
+ u_len.append(uZ_len)
+
+ grad_u_vec = sp.zeros(dim,dim)
+ for component_u in range(dim):
+ for j in range(dim):
+ # be aware that the jacobian is transposed here to get the gradient (i.e. j and component_u were swaped in comp to the jacobian)
+ grad_u_vec[j,component_u] = sum([u_values[component_u][i]*u_derivatives[component_u][i][j] for i in range(u_len[component_u])])
+ else:
+ raise HOGException(
+ "shear heating only defined for FEM function coefficients"
+ )
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+
+ it = element_matrix_iterator(trial, test, geometry)
+
+ # remove surface
+ if isinstance(blending, IdentityMap):
+ coords_on_element = trafo_ref_to_affine(geometry, symbolizer)
+ else:
+ coords_on_element = trafo_ref_to_physical(geometry, symbolizer, blending)
+
+ dim = geometry.dimensions
+ v = [coords_on_element[i] for i in range(dim)]
+ norm = sp.sqrt(sum([v[i]*v[i] for i in range(dim)]))
+
+ radiusSurface = sp.symbols("radiusSurface")
+ cutoff = sp.symbols("cutoff")
+ pos = radiusSurface - norm
+
+ posScaling = sp.Piecewise(
+ (0.0, pos < cutoff ),
+ (1.0, sp.sympify(True) )
+ )
+ # /remove surface
+
+ for data in it:
+ phi = data.trial_shape
+ psi = data.test_shape
+ grad_phi = data.trial_shape_grad
+ grad_psi = data.test_shape_grad
+
+ grad_u_vec_trafo = jac_blending_inv.T * jac_affine_inv.T * grad_u_vec
+ div_u = sum([grad_u_vec_trafo[component_u, component_u] for component_u in range(dim)])
+
+ sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
+ xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
+ epsilon_u = sigma_u + xi_u
+ tau_u = 2 * eta[0] * epsilon_u
+
+ # double contraction
+ double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
+
+ # scaling
+ scaling = invRho[0]
+
+ affine_factor = tabulation.register_factor("affine_factor_symbol", sp.Matrix([phi * psi * jac_affine_det]))[0]
+
+ form = (
+ posScaling * scaling * double_contrac * affine_factor * jac_blending_det
+ )
+
+ mat[data.row, data.col] = form
+
+ return Form(mat, tabulation, symmetric=False, docstring=docstring)
+
+def shear_heating_supg_no_surface(
+ trial: FunctionSpace,
+ test: FunctionSpace,
+ geometry: ElementGeometry,
+ symbolizer: Symbolizer,
+ blending: GeometryMap = IdentityMap(),
+ coefficient_function_space: Optional[FunctionSpace] = None,
+ eta_function_space: Optional[FunctionSpace] = None,
+ rho_function_space: Optional[FunctionSpace] = None,
+) -> Form:
+ docstring = f"""
+ Shear heating operator for the TALA
+
+ Intended for RHS use.
+
+ For the two dimensional case this is the pseudo-3D form of the operator.
+ Otherwise eps(u) would need to be := 1/2 ∇u + 1/2 (∇u)^T - 1/2 (∇ · u) I.
+
+ Typical usage sets T = 1, i.e. applying the operator to a function containing only ones.
+ TODO: Linear forms in hyteg
+
+ The scaling function for the supg stabilisation is hard coded into the form.
+
+ delta( u_abs ) := hMax / ( 2 * u_abs ) * xi( Pe )
+
+ with:
+ xi( Pe ) := ( 1 + 2 / ( exp( 2* Pe ) - 1 ) - 1 / Pe ) = coth( Pe ) - 1 / Pe
+ Pe := u_abs * h / ( 2 * k ) as the local Peclet number
+ h as the element diameter
+ k as the thermal conductivity coefficient
+ u_abs as the norm of the velocity vector at the element centroid
+
+ Note: h is calculated from the affine element geometry. If your blending map changes the element diameter too drastically, this will no longer work.
+
+ Suggested quad degrees:
+ Maximum: {(((1 if coefficient_function_space is None else coefficient_function_space.degree))-1)*2 + trial.degree + ((1 if coefficient_function_space is None else coefficient_function_space.degree)) + (test.degree-1) + 1 + (0 if isinstance(blending, IdentityMap) else 1)}
+ Reduced: {trial.degree + (test.degree-1) + 1 + (0 if isinstance(blending, IdentityMap) else 1)}
+ Minimum: {trial.degree + (test.degree-1) + 1}
+
+ Geometry map: {blending}
+
+ Weak formulation
+
+ T: trial function (scalar space: {trial})
+ w: test function (scalar space: {test})
+ u: coefficient (vector space: {coefficient_function_space})
+ rho: coefficient (scalar space: {rho_function_space})
+ eta: coefficient (scalar space: {eta_function_space})
+
+ ∫ delta 1/rho ( tau(u,eta) : eps(u) ) T * ( u · ∇w )
+
+ with
+
+ tau(u,eta) = 2 eta eps(u)
+ eps(u) := 1/2 ∇u + 1/2 (∇u)^T - 1/3 (∇ · u) I
+ I := Identity Matrix
+ """
+
+ with TimedLogger(
+ f"assembling shear_heating_supg 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):
+ jac_blending = jac_affine_to_physical(geometry, symbolizer)
+ else:
+ affine_coords = trafo_ref_to_affine(geometry, symbolizer)
+ jac_blending = blending.jacobian(affine_coords)
+
+ jac_blending_det = abs(det(jac_blending))
+ with TimedLogger("inverting blending Jacobian", level=logging.DEBUG):
+ jac_blending_inv = inv(jac_blending)
+
+ if coefficient_function_space and eta_function_space and rho_function_space:
+ vRef = symbolizer.ref_coords_as_list(dimensions=geometry.dimensions)
+ dim = geometry.dimensions
+
+ # invRho (density)
+ invRho, _ = fem_function_on_element(
+ rho_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="invRho",
+ )
+
+ # eta (viscosity)
+ eta, _ = fem_function_on_element(
+ eta_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="eta",
+ )
+
+ # uX (velocity x direction)
+ uX,uX_values,uX_basis,uX_len,uX_derivatives,uX_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uX",
+ )
+
+ # uY (velocity y direction)
+ uY,uY_values,uY_basis,uY_len,uY_derivatives,uY_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uY",
+ )
+
+ # lists
+ u = [uX, uY]
+ u_values = [uX_values,uY_values]
+ u_derivatives = [uX_derivatives,uY_derivatives]
+ u_len = [uX_len, uY_len]
+
+ if dim > 2:
+ # uZ (velocity z direction)
+ uZ,uZ_values,uZ_basis,uZ_len,uZ_derivatives,uZ_grad = fem_function_on_element_AHFC_compatibility_helper(
+ coefficient_function_space,
+ geometry,
+ symbolizer,
+ domain="reference",
+ function_id="uZ",
+ )
+
+ u.append(uZ)
+ u_values.append(uZ_values)
+ u_derivatives.append(uZ_derivatives)
+ u_len.append(uZ_len)
+
+ u_vec = sp.Matrix([[u[i]] for i in range(dim)])
+ grad_u_vec = sp.zeros(dim,dim)
+ for component_u in range(dim):
+ for j in range(dim):
+ # be aware that the jacobian is transposed here to get the gradient (i.e. j and component_u were swaped in comp to the jacobian)
+ grad_u_vec[j,component_u] = sum([u_values[component_u][i]*u_derivatives[component_u][i][j] for i in range(u_len[component_u])])
+ else:
+ raise HOGException(
+ "shear heating supg only defined for FEM function coefficients"
+ )
+
+ # init hard coded form
+ if isinstance(blending, IdentityMap):
+ coords_on_element = trafo_ref_to_affine(geometry, symbolizer)
+ else:
+ coords_on_element = trafo_ref_to_physical(geometry, symbolizer, blending)
+
+ v = [coords_on_element[i] for i in range(dim)]
+
+ # u_abs
+ centroid = {vRef[i]: sp.Rational(1,dim+1) for i in range(dim)}
+ uCentroid = [sp.simplify(u[i].subs(centroid)[0]) for i in range(dim)]
+ abs_u_centroid = sp.sqrt(sum([uCentroid[i]*uCentroid[i] for i in range(dim)]))
+
+ # h
+ affine_points = symbolizer.affine_vertices_as_vectors(geometry.dimensions, geometry.num_vertices)
+
+ h = diameter(affine_points)
+ # Alternative: h = sp.symbols("hMax")
+
+ # k
+ k = sp.symbols("thermalConductivity")
+
+ # delta function
+ delta = deltaSUPG(abs_u_centroid, h, k, True)
+
+ # remove surface
+ norm = sp.sqrt(sum([v[i]*v[i] for i in range(dim)]))
+
+ radiusSurface = sp.symbols("radiusSurface")
+ cutoff = sp.symbols("cutoff")
+ pos = radiusSurface - norm
+
+ posScaling = sp.Piecewise(
+ (0.0, pos < cutoff ),
+ (1.0, sp.sympify(True) )
+ )
+ # /remove surface
+
+ mat = create_empty_element_matrix(trial, test, geometry)
+
+ it = element_matrix_iterator(trial, test, geometry)
+
+ for data in it:
+ phi = data.trial_shape
+ psi = data.test_shape
+ grad_phi = data.trial_shape_grad
+ grad_psi = data.test_shape_grad
+
+ grad_u_vec_trafo = jac_blending_inv.T * jac_affine_inv.T * grad_u_vec
+ div_u = sum([grad_u_vec_trafo[component_u, component_u] for component_u in range(dim)])
+
+ sigma_u = sp.Rational(1,2) * (grad_u_vec_trafo + grad_u_vec_trafo.T)
+ xi_u = -sp.Rational(1,dim) * div_u * sp.eye(dim)
+ epsilon_u = sigma_u + xi_u
+ tau_u = 2 * eta[0] * epsilon_u
+
+ # double contraction
+ double_contrac = sum([tau_u[i, j] * epsilon_u[i, j] for i in range(dim) for j in range(dim)])
+
+ # scaling
+ scaling = invRho[0]
+
+ jac_affine_inv_T_grad_psi = tabulation.register_factor("jac_affine_inv_T_grad_psi", sp.Matrix([jac_affine_inv.T * grad_psi]))
+ affine_factor = tabulation.register_factor("affine_factor_symbol", sp.Matrix([phi * jac_affine_det]))[0]
+
+ grad_test_trafo = jac_blending_inv.T * jac_affine_inv_T_grad_psi
+ supg = (u_vec.T * grad_test_trafo)[0]
+
+ form = (
+ posScaling * delta * scaling * double_contrac * supg * affine_factor * jac_blending_det
+ )
+
+ mat[data.row, data.col] = form
+
return Form(mat, tabulation, symmetric=False, docstring=docstring)
\ No newline at end of file
diff --git a/hog/math_helpers.py b/hog/math_helpers.py
index 8b4d5791ad5d963aef5085cd7e78eaedfb736461..eb096be812ef0199e4d582ac0078761a47ce0c39 100644
--- a/hog/math_helpers.py
+++ b/hog/math_helpers.py
@@ -204,7 +204,7 @@ def vol(vertices: List[sp.Matrix]) -> sp.Expr:
bd = vertices[1]-vertices[3]
cd = vertices[2]-vertices[3]
- return sp.Abs(ad.dot(bd.cross(cd))) / sp.S(6)
+ return sp.Abs(ad.dot(bd.cross(cd))) * sp.Rational(1,6)
else:
raise HOGException(f"Not implemented for {len(vertices)} vertices")
@@ -219,7 +219,7 @@ def diameter(vertices: List[sp.Matrix]) -> sp.Expr:
b = norm(vertices[2]-vertices[0])
c = norm(vertices[0]-vertices[1])
- return sp.S(2) * ( a * b * c / (sp.S(4) * vol(vertices)) )
+ return 2 * ( a * b * c / (4 * vol(vertices)) )
elif len(vertices) == 4:
# a tetrahedron
a = norm(vertices[0]-vertices[3])
@@ -230,9 +230,9 @@ def diameter(vertices: List[sp.Matrix]) -> sp.Expr:
c_tilde = norm(vertices[0]-vertices[1])
# heron type formula (from Cayley Menger determinant)
- s = (a*a_tilde + b*b_tilde + c*c_tilde) / sp.S(2)
+ s = (a*a_tilde + b*b_tilde + c*c_tilde) * sp.Rational(1,2)
- return sp.S(2) * ( sp.sqrt(s * (s - a*a_tilde) * (s- b*b_tilde) * (s - c*c_tilde)) / (sp.S(6) * vol(vertices)) )
+ return 2 * ( sp.sqrt(s * (s - a*a_tilde) * (s- b*b_tilde) * (s - c*c_tilde)) / (6 * vol(vertices)) )
else:
raise HOGException(f"Not implemented for {len(vertices)} vertices")
@@ -242,7 +242,7 @@ def deltaSUPG(abs_u: sp.Expr, h: sp.Expr, k: sp. Expr, approximateXi: bool = Tru
The exp function cannot be vectorized in pystencils yet.
If approximateXi is True this returns an approximation that can be vectorized.
"""
- Pe = abs_u * h / ( sp.S(2) * k )
+ Pe = abs_u * h / ( 2 * k )
if approximateXi:
# Taylor near 0, Cubic spline in a second region, 1 - 1/Pe for the rest
@@ -254,7 +254,7 @@ def deltaSUPG(abs_u: sp.Expr, h: sp.Expr, k: sp. Expr, approximateXi: bool = Tru
((sp.S(1) - sp.S(1)/Pe) , sp.sympify(True) )
)
else:
- xi = sp.S(1) + sp.S(2) / ( sp.exp( sp.S(2) * Pe ) - sp.S(1) ) - sp.S(1) / Pe
+ xi = sp.S(1) + sp.S(2) / ( sp.exp( 2 * Pe ) - sp.S(1) ) - sp.S(1) / Pe
return h / ( 2 * abs_u ) * xi
@@ -276,12 +276,12 @@ def simpleViscosityProfile(x: sp.Expr) -> sp.Expr:
etaSimple = sp.Piecewise(
(1.25e+23, x <= 0.75),
- (3.54619628745821e+26*x**2 - 5.43063135866058e+26*x + 2.07948810730019e+26,x <= 0.758405172413793),
+ (3.54619628745821e+26*x**2 - 5.43063135866058e+26*x + 2.07948810730019e+26, x <= 0.758405172413793),
(1.86568634321297e+26*x**2 - 2.88161649064377e+26*x + 1.11289507706839e+26, x <= 0.764008620689655),
(1.06929133385375e+26*x**2 - 1.66471118539444e+26*x + 6.48032005181657e+25, x <= 0.772413793103448),
(1.5e+22 , x <= 0.85172414),
(3.85797403183421e+25*x**2 - 6.69652358427785e+25*x + 2.90638521562008e+25, x <= 0.85948275862069),
- (2.2385551110116e+25*x**2 - 3.91279830141555e+25*x + 1.71010327294176e+25, x <= 0.864655172413793),
+ (2.2385551110116e+25*x**2 - 3.91279830141555e+25*x + 1.71010327294176e+25 , x <= 0.864655172413793),
(1.38463307937214e+25*x**2 - 2.43610209842523e+25*x + 1.07168676794207e+25, x <= 0.872413793103448),
(2.5e+21, x <= 0.93793103),
(1.95514620595374e+25*x**2 - 3.65353225743697e+25*x + 1.70704057747143e+25, x <= 0.953448275862069),