Added generator for P2 epsilon operator with blending

data/codegen/forms/p2/P2Form_eps.py

+#!/usr/bin/env python3
+
+import numpy as np
+from sympy import *
+import p2_reference as p2ref
+import quadpy
+
+eta, xi, nu = symbols('eta xi nu')
+
+x1, x2, x3, x4 = symbols('coords[0][0] coords[1][0] coords[2][0] coords[3][0]', real=True)
+y1, y2, y3, y4 = symbols('coords[0][1] coords[1][1] coords[2][1] coords[3][1]', real=True)
+z1, z2, z3, z4 = symbols('coords[0][2] coords[1][2] coords[2][2] coords[3][2]', real=True)
+
+DF_11, DF_12, DF_13, DF_21, DF_22, DF_23, DF_31, DF_32, DF_33 = symbols('DF_11, DF_12, DF_13, DF_21, DF_22, DF_23, DF_31, DF_32, DF_33', real=True)
+DF_11.name = 'DF(0,0)'
+DF_12.name = 'DF(0,1)'
+DF_13.name = 'DF(0,2)'
+DF_21.name = 'DF(1,0)'
+DF_22.name = 'DF(1,1)'
+DF_23.name = 'DF(1,2)'
+DF_31.name = 'DF(2,0)'
+DF_32.name = 'DF(2,1)'
+DF_33.name = 'DF(2,2)'
+
+x_hat, y_hat, z_hat = symbols('x_hat[0], x_hat[1], x_hat[2]', real=True)
+x_tilde, y_tilde, z_tilde = symbols('x_tilde[0], x_tilde[1], x_tilde[2]', real=True)
+
+K = symbols('K', real=True)
+K.name = 'callback( x_tilde, geometryMap_ )'
+
+N_2d = 6
+N_3d = 10
+
+def generate_2d_integration_rule():
+    scheme_tri = quadpy.triangle.strang_fix_cowper_02()
+    reference_tri = np.asarray([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
+    vol = quadpy.nsimplex.get_vol(reference_tri)
+    x_q_tri = quadpy.nsimplex.transform(scheme_tri.points.T, reference_tri.T).T
+    w_tri = vol * scheme_tri.weights
+    return (w_tri, x_q_tri)
+
+def generate_2d_forms():
+    dim = 2
+    R = simplify(Matrix([[x2-x1, x3-x1], [y2-y1, y3-y1]]) * Matrix([[eta], [xi]]) + Matrix([[x1], [y1]]))
+    DR = R.jacobian([eta, xi])
+    DRinv = DR**(-1)
+
+    phi_hat, _ = p2ref.create_2d()
+
+    DF = Matrix([[DF_11, DF_12], [DF_21, DF_22]])
+    DFinv = DF**(-1)
+
+    dx_tilde = simplify(abs(det(DR)))
+    dx = simplify(abs(det(DF)) * dx_tilde)
+
+    def grad(u):
+      return Matrix([[diff(u, eta)], [diff(u, xi)]])
+
+    def vecgrad(u):
+      grad_0 = DFinv.T * DRinv.T * grad(u[0])
+      grad_1 = DFinv.T * DRinv.T * grad(u[1])
+
+      vgrad = Matrix([grad_0.T, grad_1.T])
+      vgrad = simplify(vgrad)
+      return vgrad
+
+    def epsilon(u):
+      symgrad = 0.5 * (vecgrad(u) + vecgrad(u).T)
+      symgrad = simplify(symgrad)
+      return symgrad
+
+    def inner(u, v):
+      sum_ = 0
+      for i in range(dim):
+        for j in range(dim):
+          sum_ += u[(i,j)] * v[(i,j)]
+      return sum_
+
+    forms_2d = []
+    for block_i in range(dim):
+      block_row = []
+      for block_j in range(dim):
+        local_forms = []
+        for i in range(N_2d):
+          for j in range(N_2d):
+            print('Processing 2D form block {},{} element {},{}'.format(block_i, block_j, i, j))
+
+            unitvec_i = Matrix([[0],[0]])
+            unitvec_j = Matrix([[0],[0]])
+
+            unitvec_i[block_i] = 1
+            unitvec_j[block_j] = 1
+
+            form = 2 * K * inner(epsilon(phi_hat[i] * unitvec_i), epsilon(phi_hat[j] * unitvec_j)) * dx
+            form = form.subs([(eta, x_hat), (xi, y_hat)])
+            # form = simplify(form)
+            local_forms.append(form)
+
+        block_row.append(local_forms)    
+      forms_2d.append(block_row)
+
+    return forms_2d, R
+
+def generate_3d_integration_rule():
+    scheme_tet = quadpy.tetrahedron.yu_1()
+    reference_tet = np.asarray([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
+    vol = quadpy.nsimplex.get_vol(reference_tet)
+    x_q_tet = quadpy.nsimplex.transform(scheme_tet.points.T, reference_tet.T).T
+    w_tet = vol * scheme_tet.weights
+    return (w_tet, x_q_tet)
+
+def generate_3d_forms():
+    dim = 3
+    R = simplify(Matrix([[x2-x1, x3-x1, x4-x1], [y2-y1, y3-y1, y4-y1], [z2-z1, z3-z1, z4-z1]]) * Matrix([[eta], [xi], [nu]]) + Matrix([[x1], [y1], [z1]]))
+    DR = R.jacobian([eta, xi, nu])
+    DRinv = DR**(-1)
+
+    print('generating 3d basis functions')
+    phi_hat, _ = p2ref.create_3d()
+
+    DF = Matrix([[DF_11, DF_12, DF_13], [DF_21, DF_22, DF_23], [DF_31, DF_32, DF_33]])
+    DFinv = DF**(-1)
+
+    dx_tilde = simplify(abs(det(DR)))
+    dx = simplify(abs(det(DF)) * dx_tilde)
+
+    def grad(u):
+      return Matrix([[diff(u, eta)], [diff(u, xi)], [diff(u, nu)]])
+
+    def vecgrad(u):
+      print('grad0')
+      grad_0 = (DFinv.T * DRinv.T * grad(u[0]))
+      print('grad1')
+      grad_1 = (DFinv.T * DRinv.T * grad(u[1]))
+      print('grad2')
+      grad_2 = (DFinv.T * DRinv.T * grad(u[2]))
+
+      vgrad = Matrix([grad_0.T, grad_1.T, grad_2.T])
+      # vgrad = simplify(vgrad)
+      return vgrad
+
+    def epsilon(u):
+      symgrad = 0.5 * (vecgrad(u) + vecgrad(u).T)
+      print('symgrad')
+      # symgrad = simplify(symgrad)
+      return symgrad
+
+    def inner(u, v):
+      sum_ = 0
+      for i in range(dim):
+        for j in range(dim):
+          sum_ += u[(i,j)] * v[(i,j)]
+      return sum_
+
+    forms_3d = []
+    for block_i in range(dim):
+      block_row = []
+      for block_j in range(dim):
+        local_forms = []
+        for i in range(N_3d):
+          for j in range(N_3d):
+            print('Processing 3D form block {},{} element {},{}'.format(block_i, block_j, i, j))
+
+            unitvec_i = Matrix([[0],[0],[0]])
+            unitvec_j = Matrix([[0],[0],[0]])
+
+            unitvec_i[block_i] = 1
+            unitvec_j[block_j] = 1
+
+            form = 2 * K * inner(epsilon(phi_hat[i] * unitvec_i), epsilon(phi_hat[j] * unitvec_j)) * dx
+            form = form.subs([(eta, x_hat), (xi, y_hat), (nu, z_hat)])
+            # form = simplify(form)
+            local_forms.append(form)
+
+        block_row.append(local_forms)    
+      forms_3d.append(block_row)
+
+    return forms_3d, R
+
+forms_2d, R_2d = generate_2d_forms()
+forms_3d, R_3d = generate_3d_forms()
+
+def sympyToC(classname, block_i, block_j):
+    c_code = "class {} : public P2FormHyTeG {{\n".format(classname)
+    c_code += "public:\n"
+
+    ### 2D
+
+    tmpsyms = numbered_symbols("tmp")
+
+    if block_i < 2 and block_j < 2:
+      symbols, simple = cse(forms_2d[block_i][block_j], symbols=tmpsyms)
+      w, x_q = generate_2d_integration_rule()
+
+    c_code += "  void evalQuadraturePoint2D(const Point3D& x_hat, const Point3D& x_tilde, const std::array<Point3D,3>& coords, const Matrix2r& DF, real_t w, Matrix6r& elMat) const\n"
+    c_code += "  {\n"
+    if block_i < 2 and block_j < 2:
+      for s in symbols:
+          #print s
+          c_code +=  "    real_t " +cxxcode(s[0]) + " = " + cxxcode(s[1]) + ";\n"
+
+      for i in range(N_2d):
+          for j in range(N_2d):
+              c_code +=  "    elMat({},{}) += w * ".format(i,j) + cxxcode(simple[N_2d*i + j])+";\n"
+    else:
+      c_code +=  "    WALBERLA_ABORT(\"Not available in 2D\");\n"
+    c_code += "  }\n\n"
+    c_code += "  void integrateAll( const std::array< Point3D, 3 >& coords, Matrix6r& elMat ) const final\n"
+    c_code += "  {\n"
+
+    if block_i < 2 and block_j < 2:
+      c_code += "    Point3D x_hat;\n"
+      c_code += "    Point3D x_tilde;\n"
+      c_code += "    Matrix2r DF;\n"
+      for i in range(N_2d):
+          for j in range(N_2d):
+              c_code += "    elMat({},{}) = 0;\n".format(i,j)
+
+      for q in range(len(w)):
+          for j in range(2):
+              c_code += "    x_hat[{}] = {};\n".format(j, x_q[q][j])
+
+          x_tilde = R_2d.subs([(eta, x_q[q][0]), (xi, x_q[q][1])])
+          for j in range(2):
+              c_code += "    x_tilde[{}] = {};\n".format(j, x_tilde[j])
+          c_code += "    geometryMap_->evalDF(x_tilde, DF);\n"
+          c_code += "    evalQuadraturePoint2D(x_hat, x_tilde, coords, DF, {}, elMat);\n".format(w[q])
+    else:
+      c_code +=  "    WALBERLA_ABORT(\"Not available in 2D\");\n"
+
+    c_code += "  }\n\n"
+
+    ### 3D
+
+    tmpsyms = numbered_symbols("tmp")
+    symbols, simple = cse(forms_3d[block_i][block_j], symbols=tmpsyms)
+    w, x_q = generate_3d_integration_rule()
+
+    c_code += "  void evalQuadraturePoint3D(const Point3D& x_hat, const Point3D& x_tilde, const std::array<Point3D,4>& coords, const Matrix3r& DF, real_t w, Matrix10r& elMat) const\n"
+    c_code += "  {\n"
+    for s in symbols:
+        c_code +=  "    real_t " +cxxcode(s[0]) + " = " + cxxcode(s[1]) + ";\n"
+
+    for i in range(N_3d):
+        for j in range(N_3d):
+            c_code +=  "    elMat({},{}) += w * ".format(i,j) + cxxcode(simple[N_3d*i + j])+";\n"
+    c_code += "  }\n\n"
+    c_code += "  void integrateAll( const std::array< Point3D, 4 >& coords, Matrix10r& elMat ) const final\n"
+    c_code += "  {\n"
+
+    c_code += "    Point3D x_hat;\n"
+    c_code += "    Point3D x_tilde;\n"
+    c_code += "    Matrix3r DF;\n"
+    for i in range(N_3d):
+        for j in range(N_3d):
+            c_code += "    elMat({},{}) = 0;\n".format(i,j)
+
+    for q in range(len(w)):
+        for j in range(3):
+            c_code += "    x_hat[{}] = {};\n".format(j, x_q[q][j])
+
+        x_tilde = R_3d.subs([(eta, x_q[q][0]), (xi, x_q[q][1]), (nu, x_q[q][2])])
+        for j in range(3):
+            c_code += "    x_tilde[{}] = {};\n".format(j, x_tilde[j])
+        c_code += "    geometryMap_->evalDF(x_tilde, DF);\n"
+        c_code += "    evalQuadraturePoint3D(x_hat, x_tilde, coords, DF, {}, elMat);\n".format(w[q])
+
+    # c_code += "    }\n"
+    c_code += "  }\n\n"
+
+    c_code += "  static std::function< real_t( const Point3D&, const std::shared_ptr< GeometryMap >& ) > callback;\n\n"
+    c_code += "};\n\n"
+    return c_code
+
+c_code = "#pragma once\n\n"
+c_code += "#include \"hyteg/forms/form_hyteg_base/P2FormHyTeG.hpp\"\n\n"
+c_code += "using walberla::real_c;\n\n"
+c_code += "namespace hyteg {\n\n"
+
+
+for block_i in range(3):
+  for block_j in range(3):
+    c_code += sympyToC('P2Form_epsilon_{}{}'.format(block_i+1, block_j+1), block_i, block_j)
+
+c_code += "}\n"
+
+print(c_code)

src/hyteg/elementwiseoperators/P2ElementwiseOperator.cpp

@@ -908,4 +908,15 @@ template class P2ElementwiseOperator< P2Form_laplace >;
 // P2ElementwiseLinearCombinationOperator
 template class P2ElementwiseOperator< P2LinearCombinationForm >;
 
+// P2EpsilonOperator
+template class P2ElementwiseOperator< P2Form_epsilon_11 >;
+template class P2ElementwiseOperator< P2Form_epsilon_12 >;
+template class P2ElementwiseOperator< P2Form_epsilon_13 >;
+template class P2ElementwiseOperator< P2Form_epsilon_21 >;
+template class P2ElementwiseOperator< P2Form_epsilon_22 >;
+template class P2ElementwiseOperator< P2Form_epsilon_23 >;
+template class P2ElementwiseOperator< P2Form_epsilon_31 >;
+template class P2ElementwiseOperator< P2Form_epsilon_32 >;
+template class P2ElementwiseOperator< P2Form_epsilon_33 >;
+
 } // namespace hyteg

src/hyteg/elementwiseoperators/P2ElementwiseOperator.hpp

@@ -24,6 +24,7 @@
 #include "hyteg/forms/form_fenics_base/P2FenicsForm.hpp"
 #include "hyteg/forms/form_fenics_generated/p2_polar_laplacian.h"
 #include "hyteg/forms/form_hyteg_generated/P2FormDivKGradBlending.hpp"
+#include "hyteg/forms/form_hyteg_generated/P2FormEpsilon.hpp"
 #include "hyteg/forms/form_hyteg_manual/P2FormDivKGrad.hpp"
 #include "hyteg/forms/form_hyteg_manual/P2FormLaplace.hpp"
 #include "hyteg/forms/form_hyteg_manual/P2FormMass.hpp"
@@ -250,4 +251,14 @@ typedef P2ElementwiseOperator< P2Form_divKgradBlending > P2ElementwiseDivKGradBl
 
 typedef P2ElementwiseOperator< P2LinearCombinationForm > P2ElementwiseLinearCombinationOperator;
 
+typedef P2ElementwiseOperator< P2Form_epsilon_11 > P2EpsilonOperator_11;
+typedef P2ElementwiseOperator< P2Form_epsilon_12 > P2EpsilonOperator_12;
+typedef P2ElementwiseOperator< P2Form_epsilon_13 > P2EpsilonOperator_13;
+typedef P2ElementwiseOperator< P2Form_epsilon_21 > P2EpsilonOperator_21;
+typedef P2ElementwiseOperator< P2Form_epsilon_22 > P2EpsilonOperator_22;
+typedef P2ElementwiseOperator< P2Form_epsilon_23 > P2EpsilonOperator_23;
+typedef P2ElementwiseOperator< P2Form_epsilon_31 > P2EpsilonOperator_31;
+typedef P2ElementwiseOperator< P2Form_epsilon_32 > P2EpsilonOperator_32;
+typedef P2ElementwiseOperator< P2Form_epsilon_33 > P2EpsilonOperator_33;
+
 } // namespace hyteg