dg work on 3D (unfinished)

Solver/TESTCASES/Advection3D.lua

+model = GModel  ()
+model:load ('box.geo')
+model:load ('box.msh')
+dg = dgSystemOfEquations (model)
+dg:setOrder(1)
+
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0.15)
+v:set(1,0,0)
+v:set(2,0,0)
+-- diffusivity
+nu=fullMatrix(1,1);
+nu:set(0,0,0.001)
+
+--law = ConservationLawAdvection(FunctionConstant(v):getName(),FunctionConstant(nu):getName())
+--law = ConservationLawAdvection('',FunctionConstant(nu):getName())
+law = ConservationLawAdvection(FunctionConstant(v):getName(),'')
+
+dg:setConservationLaw(law)
+
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,1)
+bndcondition=law:newOutsideValueBoundary(FunctionConstant(outside):getName())
+--[[law:addBoundaryCondition('Left',bndcondition)
+law:addBoundaryCondition('Right',bndcondition)
+--]]
+
+dg:setup()
+
+-- initial condition
+function initial_condition( xyz , f )
+  for i=0,xyz:size1()-1 do
+    x = xyz:get(i,0)
+    --y = xyz:get(i,1)
+    --z = xyz:get(i,2)
+    --f:set (i, 0, math.exp(-100*((x-0.5)^2)))
+    f:set (i, 0, x )
+  end
+end
+-- dg:L2Projection(FunctionLua(1,'initial_condition',{'XYZ'}):getName())
+
+law:addBoundaryCondition('boundary',law:new0FluxBoundary())
+
+dg:exportSolution('output/Adv3D_00000')
+
+-- main loop
+n = 5
+for i=1,100*n do
+  norm = dg:RK44(0.0003)
+  if (i % n == 0) then 
+    print('iter',i,norm)
+    dg:exportSolution(string.format("output/Adv3D-%05d", i)) 
+  end
+end

Solver/TESTCASES/AdvectionDiffusion.lua

@@ -10,7 +10,7 @@ vmodel[1]:load('')
 vmodel[2]:load('')
 
 dg = dgSystemOfEquations (model)
-dg:setOrder(5)
+dg:setOrder(1)
 
 -- conservation law
 
@@ -49,7 +49,7 @@ dg:exportSolution('output/Advection_00000')
 
 -- main loop
 for i=1,10000 do
-  norm = dg:RK44(0.001)
+  norm = dg:RK44(0.01)
   if (i % 1 == 0) then 
     print('iter',i,norm)
     dg:exportSolution(string.format("output/Advection-%05d", i)) 

Solver/TESTCASES/WavePulse.lua

@@ -22,7 +22,7 @@ myModel:load('square.msh')
 
 print'*** Create a dg solver ***'
 DG = dgSystemOfEquations (myModel)
-DG:setOrder(1)
+DG:setOrder(3)
 law=ConservationLawWaveEquation(2)
 DG:setConservationLaw(law)
 law:addBoundaryCondition('Border',law:newWallBoundary())

Solver/TESTCASES/box.geo

@@ -17,9 +17,9 @@ Line(4) = {4,1};
 Line Loop(4) = {1,2,3,4};
 Plane Surface(5) = {4};
 
-num[]  = Extrude {0, 0, 0.1} {
+Extrude {0, 0, 1} {
   Surface{5};
 }
 
-Physical Surface(30) = {22, 14, 18, 27, 5, 26};
-Physical Volume(31) = {1};
+Physical Surface("boundary") = {27, 22, 18, 14, 5, 26};
+Physical Volume("volume") = {1};

Solver/dgAlgorithm.cpp

@@ -233,6 +233,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
   }
 
   // C ) redistribute the flux to the residual (at Faces nodes)
+
   if(riemannSolver || diffusiveFluxLeft)
     residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
 
@@ -387,10 +388,10 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
               boundaryVertices[physicalName].insert( element->getVertex(0) ); 
               break;
             case 2:
-              boundaryEdges[physicalName].insert( MEdge(element->getVertex(0), element->getVertex(1)) );
+              boundaryEdges[physicalName].insert( element->getEdge(0) );
             break;
             case 3:
-              boundaryFaces[physicalName].insert( MFace(element->getVertex(0), element->getVertex(1),element->getVertex(2)) );
+              boundaryFaces[physicalName].insert( element->getFace(0));
             break;
             default :
             throw;
@@ -403,7 +404,6 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
     }
   }
   eGroups.push_back(new dgGroupOfElements(allElements,order));
-  fGroups.push_back(new dgGroupOfFaces(*eGroups[0],order));
   switch(dim) {
     case 1 : {
       std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
@@ -427,6 +427,7 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
       break;
     }
   }
+  fGroups.push_back(new dgGroupOfFaces(*eGroups[0],order));
 }
 
 // works for any number of groups 

Solver/dgGroupOfElements.cpp

@@ -136,7 +136,8 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
   MElement &elLeft = *_groupLeft.getElement(iElLeft);
   int ithFace, sign, rot;
   elLeft.getFaceInfo (topoFace, ithFace, sign, rot);
-  _closuresIdLeft.push_back(_fsLeft->getFaceClosureId(ithFace, sign, rot));
+  int closureIdLeft = _fsLeft->getFaceClosureId(ithFace, sign, rot);
+  _closuresIdLeft.push_back(closureIdLeft);
   if(iElRight>=0){
     _right.push_back(iElRight);
     MElement &elRight = *_groupRight.getElement(iElRight);
@@ -146,7 +147,7 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
   // compute the face element that correspond to the geometrical closure
   // get the vertices of the face
   std::vector<MVertex*> vertices;
-  const std::vector<int> & geomClosure = elLeft.getFunctionSpace()->getFaceClosure(_closuresIdLeft.back());
+  const std::vector<int> & geomClosure = elLeft.getFunctionSpace()->getFaceClosure(closureIdLeft);
   for (int j=0; j<geomClosure.size() ; j++)
     vertices.push_back( elLeft.getVertex(geomClosure[j]) );
   // triangular face
@@ -234,7 +235,7 @@ void dgGroupOfFaces::init(int pOrder) {
 
   // compute data on quadrature points : normals and dPsidX
   double g[256][3];
-  _normals = new fullMatrix<double> (3,_fsFace->points.size1()*_faces.size());
+  _normals = new fullMatrix<double> (3,_integration->size1()*_faces.size());
   _dPsiLeftDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsLeft->points.size1()*_faces.size());
   if(_fsRight){
     _dPsiRightDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsRight->points.size1()*_faces.size());
@@ -262,9 +263,10 @@ void dgGroupOfFaces::init(int pOrder) {
       double &nz=(*_normals)(2,index);
       double nu=0,nv=0,nw=0;
       for (size_t k=0; k<closureLeft.size(); k++){
-        nu += g[closureLeft[k]][0];
-        nv += g[closureLeft[k]][1];
-        nw += g[closureLeft[k]][2];
+        int inode=closureLeft[k];
+        nu += g[inode][0];
+        nv += g[inode][1];
+        nw += g[inode][2];
       }
       nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
       ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];

Solver/dgGroupOfElements.h

@@ -165,7 +165,7 @@ public:
   inline fullMatrix<double> & getDPsiRightDxMatrix() const { return *_dPsiRightDxOnQP;}
   //this part is common with dgGroupOfElements, we should try polymorphism
   inline int getNbElements() const {return _faces.size();}
-  inline int getNbNodes() const {return _collocation->size1();}
+  inline int getNbNodes() const {return _collocation->size2();}
   inline int getNbIntegrationPoints() const {return _collocation->size1();}
   inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
   inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}