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

rewrite the construction of dgGroupOfInterface

parent f43d2f41
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Numeric/polynomialBasis.cpp
+ 56
− 32
View file @ d4797a29
......@@ -875,179 +875,203 @@ const polynomialBasis &polynomialBases::find(int tag)
switch (tag){
case MSH_PNT:
F.numFaces = 1;
F.monomials = generate1DMonomials(0);
F.points = generate1DPoints(0);
break;
case MSH_LIN_2 :
F.numFaces = 2;
F.monomials = generate1DMonomials(1);
F.points = generate1DPoints(1);
generate1dVertexClosure(F.vertexClosure);
generate1dVertexClosure(F.closures);
break;
case MSH_LIN_3 :
F.numFaces = 2;
F.monomials = generate1DMonomials(2);
F.points = generate1DPoints(2);
generate1dVertexClosure(F.vertexClosure);
generate1dVertexClosure(F.closures);
break;
case MSH_LIN_4:
F.numFaces = 2;
F.monomials = generate1DMonomials(3);
F.points = generate1DPoints(3);
generate1dVertexClosure(F.vertexClosure);
generate1dVertexClosure(F.closures);
break;
case MSH_LIN_5:
F.numFaces = 2;
F.monomials = generate1DMonomials(4);
F.points = generate1DPoints(4);
generate1dVertexClosure(F.vertexClosure);
generate1dVertexClosure(F.closures);
break;
case MSH_LIN_6:
F.numFaces = 2;
F.monomials = generate1DMonomials(5);
F.points = generate1DPoints(5);
generate1dVertexClosure(F.vertexClosure);
generate1dVertexClosure(F.closures);
break;
case MSH_TRI_3 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(1);
F.points = gmshGeneratePointsTriangle(1, false);
generate2dEdgeClosure(F.edgeClosure, 1);
generate2dEdgeClosure(F.closures, 1);
break;
case MSH_TRI_6 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(2);
F.points = gmshGeneratePointsTriangle(2, false);
generate2dEdgeClosure(F.edgeClosure, 2);
generate2dEdgeClosure(F.closures, 2);
break;
case MSH_TRI_9 :
F.numFaces = 3;
F.monomials = generatePascalSerendipityTriangle(3);
F.points = gmshGeneratePointsTriangle(3, true);
generate2dEdgeClosure(F.edgeClosure, 3);
generate2dEdgeClosure(F.closures, 3);
break;
case MSH_TRI_10 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(3);
F.points = gmshGeneratePointsTriangle(3, false);
generate2dEdgeClosure(F.edgeClosure, 3);
generate2dEdgeClosure(F.closures, 3);
break;
case MSH_TRI_12 :
F.numFaces = 3;
F.monomials = generatePascalSerendipityTriangle(4);
F.points = gmshGeneratePointsTriangle(4, true);
generate2dEdgeClosure(F.edgeClosure, 4);
generate2dEdgeClosure(F.closures, 4);
break;
case MSH_TRI_15 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(4);
F.points = gmshGeneratePointsTriangle(4, false);
generate2dEdgeClosure(F.edgeClosure, 4);
generate2dEdgeClosure(F.closures, 4);
break;
case MSH_TRI_15I :
F.numFaces = 3;
F.monomials = generatePascalSerendipityTriangle(5);
F.points = gmshGeneratePointsTriangle(5, true);
generate2dEdgeClosure(F.edgeClosure, 5);
generate2dEdgeClosure(F.closures, 5);
break;
case MSH_TRI_21 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(5);
F.points = gmshGeneratePointsTriangle(5, false);
generate2dEdgeClosure(F.edgeClosure, 5);
generate2dEdgeClosure(F.closures, 5);
break;
case MSH_TET_4 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
generate3dFaceClosure(F.faceClosure, 1);
generate3dFaceClosure(F.closures, 1);
break;
case MSH_TET_10 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(2);
F.points = gmshGeneratePointsTetrahedron(2, false);
generate3dFaceClosure(F.faceClosure, 2);
generate3dFaceClosure(F.closures, 2);
break;
case MSH_TET_20 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(3);
F.points = gmshGeneratePointsTetrahedron(3, false);
generate3dFaceClosure(F.faceClosure, 3);
generate3dFaceClosure(F.closures, 3);
break;
case MSH_TET_35 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, false);
generate3dFaceClosure(F.faceClosure, 4);
generate3dFaceClosure(F.closures, 4);
break;
case MSH_TET_34 :
F.numFaces = 4;
F.monomials = generatePascalSerendipityTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, true);
generate3dFaceClosure(F.faceClosure, 4);
generate3dFaceClosure(F.closures, 4);
break;
case MSH_TET_52 :
F.numFaces = 4;
F.monomials = generatePascalSerendipityTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, true);
generate3dFaceClosure(F.faceClosure, 5);
generate3dFaceClosure(F.closures, 5);
break;
case MSH_TET_56 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, false);
generate3dFaceClosure(F.faceClosure, 5);
generate3dFaceClosure(F.closures, 5);
break;
case MSH_QUA_4 :
F.numFaces = 4;
F.monomials = generatePascalQuad(1);
F.points = gmshGeneratePointsQuad(1,false);
generate2dEdgeClosure(F.edgeClosure, 1, 4);
generate2dEdgeClosure(F.closures, 1, 4);
break;
case MSH_QUA_9 :
F.numFaces = 4;
F.monomials = generatePascalQuad(2);
F.points = gmshGeneratePointsQuad(2,false);
generate2dEdgeClosure(F.edgeClosure, 2, 4);
generate2dEdgeClosure(F.closures, 2, 4);
break;
case MSH_QUA_16 :
F.numFaces = 4;
F.monomials = generatePascalQuad(3);
F.points = gmshGeneratePointsQuad(3,false);
generate2dEdgeClosure(F.edgeClosure, 3, 4);
generate2dEdgeClosure(F.closures, 3, 4);
break;
case MSH_QUA_25 :
F.numFaces = 4;
F.monomials = generatePascalQuad(4);
F.points = gmshGeneratePointsQuad(4,false);
generate2dEdgeClosure(F.edgeClosure, 4, 4);
generate2dEdgeClosure(F.closures, 4, 4);
break;
case MSH_QUA_36 :
F.numFaces = 4;
F.monomials = generatePascalQuad(5);
F.points = gmshGeneratePointsQuad(5,false);
generate2dEdgeClosure(F.edgeClosure, 5, 4);
generate2dEdgeClosure(F.closures, 5, 4);
break;
case MSH_QUA_8 :
F.numFaces = 4;
F.monomials = generatePascalQuadSerendip(2);
F.points = gmshGeneratePointsQuad(2,true);
generate2dEdgeClosure(F.edgeClosure, 2, 4);
generate2dEdgeClosure(F.closures, 2, 4);
break;
case MSH_QUA_12 :
F.numFaces = 4;
F.monomials = generatePascalQuadSerendip(3);
F.points = gmshGeneratePointsQuad(3,true);
generate2dEdgeClosure(F.edgeClosure, 3, 4);
generate2dEdgeClosure(F.closures, 3, 4);
break;
case MSH_QUA_16I :
F.numFaces = 4;
F.monomials = generatePascalQuadSerendip(4);
F.points = gmshGeneratePointsQuad(4,true);
generate2dEdgeClosure(F.edgeClosure, 4, 4);
generate2dEdgeClosure(F.closures, 4, 4);
break;
case MSH_QUA_20 :
F.numFaces = 4;
F.monomials = generatePascalQuadSerendip(5);
F.points = gmshGeneratePointsQuad(5,true);
generate2dEdgeClosure(F.edgeClosure, 5, 4);
generate2dEdgeClosure(F.closures, 5, 4);
break;
case MSH_PRI_6 : // first order
F.numFaces = 5;
F.monomials = generatePascalPrism(1);
F.points = gmshGeneratePointsPrism(1, false);
generate3dFaceClosurePrism(F.faceClosure, 1);
generate3dFaceClosurePrism(F.closures, 1);
break;
case MSH_PRI_18 : // second order
F.numFaces = 5;
F.monomials = generatePascalPrism(2);
F.points = gmshGeneratePointsPrism(2, false);
generate3dFaceClosurePrism(F.faceClosure, 2);
generate3dFaceClosurePrism(F.closures, 2);
break;
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(1);
F.points = gmshGeneratePointsTetrahedron(1, false);
generate3dFaceClosure(F.faceClosure, 1);
generate3dFaceClosure(F.closures, 1);
break;
}
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
......
Numeric/polynomialBasis.h
+ 6
− 21
View file @ d4797a29
......@@ -17,34 +17,19 @@ class polynomialBasis
{
public:
typedef std::vector<std::vector<int> > clCont;
clCont faceClosure;
clCont edgeClosure;
clCont vertexClosure;
clCont closures;
fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
int numFaces;
// for a given face/edge, with both a sign and a rotation,
// give an ordered list of nodes on this face/edge
inline int getFaceClosureId(int iFace, int iSign, int iRot) const {
// return iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot; // only for tetrahedra
// return iFace + 5*(iSign == 1 ? 0 : 1) + 10*iRot; // only for prisms
return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot; // both tetrahedra and prisms
}
inline const std::vector<int> &getFaceClosure(int id) const
inline const std::vector<int> &getClosure (int id) const // return the closure of dimension dim
{
return faceClosure[id];
}
inline int getEdgeClosureId(int iEdge, int iSign) const {
return iSign == 1 ? iEdge : edgeClosure.size()/2 + iEdge;
return closures[id];
}
inline const std::vector<int> &getEdgeClosure(int id) const
{
return edgeClosure[id];
}
inline const std::vector<int> &getVertexClosure(int iVertex) const
{
return vertexClosure[iVertex];
inline int getClosureId(int iEl, int iSign=1, int iRot=0) const {
return iEl + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
}
inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
......
Solver/dgGroupOfElements.cpp
+ 219
− 558
View file @ d4797a29
This diff is collapsed.
Solver/dgGroupOfElements.h
+ 11
− 20
View file @ d4797a29
......@@ -25,6 +25,8 @@ class dgConservationLaw;
class dgDofContainer;
class dgMiniInterface;
class dgGroupCollection;
class dgElement {
MElement *_element;
// solution at points
......@@ -129,6 +131,11 @@ public:
class dgGroupOfFaces;
class dgGroupOfConnections {
// there is a finite number of combinations of orientations, senses
// and rotations of the faces (typically 24 for tets). Any pair
// is characterized by a single integer which is the combination
// this closure is for the interpolation that MAY BE DIFFERENT THAN THE
// GEOMETRICAL CLOSURE !!!
std::vector<std::vector<int> > _closures;
std::vector<int> _closuresId;
// face integration point in the coordinate of the left and right element (one fullMatrix per closure)
......@@ -168,11 +175,6 @@ class dgGroupOfFaces {
std::vector<MElement *>_faces;
// Ni integration points, matrix of size Ni x 3 (u,v,weight)
fullMatrix<double> *_integration;
// there is a finite number of combinations of orientations, senses
// and rotations of the faces (typically 24 for tets). Any pair
// is characterized by a single integer which is the combination
// this closure is for the interpolation that MAY BE DIFFERENT THAN THE
// GEOMETRICAL CLOSURE !!!
// detJac at integration points (N*Ni) x 1
fullMatrix<double> *_detJac;
// collocation matrices \psi_i (GP_j)
......@@ -182,16 +184,9 @@ class dgGroupOfFaces {
fullMatrix<double> *_redistribution;
// surface/length/1 of the interface element (sum_qp w_i detJ_i)
fullMatrix<double> *_interfaceSurface;
//common part of the 3 constructors
void init(int pOrder);
public:
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder, int numVertices = -1);
dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder, int numVertices = -1);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);
dgGroupOfFaces (dgGroupCollection &groups, std::vector<dgMiniInterface> &interfaces, int pOrder);
virtual ~dgGroupOfFaces ();
inline bool isBoundary() const {return !_boundaryTag.empty();}
inline const std::string getBoundaryTag() const {return _boundaryTag;}
//this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();}
......@@ -204,12 +199,9 @@ public:
inline double getInterfaceSurface (int iFace)const {return (*_interfaceSurface)(iFace,0);}
const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
private:
void addFace(const MFace &topoFace, const std::vector<int> &iEls);
void addEdge(const MEdge &topoEdge, const std::vector<int> &iEls);
void addVertex(MVertex *topoVertex, const std::vector<int> &iEls);
public:
// duplicate
bool isBoundary() {return _connections.size()==1;}
inline fullMatrix<double> &getNormals () const {return _connections[0]->getNormals();}
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);
......@@ -237,8 +229,6 @@ class dgGroupCollection {
std::vector<dgGroupOfFaces*> _faceGroups; //interface
std::vector<dgGroupOfFaces*> _boundaryGroups; //boundary
std::vector<dgGroupOfElements*> _ghostGroups; //ghost volume
// container of different face types (identified by the number of vertices)
std::set<int> _interfaceTypes;
//{group,id} of the elements to send to each partition for a scatter operation
std::vector< std::vector<std::pair<int,int> > >_elementsToSend;
......@@ -248,6 +238,7 @@ class dgGroupCollection {
// Multirate stuff
int _dtMaxExponent;
void buildParallelStructure();
public:
inline int getDtMaxExponent() const {return _dtMaxExponent;}
inline GModel* getModel() {return _model;}
......@@ -255,7 +246,7 @@ class dgGroupCollection {
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 dgGroupOfElements *getElementGroup(int i) const {return i<getNbElementGroups()?_elementGroups[i]:_ghostGroups[i-getNbElementGroups()];}
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];}
......
Solver/dgResidual.cpp
+ 0
− 11
View file @ d4797a29
......@@ -227,17 +227,6 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
/* Msg::Info("Left and Right for %p : gL %p %s %s %d, gR %p %s %s %d",
&faces,
&faces.getGroupLeft(),
faces.getGroupLeft().getIsInnerMultirateBuffer()?"Inner":"",
faces.getGroupLeft().getIsOuterMultirateBuffer()?"Outer":"",
faces.getGroupLeft().getMultirateExponent(),
&faces.getGroupRight(),
faces.getGroupRight().getIsInnerMultirateBuffer()?"Inner":"",
faces.getGroupRight().getIsOuterMultirateBuffer()?"Outer":"",
faces.getGroupRight().getMultirateExponent()
);*/
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
faces.mapToInterface(_nbFields, solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupRight()), solInterface);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment