dg: added 2nd order Prism18

Geo/GModel.cpp

@@ -1334,7 +1334,7 @@ void GModel::load(std::string fileName)
 {
   GModel *temp = GModel::current();
   GModel::setCurrent(this);
-  MergeFile(fileName.c_str());
+  MergeFile(fileName.c_str(),true);
   GModel::setCurrent(temp);
 }
 

Geo/MPrism.cpp

@@ -37,14 +37,14 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
     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);
+    default: Msg::Error("Order %d prism 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);
+    default: Msg::Error("Order %d prism function space not implemented", order);
     }
   }
   return 0;

Geo/MPrism.h

@@ -131,7 +131,7 @@ class MPrism : public MElement {
   }
   virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual int getVolumeSign();
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) 
+/*  virtual void getShapeFunctions(double u, double v, double w, double s[], int o) 
   {
     s[0] = (1. - u - v) * (1. - w) * 0.5;
     s[1] =       u      * (1. - w) * 0.5;
@@ -139,9 +139,10 @@ class MPrism : public MElement {
     s[3] = (1. - u - v) * (1. + w) * 0.5;
     s[4] =       u      * (1. + w) * 0.5;
     s[5] =           v  * (1. + w) * 0.5;
-  }
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) 
+  }*/
+/*  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) 
   {
+    printf("here\n");
     s[0][0] = -0.5 * (1. - w)    ;
     s[0][1] = -0.5 * (1. - w)    ;
     s[0][2] = -0.5 * (1. - u - v);
@@ -160,7 +161,7 @@ class MPrism : public MElement {
     s[5][0] =  0.                ;
     s[5][1] =  0.5 * (1. + w)    ;
     s[5][2] =  0.5 * v           ;
-  }
+  }*/
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;

Numeric/GaussQuadraturePri.cpp

@@ -10,7 +10,7 @@
 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 * GQP[28] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 
 IntPt *getGQPriPts(int order)
 { 

Numeric/polynomialBasis.cpp

@@ -244,17 +244,15 @@ static fullMatrix<double> generatePascalPrism(int order)
   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();
+  monomials.print("Pri monoms");
   return monomials;
 }  
 
@@ -625,26 +623,50 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
 
 static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip) 
 {
+  const double prism18Pts[18][3] = {
+    {0, 0, -1}, // 0
+    {1, 0, -1}, // 1
+    {0, 1, -1}, // 2
+    {0, 0, 1},  // 3
+    {1, 0, 1},  // 4
+    {0, 1, 1},  // 5
+    {0.5, 0, -1},  // 6
+    {0, 0.5, -1},  // 7
+    {0, 0, 0},  // 8
+    {0.5, 0.5, -1},  // 9
+    {1, 0, 0},  // 10
+    {0, 1, 0},  // 11
+    {0.5, 0, 1},  // 12
+    {0, 0.5, 1},  // 13
+    {0.5, 0.5, 1},  // 14
+    {0.5, 0, 0},  // 15
+    {0, 0.5, 0},  // 16
+    {0.5, 0.5, 0},  // 17
+  };
+
   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);
+
+  if (order == 2)
+    for (int i =0; i<18; i++)
+      for (int j=0; j<3;j++)
+        point(i,j) = prism18Pts[i][j];
+  else    
     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);  
+//   point.print("Pri ipts");
+
   return point;
 }
 
@@ -798,26 +820,31 @@ 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)
+  if (order > 2)
     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:
+  if (order==0) {
     closure[0] = 0;
-    break;
-  default:
-//     int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
+    return;
+  }
   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 order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15}, {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
+//   int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8}, {8, 13, 11, 7}, {9, 11, 14, 10}};
   int nVertex = isTriangle ? 3 : 4;
   for (int i = 0; i < nVertex; ++i){
-      int k;
-        k = (nVertex + (iSign * i) + iRotate) % nVertex;  //- iSign * iRotate
+    int k = (nVertex + (iSign * i) + iRotate) % nVertex;  //- iSign * iRotate
     closure[i] = order1node[iFace][k];
   }
-    break;
+  if (order==2) {
+    for (int i = 0; i < nVertex; ++i){
+      int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;  //- iSign * iRotate
+      closure[nVertex+i] = order2node[iFace][k];
+    }
+    if (!isTriangle)
+      closure[nNodes-1] = order2node[iFace][4]; // center
   }
 }
 
@@ -831,8 +858,6 @@ static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int ord
         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");
       }
     }
   }
@@ -1075,6 +1100,14 @@ const polynomialBasis &polynomialBases::find(int tag)
     break;
   }  
   F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
+//   printf("Case: %d coeffs:\n",tag);
+//   for (int i = 0; i<F.coefficients.size1(); i++) {
+//     for (int j = 0; j<F.coefficients.size2(); j++) {
+//       printf("%4.1f ",F.coefficients(i,j));
+//     }
+//     printf("\n");
+//   }
+    
   fs.insert(std::make_pair(tag, F));
   return fs[tag];
 }

Solver/TESTCASES/Diffusion3D.lua

@@ -6,7 +6,7 @@ model = GModel()
 model:load('mixed.geo')
 model:load('mixed.msh')
 dg = dgSystemOfEquations(model)
-dg:setOrder(1)
+dg:setOrder(2)
 
 -- initial condition
 function initial_condition( xyz , f )
@@ -37,17 +37,6 @@ dg:setConservationLaw(law)
 -- boundary condition
 outside=fullMatrix(1,1)
 outside:set(0,0,0.15)
--- freestream=law:newOutsideValueBoundary(functionLua(1,'initial_condition',{'XYZ'}):getName())
-
--- law:addBoundaryCondition('wall',law:new0FluxBoundary())
--- law:addBoundaryCondition('inlet',freestream)
--- law:addBoundaryCondition('outlet',law:newSymmetryBoundary())
--- law:addBoundaryCondition('symmetry',law:newSymmetryBoundary())
-
--- law:addBoundaryCondition('wall',law:new0FluxBoundary())
--- law:addBoundaryCondition('inlet',law:new0FluxBoundary())
--- law:addBoundaryCondition('outlet',law:new0FluxBoundary())
--- law:addBoundaryCondition('symmetry',law:new0FluxBoundary())
 
 law:addBoundaryCondition('side1',law:new0FluxBoundary())
 law:addBoundaryCondition('side2',law:new0FluxBoundary())
@@ -60,10 +49,10 @@ dg:setup()
 
 dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
 
-dt = 0.4
+dt = 0.0001
 --CFL = 20 -- good for lc=0.1 mesh
-CFL = 10
-dt = CFL*dg:computeInvSpectralRadius()
+CFL = 1
+-- dt = CFL*dg:computeInvSpectralRadius()
 -- export_interval = 5*dt
 export_interval = 0.05
 end_time = 3
@@ -78,11 +67,13 @@ export_count = 0
 io.write(string.format('Export: %d iter: %4d t: %.6f n: %.8f dt:%.8f\n',export_count,0,0,0,dt))
 dg:exportSolution(output_pattern .. string.format("%05d",export_count))
 for i=1,10000 do
-  dt = CFL*dg:computeInvSpectralRadius()
+--   dt = CFL*dg:computeInvSpectralRadius()
   norm = dg:RK44(dt)
   time = i*dt
   io.write('.')
   io.flush()
+--     io.write('\n')
+--     io.write(string.format('Export: %d iter: %4d t: %.6f n: %.8f dt:%.8f\n',export_count,i,time,norm,dt))
   if (time >= next_export_time) then
     export_count = export_count +1
     io.write('\n')

Solver/TESTCASES/mixed.geo

@@ -2,7 +2,7 @@
 C = 1;
 Lup = 1;
 L = 1.;
-lc = .3;
+lc = 0.3;
 
 Point(1) = {0.0, 0.0, -Lup, lc};
 Point(2) = {C  , 0.0, -Lup, lc};
@@ -23,7 +23,7 @@ outtet[] = Extrude {0,0,0.7*L} { Surface{6};};
 // Printf("volume = %g", outtet[1]);
 // Printf("side surfaces = %g %g %g %g", outtet[2], outtet[3], outtet[4], outtet[5]);
 
-outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{4};Recombine;};
+outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{Ceil(0.5*L/lc)};Recombine;};
 // outv[]= Extrude {0,0,0.5*L}{ Surface{7}; Layers{5};Recombine;};
 // Printf("top surface = %g", outpri[0]);
 // Printf("volume = %g", outpri[1]);
@@ -32,7 +32,7 @@ outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{4};Recombine;};
 Mesh.Algorithm3D=4; // frontal [lobes]
 
 Physical Surface("top") = {outpri[0]};
-Physical Surface("bottom") = {5};
+Physical Surface("bottom") = {6};
 Physical Surface("side1") = {outpri[2],outtet[2]};
 Physical Surface("side2") = {outpri[3],outtet[3]};
 Physical Surface("side3") = {outpri[4],outtet[4]};

Solver/dgDofContainer.cpp

@@ -304,6 +304,10 @@ void dgDofContainer::exportMsh(const std::string name)
     if(Msg::GetCommSize()>1)
       name_oss<<"_"<<Msg::GetCommRank();
     FILE *f = fopen (name_oss.str().c_str(),"w");
+    if(!f){
+      Msg::Error("Unable to open export file '%s'", name.c_str());
+    }
+
     int COUNT = 0;
     for (int i=0;i < _groups.getNbElementGroups() ;i++){
       COUNT += _groups.getElementGroup(i)->getNbElements();