Skip to content
Snippets Groups Projects
Commit 9cf1dfbb authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg interface term

parent bed7021b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/dgAlgorithm.cpp
+ 1
− 0
View file @ 9cf1dfbb
...@@ -126,3 +126,4 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may ...@@ -126,3 +126,4 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
// ----- 3 ---- do the redistribution at face nodes using BLAS3 // ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP); residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
} }
Solver/dgGroupOfElements.cpp
+ 125
− 75
View file @ 9cf1dfbb
...@@ -117,44 +117,7 @@ dgGroupOfElements::~dgGroupOfElements(){ ...@@ -117,44 +117,7 @@ dgGroupOfElements::~dgGroupOfElements(){
delete _imass; 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 () { void dgGroupOfFaces::computeFaceNormals () {
double g[256][3]; double g[256][3];
_normals = new fullMatrix<double> (_fsFace->points.size1()*_faces.size(),3); _normals = new fullMatrix<double> (_fsFace->points.size1()*_faces.size(),3);
...@@ -182,31 +145,64 @@ void dgGroupOfFaces::computeFaceNormals () { ...@@ -182,31 +145,64 @@ void dgGroupOfFaces::computeFaceNormals () {
} }
} }
void dgGroupOfFaces::createEdgeElements (const std::vector<MEdge> &topo_edges){ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
_faces.resize(topo_edges.size());
// compute all closures // compute all closures
for (int i=0;i<topo_edges.size();i++){ // 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; int ithEdge, sign;
getElementLeft(i)->getEdgeInfo (topo_edges[i], ithEdge, sign); elLeft.getEdgeInfo (topoEdge, ithEdge, sign);
_closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign))); _closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
getElementRight(i)->getEdgeInfo (topo_edges[i], ithEdge, sign); elRight.getEdgeInfo (topoEdge, ithEdge, sign);
_closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign))); _closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));
// get the vertices of the edge const std::vector<int> &geomClosure = elRight.getFunctionSpace()->getEdgeClosure(ithEdge, sign);
const std::vector<int> & geomClosure = getElementRight(i)->getFunctionSpace()->getEdgeClosure(ithEdge, sign); std::vector<MVertex*> vertices;
std::vector<MVertex*> _vertices; for(int j=0;j<topoEdge.getNumVertices();j++)
for (int j=0; j<geomClosure.size() ; j++){ vertices.push_back(topoEdge.getVertex(j));
int iNod = geomClosure[j]; for (int j=0; j<geomClosure.size() ; j++)
MVertex *v = getElementLeft(i)->getVertex(iNod); vertices.push_back( elRight.getVertex(geomClosure[j]) );
_vertices.push_back(v); switch(vertices.size()){
} case 2 : _faces.push_back(new MLine (vertices) ); break;
switch(_vertices.size()){ case 3 : _faces.push_back(new MLine3 (vertices) ); break;
case 2 : _faces.push_back(new MLine (_vertices) ); break; default : _faces.push_back(new MLineN (vertices) ); break;
case 3 : _faces.push_back(new MLine3 (_vertices) ); break;
default : _faces.push_back(new MLineN (_vertices) ); break;
}
} }
} }
void dgGroupOfFaces::init(int pOrder) { void dgGroupOfFaces::init(int pOrder) {
_fsFace = _faces[0]->getFunctionSpace (pOrder); _fsFace = _faces[0]->getFunctionSpace (pOrder);
_integration=dgGetIntegrationRule (_faces[0],pOrder); _integration=dgGetIntegrationRule (_faces[0],pOrder);
...@@ -223,29 +219,83 @@ void dgGroupOfFaces::init(int pOrder) { ...@@ -223,29 +219,83 @@ void dgGroupOfFaces::init(int pOrder) {
} }
} }
dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces, dgGroupOfFaces::~dgGroupOfFaces()
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))
{ {
createFaceElements (topo_faces);
init(pOrder);
} }
dgGroupOfFaces::dgGroupOfFaces (const std::vector<MEdge> &topo_edges, dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
const std::vector<dgGroupOfElements*> group_list, _groupLeft(elGroup),_groupRight(elGroup)
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))
{ {
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); 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);
}
}
} }
Solver/dgGroupOfElements.h
+ 10
− 16
View file @ 9cf1dfbb
...@@ -104,15 +104,15 @@ public: ...@@ -104,15 +104,15 @@ public:
}; };
class dgGroupOfFaces { class dgGroupOfFaces {
std::vector<dgGroupOfElements*> _group_list; const dgGroupOfElements &_groupLeft,&_groupRight;
void createFaceElements (const std::vector<MFace> &topo_faces); void addFace(const MFace &topoFace, int iElLeft, int iElRight);
void createEdgeElements (const std::vector<MEdge> &topo_faces); void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
void computeFaceNormals(); void computeFaceNormals();
// Two functionSpaces for left and right elements // Two functionSpaces for left and right elements
// the group has always the same types for left and right // the group has always the same types for left and right
const functionSpace *_fsLeft,*_fsRight, *_fsFace; const functionSpace *_fsLeft,*_fsRight, *_fsFace;
// N elements in the group // N elements in the group
std::vector<std::pair<int,int> >_left, _right; std::vector<int>_left, _right;
std::vector<MElement *>_faces; std::vector<MElement *>_faces;
// Ni integration points, matrix of size Ni x 3 (u,v,weight) // Ni integration points, matrix of size Ni x 3 (u,v,weight)
fullMatrix<double> *_integration; fullMatrix<double> *_integration;
...@@ -135,21 +135,12 @@ class dgGroupOfFaces { ...@@ -135,21 +135,12 @@ class dgGroupOfFaces {
//common part of the 3 constructors //common part of the 3 constructors
void init(int pOrder); void init(int pOrder);
public: public:
inline MElement* getElementLeft (int i) const {return _group_list[_left[i].first]->getElement(_left[i].second);} inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
inline MElement* getElementRight (int i) const {return _group_list[_right[i].first]->getElement(_right[i].second);} inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
inline MElement* getFace (int iElement) const {return _faces[iElement];} inline MElement* getFace (int iElement) const {return _faces[iElement];}
const std::vector<int> * getClosureLeft(int iFace) const{ return _closuresLeft[iFace];} const std::vector<int> * getClosureLeft(int iFace) const{ return _closuresLeft[iFace];}
const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];} const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
dgGroupOfFaces (const std::vector<MFace> &faces, dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
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);
virtual ~dgGroupOfFaces (); virtual ~dgGroupOfFaces ();
//this part is common with dgGroupOfElements, we should try polymorphism //this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();} inline int getNbElements() const {return _faces.size();}
...@@ -159,6 +150,9 @@ public: ...@@ -159,6 +150,9 @@ public:
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;} inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;} inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iElement, iGaussPoint);} 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);
}; };
......
Solver/dgMainTest.cc
+ 19
− 17
View file @ 9cf1dfbb
...@@ -9,32 +9,34 @@ ...@@ -9,32 +9,34 @@
#include "MElement.h" #include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v); void print (const char *filename,const dgGroupOfElements &els, double *v);
void buildGroupAllTri(GModel *model, int order, //in std::vector<MElement *> getAllTri(GModel *model);
std::vector<dgGroupOfElements*> &elements, //out
std::vector<dgGroupOfFaces*> &faces); //out
int main(int argc, char **argv){ int main(int argc, char **argv){
GmshMergeFile("input/mesh1.msh"); GmshMergeFile("input/mesh1.msh");
std::vector<dgGroupOfElements*> elements; std::vector<MElement *> allTri=getAllTri(GModel::current());
std::vector<dgGroupOfFaces*> faces; int order=1;
buildGroupAllTri(GModel::current(),1,elements,faces); dgGroupOfElements elements(allTri,order);
int nbNodes=elements[0]->getNbNodes(); dgGroupOfFaces faces(elements,order);
fullMatrix<double> sol(nbNodes,elements[0]->getNbElements()); int nbNodes=elements.getNbNodes();
fullMatrix<double> residu(nbNodes,elements[0]->getNbElements()); 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; dgAlgorithm algo;
dgConservationLaw *law = dgNewConservationLawAdvection(); dgConservationLaw *law = dgNewConservationLawAdvection();
algo.residualVolume(*law,*elements[0],sol,residu); algo.residualVolume(*law,elements,sol,residu);
for(int i=0;i<elements[0]->getNbElements();i++) { 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> residuEl(residu,i,1);
fullMatrix<double> solEl(sol,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<MElement *> getAllTri(GModel *model){
std::vector<dgGroupOfElements*> &elements, //out
std::vector<dgGroupOfFaces*> &faces){ //out
std::vector<GEntity*> entities; std::vector<GEntity*> entities;
model->getEntities(entities); model->getEntities(entities);
std::vector<MElement *> all_tri; std::vector<MElement *> all_tri;
...@@ -43,7 +45,7 @@ void buildGroupAllTri(GModel *model, int order, //in ...@@ -43,7 +45,7 @@ void buildGroupAllTri(GModel *model, int order, //in
for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++) for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
all_tri.push_back((*itent)->getMeshElement(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) { void print (const char *filename,const dgGroupOfElements &els, double *v) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment