Exawind · sayerhs · Jun 11, 2018 · May 31, 2018 · Jun 6, 2018 · Jun 11, 2018
diff --git a/include/AlgTraits.h b/include/AlgTraits.h
@@ -45,7 +45,7 @@ struct AlgTraitsTet4 {
 struct AlgTraitsPyr5 {
   static constexpr int nDim_ = 3;
   static constexpr int nodesPerElement_ = 5;
-  static constexpr int numScsIp_ = 8;
+  static constexpr int numScsIp_ = 12;
   static constexpr int numScvIp_ = 5;
   static constexpr int numGp_ = 5; // for FEM (not supported)
   static constexpr stk::topology::topology_t topo_ = stk::topology::PYRAMID_5;

diff --git a/include/master_element/Hex8CVFEM.h b/include/master_element/Hex8CVFEM.h
@@ -177,6 +177,8 @@ class HexSCS : public MasterElement
 
   const int * adjacentNodes();
 
+  const int * scsIpEdgeOrd() override;
+
   void shape_fcn(
     double *shpfc);
 

diff --git a/include/master_element/Hex8GeometryFunctions.h b/include/master_element/Hex8GeometryFunctions.h
@@ -162,6 +162,104 @@ namespace nalu {
     return volume;
   }
 
+  template <typename RealType>
+  RealType bhex_volume_grandy(RealType scvcoords[8][3])
+  {
+    /**
+     * The Grandy algorithm for computing the volume of a multilinear box
+     *
+     * "Efficient computation of volume ofl
+     * Hexahedral Cells", Jeffrey Grandy, LLNL, UCRL-ID-128886,
+     *  October 30, 1997.
+     * modified for non-planar bent top face which we have broken into two triangles
+     */
+    constexpr int nTri = 24;
+    constexpr int dim = 3;
+
+    constexpr int nNodes = 8;
+    constexpr int nFaces = 6;
+    constexpr int npv = nNodes + nFaces;
+
+    RealType coordv[npv][dim];
+
+    // copy coordinates
+    for (int n = 0; n < nNodes; ++n) {
+      coordv[n][0] = scvcoords[n][0];
+      coordv[n][1] = scvcoords[n][1];
+      coordv[n][2] = scvcoords[n][2];
+    }
+
+    // now we add the face midpoints
+    for (int k = 0; k < 3; ++k ) {
+      coordv[8][k] = 0.25*( scvcoords[0][k] + scvcoords[1][k]
+                            + scvcoords[2][k] + scvcoords[3][k] );
+    }
+
+    for (int k = 0; k < 3; ++k ) {
+      coordv[9][k] = 0.5*( scvcoords[5][k] + scvcoords[7][k] );
+    }
+
+    for (int k = 0; k < 3; ++k ) {
+      coordv[10][k] = 0.25*( scvcoords[0][k] + scvcoords[1][k]
+                             + scvcoords[5][k] + scvcoords[4][k] );
+    }   
+
+    for (int k = 0; k < 3; ++k ) {
+      coordv[11][k] = 0.25*( scvcoords[3][k] + scvcoords[2][k]
+                             + scvcoords[6][k] + scvcoords[7][k] );
+    }
+
+    for (int k = 0; k < 3; ++k ) {
+      coordv[12][k] = 0.25*( scvcoords[1][k] + scvcoords[2][k]
+                             + scvcoords[6][k] + scvcoords[5][k] );
+    }
+
+    for (int k = 0; k < 3; ++k ) {
+      coordv[13][k] = 0.25*( scvcoords[0][k] + scvcoords[3][k]
+                             + scvcoords[7][k] + scvcoords[4][k] );
+    }
+
+    constexpr int triangular_facets[nTri][3] = {
+      { 0,  8,  1}, { 8,  2,  1}, { 3,  2,  8},
+      { 3,  8,  0}, { 6,  9,  5}, { 7,  9,  6},
+      { 4,  9,  7}, { 4,  5,  9}, {10,  0,  1},
+      { 5, 10,  1}, { 4, 10,  5}, { 4,  0, 10},
+      { 7,  6, 11}, { 6,  2, 11}, { 2,  3, 11},
+      { 3,  7, 11}, { 6, 12,  2}, { 5, 12,  6},
+      { 5,  1, 12}, { 1,  2, 12}, { 0,  4, 13},
+      { 4,  7, 13}, { 7,  3, 13}, { 3,  0, 13}
+    };
+
+    RealType volume = 0.0;
+    for (int k = 0; k < nTri; ++k) {
+      const int p = triangular_facets[k][0];
+      const int q = triangular_facets[k][1];
+      const int r = triangular_facets[k][2];
+
+      const RealType triFaceMid[3] = {
+          coordv[p][0] + coordv[q][0] + coordv[r][0],
+          coordv[p][1] + coordv[q][1] + coordv[r][1],
+          coordv[p][2] + coordv[q][2] + coordv[r][2]
+      };
+
+      enum {XC = 0, YC = 1, ZC = 2};
+      RealType dxv[3];
+
+      dxv[0] = ( coordv[q][YC] - coordv[p][YC] ) * ( coordv[r][ZC] - coordv[p][ZC] )
+             - ( coordv[r][YC] - coordv[p][YC] ) * ( coordv[q][ZC] - coordv[p][ZC] );
+
+      dxv[1] = ( coordv[r][XC] - coordv[p][XC] ) * ( coordv[q][ZC] - coordv[p][ZC] )
+             - ( coordv[q][XC] - coordv[p][XC] ) * ( coordv[r][ZC] - coordv[p][ZC] );
+
+      dxv[2] = ( coordv[q][XC] - coordv[p][XC] ) * ( coordv[r][YC] - coordv[p][YC] )
+             - ( coordv[r][XC] - coordv[p][XC] ) * ( coordv[q][YC] - coordv[p][YC] );
+
+      volume += triFaceMid[0] * dxv[0] + triFaceMid[1] * dxv[1] + triFaceMid[2] * dxv[2];
+    }
+    volume /= RealType(18.0);
+    return volume;
+  }
+
   template <typename CoordViewType>
   void subdivide_hex_8(CoordViewType coords, typename CoordViewType::value_type coordv[27][3])
   {

diff --git a/include/master_element/MasterElement.h b/include/master_element/MasterElement.h
@@ -176,7 +176,11 @@ class MasterElement
      throw std::runtime_error("shifted_face_grad_op not implemented");}
 
   virtual const int * adjacentNodes() {
-    throw std::runtime_error("adjacentNodes not implementedunknown bc");
+    throw std::runtime_error("adjacentNodes not implemented");
+    return NULL;}
+
+  virtual const int * scsIpEdgeOrd() {
+    throw std::runtime_error("scsIpEdgeOrd not implemented");
     return NULL;}
 
   virtual const int * ipNodeMap(int ordinal = 0) {
@@ -193,7 +197,7 @@ class MasterElement
 
   virtual int opposingNodes(
     const int ordinal, const int node) {
-    throw std::runtime_error("adjacentNodes not implemented"); }
+    throw std::runtime_error("opposingNodes not implemented"); }
 
   virtual int opposingFace(
     const int ordinal, const int node) {
@@ -266,6 +270,7 @@ class MasterElement
   std::vector<double> nodeLoc_;
   std::vector<int> sideNodeOrdinals_;
   std::vector<int> sideOffset_;
+  std::vector<int> scsIpEdgeOrd_;
 
   // FEM
   std::vector<double>weights_;

diff --git a/include/master_element/Pyr5CVFEM.h b/include/master_element/Pyr5CVFEM.h
@@ -76,6 +76,11 @@ class PyrSCV : public MasterElement
     const int &npts,
     const double *par_coord, 
     double* shape_fcn);
+
+  void shifted_pyr_shape_fcn(
+    const int &npts,
+    const double *par_coord, 
+    double* shape_fcn);
 };
 
 // Pyramid 5 subcontrol surface
@@ -130,6 +135,11 @@ class PyrSCS : public MasterElement
     const double *intLoc,
     double *deriv);
 
