added no surface shear heating

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),