diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index a0346fb8a74eacc80d33463be5c01eb7dfac23e7..085b5d20dcf58427181afd46d52bbb98ac140624 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -126,3 +126,4 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
+
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 420b6b667fcde6db3e579738ae4ea91815fda80c..30c269ab42fc75f257b5f0ad4a11c6a12c0463ef 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -117,44 +117,7 @@ dgGroupOfElements::~dgGroupOfElements(){
delete _imass;
}
-// dgGroupOfFaces
-void dgGroupOfFaces::createFaceElements (const std::vector<MFace> &topo_faces){
-
- _faces.resize(topo_faces.size());
- // compute all closures
- for (int i=0;i<topo_faces.size();i++){
- // compute closures for the interpolation
- int ithFace, sign, rot;
- getElementLeft(i)->getFaceInfo (topo_faces[i], ithFace, sign, rot);
- _closuresLeft.push_back(&(_fsLeft->getFaceClosure(ithFace, sign, rot)));
- getElementRight(i)->getFaceInfo (topo_faces[i], ithFace, sign, rot);
- _closuresRight.push_back(&(_fsRight->getFaceClosure(ithFace, sign, rot)));
- // 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 = getElementRight(i)->getFunctionSpace()->getFaceClosure(ithFace, sign, rot);
- for (int j=0; j<geomClosure.size() ; j++){
- int iNod = geomClosure[j];
- MVertex *v = getElementLeft(i)->getVertex(iNod);
- _vertices.push_back(v);
- }
- // triangular face
- if (topo_faces[i].getNumVertices() == 3){
- switch(_vertices.size()){
- case 3 : _faces.push_back(new MTriangle (_vertices) ); break;
- case 6 : _faces.push_back(new MTriangle6 (_vertices) ); break;
- case 10 : _faces.push_back(new MTriangleN (_vertices, 3) ); break;
- case 15 : _faces.push_back(new MTriangleN (_vertices, 4) ); break;
- case 21 : _faces.push_back(new MTriangleN (_vertices, 5) ); break;
- default : throw;
- }
- }
- // quad face 2 do
- else throw;
- }
-}
-// dgGroupOfFaces
void dgGroupOfFaces::computeFaceNormals () {
double g[256][3];
_normals = new fullMatrix<double> (_fsFace->points.size1()*_faces.size(),3);
@@ -182,31 +145,64 @@ void dgGroupOfFaces::computeFaceNormals () {
}
}
-void dgGroupOfFaces::createEdgeElements (const std::vector<MEdge> &topo_edges){
-
- _faces.resize(topo_edges.size());
+void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
// compute all closures
- for (int i=0;i<topo_edges.size();i++){
- int ithEdge, sign;
- getElementLeft(i)->getEdgeInfo (topo_edges[i], ithEdge, sign);
- _closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
- getElementRight(i)->getEdgeInfo (topo_edges[i], ithEdge, sign);
- _closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));
- // get the vertices of the edge
- const std::vector<int> & geomClosure = getElementRight(i)->getFunctionSpace()->getEdgeClosure(ithEdge, sign);
- std::vector<MVertex*> _vertices;
- for (int j=0; j<geomClosure.size() ; j++){
- int iNod = geomClosure[j];
- MVertex *v = getElementLeft(i)->getVertex(iNod);
- _vertices.push_back(v);
- }
- switch(_vertices.size()){
- case 2 : _faces.push_back(new MLine (_vertices) ); break;
- case 3 : _faces.push_back(new MLine3 (_vertices) ); break;
- default : _faces.push_back(new MLineN (_vertices) ); break;
+ // compute closures for the interpolation
+ _left.push_back(iElLeft);
+ _right.push_back(iElRight);
+ MElement &elRight = *_groupRight.getElement(iElRight);
+ MElement &elLeft = *_groupLeft.getElement(iElLeft);
+ int ithFace, sign, rot;
+ elLeft.getFaceInfo (topoFace, ithFace, sign, rot);
+ _closuresLeft.push_back(&(_fsLeft->getFaceClosure(ithFace, sign, rot)));
+ elRight.getFaceInfo (topoFace, ithFace, sign, rot);
+ _closuresRight.push_back(&(_fsRight->getFaceClosure(ithFace, sign, rot)));
+ // compute the face element that correspond to the geometrical closure
+ // get the vertices of the face
+ std::vector<MVertex*> vertices;
+ for(int j=0;j<topoFace.getNumVertices();j++)
+ vertices.push_back(topoFace.getVertex(j));
+ const std::vector<int> & geomClosure = elRight.getFunctionSpace()->getFaceClosure(ithFace, sign, rot);
+ for (int j=0; j<geomClosure.size() ; j++)
+ vertices.push_back( elRight.getVertex(geomClosure[j]) );
+ // triangular face
+ if (topoFace.getNumVertices() == 3){
+ switch(vertices.size()){
+ case 3 : _faces.push_back(new MTriangle (vertices) ); break;
+ case 6 : _faces.push_back(new MTriangle6 (vertices) ); break;
+ case 10 : _faces.push_back(new MTriangleN (vertices, 3) ); break;
+ case 15 : _faces.push_back(new MTriangleN (vertices, 4) ); break;
+ case 21 : _faces.push_back(new MTriangleN (vertices, 5) ); break;
+ default : throw;
}
}
+ // quad face 2 do
+ else throw;
+}
+
+void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight){
+ _left.push_back(iElLeft);
+ _right.push_back(iElRight);
+ MElement &elRight = *_groupRight.getElement(iElRight);
+ MElement &elLeft = *_groupLeft.getElement(iElLeft);
+ int ithEdge, sign;
+ elLeft.getEdgeInfo (topoEdge, ithEdge, sign);
+ _closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
+ elRight.getEdgeInfo (topoEdge, ithEdge, sign);
+ _closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));
+ const std::vector<int> &geomClosure = elRight.getFunctionSpace()->getEdgeClosure(ithEdge, sign);
+ std::vector<MVertex*> vertices;
+ for(int j=0;j<topoEdge.getNumVertices();j++)
+ vertices.push_back(topoEdge.getVertex(j));
+ for (int j=0; j<geomClosure.size() ; j++)
+ vertices.push_back( elRight.getVertex(geomClosure[j]) );
+ switch(vertices.size()){
+ case 2 : _faces.push_back(new MLine (vertices) ); break;
+ case 3 : _faces.push_back(new MLine3 (vertices) ); break;
+ default : _faces.push_back(new MLineN (vertices) ); break;
+ }
}
+
void dgGroupOfFaces::init(int pOrder) {
_fsFace = _faces[0]->getFunctionSpace (pOrder);
_integration=dgGetIntegrationRule (_faces[0],pOrder);
@@ -223,29 +219,83 @@ void dgGroupOfFaces::init(int pOrder) {
}
}
-dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces,
- const std::vector<dgGroupOfElements*> group_list,
- const std::vector<std::pair<int,int> > &l,
- const std::vector<std::pair<int,int> > &r,
- int pOrder) : _group_list(group_list),_left(l), _right(r),
- _fsLeft(getElementLeft(0)->getFunctionSpace (pOrder)), _fsRight (getElementRight(0)->getFunctionSpace (pOrder))
+dgGroupOfFaces::~dgGroupOfFaces()
{
- createFaceElements (topo_faces);
- init(pOrder);
+
}
-dgGroupOfFaces::dgGroupOfFaces (const std::vector<MEdge> &topo_edges,
- const std::vector<dgGroupOfElements*> group_list,
- const std::vector<std::pair<int,int> > &l,
- const std::vector<std::pair<int,int> > &r,
- int pOrder) : _group_list(group_list),_left(l), _right(r),
- _fsLeft(getElementLeft(0)->getFunctionSpace (pOrder)), _fsRight (getElementRight(0)->getFunctionSpace (pOrder))
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
+ _groupLeft(elGroup),_groupRight(elGroup)
{
- createEdgeElements (topo_edges);
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
+ switch (_groupLeft.getElement(0)->getDim()) {
+ case 2 : {
+ std::map<MEdge,int,Less_Edge> edgeMap;
+ for(int i=0; i<elGroup.getNbElements(); i++){
+ MElement &el = *elGroup.getElement(i);
+ for (int j=0; j<el.getNumEdges(); j++){
+ MEdge edge = el.getEdge(j);
+ if(edgeMap.find(edge) == edgeMap.end()){
+ edgeMap[edge] = i;
+ }else{
+ addEdge(edge,edgeMap[edge],i);
+ }
+ }
+ }
+ break;
+ }
+ case 3 : {
+ std::map<MFace,int,Less_Face> faceMap;
+ for(int i=0; i<elGroup.getNbElements(); i++){
+ MElement &el = *elGroup.getElement(i);
+ for (int j=0; j<el.getNumFaces(); j++){
+ MFace face = el.getFace(j);
+ if(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = i;
+ }else{
+ addFace(face,faceMap[face],i);
+ }
+ }
+ }
+ break;
+ }
+ default : throw;
+ }
init(pOrder);
}
-dgGroupOfFaces::~dgGroupOfFaces()
+void dgGroupOfFaces::mapToInterface ( int nFields,
+ const fullMatrix<double> &vLeft,
+ const fullMatrix<double> &vRight,
+ fullMatrix<double> &v)
{
-
+ for(int i=0; i<getNbElements(); i++) {
+ const std::vector<int> &closureRight = *getClosureRight(i);
+ const std::vector<int> &closureLeft = *getClosureLeft(i);
+ for (int iField=0; iField<nFields; iField++){
+ for(size_t j =0; j < closureLeft.size(); j++)
+ v(j, i*2*nFields + iField) = vLeft(closureLeft[j], _left[i]*nFields + iField);
+ for(size_t j =0; j < closureRight.size(); j++)
+ v(j, (i*2+1)*nFields + iField) = vRight(closureRight[j], _right[i]*nFields + iField);
+ }
+ }
+}
+
+void dgGroupOfFaces::mapFromInterface ( int nFields,
+ const fullMatrix<double> &v,
+ fullMatrix<double> &vLeft,
+ fullMatrix<double> &vRight
+ )
+{
+ for(int i=0; i<getNbElements(); i++) {
+ const std::vector<int> &closureRight = *getClosureRight(i);
+ const std::vector<int> &closureLeft = *getClosureLeft(i);
+ for (int iField=0; iField<nFields; iField++){
+ for(size_t j =0; j < closureLeft.size(); j++)
+ vLeft(closureLeft[j], _left[i]*nFields + iField) = v(j, i*2*nFields + iField);
+ for(size_t j =0; j < closureRight.size(); j++)
+ vRight(closureRight[j], _right[i]*nFields + iField) = v(j, (i*2+1)*nFields + iField);
+ }
+ }
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 2b34680d6c53ce28713439131345d9230f42942c..6277359e3412fb8dc8addcff47660360c252c59a 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -104,15 +104,15 @@ public:
};
class dgGroupOfFaces {
- std::vector<dgGroupOfElements*> _group_list;
- void createFaceElements (const std::vector<MFace> &topo_faces);
- void createEdgeElements (const std::vector<MEdge> &topo_faces);
+ const dgGroupOfElements &_groupLeft,&_groupRight;
+ void addFace(const MFace &topoFace, int iElLeft, int iElRight);
+ void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
void computeFaceNormals();
// Two functionSpaces for left and right elements
// the group has always the same types for left and right
const functionSpace *_fsLeft,*_fsRight, *_fsFace;
// N elements in the group
- std::vector<std::pair<int,int> >_left, _right;
+ std::vector<int>_left, _right;
std::vector<MElement *>_faces;
// Ni integration points, matrix of size Ni x 3 (u,v,weight)
fullMatrix<double> *_integration;
@@ -135,21 +135,12 @@ class dgGroupOfFaces {
//common part of the 3 constructors
void init(int pOrder);
public:
- inline MElement* getElementLeft (int i) const {return _group_list[_left[i].first]->getElement(_left[i].second);}
- inline MElement* getElementRight (int i) const {return _group_list[_right[i].first]->getElement(_right[i].second);}
+ inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
+ inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
const std::vector<int> * getClosureLeft(int iFace) const{ return _closuresLeft[iFace];}
const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
- dgGroupOfFaces (const std::vector<MFace> &faces,
- const std::vector<dgGroupOfElements*> group_list,
- const std::vector<std::pair<int,int> > &l,
- const std::vector<std::pair<int,int> > &r,
- int pOrder);
- dgGroupOfFaces (const std::vector<MEdge> &edges,
- const std::vector<dgGroupOfElements*> group_list,
- const std::vector<std::pair<int,int> > &l,
- const std::vector<std::pair<int,int> > &r,
- int pOrder);
+ dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
virtual ~dgGroupOfFaces ();
//this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();}
@@ -159,6 +150,9 @@ public:
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iElement, iGaussPoint);}
+ //keep this outside the Algorithm because this is the only place where data overlap
+ void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
+ void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
};
diff --git a/Solver/dgMainTest.cc b/Solver/dgMainTest.cc
index cc0ea70b0f7c678f7c18e069dfc7577b70d902a6..a0ca1d73709dea69ec2b6f8142d5a15d90c56712 100644
--- a/Solver/dgMainTest.cc
+++ b/Solver/dgMainTest.cc
@@ -9,32 +9,34 @@
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v);
-void buildGroupAllTri(GModel *model, int order, //in
- std::vector<dgGroupOfElements*> &elements, //out
- std::vector<dgGroupOfFaces*> &faces); //out
+std::vector<MElement *> getAllTri(GModel *model);
int main(int argc, char **argv){
GmshMergeFile("input/mesh1.msh");
- std::vector<dgGroupOfElements*> elements;
- std::vector<dgGroupOfFaces*> faces;
- buildGroupAllTri(GModel::current(),1,elements,faces);
- int nbNodes=elements[0]->getNbNodes();
- fullMatrix<double> sol(nbNodes,elements[0]->getNbElements());
- fullMatrix<double> residu(nbNodes,elements[0]->getNbElements());
+ std::vector<MElement *> allTri=getAllTri(GModel::current());
+ int order=1;
+ dgGroupOfElements elements(allTri,order);
+ dgGroupOfFaces faces(elements,order);
+ int nbNodes=elements.getNbNodes();
+ fullMatrix<double> sol(nbNodes,elements.getNbElements());
+ fullMatrix<double> solInterface(nbNodes,faces.getNbElements()*2);
+ fullMatrix<double> residu(nbNodes,elements.getNbElements());
+ fullMatrix<double> residuInterface(nbNodes,faces.getNbElements()*2);
dgAlgorithm algo;
dgConservationLaw *law = dgNewConservationLawAdvection();
- algo.residualVolume(*law,*elements[0],sol,residu);
- for(int i=0;i<elements[0]->getNbElements();i++) {
+ algo.residualVolume(*law,elements,sol,residu);
+ faces.mapToInterface(1, sol, sol, solInterface);
+ algo.residualInterface(*law,faces,solInterface,residuInterface);
+ faces.mapFromInterface(1, residuInterface, residu, residu);
+ for(int i=0;i<elements.getNbElements();i++) {
fullMatrix<double> residuEl(residu,i,1);
fullMatrix<double> solEl(sol,i,1);
- solEl.gemm(elements[0]->getInverseMassMatrix(i),residuEl);
+ solEl.gemm(elements.getInverseMassMatrix(i),residuEl);
}
- print("test.pos",*elements[0],&sol(0,0));
+ print("test.pos",elements,&sol(0,0));
}
-void buildGroupAllTri(GModel *model, int order, //in
- std::vector<dgGroupOfElements*> &elements, //out
- std::vector<dgGroupOfFaces*> &faces){ //out
+std::vector<MElement *> getAllTri(GModel *model){
std::vector<GEntity*> entities;
model->getEntities(entities);
std::vector<MElement *> all_tri;
@@ -43,7 +45,7 @@ void buildGroupAllTri(GModel *model, int order, //in
for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
all_tri.push_back((*itent)->getMeshElement(iel));
}
- elements.push_back(new dgGroupOfElements(all_tri,order));
+ return all_tri;
}
void print (const char *filename,const dgGroupOfElements &els, double *v) {