+  void shifted_pyr_derivative(
+    const int npts,
+    const double *intLoc,
+    double *deriv);
+
   void gij( 
     SharedMemView<DoubleType**>& coords,
     SharedMemView<DoubleType***>& gupper,
@@ -144,6 +154,8 @@ class PyrSCS : public MasterElement
 
   const int * adjacentNodes();
 
+  const int * scsIpEdgeOrd() override;
+
   void shape_fcn(
     double *shpfc);
 
@@ -155,6 +167,11 @@ class PyrSCS : public MasterElement
     const double *par_coord, 
     double* shape_fcn);
 
+  void shifted_pyr_shape_fcn(
+    const int &npts,
+    const double *par_coord, 
+    double* shape_fcn);
+
   void
   general_shape_fcn(const int numIp, const double* isoParCoord, double* shpfc)
   {
@@ -222,13 +239,6 @@ class PyrSCS : public MasterElement
     const double *field,
     double *result);
 
-private:
-  using QuadFaceGradType = SharedMemView<DoubleType***>;
-  using TriFaceGradType = SharedMemView<DoubleType***>;
-
-  void face_grad_op(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, QuadFaceGradType& gradop);
-  void face_grad_op_quad(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, QuadFaceGradType& gradop);
-  void face_grad_op_tri(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, TriFaceGradType& gradop);
 };
 
 } // namespace nalu

