diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 6a82c3809b1086887eb129742460b5d302665174..46f7c9e38e88ea4883b0a4ab7948ea879abaae58 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -749,19 +749,17 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order)
{
closure.clear();
- // first give edge nodes of the three edges in direct order
- int index = 3;
- std::vector<int> c0,c1,c2;
- for (int i = 0; i < order - 1; i++, index++) c0.push_back(index);
- for (int i = 0; i < order - 1; i++, index++) c1.push_back(index);
- for (int i = 0; i < order - 1; i++, index++) c2.push_back(index);
- closure.push_back(c0);closure.push_back(c1);closure.push_back(c2);
- // then give edge nodes in reverse order
- std::vector<int> c3,c4,c5;
- for (int i = c0.size() - 1; i >= 0; i--) c3.push_back(c0[i]);
- for (int i = c1.size() - 1; i >= 0; i--) c4.push_back(c1[i]);
- for (int i = c2.size() - 1; i >= 0; i--) c5.push_back(c2[i]);
- closure.push_back(c3); closure.push_back(c4); closure.push_back(c5);
+ closure.resize(6);
+ for (int j = 0; j < 3 ; j++){
+ closure[j].push_back(j);
+ closure[j].push_back((j+1)%3);
+ closure[3+j].push_back((j+1)%3);
+ closure[3+j].push_back(j);
+ for (int i=0; i < order-1; i++){
+ closure[j].push_back( 3 + (order-1)*j + i );
+ closure[3+j].push_back(3 + (order-1)*(j+1) -i -1);
+ }
+ }
}
std::map<int, polynomialBasis> polynomialBases::fs;
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 8025824daed907c7fb0ebbe45d2daa5fd3ac1854..6e119af54709816b96470916f81004cd8ba8241a 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -1,7 +1,11 @@
+#include <set>
#include "dgAlgorithm.h"
#include "dgGroupOfElements.h"
#include "dgConservationLaw.h"
-
+#include "GEntity.h"
+#include "MElement.h"
+#include "GModel.h"
+#include "MEdge.h"
/*
compute
\int \vec{f} \cdot \grad \phi dv
@@ -136,3 +140,124 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
solEl.gemm(group.getInverseMassMatrix(i),residuEl);
}
}
+
+void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &residual // residual !! at faces nodes
+ )
+{
+ int nbFields = claw.nbFields();
+ const dgBoundaryCondition *boundaryCondition = claw.boundaryCondition(group.getBoundaryTag());
+ // ----- 1 ---- get the solution at quadrature points
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
+ group.getCollocationMatrix().mult(solution, solutionQP);
+ // ----- 2 ---- compute normal fluxes at integration points
+ fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
+ fullMatrix<double> solutionQPRight(group.getNbIntegrationPoints(), nbFields);
+ for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
+ // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
+ fullMatrix<double> solutionQPLeft (solutionQP, iFace*nbFields, nbFields);
+ fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields, nbFields);
+ dgFace DGF( group.getFace(iFace), group.getElementLeft(iFace), NULL,
+ solutionQPLeft, solutionQPRight, group.getIntegrationPointsMatrix(),group.getNormals(iFace));
+ // ----- 2.3.2 --- compute fluxes
+ switch (boundaryCondition->type()){
+ case dgBoundaryCondition::EXTERNAL_VALUES : {
+ (*boundaryCondition)(DGF, &solutionQPRight);
+ (*claw.riemannSolver())(DGF, &normalFluxQP);
+ break;
+ }
+ case dgBoundaryCondition::FLUX :{
+ (*boundaryCondition)(DGF, &normalFluxQP);
+ }
+ default :
+ throw;
+ }
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ for (int k=0;k<nbFields;k++)
+ normalFluxQP(iPt,k) *= detJ;
+ }
+ }
+ // ----- 3 ---- do the redistribution at face nodes using BLAS3
+ residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+}
+
+void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups,
+ std::vector<dgGroupOfFaces*> &bGroups)
+{
+ 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::vector<MElement *> allElements;
+ 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);
+ if(dim==2)
+ boundaryEdges[physicalName].insert( MEdge(element->getVertex(0), element->getVertex(1)) );
+ else
+ boundaryFaces[physicalName].insert( MFace(element->getVertex(0), element->getVertex(1),element->getVertex(2)) );
+ }
+ }
+ }else if(entity->dim() == dim){
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++)
+ allElements.push_back(entity->getMeshElement(iel));
+ }
+ }
+ eGroups.push_back(new dgGroupOfElements(allElements,order));
+ fGroups.push_back(new dgGroupOfFaces(*eGroups[0],order));
+ if(dim==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[0],mapIt->first,order,mapIt->second));
+ }
+ }else if(dim=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[0],mapIt->first,order,mapIt->second));
+ }
+ }else throw;
+}
+
+// works only if there is only 1 group of element
+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
+ const fullMatrix<double> &solution, // solution
+ fullMatrix<double> &residu) // residual
+{
+ dgAlgorithm algo;
+ residu.scale(0);
+ //volume term
+ for(size_t i=0;i<eGroups.size() ; i++) {
+ algo.residualVolume(claw,*eGroups[i],solution,residu);
+ }
+ //interface term
+ for(size_t i=0;i<fGroups.size() ; i++) {
+ dgGroupOfFaces &faces = *fGroups[i];
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2);
+ faces.mapToInterface(1, solution, solution, solInterface);
+ algo.residualInterface(claw,faces,solInterface,residuInterface);
+ faces.mapFromInterface(1, residuInterface, residu, residu);
+ }
+ //boundaries
+ for(size_t i=0;i<bGroups.size() ; i++) {
+ dgGroupOfFaces &faces = *bGroups[i];
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements());
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements());
+ faces.mapToInterface(1, solution, solution, solInterface);
+ algo.residualBoundary(claw,faces,solInterface,residuInterface);
+ faces.mapFromInterface(1, residuInterface, residu, residu);
+ }
+}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 2af608c6df2d32098b24556d7dde4194e05a508c..9a0c3c96b123ec03ae4d904e0e07ac45809fda52 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -1,8 +1,9 @@
#ifndef _DG_ALGORITHM_H_
#define _DG_ALGORITHM_H_
-#include "dofManager.h"
-
+#include "fullMatrix.h"
+#include <vector>
+class GModel;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
@@ -20,13 +21,28 @@ class dgAlgorithm {
const dgGroupOfFaces & group,
const fullMatrix<double> &solution, // ! at face nodes
fullMatrix<double> &residual); // ! at face nodes
- void residualBoundaryConditions ( /*dofManager &dof,*/
- const dgConservationLaw &law,
- const dgGroupOfElements & group);
+ void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &residual // residual !! at faces nodes
+ );
+ // works only if there is only 1 group of element
+ 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
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &residual // residual !! at faces nodes
+ );
void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
fullMatrix<double> &residu,
fullMatrix<double> &sol);
+ void buildGroups(GModel *model, int dim, int order,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups,
+ std::vector<dgGroupOfFaces*> &bGroups);
};
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 05eb28348a0b8fe5ca81e921b3e58551dd08240c..fd1441f14e9b1dab4a6f5a3a893803d88e91b382 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -22,11 +22,20 @@ class dgFaceTerm{
virtual void operator () (const dgFace &, fullMatrix<double> fcx[]) const = 0;
};
+class dgBoundaryCondition {
+public:
+ virtual ~dgBoundaryCondition () {}
+ typedef enum {FLUX=0,EXTERNAL_VALUES}boundaryType;
+ virtual void operator () (const dgFace&, fullMatrix<double> fcx[]) const = 0;
+ virtual boundaryType type() const =0;
+};
+
class dgConservationLaw {
protected :
int _nbf;
dgTerm *_diffusive, *_convective, *_source, *_maxConvectiveSpeed;
dgFaceTerm *_riemannSolver;
+ std::map<const std::string,dgBoundaryCondition*> _boundaryConditions;
public:
dgConservationLaw () : _diffusive(0), _convective(0), _source (0),
_riemannSolver(0),_maxConvectiveSpeed (0) {}
@@ -39,10 +48,15 @@ public:
inline const dgTerm * sourceTerm () const {return _source;}
inline const dgFaceTerm * riemannSolver () const {return _riemannSolver;}
inline const dgTerm * maxConvectiveSpeed () const {return _maxConvectiveSpeed;}
+ inline const dgBoundaryCondition *boundaryCondition(const std::string tag) const {
+ std::map<const std::string,dgBoundaryCondition*>::const_iterator it = _boundaryConditions.find(tag);
+ if(it==_boundaryConditions.end())
+ throw;
+ return it->second;
+ }
};
-
dgConservationLaw *dgNewConservationLawAdvection();
#endif
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index 269d43538a685a3740f2dc32c5b62909546815ab..8f4677c9bc248eba70297de1fae432bcf19b7531 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -4,10 +4,13 @@
#include "SPoint3.h"
#include "MElement.h"
void getV(SPoint3 &p, double (&v)[3]){
- double r=hypot(p.x(),p.y());
+ /*double r=hypot(p.x(),p.y());
double alpha=atan2(p.y(),p.x());
v[0]=r*sin(alpha);
v[1]=-r*cos(alpha);
+ v[2]=0;*/
+ v[0]=0.15;
+ v[1]=0;
v[2]=0;
}
@@ -44,14 +47,27 @@ class dgAdvectionRiemanTerm : public dgFaceTerm{
double un=v[0]*normals(0,i)+v[1]*normals(1,i)+v[2]*normals(2,i);
if(un>0){
fcx[0](i,0) = -solLeft(i,0)*un;
- fcx[0](i,1) = solLeft(i,0)*un;
+ if(fcx[0].size2()==2)
+ fcx[0](i,1) = solLeft(i,0)*un;
}else{
fcx[0](i,0) = -solRight(i,0)*un;
- fcx[0](i,1) = solRight(i,0)*un;
+ if(fcx[0].size2()==2)
+ fcx[0](i,1) = solRight(i,0)*un;
}
}
}
};
+class dgAdvectionBoundary0External : public dgBoundaryCondition {
+ boundaryType type() const {
+ return EXTERNAL_VALUES;
+ }
+ void operator () (const dgFace &face, fullMatrix<double> fcx[]) const{
+ const fullMatrix<double> &solLeft = face.solutionLeft();
+ for(int i=0; i< solLeft.size1(); i++) {
+ fcx[0](i,0) = 0;
+ }
+ }
+};
class dgConservationLawAdvection : public dgConservationLaw {
public:
@@ -60,11 +76,19 @@ class dgConservationLawAdvection : public dgConservationLaw {
_nbf = 1;
_convective = new dgAdvectionFluxTerm;
_riemannSolver = new dgAdvectionRiemanTerm;
+ _boundaryConditions["Left"] = new dgAdvectionBoundary0External();
+ _boundaryConditions["Right"] = new dgAdvectionBoundary0External();
+ _boundaryConditions["Top"] = new dgAdvectionBoundary0External();
+ _boundaryConditions["Bottom"] = new dgAdvectionBoundary0External();
}
~dgConservationLawAdvection()
{
delete _convective;
delete _riemannSolver;
+ delete _boundaryConditions["Left"];
+ delete _boundaryConditions["Right"];
+ delete _boundaryConditions["Top"];
+ delete _boundaryConditions["Bottom"];
}
};
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 6d10a34e4ebbc1a8a716ca4af83cb84355afaea4..abbb35e3263bfcbad7c79518962598707147abf7 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -4,6 +4,7 @@
#include "Numeric.h"
#include "MTriangle.h"
#include "MLine.h"
+#include "GModel.h"
static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
int npts;
@@ -156,22 +157,24 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
// compute all closures
// 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)));
+ const std::vector<int> & geomClosure = elLeft.getFunctionSpace()->getFaceClosure(ithFace, sign, rot);
+ if(iElRight>=0){
+ MElement &elRight = *_groupRight.getElement(iElRight);
+ _right.push_back(iElRight);
+ 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]) );
+ vertices.push_back( elLeft.getVertex(geomClosure[j]) );
// triangular face
if (topoFace.getNumVertices() == 3){
switch(vertices.size()){
@@ -186,31 +189,20 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
// quad face 2 do
else throw;
}
-static std::vector<int> *fakeClosure2d(const polynomialBasis *fs, int ithEdge, int sign){
- std::vector<int> closure;
- if(sign==1){
- closure.push_back(ithEdge);
- closure.push_back((ithEdge+1)%3);
- }else{
- closure.push_back((ithEdge+1)%3);
- closure.push_back(ithEdge);
- }
- std::vector<int> closureHO = fs->getEdgeClosure(ithEdge, sign);
- closure.insert(closure.end(),closureHO.begin(),closureHO.end());
- return new std::vector<int>(closure);
-}
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(fakeClosure2d(_fsLeft, ithEdge, sign));
- elRight.getEdgeInfo (topoEdge, ithEdge, sign);
- _closuresRight.push_back(fakeClosure2d(_fsRight, ithEdge, sign));
- const std::vector<int> &geomClosure = elRight.getFunctionSpace()->getEdgeClosure(ithEdge, sign);
+ _closuresLeft.push_back(&_fsLeft->getEdgeClosure(ithEdge, sign));
+ const std::vector<int> &geomClosure = elLeft.getFunctionSpace()->getEdgeClosure(ithEdge, sign);
+ if(iElRight>=0){
+ _right.push_back(iElRight);
+ MElement &elRight = *_groupRight.getElement(iElRight);
+ elRight.getEdgeInfo (topoEdge, ithEdge, sign);
+ _closuresRight.push_back(&_fsRight->getEdgeClosure(ithEdge, sign));
+ }
std::vector<MVertex*> vertices;
if(sign==1)
for(int j=0;j<topoEdge.getNumVertices();j++)
@@ -219,7 +211,7 @@ void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight){
for(int j=topoEdge.getNumVertices()-1;j>=0;j--)
vertices.push_back(topoEdge.getVertex(j));
for (int j=0; j<geomClosure.size() ; j++)
- vertices.push_back( elRight.getVertex(geomClosure[j]) );
+ vertices.push_back( elLeft.getVertex(geomClosure[j]) );
switch(vertices.size()){
case 2 : _faces.push_back(new MLine (vertices) ); break;
case 3 : _faces.push_back(new MLine3 (vertices) ); break;
@@ -258,6 +250,42 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _collocation;
delete _detJac;
}
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges):
+ _groupLeft(elGroup),_groupRight(elGroup)
+{
+ _boundaryTag=boundaryTag;
+ if(boundaryTag=="")
+ throw;
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=NULL;
+ 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(boundaryEdges.find(edge)!=boundaryEdges.end())
+ addEdge(edge,i,-1);
+ }
+ }
+ init(pOrder);
+}
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces):
+ _groupLeft(elGroup),_groupRight(elGroup)
+{
+ _boundaryTag=boundaryTag;
+ if(boundaryTag=="")
+ throw;
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=NULL;
+ 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(boundaryFaces.find(face)!=boundaryFaces.end())
+ addFace(face,i,-1);
+ }
+ }
+ init(pOrder);
+}
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
_groupLeft(elGroup),_groupRight(elGroup)
@@ -305,15 +333,27 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
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);
+ if(isBoundary()){
+ for(int i=0; i<getNbElements(); 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*nFields + iField) = vLeft(closureLeft[j], _left[i]*nFields + iField);
+ }
+ }
+ }
+
+ }else{
+ 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);
}
- for(size_t j =0; j < closureRight.size(); j++)
- v(j, (i*2+1)*nFields + iField) = vRight(closureRight[j], _right[i]*nFields + iField);
}
}
}
@@ -324,14 +364,25 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
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);
+ if(isBoundary()){
+ for(int i=0; i<getNbElements(); 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*nFields + iField);
+ }
+ }
+ }else{
+ 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 a2f5035d5a72d4eb7074ad7fac40f8d547a9e708..da35c02100acbff7c886136eeb5ada8f5d07b68e 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -13,10 +13,11 @@
we DO NOT store N, Ni or Nv (matrices and vectors
have the right sizes)
*/
+#include "MFace.h"
+#include "MEdge.h"
class MElement;
-class MFace;
-class MEdge;
class polynomialBasis;
+class GEntity;
class dgElement {
MElement *_element;
@@ -105,6 +106,8 @@ public:
};
class dgGroupOfFaces {
+ // only used if this is a group of boundary faces
+ std::string _boundaryTag;
const dgGroupOfElements &_groupLeft,&_groupRight;
void addFace(const MFace &topoFace, int iElLeft, int iElRight);
void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
@@ -143,7 +146,11 @@ public:
const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
inline fullMatrix<double> getNormals (int iFace) const {return fullMatrix<double>(*_normals,iFace*getNbIntegrationPoints(),getNbIntegrationPoints());}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
+ 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);
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();}
inline int getNbNodes() const {return _collocation->size1();}
@@ -158,5 +165,4 @@ public:
};
-
#endif
diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp
index d162aa98ff196865da2d6f5c484698c54a77d47b..cd53ba9c7ac9c1caf47a3d25a6102de78d95c421 100644
--- a/Solver/dgMainTest.cpp
+++ b/Solver/dgMainTest.cpp
@@ -9,7 +9,6 @@
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v);
-std::vector<MElement *> getAllTri(GModel *model);
class testSourceTerm : public dgTerm {
void operator () (const dgElement &el, fullMatrix<double> fcx[]) const{
@@ -18,7 +17,7 @@ class testSourceTerm : public dgTerm {
SPoint3 p;
for(int i=0; i< sol.size1(); i++) {
el.element()->pnt(qp(i,0),qp(i,1),qp(i,2),p);
- fcx[0](i,0)=exp(-(pow(p.x()-0.2,2)+pow(p.y(),2))*100);
+ fcx[0](i,0)=exp(-(pow(p.x()-0.2,2)+pow(p.y()-0.3,2))*100);
}
}
};
@@ -35,60 +34,43 @@ class dgConservationLawInitialCondition : public dgConservationLaw {
};
int main(int argc, char **argv){
- GmshMergeFile("mesh1.msh");
- //GmshMergeFile("square2.msh");
- std::vector<MElement *> allTri=getAllTri(GModel::current());
+ GmshMergeFile("square1.msh");
+ //we probably need a class to group those three ones
+ std::vector<dgGroupOfElements*> elementGroups;
+ std::vector<dgGroupOfFaces*> faceGroups;
+ std::vector<dgGroupOfFaces*> boundaryGroups;
int order=1;
- dgGroupOfElements elements(allTri,order);
- dgGroupOfFaces faces(elements,order);
- //dgGroupOfFaces faces(elements0,elements1,order);
- int nbNodes = elements.getNbNodes();
- fullMatrix<double> sol(nbNodes,elements.getNbElements());
- fullMatrix<double> residu(nbNodes,elements.getNbElements());
+ int dimension=2;
dgAlgorithm algo;
+ algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
+
+ //for now, we suppose there is only one group of elements
+ int nbNodes = elementGroups[0]->getNbNodes();
+ fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements());
+ fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements());
// initial condition
{
dgConservationLawInitialCondition initLaw;
- algo.residualVolume(initLaw,elements,sol,residu);
- algo.multAddInverseMassMatrix(elements,residu,sol);
+ algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
+ algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
}
- print("init.pos",elements,&sol(0,0));
+ print("init.pos",*elementGroups[0],&sol(0,0));
//advection
dgConservationLaw *law = dgNewConservationLawAdvection();
for(int iT=0; iT<1000; iT++) {
- residu.scale(0);
- algo.residualVolume(*law,elements,sol,residu);
- { //interface term
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2);
- faces.mapToInterface(1, sol, sol, solInterface);
- algo.residualInterface(*law,faces,solInterface,residuInterface);
- faces.mapFromInterface(1, residuInterface, residu, residu);
- }
+ algo.residual(*law,elementGroups,faceGroups,boundaryGroups,sol,residu);
residu.scale(0.01);
- algo.multAddInverseMassMatrix(elements,residu,sol);
+ algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
if(iT%10==0){
char name[100];
sprintf(name,"test_%05i.pos",iT/10);
printf("%i\n",iT);
- print(name,elements,&sol(0,0));
+ print(name,*elementGroups[0],&sol(0,0));
}
}
delete law;
}
-std::vector<MElement *> getAllTri(GModel *model){
- std::vector<GEntity*> entities;
- model->getEntities(entities);
- std::vector<MElement *> all_tri;
- for(std::vector<GEntity *>::iterator itent=entities.begin(); itent!=entities.end(); itent++){
- if ((*itent)->dim()!=2) continue;
- for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
- all_tri.push_back((*itent)->getMeshElement(iel));
- }
- return all_tri;
-}
-
void print (const char *filename,const dgGroupOfElements &els, double *v) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);