dg : added 3D prisms, supports prism/tetrahedron mixed mesh

Geo/MElement.h

@@ -237,7 +237,7 @@ class MElement
   // integration routines
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   {
-    Msg::Error("No integration points defined for this type of element");
+    Msg::Error("No integration points defined for this type of element: %d",this->getType());
   }
 
   // IO routines

Geo/MPrism.cpp

@@ -20,3 +20,91 @@ int MPrism::getVolumeSign()
   mat[2][2] = _v[3]->z() - _v[0]->z();
   return sign(det3x3(mat));
 }
+
+void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+  *npts = getNGQPriPts(pOrder);
+  *pts = getGQPriPts(pOrder);
+}
+
+const polynomialBasis* MPrism::getFunctionSpace(int o) const
+{
+  int order = (o == -1) ? getPolynomialOrder() : o;
+
+  int nv = getNumVolumeVertices();
+  
+  if ((nv == 0) && (o == -1)) {
+    switch (order) {
+    case 1: return &polynomialBases::find(MSH_PRI_6);
+    case 2: return &polynomialBases::find(MSH_PRI_18);
+    default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+    }
+  }
+  else { 
+    switch (order) {
+    case 1: return &polynomialBases::find(MSH_PRI_6);
+    case 2: return &polynomialBases::find(MSH_PRI_18);
+    default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+    }
+  }
+  return 0;
+}
+
+void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) const
+{
+  for (ithFace = 0; ithFace < 5; ithFace++){
+    MVertex *v0 = _v[faces_prism(ithFace, 0)];
+    MVertex *v1 = _v[faces_prism(ithFace, 1)];
+    MVertex *v2 = _v[faces_prism(ithFace, 2)];
+
+    if (face.getNumVertices()==3) {
+      if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2)){
+        sign = 1; rot = 0; return;
+      }
+      if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(0)){
+        sign = 1; rot = 1; return;
+      }
+      if (v0 == face.getVertex(2) && v1 == face.getVertex(0) && v2 == face.getVertex(1)){
+        sign = 1; rot = 2; return;
+      }
+      if (v0 == face.getVertex(0) && v1 == face.getVertex(2) && v2 == face.getVertex(1)){
+        sign = -1; rot = 0; return;
+      }
+      if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(2)){
+        sign = -1; rot = 1; return;
+      }
+      if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0)){
+        sign = -1; rot = 2; return;
+      }
+    }
+    else {
+      MVertex *v3 = _v[faces_prism(ithFace, 3)];
+      if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && v2 == face.getVertex(2) && v3 == face.getVertex(3)){
+        sign = 1; rot = 0; return;
+      }
+      if (v0 == face.getVertex(1) && v1 == face.getVertex(2) && v2 == face.getVertex(3) && v3 == face.getVertex(0)){
+        sign = 1; rot = 1; return;
+      }
+      if (v0 == face.getVertex(2) && v1 == face.getVertex(3) && v2 == face.getVertex(0) && v3 == face.getVertex(1)){
+        sign = 1; rot = 2; return;
+      }
+      if (v0 == face.getVertex(3) && v1 == face.getVertex(0) && v2 == face.getVertex(1) && v3 == face.getVertex(2)){
+        sign = 1; rot = 3; return;
+      }
+      if (v0 == face.getVertex(0) && v1 == face.getVertex(3) && v2 == face.getVertex(2) && v3 == face.getVertex(1)){
+        sign = -1; rot = 0; return;
+      }
+      if (v0 == face.getVertex(1) && v1 == face.getVertex(0) && v2 == face.getVertex(3) && v3 == face.getVertex(2)){
+        sign = -1; rot = 1; return;
+      }
+      if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && v2 == face.getVertex(0) && v3 == face.getVertex(3)){
+        sign = -1; rot = 2; return;
+      }
+      if (v0 == face.getVertex(3) && v1 == face.getVertex(2) && v2 == face.getVertex(1) && v3 == face.getVertex(0)){
+        sign = -1; rot = 3; return;
+      }
+    }
+      
+  }
+  Msg::Error("Could not get face information for prism %d", getNum());
+}

Geo/MPrism.h

@@ -86,6 +86,7 @@ class MPrism : public MElement {
     _getEdgeVertices(num, v);
   }
   virtual int getNumFaces(){ return 5; }
+  virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign, int &rot) const; 
   virtual MFace getFace(int num)
   {
     if(num < 2)
@@ -128,6 +129,7 @@ class MPrism : public MElement {
     tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
     tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
   }
+  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual int getVolumeSign();
   virtual void getShapeFunctions(double u, double v, double w, double s[], int o) 
   {
@@ -167,6 +169,7 @@ class MPrism : public MElement {
       return false;
     return true;
   }
+  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
  private:
   int edges_prism(const int edge, const int vert) const
   {

Numeric/CMakeLists.txt

@@ -12,6 +12,7 @@ set(SRC
     GaussQuadratureQuad.cpp
     GaussQuadratureTet.cpp
     GaussQuadratureHex.cpp
+    GaussQuadraturePri.cpp
     GaussLegendreSimplex.cpp
   robustPredicates.cpp
   mathEvaluator.cpp

Numeric/Gauss.h

@@ -14,6 +14,7 @@ struct IntPt{
 int GaussLegendreTri(int n1, int n2, IntPt *pts);
 int GaussLegendreTet(int n1, int n2, int n3, IntPt *pts); 
 int GaussLegendreHex(int n1, int n2, int n3, IntPt *pts);
+int GaussLegendrePri(int n1, int n2, int n3, IntPt *pts); 
 
 int getNGQLPts (int order);
 IntPt *getGQLPts (int order);
@@ -27,6 +28,9 @@ IntPt *getGQQPts(int order);
 int getNGQTetPts(int order);
 IntPt *getGQTetPts(int order);
 
+int getNGQPriPts(int order);
+IntPt *getGQPriPts(int order);
+
 int getNGQHPts(int order);
 IntPt *getGQHPts(int order);
 

Numeric/GaussQuadraturePri.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "Gauss.h"
+#include "GaussLegendre1D.h"
+#include "stdio.h"
+
+IntPt *getGQPriPts(int order);
+int getNGQPriPts(int order);
+
+IntPt * GQP[20] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+
+IntPt *getGQPriPts(int order)
+{ 
+  int nLin = (order+3)/2;
+  int nTri = getNGQTPts(order);
+  int n = nLin*nTri;
+  int index = n;
+  if(!GQP[index])
+    {
+      double *linPt,*linWt;
+      IntPt *triPts = getGQTPts(order);
+//            printf("o = %d n = %d nT = %d nL = %d i = %d\n",order,n,nTri,nLin,index);
+      gmshGaussLegendre1D(nLin,&linPt,&linWt);
+      GQP[index] = new IntPt[n];
+      int l = 0;
+      for(int i=0; i < nTri; i++) {
+        for(int j=0; j < nLin; j++) {
+            GQP[index][l].pt[0] = triPts[i].pt[0];
+            GQP[index][l].pt[1] = triPts[i].pt[1];
+            GQP[index][l].pt[2] = linPt[j];
+            GQP[index][l++].weight = triPts[i].weight*linWt[j];
+//                  printf ("%d: %f %f %f %f\n",l-1,triPts[i].pt[0],triPts[i].pt[1],linPt[j],triPts[i].weight*linWt[j]);
+        }
+      }
+    }
+  return GQP[index];
+}
+
+int getNGQPriPts(int order)
+{ 
+  int nLin = (order+3)/2;
+  int nTri = getNGQTPts(order);
+  return nLin*nTri;
+
+//   if(order == 3)return 8;
+//   if(order == 2)return 8;
+//   if(order < 2)
+//     return GQPnPt[order]; 
+//   return ((order+3)/2)*((order+3)/2)*((order+3)/2);
+}
+
+
+
+
+

Numeric/polynomialBasis.cpp

@@ -233,6 +233,32 @@ static fullMatrix<double> generatePascalTetrahedron(int order)
   return monomials;
 }  
 
+// generate Pascal prism
+
+static fullMatrix<double> generatePascalPrism(int order)
+{
+  int nbMonomials = (order + 1) * (order + 1) * (order + 2) / 2;
+
+  fullMatrix<double> monomials(nbMonomials, 3);
+
+  int index = 0;
+  fullMatrix<double> lineMonoms = generate1DMonomials(order);
+  fullMatrix<double> triMonoms = generatePascalTriangle(order);
+//   printf("nb: %d nT: %d %d nL: %d %d\n",nbMonomials,triMonoms.size1(),(order+1)*(order+2)/2,lineMonoms.size1(),order+1);
+  for (int j = 0; j < lineMonoms.size1(); j++) {
+    for (int i = 0; i < triMonoms.size1(); i++) {
+      monomials(index,0) = triMonoms(i,0);
+      monomials(index,1) = triMonoms(i,1);
+      monomials(index,2) = lineMonoms(j,0);
+      index ++;
+//       printf("%d: %f %f %f\n",index,triMonoms(i,0),triMonoms(i,1),lineMonoms(j,0));
+    }
+  }
+//   monomials.print();
+  return monomials;
+}  
+
+
 static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
 //static int nbdoftriangleserendip(int order) { return 3 * order; }
 
@@ -597,6 +623,31 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
   return point;  
 }
 
+static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip) 
+{
+  int nbPoints = (order + 1)*(order + 1)*(order + 2)/2; 
+  fullMatrix<double> point(nbPoints, 3);
+  
+  double overOrder = (order == 0 ? 1. : 1. / order);
+  int index = 0;
+  fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
+  fullMatrix<double> linePoints = generate1DPoints(order);
+//   printf("nb: %d nT: %d %d nL: %d %d\n",nbPoints,triPoints.size1(),(order+1)*(order+2)/2,linePoints.size1(),order+1);
+  for (int j = 0; j < linePoints.size1(); j++) {
+    for (int i = 0; i < triPoints.size1(); i++) {
+      point(index,0) = triPoints(i,0);
+      point(index,1) = triPoints(i,1);
+      point(index,2) = linePoints(j,0);
+      index ++;
+//       printf("%d: %f %f %f\n",index,triPoints(i,0),triPoints(i,1),linePoints(j,0));
+    }
+  }
+//   point.print();
+
+  point.scale(overOrder);  
+  return point;
+}
+
 static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip) 
 {
   int nbPoints = serendip ? order*4 : (order+1)*(order+1);
@@ -744,6 +795,50 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
   }
 }
 