diff --git a/include/master_element/Quad42DCVFEM.h b/include/master_element/Quad42DCVFEM.h
@@ -154,6 +154,8 @@ class Quad42DSCS : public MasterElement
 
   const int * adjacentNodes() override;
 
+  const int * scsIpEdgeOrd() override;
+
   int opposingNodes(
     const int ordinal, const int node) override;
 

diff --git a/include/master_element/Tet4CVFEM.h b/include/master_element/Tet4CVFEM.h
@@ -142,6 +142,8 @@ class TetSCS : public MasterElement
 
   const int * adjacentNodes();
 
+  const int * scsIpEdgeOrd() override;
+
   void shape_fcn(
     double *shpfc);
 

diff --git a/include/master_element/Tri32DCVFEM.h b/include/master_element/Tri32DCVFEM.h
@@ -153,6 +153,8 @@ class Tri32DSCS : public MasterElement
 
   const int * adjacentNodes() override;
 
+  const int * scsIpEdgeOrd() override;
+
   void shape_fcn(
     double *shpfc) override;
 

diff --git a/include/master_element/Wed6CVFEM.h b/include/master_element/Wed6CVFEM.h
@@ -143,6 +143,8 @@ class WedSCS : public MasterElement
     double *deriv);
 
   const int * adjacentNodes();
+
+  const int * scsIpEdgeOrd() override;
 
   int opposingNodes(
     const int ordinal, const int node);
@@ -197,15 +199,6 @@ class WedSCS : public MasterElement
   double parametric_distance( const std::vector<double> &x);
 
   const int* side_node_ordinals(int sideOrdinal) final;
-
-private:
-  using QuadFaceGradType = SharedMemView<DoubleType***>;
-  using TriFaceGradType = SharedMemView<DoubleType***>;
-
-  void face_grad_op(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, TriFaceGradType& gradop);
-  void face_grad_op_tri(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, TriFaceGradType& gradop);
-  void face_grad_op_quad(const int face_ordinal, const bool shifted, SharedMemView<DoubleType**>& coords, QuadFaceGradType& gradop);
-
 };
 
 

