diff --git a/Solver/TESTCASES/Advection3D.lua b/Solver/TESTCASES/Advection3D.lua
new file mode 100644
index 0000000000000000000000000000000000000000..14de45f995fb91c58f892493fd19132f454d8c1c
--- /dev/null
+++ b/Solver/TESTCASES/Advection3D.lua
@@ -0,0 +1,58 @@
+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
diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
index 36e1f84f688cb5b81131da72af52f8e9d4566dc6..3be0149715359c7320b28a4af8f7b8558ced965e 100644
--- a/Solver/TESTCASES/AdvectionDiffusion.lua
+++ b/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))
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 0ef157ba57622cc354f45d21315b79bb459d95cd..68bd5f0dc60430c8a556650b6542f96a7f44802f 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/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())
diff --git a/Solver/TESTCASES/box.geo b/Solver/TESTCASES/box.geo
index 1f012d2f05c33e1a78d90ac98d71232cdaed2f6a..8dd3a88c9b11c06f6fc73674316514e0c617eeff 100644
--- a/Solver/TESTCASES/box.geo
+++ b/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};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 3c32bb7d7029d5dd478d977dce0abe7c1a9b5188..af63131bc0334574e59c71605440b37155cb4957 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/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
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c1cea4499dab60bd88bca5d3d72c90a7320472b5..182e6652cb155947fe83f441f4529147a46a3323 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/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];
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index bbea72d73d577ad0b8669351304b13d53cb64bb6..224e71983fb9dae32c5d50e1b138648786b79cda 100644
--- a/Solver/dgGroupOfElements.h
+++ b/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;}