+static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<int> &closure,
+                           int order)
+{
+  if (order > 1)
+    Msg::Error("FaceClosure not implemented for prisms of order %d",order);
+  bool isTriangle = iFace<2;
+  int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
+  closure.clear();
+  closure.resize(nNodes);
+  switch (order){
+  case 0:
+    closure[0] = 0;
+    break;
+  default:
+//     int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
+    int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
+    int nVertex = isTriangle ? 3 : 4;
+    for (int i = 0; i < nVertex; ++i){
+      int k;
+        k = (nVertex + (iSign * i) + iRotate) % nVertex;  //- iSign * iRotate
+      closure[i] = order1node[iFace][k];
+    }
+    break;
+  }
+}
+
+static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int order)
+{
+
+  closure.clear();
+  for (int iRotate = 0; iRotate < 4; iRotate++){
+    for (int iSign = 1; iSign >= -1; iSign -= 2){
+      for (int iFace = 0; iFace < 5; iFace++){
+    std::vector<int> closure_face;
+    getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order); 
+    closure.push_back(closure_face);
+//       for(int i=0;i < closure_face.size(); i++) { printf("%d ",closure_face.at(i)); }
+//       printf("\n");
+      }
+    }
+  }
+}
+
+
 static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
 {
   closure.clear();
@@ -776,6 +871,7 @@ const polynomialBasis &polynomialBases::find(int tag)
   std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
   if (it != fs.end())     return it->second;
   polynomialBasis F;
+  F.numFaces = -1;
   
   switch (tag){
   case MSH_PNT:
@@ -848,36 +944,43 @@ const polynomialBasis &polynomialBases::find(int tag)
     generate2dEdgeClosure(F.edgeClosure, 5);
     break;
   case MSH_TET_4 :
+    F.numFaces = 4;
     F.monomials = generatePascalTetrahedron(1);
     F.points =    gmshGeneratePointsTetrahedron(1, false);
     generate3dFaceClosure(F.faceClosure, 1);
     break;
   case MSH_TET_10 :
+    F.numFaces = 4;
     F.monomials = generatePascalTetrahedron(2);
     F.points =    gmshGeneratePointsTetrahedron(2, false);
     generate3dFaceClosure(F.faceClosure, 2);
     break;
   case MSH_TET_20 :
+    F.numFaces = 4;
     F.monomials = generatePascalTetrahedron(3);
     F.points =    gmshGeneratePointsTetrahedron(3, false);
     generate3dFaceClosure(F.faceClosure, 3);
     break;
   case MSH_TET_35 :
+    F.numFaces = 4;
     F.monomials = generatePascalTetrahedron(4);
     F.points =    gmshGeneratePointsTetrahedron(4, false);
     generate3dFaceClosure(F.faceClosure, 4);
     break;
   case MSH_TET_34 :
+    F.numFaces = 4;
     F.monomials = generatePascalSerendipityTetrahedron(4);
     F.points =    gmshGeneratePointsTetrahedron(4, true);
     generate3dFaceClosure(F.faceClosure, 4);
     break;
   case MSH_TET_52 :
+    F.numFaces = 4;
     F.monomials = generatePascalSerendipityTetrahedron(5);
     F.points =    gmshGeneratePointsTetrahedron(5, true);
     generate3dFaceClosure(F.faceClosure, 5);
     break;
   case MSH_TET_56 :
+    F.numFaces = 4;
     F.monomials = generatePascalTetrahedron(5);
     F.points =    gmshGeneratePointsTetrahedron(5, false);
     generate3dFaceClosure(F.faceClosure, 5);
@@ -927,6 +1030,19 @@ const polynomialBasis &polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(5,true);
     generate2dEdgeClosure(F.edgeClosure, 5, 4);
     break;
+  case MSH_PRI_6 : // first order
+    F.numFaces = 5;
+    F.monomials = generatePascalPrism(1);
+    F.points =    gmshGeneratePointsPrism(1, false);
+    generate3dFaceClosurePrism(F.faceClosure, 1);
+    break;
+  case MSH_PRI_18 : // second order
+    F.numFaces = 5;
+    F.monomials = generatePascalPrism(2);
+    F.points =    gmshGeneratePointsPrism(2, false);
+    generate3dFaceClosurePrism(F.faceClosure, 2);
+    break;
+    
   default :
     Msg::Error("Unknown function space %d: reverting to TET_4", tag);
     F.monomials = generatePascalTetrahedron(1);

Numeric/polynomialBasis.h

@@ -23,10 +23,13 @@ class polynomialBasis
   fullMatrix<double> points;
   fullMatrix<double> monomials;
   fullMatrix<double> coefficients;
+  int numFaces;
   // for a given face/edge, with both a sign and a rotation,
   // give an ordered list of nodes on this face/edge
   inline int getFaceClosureId(int iFace, int iSign, int iRot) const {
-    return iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot;
+//     return iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot; // only for tetrahedra
+//     return iFace + 5*(iSign == 1 ? 0 : 1) + 10*iRot; // only for prisms
+    return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot; // both tetrahedra and prisms
   }
   inline const std::vector<int> &getFaceClosure(int id) const
   {

Solver/TESTCASES/Advection1D.lua

@@ -16,7 +16,7 @@ v:set(2,0,0)
 nu=fullMatrix(1,1);
 nu:set(0,0,0)
 
-law = dgConservationLawAdvection(functionConstant(v):getName(), '') 
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '') 
 --FunctionConstant(nu):getName())
 
 dg:setConservationLaw(law)

Solver/TESTCASES/AdvectionDiffusion.lua

@@ -24,7 +24,7 @@ v:set(2,0,0)
 nu=fullMatrix(1,1);
 nu:set(0,0,0.001)
 
-law = dgConservationLawAdvection(functionConstant(v):getName(),functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),functionConstant(nu):getName())
 dg:setConservationLaw(law)
 
 -- boundary condition

Solver/TESTCASES/edge.msh

 $MeshFormat
-2.1 0 8
+2.2 0 8
 $EndMeshFormat
 $PhysicalNames
 3
@@ -212,205 +212,205 @@ $Nodes
 $EndNodes
 $Elements
 201
-1 15 3 1 1 0 1
-2 15 3 2 2 0 2
-3 1 3 3 1 0 1 3
-4 1 3 3 1 0 3 4
-5 1 3 3 1 0 4 5
-6 1 3 3 1 0 5 6
-7 1 3 3 1 0 6 7
-8 1 3 3 1 0 7 8
-9 1 3 3 1 0 8 9
-10 1 3 3 1 0 9 10
-11 1 3 3 1 0 10 11
-12 1 3 3 1 0 11 12
-13 1 3 3 1 0 12 13
-14 1 3 3 1 0 13 14
-15 1 3 3 1 0 14 15
-16 1 3 3 1 0 15 16
-17 1 3 3 1 0 16 17
-18 1 3 3 1 0 17 18
-19 1 3 3 1 0 18 19
-20 1 3 3 1 0 19 20
-21 1 3 3 1 0 20 21
-22 1 3 3 1 0 21 22
-23 1 3 3 1 0 22 23
-24 1 3 3 1 0 23 24
-25 1 3 3 1 0 24 25
-26 1 3 3 1 0 25 26
-27 1 3 3 1 0 26 27
-28 1 3 3 1 0 27 28
-29 1 3 3 1 0 28 29
-30 1 3 3 1 0 29 30
-31 1 3 3 1 0 30 31
-32 1 3 3 1 0 31 32
-33 1 3 3 1 0 32 33
-34 1 3 3 1 0 33 34
-35 1 3 3 1 0 34 35
-36 1 3 3 1 0 35 36
-37 1 3 3 1 0 36 37
-38 1 3 3 1 0 37 38
-39 1 3 3 1 0 38 39
-40 1 3 3 1 0 39 40
-41 1 3 3 1 0 40 41
-42 1 3 3 1 0 41 42
-43 1 3 3 1 0 42 43
-44 1 3 3 1 0 43 44
-45 1 3 3 1 0 44 45
-46 1 3 3 1 0 45 46
-47 1 3 3 1 0 46 47
-48 1 3 3 1 0 47 48
-49 1 3 3 1 0 48 49
-50 1 3 3 1 0 49 50
-51 1 3 3 1 0 50 51
-52 1 3 3 1 0 51 52
-53 1 3 3 1 0 52 53
-54 1 3 3 1 0 53 54
-55 1 3 3 1 0 54 55
-56 1 3 3 1 0 55 56
-57 1 3 3 1 0 56 57
-58 1 3 3 1 0 57 58
-59 1 3 3 1 0 58 59
-60 1 3 3 1 0 59 60
-61 1 3 3 1 0 60 61
-62 1 3 3 1 0 61 62
-63 1 3 3 1 0 62 63
-64 1 3 3 1 0 63 64
-65 1 3 3 1 0 64 65
-66 1 3 3 1 0 65 66
-67 1 3 3 1 0 66 67
-68 1 3 3 1 0 67 68
-69 1 3 3 1 0 68 69
-70 1 3 3 1 0 69 70
-71 1 3 3 1 0 70 71
-72 1 3 3 1 0 71 72
-73 1 3 3 1 0 72 73
-74 1 3 3 1 0 73 74
-75 1 3 3 1 0 74 75
-76 1 3 3 1 0 75 76
-77 1 3 3 1 0 76 77
-78 1 3 3 1 0 77 78
-79 1 3 3 1 0 78 79
-80 1 3 3 1 0 79 80
-81 1 3 3 1 0 80 81
-82 1 3 3 1 0 81 82
-83 1 3 3 1 0 82 83
-84 1 3 3 1 0 83 84
-85 1 3 3 1 0 84 85
-86 1 3 3 1 0 85 86
-87 1 3 3 1 0 86 87
-88 1 3 3 1 0 87 88
-89 1 3 3 1 0 88 89
-90 1 3 3 1 0 89 90
-91 1 3 3 1 0 90 91
-92 1 3 3 1 0 91 92
-93 1 3 3 1 0 92 93
-94 1 3 3 1 0 93 94
-95 1 3 3 1 0 94 95
-96 1 3 3 1 0 95 96
-97 1 3 3 1 0 96 97
-98 1 3 3 1 0 97 98
-99 1 3 3 1 0 98 99
-100 1 3 3 1 0 99 100
-101 1 3 3 1 0 100 101
-102 1 3 3 1 0 101 102
-103 1 3 3 1 0 102 103
-104 1 3 3 1 0 103 104
-105 1 3 3 1 0 104 105
-106 1 3 3 1 0 105 106
-107 1 3 3 1 0 106 107
-108 1 3 3 1 0 107 108
-109 1 3 3 1 0 108 109
-110 1 3 3 1 0 109 110
-111 1 3 3 1 0 110 111
-112 1 3 3 1 0 111 112
-113 1 3 3 1 0 112 113
-114 1 3 3 1 0 113 114
-115 1 3 3 1 0 114 115
-116 1 3 3 1 0 115 116
-117 1 3 3 1 0 116 117
-118 1 3 3 1 0 117 118
-119 1 3 3 1 0 118 119
-120 1 3 3 1 0 119 120
-121 1 3 3 1 0 120 121
-122 1 3 3 1 0 121 122
-123 1 3 3 1 0 122 123
-124 1 3 3 1 0 123 124
-125 1 3 3 1 0 124 125
-126 1 3 3 1 0 125 126
-127 1 3 3 1 0 126 127
-128 1 3 3 1 0 127 128
-129 1 3 3 1 0 128 129
-130 1 3 3 1 0 129 130
-131 1 3 3 1 0 130 131
-132 1 3 3 1 0 131 132
-133 1 3 3 1 0 132 133
-134 1 3 3 1 0 133 134
-135 1 3 3 1 0 134 135
-136 1 3 3 1 0 135 136
-137 1 3 3 1 0 136 137
-138 1 3 3 1 0 137 138
-139 1 3 3 1 0 138 139
-140 1 3 3 1 0 139 140
-141 1 3 3 1 0 140 141
-142 1 3 3 1 0 141 142
-143 1 3 3 1 0 142 143
-144 1 3 3 1 0 143 144
-145 1 3 3 1 0 144 145
-146 1 3 3 1 0 145 146
-147 1 3 3 1 0 146 147
-148 1 3 3 1 0 147 148
-149 1 3 3 1 0 148 149
-150 1 3 3 1 0 149 150
-151 1 3 3 1 0 150 151
-152 1 3 3 1 0 151 152
-153 1 3 3 1 0 152 153
-154 1 3 3 1 0 153 154
-155 1 3 3 1 0 154 155
-156 1 3 3 1 0 155 156
-157 1 3 3 1 0 156 157
-158 1 3 3 1 0 157 158
-159 1 3 3 1 0 158 159
-160 1 3 3 1 0 159 160
-161 1 3 3 1 0 160 161
-162 1 3 3 1 0 161 162
-163 1 3 3 1 0 162 163
-164 1 3 3 1 0 163 164
-165 1 3 3 1 0 164 165
-166 1 3 3 1 0 165 166
-167 1 3 3 1 0 166 167
-168 1 3 3 1 0 167 168
-169 1 3 3 1 0 168 169
-170 1 3 3 1 0 169 170
-171 1 3 3 1 0 170 171
-172 1 3 3 1 0 171 172
-173 1 3 3 1 0 172 173
-174 1 3 3 1 0 173 174
-175 1 3 3 1 0 174 175
-176 1 3 3 1 0 175 176
-177 1 3 3 1 0 176 177
-178 1 3 3 1 0 177 178
-179 1 3 3 1 0 178 179
-180 1 3 3 1 0 179 180
-181 1 3 3 1 0 180 181
-182 1 3 3 1 0 181 182
-183 1 3 3 1 0 182 183
-184 1 3 3 1 0 183 184
-185 1 3 3 1 0 184 185
-186 1 3 3 1 0 185 186
-187 1 3 3 1 0 186 187
-188 1 3 3 1 0 187 188
-189 1 3 3 1 0 188 189
-190 1 3 3 1 0 189 190
-191 1 3 3 1 0 190 191
-192 1 3 3 1 0 191 192
-193 1 3 3 1 0 192 193
-194 1 3 3 1 0 193 194
-195 1 3 3 1 0 194 195
-196 1 3 3 1 0 195 196
-197 1 3 3 1 0 196 197
-198 1 3 3 1 0 197 198
-199 1 3 3 1 0 198 199
-200 1 3 3 1 0 199 200
-201 1 3 3 1 0 200 2
+1 15 2 1 1 1
+2 15 2 2 2 2
+3 1 2 3 1 1 3
+4 1 2 3 1 3 4
+5 1 2 3 1 4 5
+6 1 2 3 1 5 6
+7 1 2 3 1 6 7
+8 1 2 3 1 7 8
+9 1 2 3 1 8 9
+10 1 2 3 1 9 10
+11 1 2 3 1 10 11
+12 1 2 3 1 11 12
+13 1 2 3 1 12 13
+14 1 2 3 1 13 14
+15 1 2 3 1 14 15
+16 1 2 3 1 15 16
+17 1 2 3 1 16 17
+18 1 2 3 1 17 18
+19 1 2 3 1 18 19
+20 1 2 3 1 19 20
+21 1 2 3 1 20 21
+22 1 2 3 1 21 22
+23 1 2 3 1 22 23
+24 1 2 3 1 23 24
+25 1 2 3 1 24 25
+26 1 2 3 1 25 26
+27 1 2 3 1 26 27
+28 1 2 3 1 27 28
+29 1 2 3 1 28 29
+30 1 2 3 1 29 30
+31 1 2 3 1 30 31
+32 1 2 3 1 31 32
+33 1 2 3 1 32 33
+34 1 2 3 1 33 34
+35 1 2 3 1 34 35
+36 1 2 3 1 35 36
+37 1 2 3 1 36 37
+38 1 2 3 1 37 38
+39 1 2 3 1 38 39
+40 1 2 3 1 39 40
+41 1 2 3 1 40 41
+42 1 2 3 1 41 42
+43 1 2 3 1 42 43
+44 1 2 3 1 43 44
+45 1 2 3 1 44 45
+46 1 2 3 1 45 46
+47 1 2 3 1 46 47
+48 1 2 3 1 47 48
+49 1 2 3 1 48 49
+50 1 2 3 1 49 50
+51 1 2 3 1 50 51
+52 1 2 3 1 51 52
+53 1 2 3 1 52 53
+54 1 2 3 1 53 54
+55 1 2 3 1 54 55
+56 1 2 3 1 55 56
+57 1 2 3 1 56 57
+58 1 2 3 1 57 58
+59 1 2 3 1 58 59
+60 1 2 3 1 59 60
+61 1 2 3 1 60 61
+62 1 2 3 1 61 62
+63 1 2 3 1 62 63
+64 1 2 3 1 63 64
+65 1 2 3 1 64 65
+66 1 2 3 1 65 66
+67 1 2 3 1 66 67
+68 1 2 3 1 67 68
+69 1 2 3 1 68 69
+70 1 2 3 1 69 70
+71 1 2 3 1 70 71
+72 1 2 3 1 71 72
+73 1 2 3 1 72 73
+74 1 2 3 1 73 74
+75 1 2 3 1 74 75
+76 1 2 3 1 75 76
+77 1 2 3 1 76 77
+78 1 2 3 1 77 78
+79 1 2 3 1 78 79
+80 1 2 3 1 79 80
+81 1 2 3 1 80 81
+82 1 2 3 1 81 82
+83 1 2 3 1 82 83
+84 1 2 3 1 83 84
+85 1 2 3 1 84 85
+86 1 2 3 1 85 86
+87 1 2 3 1 86 87
+88 1 2 3 1 87 88
+89 1 2 3 1 88 89
+90 1 2 3 1 89 90
+91 1 2 3 1 90 91
+92 1 2 3 1 91 92
+93 1 2 3 1 92 93
+94 1 2 3 1 93 94
+95 1 2 3 1 94 95
+96 1 2 3 1 95 96
+97 1 2 3 1 96 97
+98 1 2 3 1 97 98
+99 1 2 3 1 98 99
+100 1 2 3 1 99 100
+101 1 2 3 1 100 101
+102 1 2 3 1 101 102
+103 1 2 3 1 102 103
+104 1 2 3 1 103 104
+105 1 2 3 1 104 105
+106 1 2 3 1 105 106
+107 1 2 3 1 106 107
+108 1 2 3 1 107 108
+109 1 2 3 1 108 109
+110 1 2 3 1 109 110
+111 1 2 3 1 110 111
+112 1 2 3 1 111 112
+113 1 2 3 1 112 113
+114 1 2 3 1 113 114
+115 1 2 3 1 114 115
+116 1 2 3 1 115 116
+117 1 2 3 1 116 117
+118 1 2 3 1 117 118
+119 1 2 3 1 118 119
+120 1 2 3 1 119 120
+121 1 2 3 1 120 121
+122 1 2 3 1 121 122
+123 1 2 3 1 122 123
+124 1 2 3 1 123 124
+125 1 2 3 1 124 125
+126 1 2 3 1 125 126
+127 1 2 3 1 126 127
+128 1 2 3 1 127 128
+129 1 2 3 1 128 129
+130 1 2 3 1 129 130
+131 1 2 3 1 130 131
+132 1 2 3 1 131 132
+133 1 2 3 1 132 133
+134 1 2 3 1 133 134
+135 1 2 3 1 134 135
+136 1 2 3 1 135 136
+137 1 2 3 1 136 137
+138 1 2 3 1 137 138
+139 1 2 3 1 138 139
+140 1 2 3 1 139 140
+141 1 2 3 1 140 141
+142 1 2 3 1 141 142
+143 1 2 3 1 142 143
+144 1 2 3 1 143 144
+145 1 2 3 1 144 145
+146 1 2 3 1 145 146
+147 1 2 3 1 146 147
+148 1 2 3 1 147 148
+149 1 2 3 1 148 149
+150 1 2 3 1 149 150
+151 1 2 3 1 150 151
+152 1 2 3 1 151 152
+153 1 2 3 1 152 153
+154 1 2 3 1 153 154
+155 1 2 3 1 154 155
+156 1 2 3 1 155 156
+157 1 2 3 1 156 157
+158 1 2 3 1 157 158
+159 1 2 3 1 158 159
+160 1 2 3 1 159 160
+161 1 2 3 1 160 161
+162 1 2 3 1 161 162
+163 1 2 3 1 162 163
+164 1 2 3 1 163 164
+165 1 2 3 1 164 165
+166 1 2 3 1 165 166
+167 1 2 3 1 166 167
+168 1 2 3 1 167 168
+169 1 2 3 1 168 169
+170 1 2 3 1 169 170
+171 1 2 3 1 170 171
+172 1 2 3 1 171 172
+173 1 2 3 1 172 173
+174 1 2 3 1 173 174
+175 1 2 3 1 174 175
+176 1 2 3 1 175 176
+177 1 2 3 1 176 177
+178 1 2 3 1 177 178
+179 1 2 3 1 178 179
+180 1 2 3 1 179 180
+181 1 2 3 1 180 181
+182 1 2 3 1 181 182
+183 1 2 3 1 182 183
+184 1 2 3 1 183 184
+185 1 2 3 1 184 185
+186 1 2 3 1 185 186
+187 1 2 3 1 186 187
+188 1 2 3 1 187 188
+189 1 2 3 1 188 189
+190 1 2 3 1 189 190
+191 1 2 3 1 190 191
+192 1 2 3 1 191 192
+193 1 2 3 1 192 193
+194 1 2 3 1 193 194
+195 1 2 3 1 194 195
+196 1 2 3 1 195 196
+197 1 2 3 1 196 197
+198 1 2 3 1 197 198
+199 1 2 3 1 198 199
+200 1 2 3 1 199 200
+201 1 2 3 1 200 2
 $EndElements

Solver/dgGroupOfElements.cpp

@@ -408,7 +408,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
   init(pOrder);
 }
 
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, int numVertices):
   _groupLeft(elGroup),_groupRight(elGroup)
 {
   _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
@@ -456,6 +456,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
         MElement &el = *elGroup.getElement(i);
         for (int j=0; j<el.getNumFaces(); j++){
           MFace face = el.getFace(j);
+          if (numVertices < 0 || face.getNumVertices() == numVertices) {
             if(faceMap.find(face) == faceMap.end()){
               faceMap[face] = i;
             }else{
@@ -463,6 +464,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
             }
           }
         }
+      }
       break;
     }
     default : throw;
@@ -470,7 +472,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
   init(pOrder);
 }
 
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder):
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder, int numVertices):
   _groupLeft(elGroup1),_groupRight(elGroup2)
 {
   _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
@@ -546,6 +548,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
         MElement &el = *elGroup1.getElement(i);
         for (int j=0; j<el.getNumFaces(); j++){
           MFace face = el.getFace(j);
+          if (numVertices < 0 || face.getNumVertices() == numVertices) {
             if(faceMap.find(face) == faceMap.end()){
               faceMap[face] = i;
             }else{
@@ -553,6 +556,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
             }
           }
         }
+      }
       for(int i=0; i<elGroup2.getNbElements(); i++){
         MElement &el = *elGroup2.getElement(i);
         for (int j=0; j<el.getNumFaces(); j++){
@@ -703,6 +707,36 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
         if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
           localElements[el->getType()].push_back(el);
           nlocal++;
+          switch(dim) {
+            case 1: {
+              int interfaceType = 1; // number of vertices
+              if(_interfaceTypes.find(interfaceType) == _interfaceTypes.end()) {
+                _interfaceTypes.insert(interfaceType);
+//                 printf("Inserted interfaceType: %d el: %d dim: %d\n",interfaceType,el->getType(),dim);
+              }
+              break;
+            }
+            case 2: {
+              int interfaceType = 2; // number of vertices
+              if(_interfaceTypes.find(interfaceType) == _interfaceTypes.end()) {
+                _interfaceTypes.insert(interfaceType);
+//                 printf("Inserted interfaceType: %d el: %d dim: %d\n",interfaceType,el->getType(),dim);
+              }
+              break;
+            }
+            case 3: {
+              for(int iFace=0; iFace<el->getNumFaces(); iFace++) {
+                MFace face = el->getFace(iFace);
+                if(_interfaceTypes.find(face.getNumVertices()) == _interfaceTypes.end()) {
+                  _interfaceTypes.insert(face.getNumVertices());
+//                   printf("Inserted interfaceType: %d el: %d dim: %d\n",face.getNumVertices(),el->getType(),dim);
+                }
+              }
+              break;
+            }
+            default :
+            throw;
+          }
         }
       }
     }
@@ -730,7 +764,7 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
 
   std::map<const std::string,std::set<MVertex*> > boundaryVertices;
   std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
-  std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
+  std::map<const int,std::map<const std::string,std::set<MFace, Less_Face> > > boundaryFaces; // [elementType][bndString]
 
   std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
   int nghosts=0;
@@ -754,9 +788,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
             case 2:
               boundaryEdges[physicalName].insert( element->getEdge(0) );
             break;
-            case 3:
-              boundaryFaces[physicalName].insert( element->getFace(0));
+            case 3: {
+              MFace face = element->getFace(0);
+              boundaryFaces[element->getType()][physicalName].insert( face );
               break;
+            }
             default :
             throw;
           }
@@ -807,32 +843,38 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
       break;
     }
     case 3 : {
+    std::map<const int,std::map<const std::string,std::set<MFace, Less_Face> > >::iterator elemTypeIt;
     std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
-      for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
+      for(elemTypeIt=boundaryFaces.begin(); elemTypeIt!=boundaryFaces.end(); elemTypeIt++) {
+        for(mapIt=elemTypeIt->second.begin(); mapIt!=elemTypeIt->second.end(); mapIt++) {
           dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
           if(gof->getNbElements())
             _boundaryGroups.push_back(gof);
           else
             delete gof;
         }
+      }
       break;
     }
     }
     // create a group of faces for every face that is common to elements of the same group