diff --git a/reg_tests/test_files/cvfemHC/cvfemHC.norm.gold b/reg_tests/test_files/cvfemHC/cvfemHC.norm.gold
@@ -1,10 +1,10 @@
-0.2628771158640655 1 0.002
-0.3740009686945356 2 0.004
-0.4158920415683592 3 0.006
-0.4266331544604634 4 0.008
-0.4230444195793962 5 0.01
-0.4128160274981106 6 0.012
-0.3996919647873065 7 0.014
-0.3855662411850521 8 0.016
-0.3713926724760221 9 0.018
-0.3576293361290214 10 0.02
+0.2627997981379542 1 0.002
+0.3739570323910922 2 0.004
+0.4158868201964299 3 0.006
+0.4266528240801 4 0.008
+0.4230776743731539 5 0.01
+0.4128556496220737 6 0.012
+0.3997337729424268 7 0.014
+0.3856080447944304 8 0.016
+0.3714334464751312 9 0.018
+0.3576686726076596 10 0.02
diff --git a/reg_tests/test_files/dgNonConformalEdgeCylinder/dgNonConformalEdgeCylinder.norm.gold b/reg_tests/test_files/dgNonConformalEdgeCylinder/dgNonConformalEdgeCylinder.norm.gold
@@ -1,10 +1,10 @@
-0.0257392190579548 1 0.002
-0.02574740104224447 2 0.004
-0.007430198932416332 3 0.006
-0.002244591729982638 4 0.008
-0.001284784350631387 5 0.01
-0.0009345657653626473 6 0.012
-0.0007968521337429102 7 0.014
-0.0007702662515903432 8 0.016
-0.0008062636390563359 9 0.018
-0.0008755843133531737 10 0.02
+0.02573816754026069 1 0.002
+0.02574926079128527 2 0.004
+0.007427623748310587 3 0.006
+0.002245139667198144 4 0.008
+0.001284820897687786 5 0.01
+0.0009345793445916078 6 0.012
+0.0007968955767543434 7 0.014
+0.0007703254790273972 8 0.016
+0.0008063558961797211 9 0.018
+0.0008757271998318686 10 0.02
diff --git a/reg_tests/test_files/dgNonConformalElemCylinder/dgNonConformalElemCylinder.norm.gold b/reg_tests/test_files/dgNonConformalElemCylinder/dgNonConformalElemCylinder.norm.gold
@@ -1,10 +1,10 @@
-0.01664684637424995 1 0.002
-0.01242786016167628 2 0.004
-0.004467175260558592 3 0.006
-0.001785630924956796 4 0.008
-0.001151317394367551 5 0.01
-0.0008223057256807056 6 0.012
-0.0006522880498152978 7 0.014
-0.00058519172715382 8 0.016
-0.0005786660921017782 9 0.018
-0.0006036230765110883 10 0.02
+0.01664646563202813 1 0.002
+0.01242849768020153 2 0.004
+0.004467313637199031 3 0.006
+0.001785640455910638 4 0.008
+0.001151299749497504 5 0.01
+0.0008223224078619498 6 0.012
+0.0006523164571680431 7 0.014
+0.0005852609531605202 8 0.016
+0.0005787519331727764 9 0.018
+0.0006037167241515856 10 0.02
diff --git a/reg_tests/test_files/edgeHybridFluids/edgeHybridFluids.norm.gold b/reg_tests/test_files/edgeHybridFluids/edgeHybridFluids.norm.gold
@@ -1,30 +1,30 @@
-0.002433239019241183 1 0.00833333
-0.00182187455877493 2 0.0152778
-0.001372110346070923 3 0.0210648
-0.001117003702144699 4 0.0258873
-0.0009775100620998849 5 0.0299061
-0.0009146440951029971 6 0.0332551
-0.0008842004618467895 7 0.0360794
-0.0008781867243786976 8 0.0389037
-0.0008685511301787267 9 0.041728
-0.0008550866333398277 10 0.0445523
-0.0008532934688223621 11 0.0473766
-0.0008605125228310587 12 0.0502009
-0.0008675541741859722 13 0.0530252
-0.0008765929330727491 14 0.0558495
-0.0008850125572882093 15 0.0586738
-0.0008941456737045214 16 0.0614981
-0.0009026417083060649 17 0.0643224
-0.0009153427466626436 18 0.0671467
-0.0009252790711106383 19 0.069971
-0.000943545567097144 20 0.0727952
-0.0009738949867587277 21 0.0756195
-0.001011128892709552 22 0.0784438
-0.001039891540612789 23 0.0812681
-0.001053094543777811 24 0.0840924
-0.001059991393360899 25 0.0869167
-0.00106232513551189 26 0.089741
-0.001060662818711306 27 0.0925653
-0.001058973881926617 28 0.0953896
-0.001061263398261228 29 0.0982139
-0.001064396874900402 30 0.101038
+0.002433225648733176 1 0.00833333
+0.001821873041579242 2 0.0152778
+0.001372111221006272 3 0.0210648
+0.001117004521992668 4 0.0258873
+0.0009775105131627547 5 0.0299061
+0.0009146442997356928 6 0.0332551
+0.0008842005406580971 7 0.0360794
+0.0008781867702797445 8 0.0389037
+0.0008685511773113951 9 0.041728
+0.0008550866698412334 10 0.0445523
+0.0008532934851332336 11 0.0473766
+0.0008605125558776005 12 0.0502009
+0.0008675542422153722 13 0.0530252
+0.000876593017054678 14 0.0558495
+0.0008850126091891584 15 0.0586738
+0.0008941457181818259 16 0.0614981
+0.0009026418008668091 17 0.0643224
+0.0009153428204842487 18 0.0671467
+0.0009252791061810112 19 0.069971
+0.0009435456228732972 20 0.0727952
+0.0009738950630714776 21 0.0756195
+0.001011128960886205 22 0.0784438
+0.001039891590172271 23 0.0812681
+0.001053094597904249 24 0.0840924
+0.001059991453247039 25 0.0869167
+0.001062325189813977 26 0.089741
+0.001060662853125896 27 0.0925653
+0.001058973923059289 28 0.0953896
+0.001061263439784936 29 0.0982139
+0.001064396924961079 30 0.101038