From ae69fdf5abe0438a00cc4cc2ca9413587a2a7f6a Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 15 Jan 2010 12:41:31 +0000
Subject: [PATCH] add dgGroupCollection class
---
Solver/dgAlgorithm.cpp | 317 +++++++++------------------------
Solver/dgAlgorithm.h | 17 +-
Solver/dgGroupOfElements.cpp | 157 ++++++++++++++++
Solver/dgGroupOfElements.h | 21 +++
Solver/dgLimiter.h | 10 +-
Solver/dgSlopeLimiter.cpp | 11 +-
Solver/dgSystemOfEquations.cpp | 165 +++++++----------
Solver/dgSystemOfEquations.h | 20 +--
8 files changed, 345 insertions(+), 373 deletions(-)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index c05c3e8571..498f37b4bb 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -255,10 +255,7 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation law
- std::vector<dgGroupOfElements*> &eGroups, // group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ dgGroupCollection &groups,
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
@@ -266,68 +263,66 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
dgLimiter *limiter,
int orderRK
) // order of RK integrator
- {
+{
- // U_{n+1}=U_n+h*(SUM_i a_i*K_i)
- // K_i=M^(-1)R(U_n+b_i*K_{i-1})
-
- double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
- double b[4] = {0.,h/2.0,h/2.0,h};
-
- if (orderRK == 1){
- a[0] = h;
- }
+ // U_{n+1}=U_n+h*(SUM_i a_i*K_i)
+ // K_i=M^(-1)R(U_n+b_i*K_{i-1})
+
+ double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
+ double b[4] = {0.,h/2.0,h/2.0,h};
+
+ if (orderRK == 1){
+ a[0] = h;
+ }
+
+ // Current updated solution
+ fullMatrix<double> residuEl, KEl;
+ fullMatrix<double> iMassEl;
+
+ int nbFields = claw.nbFields();
+
+ dgDofContainer K (groups,claw);
+ dgDofContainer Unp (groups,claw);
+
+ K._data->scale(0.0);
+ K._data->axpy(*(sol._data));
+ Unp._data->scale(0.0);
+ Unp._data->axpy(*(sol._data));
+
+ for(int j=0; j<orderRK;j++){
+ if(j){
+ K._data->scale(b[j]);
+ K._data->axpy(*(sol._data));
+ }
- // Current updated solution
- fullMatrix<double> residuEl, KEl;
- fullMatrix<double> iMassEl;
-
- int nbFields = claw.nbFields();
-
- dgDofContainer K (eGroups,ghostGroups,claw);
- dgDofContainer Unp (eGroups,ghostGroups,claw);
-
- K._data->scale(0.0);
- K._data->axpy(*(sol._data));
- Unp._data->scale(0.0);
- Unp._data->axpy(*(sol._data));
-
- for(int j=0; j<orderRK;j++){
- if(j){
- K._data->scale(b[j]);
- K._data->axpy(*(sol._data));
- }
-
if (limiter){
- limiter->apply(K, eGroups, fGroups);
+ limiter->apply(K, groups);
}
syst->scatter(&K);
- this->residual(claw,eGroups,fGroups,bGroups,ghostGroups,K._dataProxys,resd._dataProxys);
- K._data->scale(0.);
- for(int k=0;k < eGroups.size();k++) {
- int nbNodes = eGroups[k]->getNbNodes();
- for(int i=0;i<eGroups[k]->getNbElements();i++) {
- residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
- KEl.setAsProxy(*(K._dataProxys[k]),i*nbFields,nbFields);
- iMassEl.setAsProxy(eGroups[k]->getInverseMassMatrix(),i*nbNodes,nbNodes);
- iMassEl.mult(residuEl,KEl);
- }
- }
- Unp._data->axpy(*(K._data),a[j]);
- }
- if (limiter) limiter->apply(Unp, eGroups, fGroups);
- for (int i=0;i<sol._dataSize;i++){
- (*sol._data)(i)=(*Unp._data)(i);
- }
- }
+ this->residual(claw,groups,K._dataProxys,resd._dataProxys);
+ K._data->scale(0.);
+ for(int k=0; k < groups.getNbElementGroups(); k++) {
+ dgGroupOfElements *group = groups.getElementGroup(k);
+ int nbNodes = group->getNbNodes();
+ for(int i=0;i<group->getNbElements();i++) {
+ residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
+ KEl.setAsProxy(*(K._dataProxys[k]),i*nbFields,nbFields);
+ iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
+ }
+ Unp._data->axpy(*(K._data),a[j]);
+ }
+ if (limiter) limiter->apply(Unp, groups);
+ for (int i=0;i<sol._dataSize;i++){
+ (*sol._data)(i)=(*Unp._data)(i);
+ }
+}
void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
- std::vector<dgGroupOfElements*> &eGroups, // group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ dgGroupCollection &groups,
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
@@ -413,11 +408,11 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
dgDofContainer **K;
K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
- K[i]=new dgDofContainer(eGroups,claw);
+ K[i]=new dgDofContainer(groups,claw);
K[i]->_data->scale(0.0);
}
- dgDofContainer Unp (eGroups,claw);
- dgDofContainer tmp (eGroups,claw);
+ dgDofContainer Unp (groups,claw);
+ dgDofContainer tmp (groups,claw);
Unp._data->scale(0.0);
Unp._data->axpy(*(sol._data));
@@ -425,20 +420,21 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
for(int j=0; j<nStages;j++){
tmp._data->scale(0.0);
tmp._data->axpy(*(sol._data));
- for(int k=0;k < eGroups.size();k++) {
+ for(int k=0;k < groups.getNbElementGroups();k++) {
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]);
}
}
}
- this->residual(claw,eGroups,fGroups,bGroups,ghostGroups,tmp._dataProxys,resd._dataProxys);
- for(int k=0;k < eGroups.size();k++) {
- int nbNodes = eGroups[k]->getNbNodes();
- for(int i=0;i<eGroups[k]->getNbElements();i++) {
+ this->residual(claw,groups,tmp._dataProxys,resd._dataProxys);
+ for(int k=0;k < groups.getNbElementGroups(); k++) {
+ dgGroupOfElements *group = groups.getElementGroup(k);
+ int nbNodes = group->getNbNodes();
+ for(int i=0;i<group->getNbElements();i++) {
residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
KEl.setAsProxy(*(K[j]->_dataProxys[k]),i*nbFields,nbFields);
- iMassEl.setAsProxy(eGroups[k]->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
Unp.getGroupProxy(k).add(K[j]->getGroupProxy(k),h*b[j]);
@@ -538,182 +534,33 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
-/*
-void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups){
-}
-
-void dgAlgorithm::buildMandatoryGroups(dgGroupOfElements &eGroup,
- std::vector<dgGroupOfElements*> &partitionedGroups){
-}
-
-void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
- std::vector<dgGroupOfElements*> &partitionedGroups){
-}
-*/
-
-//static void partitionGroup (std::vector<MElement *>,){
-//}
-
-// this function creates groups of elements
-// First, groups of faces are created on every physical group
-// of dimension equal to the dimension of the problem minus one
-// Then, groups of elements are created
-// -) Elements of different types are separated
-// -) Then those groups are split into partitions
-// -) Finally, those groups are re-numbered
-// Finally, group of interfaces are created
-// -) Groups of faces internal to a given group
-// -) Groups of faces between groups.
-
-void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups,
- std::vector<dgGroupOfElements*> &ghostGroups
- )
-{
- 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){
- for (int iel=0; iel<entity->getNumMeshElements(); iel++){
- MElement *el=entity->getMeshElement(iel);
- if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
- 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){
- nghosts+=(abs(ghost->second)-1==Msg::GetCommRank());
- ghostElements[el->getPartition()-1][el->getType()].push_back(el);
- ghost++;
- }
- }
- }
- }
- }
-
-
- 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){
- eGroups.push_back(new dgGroupOfElements(it->second,order));
- }
-
-
- for (int i=0;i<eGroups.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++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
- }
- break;
- }
- case 2 : {
- std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
- for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
- }
- break;
- }
- case 3 : {
- std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
- for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
- }
- break;
- }
- }
- // create a group of faces for every face that is common to elements of the same group
- fGroups.push_back(new dgGroupOfFaces(*eGroups[i],order));
- // create a groupe of d
- for (int j=i+1;j<eGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*eGroups[i],*eGroups[j],order);
- if (gof->getNbElements())
- fGroups.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<eGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*ghostGroups[i],*eGroups[j],order);
- if (gof->getNbElements())
- fGroups.push_back(gof);
- else
- delete gof;
- }
- }
-}
// works for any number of groups
void dgAlgorithm::residual( const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ dgGroupCollection &groups,
std::vector<fullMatrix<double> *> &solution, // solution
std::vector<fullMatrix<double> *> &residu) // residual
{
int nbFields=claw.nbFields();
//volume term
- for(size_t i=0;i<eGroups.size() ; i++) {
+ for(size_t i=0;i<groups.getNbElementGroups() ; i++) {
residu[i]->scale(0);
-
- residualVolume(claw,*eGroups[i],*solution[i],*residu[i]);
+ residualVolume(claw,*groups.getElementGroup(i),*solution[i],*residu[i]);
}
// residu[0]->print("Volume");
//interface term
- for(size_t i=0;i<fGroups.size() ; i++) {
- dgGroupOfFaces &faces = *fGroups[i];
+ for(size_t i=0;i<groups.getNbFaceGroups() ; i++) {
+ dgGroupOfFaces &faces = *groups.getFaceGroup(i);
int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<eGroups.size() ; j++) {
- if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
- if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+ for(size_t j=0;j<groups.getNbElementGroups(); j++) {
+ dgGroupOfElements *groupj = groups.getElementGroup(j);
+ if (groupj == &faces.getGroupLeft()) iGroupLeft = j;
+ if (groupj == &faces.getGroupRight()) iGroupRight= j;
}
- for(size_t j=0;j<ghostGroups.size() ; j++) {
- if (ghostGroups[j] == &faces.getGroupLeft())iGroupLeft = j+eGroups.size();
- if (ghostGroups[j] == &faces.getGroupRight())iGroupRight= j+eGroups.size();
+ for(size_t j=0;j<groups.getNbGhostGroups(); j++) {
+ dgGroupOfElements *groupj = groups.getGhostGroup(j);
+ if (groupj == &faces.getGroupLeft()) iGroupLeft = j;
+ if (groupj == &faces.getGroupLeft())iGroupLeft = j + groups.getNbElementGroups();
+ if (groupj == &faces.getGroupRight())iGroupRight= j + groups.getNbElementGroups();
}
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
@@ -723,12 +570,13 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
}
// residu[0]->print("Interfaces");
//boundaries
- for(size_t i=0;i<bGroups.size() ; i++) {
- dgGroupOfFaces &faces = *bGroups[i];
+ for(size_t i=0;i<groups.getNbBoundaryGroups() ; i++) {
+ dgGroupOfFaces &faces = *groups.getBoundaryGroup(i);
int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<eGroups.size() ; j++) {
- if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
- if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+ for(size_t j=0;j<groups.getNbElementGroups(); j++) {
+ dgGroupOfElements *groupj = groups.getElementGroup(j);
+ if (groupj == &faces.getGroupLeft())iGroupLeft = j;
+ if (groupj == &faces.getGroupRight())iGroupRight= j;
}
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
@@ -783,4 +631,3 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
DT[iElement] = 1./spectralRadius;
}
}
-
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 2ca505a71d..c3dd4aa656 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -4,6 +4,7 @@
#include "fullMatrix.h"
#include <vector>
class GModel;
+class dgGroupCollection;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
@@ -34,20 +35,13 @@ class dgAlgorithm {
fullMatrix<double> &residual // residual !! at faces nodes
);
// works only if there is only 1 group of element
- static void residual( const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ static void residual( const dgConservationLaw &claw, dgGroupCollection &groups,
std::vector<fullMatrix<double> *> &solution, // solution !! at faces nodes
std::vector<fullMatrix<double> *> &residual // residual !! at faces nodes
);
void rungeKutta (
const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ dgGroupCollection &groups,
double h,
dgDofContainer &residu,
dgDofContainer &sol,
@@ -62,10 +56,7 @@ class dgAlgorithm {
);
void multirateRungeKutta (
const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
+ dgGroupCollection &groups,
double h,
dgDofContainer &residu,
dgDofContainer &sol,
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index fc5080aad9..a1867f108e 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -606,3 +606,160 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
}
}
}
+
+/*
+void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups,
+ std::vector<dgGroupOfFaces*> &bGroups){
+}
+
+void dgAlgorithm::buildMandatoryGroups(dgGroupOfElements &eGroup,
+ std::vector<dgGroupOfElements*> &partitionedGroups){
+}
+
+void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
+ std::vector<dgGroupOfElements*> &partitionedGroups){
+}
+*/
+
+//static void partitionGroup (std::vector<MElement *>,){
+//}
+
+// this function creates groups of elements
+// First, groups of faces are created on every physical group
+// of dimension equal to the dimension of the problem minus one
+// Then, groups of elements are created
+// -) Elements of different types are separated
+// -) Then those groups are split into partitions
+// -) Finally, those groups are re-numbered
+// Finally, group of interfaces are created
+// -) Groups of faces internal to a given group
+// -) Groups of faces between groups.
+
+void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
+{
+ _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){
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+ MElement *el=entity->getMeshElement(iel);
+ if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
+ 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){
+ nghosts+=(abs(ghost->second)-1==Msg::GetCommRank());
+ ghostElements[el->getPartition()-1][el->getType()].push_back(el);
+ ghost++;
+ }
+ }
+ }
+ }
+ }
+
+
+ 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){
+ _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++) {
+ _boundaryGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second));
+ }
+ break;
+ }
+ case 2 : {
+ std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
+ for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
+ _boundaryGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second));
+ }
+ break;
+ }
+ case 3 : {
+ std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
+ for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
+ _boundaryGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second));
+ }
+ 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;
+ }
+ }
+ printf("nbGhostGroups : %i nbElementsGroups : %i\n",getNbGhostGroups(), getNbElementGroups());
+}
+
+dgGroupCollection::~dgGroupCollection() {
+ for (int i=0; i< _elementGroups.size(); i++)
+ delete _elementGroups[i];
+ for (int i=0; i< _faceGroups.size(); i++)
+ delete _faceGroups[i];
+ for (int i=0; i< _boundaryGroups.size(); i++)
+ delete _boundaryGroups[i];
+ for (int i=0; i< _ghostGroups.size(); i++)
+ delete _ghostGroups[i];
+}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 61c9925384..7c0f177772 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -18,6 +18,7 @@
class MElement;
class polynomialBasis;
class GEntity;
+class GModel;
class dgElement {
MElement *_element;
@@ -180,4 +181,24 @@ public:
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
};
+
+class dgGroupCollection {
+ GModel *_model;
+ std::vector<dgGroupOfElements*> _elementGroups; //volume
+ std::vector<dgGroupOfFaces*> _faceGroups; //interface
+ std::vector<dgGroupOfFaces*> _boundaryGroups; //boundary
+ std::vector<dgGroupOfElements*> _ghostGroups; //ghost volume
+ public:
+ inline int getNbElementGroups() const {return _elementGroups.size();}
+ inline int getNbFaceGroups() const {return _faceGroups.size();}
+ inline int getNbBoundaryGroups() const {return _boundaryGroups.size();}
+ inline int getNbGhostGroups() const {return _ghostGroups.size();}
+ inline dgGroupOfElements *getElementGroup(int i) const {return _elementGroups[i];}
+ inline dgGroupOfFaces *getFaceGroup(int i) const {return _faceGroups[i];}
+ inline dgGroupOfFaces *getBoundaryGroup(int i) const {return _boundaryGroups[i];}
+ inline dgGroupOfElements *getGhostGroup(int i) const {return _ghostGroups[i];}
+
+ void buildGroups (GModel *model,int dimension, int order);
+ ~dgGroupCollection();
+};
#endif
diff --git a/Solver/dgLimiter.h b/Solver/dgLimiter.h
index 860dd42822..93951c62a2 100644
--- a/Solver/dgLimiter.h
+++ b/Solver/dgLimiter.h
@@ -4,8 +4,7 @@
#include "fullMatrix.h"
#include <vector>
struct dgDofContainer;
-class dgGroupOfElements;
-class dgGroupOfFaces;
+class dgGroupCollection;
class dgConservationLaw;
class dgLimiter{
@@ -13,16 +12,13 @@ protected:
dgConservationLaw *_claw;
public:
dgLimiter (dgConservationLaw *claw) : _claw(claw) {}
- virtual bool apply ( dgDofContainer &sol, std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroup ) = 0;
+ virtual bool apply ( dgDofContainer &sol, dgGroupCollection &groups)=0;
};
class dgSlopeLimiter : public dgLimiter{
public :
dgSlopeLimiter (dgConservationLaw *claw) : dgLimiter (claw) {}
- virtual bool apply ( dgDofContainer &solution,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroup);
+ virtual bool apply ( dgDofContainer &solution, dgGroupCollection &groups);
};
#endif
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index 0be45a3979..814cb566b2 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -4,15 +4,12 @@
#include "function.h"
//----------------------------------------------------------------------------------
-bool dgSlopeLimiter::apply ( dgDofContainer &solution,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups)
+bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups)
{
-
//WARNING: ONLY FOR 1 GROUP OF FACES
//TODO: make this more general
- dgGroupOfFaces* group = fGroups[0];
+ dgGroupOfFaces* group = groups.getFaceGroup(0);
fullMatrix<double> &solleft = solution.getGroupProxy(0);
fullMatrix<double> &solright = solution.getGroupProxy(0);
int nbFields =_claw->nbFields();
@@ -93,8 +90,8 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
}
}
// --- CLIPPING: check unphysical values
- for (int iG = 0; iG < eGroups.size(); iG++){
- dgGroupOfElements* egroup = eGroups[iG];
+ for (int iG = 0; iG < groups.getNbElementGroups(); iG++){
+ dgGroupOfElements* egroup = groups.getElementGroup(iG);
fullMatrix<double> &solGroup = solution.getGroupProxy(iG);
dataCacheMap cacheMap(egroup->getNbNodes());//nbdofs for each element
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 45f7ec2b17..ddd220d488 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -89,31 +89,26 @@ void dgSystemOfEquations::registerBindings(binding *b) {
// do a L2 projection
void dgSystemOfEquations::L2Projection (std::string functionName){
dgConservationLawL2Projection Law(functionName,*_claw);
- for (int i=0;i<_elementGroups.size();i++){
- _algo->residualVolume(Law,*_elementGroups[i],*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
- _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
+ for (int i=0;i<_groups.getNbElementGroups();i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ _algo->residualVolume(Law,*group,*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
+ _algo->multAddInverseMassMatrix(*group,*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
}
}
// ok, we can setup the groups and create solution vectors
void dgSystemOfEquations::setup(){
if (!_claw) throw;
- _algo->buildGroups(_gm,
- _dimension,
- _order,
- _elementGroups,
- _faceGroups,
- _boundaryGroups,
- _ghostGroups
- );
+ _groups.buildGroups(_gm,_dimension,_order);
// 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< _ghostGroups.size(); i++){
- nGhostElements[_ghostGroups[i]->getGhostPartition()] += _ghostGroups[i]->getNbElements();
+ for (size_t i = 0; i< _groups.getNbGhostGroups(); i++){
+ dgGroupOfElements *group=_groups.getGhostGroup(i);
+ nGhostElements[group->getGhostPartition()] += group->getNbElements();
}
#ifdef HAVE_MPI
MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD);
@@ -134,18 +129,20 @@ void dgSystemOfEquations::setup(){
int *idSend = new int[totalSend];
int *idRecv = new int[totalRecv];
- for (size_t i = 0; i< _ghostGroups.size(); i++){
- int part = _ghostGroups[i]->getGhostPartition();
- for (int j=0; j< _ghostGroups[i]->getNbElements() ; j++){
- idSend[curShiftSend[part]++] = _ghostGroups[i]->getElement(j)->getNum();
+ for (int i = 0; i< _groups.getNbGhostGroups(); i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
+ int part = group->getGhostPartition();
+ for (int j=0; j< group->getNbElements() ; j++){
+ idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
}
}
MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
//create a Map elementNum :: group, position in group
std::map<int, std::pair<int,int> > elementMap;
- for(size_t i = 0; i<_elementGroups.size(); i++) {
- for(int j=0; j<_elementGroups[i]->getNbElements(); j++){
- elementMap[_elementGroups[i]->getElement(j)->getNum()]=std::pair<int,int>(i,j);
+ for(size_t i = 0; i< _groups.getNbElementGroups(); i++) {
+ dgGroupOfElements *group = _groups.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());
@@ -158,21 +155,21 @@ void dgSystemOfEquations::setup(){
}
// compute the total size of the soution
delete curShiftSend;
- _solution = new dgDofContainer(_elementGroups,_ghostGroups,*_claw);
- _rightHandSide = new dgDofContainer(_elementGroups,_ghostGroups,*_claw);
+ _solution = new dgDofContainer(_groups,*_claw);
+ _rightHandSide = new dgDofContainer(_groups,*_claw);
}
double dgSystemOfEquations::RK44(double dt){
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide, this);
+ _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide, this);
return _solution->_data->norm();
}
double dgSystemOfEquations::computeInvSpectralRadius(){
double sr = 1.e22;
- for (int i=0;i<_elementGroups.size();i++){
+ for (int i=0;i<_groups.getNbElementGroups();i++){
std::vector<double> DTS;
- _algo->computeElementaryTimeSteps(*_claw, *_elementGroups[i], *_solution->_dataProxys[i], DTS);
+ _algo->computeElementaryTimeSteps(*_claw, *_groups.getElementGroup(i), *_solution->_dataProxys[i], DTS);
for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
}
#ifdef HAVE_MPI
@@ -187,17 +184,17 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter(_claw);
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide,this, sl);
+ _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide,this, sl);
delete sl;
return _solution->_data->norm();
}
double dgSystemOfEquations::ForwardEuler(double dt){
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide,this, NULL,1);
+ _algo->rungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide,this, NULL,1);
return _solution->_data->norm();
}
double dgSystemOfEquations::multirateRK43(double dt){
- _algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide);
+ _algo->multirateRungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
}
@@ -207,15 +204,12 @@ void dgSystemOfEquations::exportSolution(std::string outputFile){
void dgSystemOfEquations::limitSolution(){
dgLimiter *sl = new dgSlopeLimiter(_claw);
- sl->apply(*_solution,_elementGroups,_faceGroups);
+ sl->apply(*_solution,_groups);
delete sl;
}
dgSystemOfEquations::~dgSystemOfEquations(){
- for (int i=0;i<_elementGroups.size();i++){
- delete _elementGroups[i];
- }
if (_solution){
delete _solution;
delete _rightHandSide;
@@ -244,8 +238,8 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
int COUNT = 0;
- for (int i=0;i < _elementGroups.size() ;i++){
- COUNT += _elementGroups[i]->getNbElements();
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ COUNT += _groups.getElementGroup(i)->getNbElements();
}
fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
fprintf(f,"$ElementNodeData\n");
@@ -255,16 +249,17 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
fprintf(f,"0.0\n");
fprintf(f,"3\n");
fprintf(f,"0\n 1\n %d\n",COUNT);
- for (int i=0;i < _elementGroups.size() ;i++){
- for (int iElement=0 ; iElement<_elementGroups[i]->getNbElements() ;++iElement) {
- MElement *e = _elementGroups[i]->getElement(iElement);
- int num = e->getNum();
- fullMatrix<double> sol = getSolutionProxy (i, iElement);
- fprintf(f,"%d %d",num,sol.size1());
- for (int k=0;k<sol.size1();++k) {
- fprintf(f," %12.5E ",sol(k,ICOMP));
- }
- fprintf(f,"\n");
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
+ MElement *e = group->getElement(iElement);
+ int num = e->getNum();
+ fullMatrix<double> sol = getSolutionProxy (i, iElement);
+ fprintf(f,"%d %d",num,sol.size1());
+ for (int k=0;k<sol.size1();++k) {
+ fprintf(f," %12.5E ",sol(k,ICOMP));
+ }
+ fprintf(f,"\n");
}
}
fprintf(f,"$EndElementNodeData\n");
@@ -276,9 +271,10 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
std::map<int, std::vector<double> > data;
- for (int i=0;i < _elementGroups.size() ;i++){
- for (int iElement=0 ; iElement<_elementGroups[i]->getNbElements() ;++iElement) {
- MElement *e = _elementGroups[i]->getElement(iElement);
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
+ MElement *e = group->getElement(iElement);
int num = e->getNum();
fullMatrix<double> sol = getSolutionProxy (i, iElement);
std::vector<double> val;
@@ -297,21 +293,23 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
delete pv;
}
-dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups,std::vector<dgGroupOfElements*> ghostGroups, const dgConservationLaw &claw){
+dgDofContainer::dgDofContainer (dgGroupCollection &groups, const dgConservationLaw &claw):
+ _groups(groups)
+{
_dataSize = 0;
_dataSizeGhost=0;
totalNbElements = 0;
int totalNbElementsGhost =0;
nbFields = claw.nbFields();
- for (int i=0;i<elementGroups.size();i++){
- int nbNodes = elementGroups[i]->getNbNodes();
- int nbElements = elementGroups[i]->getNbElements();
+ for (int i=0; i< _groups.getNbElementGroups();i++){
+ int nbNodes = _groups.getElementGroup(i)->getNbNodes();
+ int nbElements = _groups.getElementGroup(i)->getNbElements();
totalNbElements +=nbElements;
_dataSize += nbNodes*nbFields*nbElements;
}
- for (int i=0;i<ghostGroups.size();i++){
- int nbNodes = ghostGroups[i]->getNbNodes();
- int nbElements = ghostGroups[i]->getNbElements();
+ for (int i=0; i<_groups.getNbGhostGroups(); i++){
+ int nbNodes = _groups.getGhostGroup(i)->getNbNodes();
+ int nbElements = _groups.getGhostGroup(i)->getNbElements();
totalNbElementsGhost +=nbElements;
_dataSizeGhost += nbNodes*nbFields*nbElements;
}
@@ -321,49 +319,22 @@ dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups,s
_ghostData = new fullVector<double>(_dataSizeGhost);
// create proxys for each group
int offset = 0;
- _dataProxys.resize(elementGroups.size()+ghostGroups.size());
- for (int i=0;i<elementGroups.size();i++){
- dgGroupOfElements *group=elementGroups[i];
+ _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
+ for (int i=0;i<_groups.getNbElementGroups();i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
int nbNodes = group->getNbNodes();
int nbElements = group->getNbElements();
_dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, nbFields*nbElements);
offset += nbNodes*nbFields*nbElements;
}
offset=0;
- for (int i=0;i<ghostGroups.size();i++){
- dgGroupOfElements *group=ghostGroups[i];
- int nbNodes = group->getNbNodes();
- int nbElements = group->getNbElements();
- _dataProxys[i+elementGroups.size()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, nbFields*nbElements);
- offset += nbNodes*nbFields*nbElements;
- }
- // printf("datasize = %d\n",_dataSize);
-}
-
-dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups, const dgConservationLaw &claw){
- _dataSize = 0;
- totalNbElements = 0;
- nbFields = claw.nbFields();
- for (int i=0;i<elementGroups.size();i++){
- int nbNodes = elementGroups[i]->getNbNodes();
- int nbElements = elementGroups[i]->getNbElements();
- totalNbElements +=nbElements;
- _dataSize += nbNodes*nbFields*nbElements;
- }
-
- // allocate the big vectors
- _data = new fullVector<double>(_dataSize);
- // create proxys for each group
- int offset = 0;
- _dataProxys.resize(elementGroups.size());
- for (int i=0;i<elementGroups.size();i++){
- dgGroupOfElements *group=elementGroups[i];
+ for (int i=0;i<_groups.getNbGhostGroups();i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
int nbNodes = group->getNbNodes();
int nbElements = group->getNbElements();
- _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, nbFields*nbElements);
+ _dataProxys[i+_groups.getNbElementGroups()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, nbFields*nbElements);
offset += nbNodes*nbFields*nbElements;
}
- // printf("datasize = %d\n",_dataSize);
}
dgDofContainer::~dgDofContainer (){
@@ -382,11 +353,11 @@ void dgSystemOfEquations::scatter(dgDofContainer *vector){
countS[i]=0;
countR[i]=0;
for(size_t j=0;j<_elementsToSend[i].size();j++){
- countS[i] += _elementGroups[_elementsToSend[i][j].first]->getNbNodes()*_claw->nbFields();
+ countS[i] += _groups.getElementGroup(_elementsToSend[i][j].first)->getNbNodes()*_claw->nbFields();
}
}
- for(size_t i=0; i<_ghostGroups.size(); i++){
- dgGroupOfElements *group=_ghostGroups[i];
+ for(size_t i=0; i<_groups.getNbGhostGroups(); i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
countR[group->getGhostPartition()]+=group->getNbElements()*group->getNbNodes()*_claw->nbFields();
}
shiftS[0]=shiftR[0]=0;
@@ -410,7 +381,7 @@ void dgSystemOfEquations::scatter(dgDofContainer *vector){
for(int l=0;l<_claw->nbFields();l++){
for(int k=0;k<sol.size1();k++){
send[index++] = sol(k,eid*_claw->nbFields()+l);
-// send[index++] = _elementGroups[gid]->getElement(eid)->getNum()*9+l*3+k;
+ //send[index++] = _groups.getElementGroup(gid)->getElement(eid)->getNum()*9+l*3+k;
}
}
}
@@ -419,16 +390,16 @@ void dgSystemOfEquations::scatter(dgDofContainer *vector){
//3) send
MPI_Alltoallv(send,countS,shiftS,MPI_DOUBLE,recv,countR,shiftR,MPI_DOUBLE,MPI_COMM_WORLD);
//4) distribute
- for(int i=0; i< _ghostGroups.size();i++){
- fullMatrix<double>&sol = vector->getGroupProxy(i+_elementGroups.size());
- int part = _ghostGroups[i]->getGhostPartition();
+ for(int i=0; i< _groups.getNbGhostGroups();i++){
+ fullMatrix<double>&sol = vector->getGroupProxy(i+_groups.getNbElementGroups());
+ int part = _groups.getGhostGroup(i)->getGhostPartition();
int ndof = sol.size1()*sol.size2();
for(int j=0;j<sol.size2();j++){
for(int k=0;k<sol.size1();k++){
sol(k,j)=recv[shiftR[part]++];
-/* int a=(int)recv[shiftR[part]++];
- int b= _ghostGroups[i]->getElement(j/3)->getNum()*9+3*(j%3)+k;
- if(a!=b)printf("%i %i\n",a,b);*/
+ //int a=(int)recv[shiftR[part]++];
+ //int b= _groups.getGhostGroup(i)->getElement(j/3)->getNum()*9+3*(j%3)+k;
+ ///*if(a!=b)*/printf("%i %i\n",a,b);
}
}
}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index e8f4bcd583..81182a3cfb 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -16,6 +16,7 @@ private:
dgDofContainer (const dgDofContainer&);
int totalNbElements;
int nbFields;
+ dgGroupCollection &_groups;
public:
int _dataSize; // the full data size i.e. concerning all groups (not ghost, see bellow)
int _dataSizeGhost;
@@ -24,16 +25,14 @@ public:
fullVector<double> * _ghostData;
inline int getDataSize(){return _dataSize;}
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
- dgDofContainer (std::vector<dgGroupOfElements*> &groups,std::vector<dgGroupOfElements*> ghostGroups, const dgConservationLaw &claw);
- dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
+ dgDofContainer (dgGroupCollection &groups, const dgConservationLaw &claw);
~dgDofContainer ();
int getNbElements() {return totalNbElements;}
- //collective, should be called on all proc
- void setGhostedGroups (std::vector<dgGroupOfElements*> &ghostGroups);
};
class binding;
+
class dgSystemOfEquations {
//////////////
//parallel section, should be moved to a new class
@@ -42,10 +41,12 @@ class dgSystemOfEquations {
int totalSend, totalRecv;
public :
void scatter(dgDofContainer *solution);
+ std::vector< std::vector<std::pair<int,int> > >_elementsToSend; //{group,id} of the elements to send to each proc
//////////////
private:
// the mesh and the model
GModel *_gm;
+ dgGroupCollection _groups;
// the algorithm that computes DG operators
dgAlgorithm *_algo;
// the conservation law
@@ -58,15 +59,6 @@ class dgSystemOfEquations {
// solution and righ hand sides
dgDofContainer *_solution;
dgDofContainer *_rightHandSide;
- // groups of elements (volume terms)
- std::vector<dgGroupOfElements*> _elementGroups;
- //ghost structure
- std::vector<dgGroupOfElements*> _ghostGroups;
- std::vector< std::vector<std::pair<int,int> > >_elementsToSend; //{group,id} of the elements to send to each proc
- // groups of faces (interface terms)
- std::vector<dgGroupOfFaces*> _faceGroups;
- // groups of faces (boundary conditions)
- std::vector<dgGroupOfFaces*> _boundaryGroups;
dgSystemOfEquations(const dgSystemOfEquations &) {}
double computeTimeStepMethodOfLines () const;
public:
@@ -93,7 +85,7 @@ public:
void saveSolution (std::string fileName) ;
void loadSolution (std::string fileName);
~dgSystemOfEquations();
-private:
};
+
#endif // _DG_SYSTEM_OF_EQUATIONS_
--
GitLab