-    dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],order);
+    // create separate groups for each face type
+    for(std::set<int>::iterator faceTypeIt=_interfaceTypes.begin(); faceTypeIt!=_interfaceTypes.end(); faceTypeIt++) {
+      dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],order,*faceTypeIt);
       if (gof->getNbElements())
         _faceGroups.push_back(gof);
       else
         delete gof;
 //       create a groupe of d
       for (int j=i+1;j<_elementGroups.size();j++){
-      dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
+        dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order,*faceTypeIt);
         if (gof->getNbElements())
           _faceGroups.push_back(gof);
         else
           delete gof;
       }
     }
+  }
   //create ghost groups
   for(int i=0;i<Msg::GetCommSize();i++){
     for (std::map<int, std::vector<MElement *> >::iterator it = ghostElements[i].begin(); it !=ghostElements[i].end() ; ++it){
@@ -846,14 +888,16 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
   //create face group for ghostGroups
   for (int i=0; i<_ghostGroups.size(); i++){
     for (int j=0;j<_elementGroups.size();j++){
-      dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order);
+      for(std::set<int>::iterator faceTypeIt=_interfaceTypes.begin(); faceTypeIt!=_interfaceTypes.end(); faceTypeIt++) {
+        dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order,*faceTypeIt);
         if (gof->getNbElements())
           _faceGroups.push_back(gof);
         else
           delete gof;
       }
     }
-
+  }
+  Msg::Info("%d groups of interfaces",_faceGroups.size());
   // build the ghosts structure
   int *nGhostElements = new int[Msg::GetCommSize()];
   int *nParentElements = new int[Msg::GetCommSize()];

Solver/dgGroupOfElements.h

@@ -184,8 +184,8 @@ public:
   inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _integrationPointsRight[_closuresIdRight[iFace]];}
   
   inline fullMatrix<double> &getNormals () const {return *_normals;}
-  dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
-  dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder);
+  dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder, int numVertices = -1);
+  dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder, int numVertices = -1);
   dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
   dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
   dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);
@@ -217,6 +217,8 @@ class dgGroupCollection {
   std::vector<dgGroupOfFaces*> _faceGroups; //interface
   std::vector<dgGroupOfFaces*> _boundaryGroups; //boundary
   std::vector<dgGroupOfElements*> _ghostGroups; //ghost volume
+  // container of different face types (identified by the number of vertices) 
+  std::set<int> _interfaceTypes;
 
   //{group,id} of the elements to send to each partition for a scatter operation
   std::vector< std::vector<std::pair<int,int> > >_elementsToSend;