diff --git a/Solver/TESTCASES/Advection_1.lua b/Solver/TESTCASES/Advection_1.lua
new file mode 100644
index 0000000000000000000000000000000000000000..57198d6827aa0ad2ac6f6a328af2217f63369554
--- /dev/null
+++ b/Solver/TESTCASES/Advection_1.lua
@@ -0,0 +1,69 @@
+model = GModel ()
+print'*** Loading the geo ***'
+model:load ('square_m.geo')
+-- model:load ('square_p.geo')
+print'*** Loading the mesh ***'
+model:load ('square_m.msh')
+-- model:load('square_p.msh')
+print'*** mesh loaded ***'
+dg = dgSystemOfEquations (model)
+order=3
+dg:setOrder(order)
+print("*** order = " .. order .. " ***")
+print'*** Model set ***'
+
+print("*** set initial condition ***")
+-- 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.2)^2 +(y+0.3)^2)))
+ end
+end
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0.01)
+v:set(1,0,0.025)
+v:set(2,0,0)
+--print("*** set advection speed to: v_x=" .. v:get(0,0) ..", v_y=" .. v:get(1,0) .. ", v_z=" .. v:get(2,0) .. " ***")
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),'')
+dg:setConservationLaw(law)
+
+print("*** set boundary conditions *** ")
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.0)
+bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+law:addBoundaryCondition('Border',bndcondition)
+
+--dg:setup()
+-- create 2 groups of elements; criterion=maxEdge
+dg:createGroups("maxEdge")
+-- dg:createGroups("minEdge")
+-- dg:createGroups("elementType")
+
+dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+
+dg:exportSolution('output/Advection_00000')
+
+dg:computeInvSpectralRadius();
+T=40
+dt=0.1
+
+print("*** default time step (dt)=" .. dt .. ", Period (T)=" ..T .." ***")
+
+N=T/dt
+print("*** Number of time steps=" .. N .. " ***")
+-- main loop
+for i=1,N do
+ norm = dg:multirateRK43(dt)
+ print('iter',i,norm)
+ if (i % 20 == 0) then
+ dg:exportSolution(string.format("output/Advection_%05d", i))
+ end
+end
+print '*** adding the mesh and the model ***'
diff --git a/Solver/TESTCASES/square_m.geo b/Solver/TESTCASES/square_m.geo
new file mode 100644
index 0000000000000000000000000000000000000000..808e014b8ad7fcf76dd04d90ecd05dae1de830eb
--- /dev/null
+++ b/Solver/TESTCASES/square_m.geo
@@ -0,0 +1,15 @@
+Point(1) = {-0.5, -0.5, 0, 0.025};
+Point(2) = {0.5, -0.5, 0, 0.025};
+Point(3) = {0.5, 0.5, 0, 0.125};
+Point(4) = {-0.5, 0.5, 0, 0.125};
+Line(1) = {4, 3};
+Line(2) = {3, 2};
+Line(3) = {2, 1};
+Line(4) = {1, 4};
+Line Loop(5) = {2, 3, 4, 1};
+Plane Surface(6) = {5};
+Physical Line("Border") = {1, 2, 3, 4};
+Physical Surface("Inside") = {6};
+//Mesh.CharacteristicLengthExtendFromBoundary=1;
+//Transfinite Line {1, 2, 4, 3} = 20 Using Progression 1;
+//Transfinite Surface {6};
diff --git a/Solver/TESTCASES/square_p.geo b/Solver/TESTCASES/square_p.geo
new file mode 100644
index 0000000000000000000000000000000000000000..5426e6d6a87cee9ebc13e9739bc9cde356d07f8d
--- /dev/null
+++ b/Solver/TESTCASES/square_p.geo
@@ -0,0 +1,31 @@
+Point(1) = {-0.5, -0.5, 0, .03};
+Point(2) = {0, -0.5, 0, .03};
+Point(3) = {0, 0.5, 0, .03};
+Point(4) = {-0.5, 0.5, 0, .03};
+
+Point(5) = {0.5, -0.5, 0, .03};
+Point(6) = {0.5, 0.5, 0, .03};
+
+Line(1) = {4, 3};
+Line(2) = {3, 2};
+Line(3) = {2, 1};
+Line(4) = {1, 4};
+
+Line(7) = {3, 6};
+Line(8) = {6, 5};
+Line(9) = {5, 2};
+
+Line Loop(5) = {2, 3, 4, 1};
+Line Loop(6) = {8, 9, -2, 7};
+Plane Surface(7) = {5};
+Plane Surface(8) = {6};
+Transfinite Line {4,8} = 20 ;
+Transfinite Line{2}=20;
+Transfinite Surface {7};
+Recombine Surface {7};
+Transfinite Line {3,1,7,9} = 10;
+Transfinite Surface{8};
+
+Physical Line("Border") = {1, 3, 4, 7, 8, 9};
+Physical Surface("Inside1") = {7};
+Physical Surface("Inside2") = {8};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 8035d305338b5cefa2c20b75e02e06c052f66306..fde2d637dbc8dc0cf3db0cadefd615cc32749f70 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -357,9 +357,11 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
// double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
// double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
-
+ double A[10][10];
+ double *b;
+ double *c;
// Small step RK43
- double A[10][10]={
+ double A1[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
@@ -371,24 +373,24 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
};
- double b[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
- double c[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+ double b1[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
+ double c1[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
// Big step RK43
-// double A[10][10]={
-//{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-//{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-//{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-//{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-//{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
-//{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
-//{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
-//{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
-//{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
-//{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
-// };
-// double b[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
-// double c[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+ double A2[10][10]={
+{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
+{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
+{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
+{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
+{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
+ };
+ double b2[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
+ double c2[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
// fullMatrix<double> K(sol);
@@ -410,11 +412,25 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
Unp.scale(0.0);
Unp.axpy(sol);
-
+ // Case of 2 groups: first with big steps, second with small steps
for(int j=0; j<nStages;j++){
tmp.scale(0.0);
tmp.axpy(sol);
for(int k=0;k < groups.getNbElementGroups();k++) {
+ if (k==0) {
+ for (int ii=0; ii<10; ii++) {
+ for (int jj=0; jj<10; jj++) {
+ A[ii][jj]=A2[ii][jj];
+ }
+ }
+ }
+ else{
+ for (int ii=0; ii<10; ii++) {
+ for (int jj=0; jj<10; jj++) {
+ A[ii][jj]=A1[ii][jj];
+ }
+ }
+ }
for(int i=0;i<j;i++){
if(fabs(A[j][i])>1e-12){
tmp.getGroupProxy(k).add(K[i]->getGroupProxy(k),h*A[j][i]);
@@ -423,6 +439,15 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
}
this->residual(claw,groups,tmp,resd);
for(int k=0;k < groups.getNbElementGroups(); k++) {
+ if (k==0) {
+ b=b2;
+ c=c2;
+ }
+ else{
+ b=b1;
+ c=c1;
+ }
+
dgGroupOfElements *group = groups.getElementGroup(k);
int nbNodes = group->getNbNodes();
for(int i=0;i<group->getNbElements();i++) {
@@ -592,7 +617,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
// provided dataCache
dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
-
+
const int nbFields = claw.nbFields();
/* This is an estimate on how lengths changes with p
It is merely the smallest distance between gauss
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 25a37420ee7a09bad5ae1b97c4ebd9b36546770f..2baaa653050e6bba98074ec8004e68ec108e306a 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -51,6 +51,8 @@ class dgAlgorithm {
fullMatrix<double> &solution, // the solution
std::vector<double> & DT // elementary time steps
);
+ // works only (for the moment) for two groups of elements, small and big step
+ // Uses RK43 from Schlegel (10 stages)
void multirateRungeKutta (
const dgConservationLaw &claw,
dgGroupCollection &groups,
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index e9a3252442bd1b2a70eba3f4d432c522b0b08072..c897f2788733c2a66fb71c4a501602235475f5e5 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -700,8 +700,8 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
}
}
}
-
-
+
+
Msg::Info("%d groups of elements",localElements.size());
// do a group of elements for every element type in the mesh
for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
@@ -838,6 +838,248 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
delete []shiftRecv;
}
+
+void dgGroupCollection::buildGroups(GModel *model, int dim, int order,std::string groupType)
+{
+ _model = model;
+ 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::vector<GEntity*> entities;
+ model->getEntities(entities);
+ std::map<int, std::vector<MElement *> >localElements;
+ std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
+ int nlocal=0;
+ int nghosts=0;
+ std::multimap<MElement*, short> &ghostsCells = model->getGhostCells();
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity *entity = entities[i];
+ if(entity->dim() == dim-1){
+ for(unsigned int j = 0; j < entity->physicals.size(); j++){
+ const std::string physicalName = model->getPhysicalName(entity->dim(), entity->physicals[j]);
+ for (int k = 0; k < entity->getNumMeshElements(); k++) {
+ MElement *element = entity->getMeshElement(k);
+ switch(dim) {
+ case 1:
+ boundaryVertices[physicalName].insert( element->getVertex(0) );
+ break;
+ case 2:
+ boundaryEdges[physicalName].insert( element->getEdge(0) );
+ break;
+ case 3:
+ boundaryFaces[physicalName].insert( element->getFace(0));
+ break;
+ default :
+ throw;
+ }
+ }
+ }
+ }else if(entity->dim() == dim){
+ double max_edgeM=-1.e22;
+ double max_edgem=1.e22;
+ double min_edgeM=-1.e22;
+ double min_edgem=1.e22;
+
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+ MElement *el=entity->getMeshElement(iel);
+ if (el->maxEdge()>max_edgeM) {
+ max_edgeM=el->maxEdge();
+ }
+ if (el->maxEdge()<max_edgem) {
+ max_edgem=el->maxEdge();
+ }
+ if (el->minEdge()>min_edgeM) {
+ min_edgeM=el->minEdge();
+ }
+ if (el->minEdge()<min_edgeM) {
+ min_edgem=el->minEdge();
+ }
+ }
+ //printf("\n a1=%f, a2=%f \n",(min_edgeM+min_edgem)/2.,(max_edgeM+max_edgem)/2.);
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+ MElement *el=entity->getMeshElement(iel);
+ if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
+
+ if (groupType=="minEdge") {
+ if (el->minEdge()> (min_edgeM+min_edgem)/2.) {
+ localElements[0].push_back(el);
+ }
+ else{
+ localElements[1].push_back(el);
+ }
+ }
+ else if(groupType=="maxEdge"){
+ if (el->maxEdge()> (max_edgeM+max_edgem)/2.) {
+ localElements[0].push_back(el);
+ }
+ else{
+ localElements[1].push_back(el);
+ }
+
+ }
+ else{
+ localElements[el->getType()].push_back(el);
+ }
+
+
+ nlocal++;
+ }else{
+ std::multimap<MElement*, short>::iterator ghost=ghostsCells.lower_bound(el);
+ while(ghost!= ghostsCells.end() && ghost->first==el){
+ if(abs(ghost->second)-1==Msg::GetCommRank()){
+ ghostElements[el->getPartition()-1][el->getType()].push_back(el);
+ nghosts+=1;
+ }
+ ghost++;
+ }
+ }
+ }
+ }
+ }
+
+
+ Msg::Info("%d groups of elements",localElements.size());
+ for (int n=0; n<localElements.size(); n++) {
+ printf("\n **** Group %d: %d elements **** \n", n,localElements[n].size());
+ }
+ // do a group of elements for every element type in the mesh
+ for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
+ _elementGroups.push_back(new dgGroupOfElements(it->second,order));
+ }
+
+
+ for (int i=0;i<_elementGroups.size();i++){
+ // create a group of faces on the boundary for every group ef elements
+ switch(dim) {
+ case 1 : {
+ std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
+ for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+ if (gof->getNbElements())
+ _boundaryGroups.push_back(gof);
+ else
+ delete gof;
+ break;
+ }
+ }
+ case 2 : {
+ std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
+ for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
+ dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+ if(gof->getNbElements())
+ _boundaryGroups.push_back(gof);
+ else
+ delete gof;
+ }
+ break;
+ }
+ case 3 : {
+ std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
+ for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.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
+ _faceGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],order));
+ // create a groupe of d
+ for (int j=i+1;j<_elementGroups.size();j++){
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
+ 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){
+ _ghostGroups.push_back(new dgGroupOfElements(it->second,order,i));
+ }
+ }
+ //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);
+ if (gof->getNbElements())
+ _faceGroups.push_back(gof);
+ else
+ delete gof;
+ }
+ }
+
+ // build the ghosts structure
+ int *nGhostElements = new int[Msg::GetCommSize()];
+ int *nParentElements = new int[Msg::GetCommSize()];
+ for (int i=0;i<Msg::GetCommSize();i++) {
+ nGhostElements[i]=0;
+ }
+ for (size_t i = 0; i< getNbGhostGroups(); i++){
+ dgGroupOfElements *group = getGhostGroup(i);
+ nGhostElements[group->getGhostPartition()] += group->getNbElements();
+ }
+#ifdef HAVE_MPI
+ MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD);
+#else
+ nParentElements[0]=nGhostElements[0];
+#endif
+ int *shiftSend = new int[Msg::GetCommSize()];
+ int *shiftRecv = new int[Msg::GetCommSize()];
+ int *curShiftSend = new int[Msg::GetCommSize()];
+ int totalSend=0,totalRecv=0;
+ for(int i= 0; i<Msg::GetCommSize();i++){
+ shiftSend[i] = (i==0 ? 0 : shiftSend[i-1]+nGhostElements[i-1]);
+ curShiftSend[i] = shiftSend[i];
+ shiftRecv[i] = (i==0 ? 0 : shiftRecv[i-1]+nParentElements[i-1]);
+ totalSend += nGhostElements[i];
+ totalRecv += nParentElements[i];
+ }
+
+ int *idSend = new int[totalSend];
+ int *idRecv = new int[totalRecv];
+ for (int i = 0; i<getNbGhostGroups(); i++){
+ dgGroupOfElements *group = getGhostGroup(i);
+ int part = group->getGhostPartition();
+ for (int j=0; j< group->getNbElements() ; j++){
+ idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
+ }
+ }
+#ifdef HAVE_MPI
+ MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
+#else
+ memcpy(idRecv,idSend,nParentElements[0]*sizeof(int));
+#endif
+ //create a Map elementNum :: group, position in group
+ std::map<int, std::pair<int,int> > elementMap;
+ for(size_t i = 0; i< getNbElementGroups(); i++) {
+ dgGroupOfElements *group = getElementGroup(i);
+ for(int j=0; j< group->getNbElements(); j++){
+ elementMap[group->getElement(j)->getNum()]=std::pair<int,int>(i,j);
+ }
+ }
+ _elementsToSend.resize(Msg::GetCommSize());
+ for(int i = 0; i<Msg::GetCommSize();i++){
+ _elementsToSend[i].resize(nParentElements[i]);
+ for(int j=0;j<nParentElements[i];j++){
+ int num = idRecv[shiftRecv[i]+j];
+ _elementsToSend[i][j]=elementMap[num];
+ }
+ }
+ delete []idRecv;
+ delete []idSend;
+ delete []curShiftSend;
+ delete []shiftSend;
+ delete []shiftRecv;
+}
+
+
+
+
dgGroupCollection::~dgGroupCollection() {
for (int i=0; i< _elementGroups.size(); i++)
delete _elementGroups[i];
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 0b9be6a2d1ac8da1478369df631f5ea7f1b9a657..3452aa31d458ac13bb2b2ba6234999f7c253231e 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -208,6 +208,7 @@ class dgGroupCollection {
inline int getImageElementPositionInGroup(int partId, int i) const {return _elementsToSend[partId][i].second;}
void buildGroups (GModel *model,int dimension, int order);
+ void buildGroups (GModel *model,int dimension, int order, std::string groupType);
~dgGroupCollection();
};
#endif
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index 814cb566b262eaaa0bb1f3fd48610e0380192e78..adcdd9355df4d6aeaf6379a1bbf25386ed61766c 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -14,7 +14,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups
fullMatrix<double> &solright = solution.getGroupProxy(0);
int nbFields =_claw->nbFields();
int totNbElems = solution.getNbElements();
-
+
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
@@ -32,7 +32,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups
fullMatrix<double> TempL, TempR;
TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
-
+
fSize = TempL.size1();
for (int k=0; k< nbFields; ++k){
double AVGL = 0;
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 68cfdb03a501d71d0c225320fb3e4cd5b52aeab8..29af76404b3f27e5afc402f11c5cb5dd14656ecf 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -58,6 +58,9 @@ void dgSystemOfEquations::registerBindings(binding *b) {
cm->setDescription("set the conservation law this system solves");
cm = cb->addMethod("setup",&dgSystemOfEquations::setup);
cm->setDescription("allocate and init internal stuff, call this function after setOrder and setLaw and before anything else on this instance");
+ cm = cb->addMethod("createGroups",&dgSystemOfEquations::createGroups);
+ cm->setArgNames("groupType",NULL);
+ cm->setDescription("allocate and init internal stuff, creates groups form criterion groupType");
cm = cb->addMethod("exportSolution",&dgSystemOfEquations::exportSolution);
cm->setArgNames("filename",NULL);
cm->setDescription("Print the solution into a file. This file does not contain the mesh. To visualize the solution in gmsh you have to open the mesh file first.");
@@ -101,6 +104,16 @@ void dgSystemOfEquations::L2Projection (std::string functionName){
}
}
+// dgSystemOfEquations::setup() + build groups with criterion:
+// - default: elementType
+// - minEdge (based on minimum edges of elements) , for the moment only two groups possible
+// - maxEdge (based on maximum edges of elements) , for the moment only two groups possible
+void dgSystemOfEquations::createGroups(std::string groupType){
+ _groups.buildGroups(_gm,_dimension,_order,groupType);
+ _solution = new dgDofContainer(_groups,_claw->nbFields());
+ _rightHandSide = new dgDofContainer(_groups,_claw->nbFields());
+ }
+
// ok, we can setup the groups and create solution vectors
void dgSystemOfEquations::setup(){
if (!_claw) throw;
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index fdb887dc933e4232f9b07dd945654a6c8c44e914..85625e8eea32ed15ae837e997ded85d6f461fe62 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -42,6 +42,7 @@ public:
void setOrder (int order); // set the polynomial order
void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
dgSystemOfEquations(GModel *_gm);
+ void createGroups(std::string groupType); // create groups from criterion: minEdge (2 groups), maxEdge (2 groups), elementType, allocate (same as setup())
void setup (); // setup the groups and allocate
void exportSolution (std::string filename); // export the solution
void limitSolution (); // apply the limiter on the solution