diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index dd6018c3209e9275bef5bd8045ac2e770c9a355e..7d3d7d1725d0bead863bbe28f972f9161ce1c5b0 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -155,7 +155,7 @@ class fullMatrix
}
return false; // no reallocation
}
- void setAsProxy(fullMatrix<scalar> &original, int c_start, int c)
+ void setAsProxy(const fullMatrix<scalar> &original, int c_start, int c)
{
if(_data && _own_data)
delete [] _data;
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 8058c567ac9dd4bc6768f19be814a2cda06ede65..37c46b7d931ffa182a8597fac4407e6399bf4469 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -48,8 +48,10 @@ class dgConservationLaw {
inline const dgBoundaryCondition *getBoundaryCondition(const std::string tag) const {
std::map<const std::string,dgBoundaryCondition*>::const_iterator it = _boundaryConditions.find(tag);
- if(it==_boundaryConditions.end())
+ if(it==_boundaryConditions.end()) {
+ Msg::Error("no boundary condition defined with tag '%s'", tag.c_str());
throw;
+ }
return it->second;
}
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index aae4f59b214f11fd3a828885cd546956d9372ef8..1836721e214d16e51de15e55b158589ca65bf800 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -53,7 +53,6 @@ static fullMatrix<double> dgGetFaceIntegrationRuleOnElement (
return intgElement;
}
-
dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
int polyOrder,
int ghostPartition)
@@ -175,9 +174,9 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
const fullMatrix<double> &vRight,
fullMatrix<double> &v)
{
- if(isBoundary()){
+ if(_connections.size()==1){
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureLeft = _connections[0]->getClosure(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], _connections[0]->getElementId(i)*nFields + iField);
@@ -187,8 +186,8 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
}else{
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureRight = getClosureRight(i);
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
+ const std::vector<int> &closureRight = _connections[1]->getClosure(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],_connections[0]->getElementId(i)*nFields + iField);
@@ -206,9 +205,9 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
fullMatrix<double> &vRight
)
{
- if(isBoundary()){
+ if(_connections.size()==1){
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*nFields + iField);
@@ -216,8 +215,8 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
}
}else{
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureRight = getClosureRight(i);
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
+ const std::vector<int> &closureRight = _connections[1]->getClosure(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*2*nFields + iField);
@@ -227,14 +226,14 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
}
}
}
+
void dgGroupOfFaces::mapLeftFromInterface ( int nFields,
const fullMatrix<double> &v,
fullMatrix<double> &vLeft
)
{
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureRight = getClosureRight(i);
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*2*nFields + iField);
@@ -246,11 +245,10 @@ void dgGroupOfFaces::mapRightFromInterface ( int nFields,
fullMatrix<double> &vRight
)
{
- if(isBoundary())
+ if(_connections.size()==1)
return;
for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureRight = getClosureRight(i);
- const std::vector<int> &closureLeft = getClosureLeft(i);
+ const std::vector<int> &closureRight = _connections[1]->getClosure(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureRight.size(); j++)
vRight(closureRight[j], _connections[1]->getElementId(i)*nFields + iField) += v(j, (i*2+1)*nFields + iField);
@@ -271,7 +269,7 @@ void dgGroupOfConnections::init() {
// compute data on quadrature points : normals and dPsidX
double g[256][3];
_normals = fullMatrix<double> (3, nQP*size);
- _dPsiDxOnQP = fullMatrix<double> ( nQP*3, nPsi*size);
+ _dPsiDx = fullMatrix<double> ( nQP*3, nPsi*size);
int index = 0;
for (size_t i=0; i<size;i++){
const std::vector<int> &closure = getClosure(i);
@@ -284,9 +282,9 @@ void dgGroupOfConnections::init() {
//compute dPsiDxOnQP
// (iQP*3+iXYZ , iFace*NPsi+iPsi)
for (int iPsi=0; iPsi< nPsi; iPsi++) {
- _dPsiDxOnQP(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
- _dPsiDxOnQP(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
- _dPsiDxOnQP(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
+ _dPsiDx(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
+ _dPsiDx(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
+ _dPsiDx(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
}
//compute face normals
double &nx=_normals(0,index);
@@ -311,14 +309,12 @@ void dgGroupOfConnections::init() {
}
}
-
// 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
void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
{
@@ -371,10 +367,10 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
}
}
-class dgMiniConnection {
+class dgMiniConnection
+{
public:
int iGroup, iElement, iClosure;
- int numVertices;
dgMiniConnection (int iGroup_, int iElement_, int iClosure_)
{
iGroup = iGroup_;
@@ -392,7 +388,8 @@ class dgMiniConnection {
}
};
-class dgMiniInterface {
+class dgMiniInterface
+{
public:
int physicalTag;
int numVertices;
@@ -411,7 +408,8 @@ class dgMiniInterface {
}
};
-static std::vector<dgMiniInterface> *_createMiniInterfaces(dgGroupCollection &groups) {
+static std::vector<dgMiniInterface> *_createMiniInterfaces(dgGroupCollection &groups)
+{
std::vector<GEntity*> entities;
groups.getModel()->getEntities(entities);
std::map<MVertex*, dgMiniInterface> vertexInterfaces;
@@ -477,18 +475,24 @@ static std::vector<dgMiniInterface> *_createMiniInterfaces(dgGroupCollection &gr
switch(dim) {
case 1:
interfaces->reserve(vertexInterfaces.size());
- for(std::map<MVertex*, dgMiniInterface>::iterator it = vertexInterfaces.begin(); it != vertexInterfaces.end(); it++)
+ for(std::map<MVertex*, dgMiniInterface>::iterator it = vertexInterfaces.begin(); it != vertexInterfaces.end(); it++) {
+ it->second.numVertices = 1;
interfaces->push_back(it->second);
+ }
break;
case 2:
interfaces->reserve(edgeInterfaces.size());
- for(std::map<MEdge, dgMiniInterface, Less_Edge>::iterator it = edgeInterfaces.begin(); it != edgeInterfaces.end(); it++)
+ for(std::map<MEdge, dgMiniInterface, Less_Edge>::iterator it = edgeInterfaces.begin(); it != edgeInterfaces.end(); it++) {
+ it->second.numVertices = it->first.getNumVertices();
interfaces->push_back(it->second);
+ }
break;
case 3:
interfaces->reserve(faceInterfaces.size());
- for(std::map<MFace, dgMiniInterface, Less_Face>::iterator it = faceInterfaces.begin(); it != faceInterfaces.end(); it++)
+ for(std::map<MFace, dgMiniInterface,Less_Face>::iterator it = faceInterfaces.begin(); it != faceInterfaces.end(); it++) {
+ it->second.numVertices = it->first.getNumVertices();
interfaces->push_back(it->second);
+ }
break;
}
for (size_t i = 0; i < interfaces->size(); i++) {
@@ -565,8 +569,13 @@ dgGroupOfFaces::dgGroupOfFaces (dgGroupCollection &groups, std::vector<dgMiniInt
dgMiniInterface &first = interfaces.front();
size_t nconnections = first.connections.size();
int dim = groups.getElementGroup(first.connections[0].iGroup)->getElement(first.connections[0].iElement)->getDim()-1;
+ _physicalTag.clear();
if(first.physicalTag>=0)
- _boundaryTag = groups.getModel()->getPhysicalName(dim, first.physicalTag);
+ _physicalTag = groups.getModel()->getPhysicalName(dim, first.physicalTag);
+ if(nconnections==1 && (_physicalTag.empty() || _physicalTag=="")) {
+ Msg::Error("boundary face without tag in the mesh");
+ throw;
+ }
for (size_t i=0; i<nconnections; i++) {
_connections.push_back (new dgGroupOfConnections(*groups.getElementGroup(first.connections[i].iGroup), *this, pOrder));
}
@@ -592,7 +601,7 @@ dgGroupOfFaces::dgGroupOfFaces (dgGroupCollection &groups, std::vector<dgMiniInt
}
break;
case 3 :
- if (connection.numVertices == 3) {
+ if (interface.numVertices == 3) {
switch(vertices.size()){
case 3 : _faces.push_back (new MTriangle (vertices)); break;
case 6 : _faces.push_back (new MTriangle6 (vertices)); break;
@@ -601,7 +610,7 @@ dgGroupOfFaces::dgGroupOfFaces (dgGroupCollection &groups, std::vector<dgMiniInt
case 21 : _faces.push_back (new MTriangleN (vertices, 5)); break;
default : throw;
}
- } else if (connection.numVertices == 4) {
+ } else if (interface.numVertices == 4) {
switch(vertices.size()){
case 4 : _faces.push_back (new MQuadrangle (vertices)); break;
case 8 : _faces.push_back (new MQuadrangle8 (vertices)); break;
@@ -712,7 +721,6 @@ void dgGroupCollection::buildParallelStructure()
delete []shiftRecv;
}
-
// Split the groups of elements depending on their local time step
double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution) {
// What are the levels/layers:
@@ -959,8 +967,6 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
lowerLevelGroupIdEnd=currentNewGroupId;
isOuterBufferLayer=!isOuterBufferLayer;
}
-
-
// Some stats
int count=0;
for(int i=0;i<newGroups.size();i++){
@@ -980,11 +986,11 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
return dtRef;
}
-
dgGroupCollection::dgGroupCollection()
{
_groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
}
+
dgGroupCollection::dgGroupCollection(GModel *model, int dimension, int order)
{
_groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
@@ -1010,6 +1016,7 @@ void dgGroupCollection::find (MElement*e, int &ig, int &index)
if (index != -1)return;
}
}
+
#include "LuaBindings.h"
void dgGroupCollection::registerBindings(binding *b)
{
@@ -1019,7 +1026,6 @@ void dgGroupCollection::registerBindings(binding *b)
cb->setDescription("a group of element sharing the same properties");
cb = b->addClass<dgGroupOfFaces>("dgGroupOfFaces");
cb->setDescription("a group of faces. All faces of this group are either interior of the same group of elements or at the boundary between two group of elements");
-
cb = b->addClass<dgGroupCollection>("dgGroupCollection");
cb->setDescription("The GroupCollection class let you access to group partitioning functions");
cm = cb->setConstructor<dgGroupCollection,GModel*,int,int>();
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 3656a7055cd0b84599083b6d49a3129fe2ed1269..987e86593a50aad9b08d9235ece254d93d44a126 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -45,7 +45,6 @@ public:
class dgGroupOfFaces;
-
// store topological and geometrical data for 1 group for 1 discretisation
class dgGroupOfElements {
int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition
@@ -142,7 +141,7 @@ class dgGroupOfConnections {
std::vector<fullMatrix<double> > _integrationPoints;
// XYZ gradient of the shape functions of both elements on the integrations points of the face
// (iQP*3+iXYZ , iFace*NPsi+iPsi)
- fullMatrix<double> _dPsiDxOnQP;
+ fullMatrix<double> _dPsiDx;
const polynomialBasis *_fs;
const dgGroupOfElements &_elementGroup;
const dgGroupOfFaces &_faceGroup;
@@ -152,14 +151,14 @@ class dgGroupOfConnections {
public:
void addElement(int iElement, int iClosure);
dgGroupOfConnections(const dgGroupOfElements &elementGroup, const dgGroupOfFaces &face, int pOrder);
- inline const polynomialBasis *getFunctionSpace() {return _fs;}
- inline const dgGroupOfElements &getGroupOfElements() {return _elementGroup;}
- inline int getElementId(int i) {return _elementId[i];}
- inline MElement *getElement(int i) {return _elementGroup.getElement(_elementId[i]);}
- inline const std::vector<int>& getClosure(int i) {return _closures[_closuresId[i]];}
- inline fullMatrix<double>& getIntegrationPointsOnElement(int i) {return _integrationPoints[_closuresId[i]];}
- inline fullMatrix<double>& getDPsiDxOnQP() {return _dPsiDxOnQP;}
- inline fullMatrix<double>& getNormals() {return _normals;}
+ inline const polynomialBasis *getFunctionSpace() const {return _fs;}
+ inline const dgGroupOfElements &getGroupOfElements() const {return _elementGroup;}
+ inline int getElementId(int i) const {return _elementId[i];}
+ inline MElement *getElement(int i) const {return _elementGroup.getElement(_elementId[i]);}
+ inline const std::vector<int>& getClosure(int i) const {return _closures[_closuresId[i]];}
+ inline const fullMatrix<double>& getIntegrationPointsOnElement(int i) const {return _integrationPoints[_closuresId[i]];}
+ inline const fullMatrix<double>& getDPsiDx() const {return _dPsiDx;}
+ inline const fullMatrix<double>& getNormals() const {return _normals;}
void init();
};
@@ -167,7 +166,7 @@ class dgGroupOfFaces {
std::vector<dgGroupOfConnections*> _connections;
// normals always point outside left to right
// only used if this is a group of boundary faces
- std::string _boundaryTag;
+ std::string _physicalTag;
// Two polynomialBases for left and right elements
// the group has always the same types for left and right
const polynomialBasis *_fsFace;
@@ -187,7 +186,9 @@ class dgGroupOfFaces {
public:
dgGroupOfFaces (dgGroupCollection &groups, std::vector<dgMiniInterface> &interfaces, int pOrder);
virtual ~dgGroupOfFaces ();
- inline const std::string getBoundaryTag() const {return _boundaryTag;}
+ inline const std::string getPhysicalTag() const {return _physicalTag;}
+ inline const dgGroupOfConnections &getGroupOfConnections(int i) const {return *_connections[i];}
+ inline int getNbGroupOfConnections() const {return _connections.size();}
//this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();}
inline int getNbNodes() const {return _collocation->size2();}
@@ -199,28 +200,11 @@ 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];}
-public:
- // duplicate
- bool isBoundary() {return _connections.size()==1;}
- inline fullMatrix<double> &getNormals () const {return _connections[0]->getNormals();}
+public: // should be more generic on the number of connections
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);
void mapLeftFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft);
void mapRightFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vRight);
- inline fullMatrix<double> & getDPsiLeftDxMatrix() const { return _connections[0]->getDPsiDxOnQP();}
- inline fullMatrix<double> & getDPsiRightDxMatrix() const { return _connections[1]->getDPsiDxOnQP();}
- inline const dgGroupOfElements &getGroupLeft()const {return _connections[0]->getGroupOfElements(); }
- inline const dgGroupOfElements &getGroupRight()const {return _connections[1]->getGroupOfElements(); }
- inline int getElementLeftId (int i) const {return _connections[0]->getElementId(i);}
- inline int getElementRightId (int i) const {return _connections[1]->getElementId(i);}
- inline MElement* getElementLeft (int i) const {return _connections[0]->getElement(i);}
- inline MElement* getElementRight (int i) const {return _connections[1]->getElement(i);}
- inline double getElementVolumeLeft(int iFace) const {return _connections[0]->getGroupOfElements().getElementVolume(getElementLeftId(iFace));}
- inline double getElementVolumeRight(int iFace) const {return _connections[1]->getGroupOfElements().getElementVolume(getElementRightId(iFace));}
- inline const std::vector<int> &getClosureLeft(int iFace) const{ return _connections[0]->getClosure(iFace);}
- inline const std::vector<int> &getClosureRight(int iFace) const{ return _connections[1]->getClosure(iFace);}
- inline fullMatrix<double> &getIntegrationOnElementLeft(int iFace) { return _connections[0]->getIntegrationPointsOnElement(iFace);}
- inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _connections[1]->getIntegrationPointsOnElement(iFace);}
};
class dgGroupCollection {
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index 6b9c9fb2947ff07bf47eb639ef81026b93bcdd54..2bce19da120043de4576bd1d0bdc6ce05b677ebf 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -22,8 +22,10 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
fullMatrix<double> TempL, TempR;
for( int iGFace=0; iGFace<groups->getNbFaceGroups(); iGFace++) {
dgGroupOfFaces* group = groups->getFaceGroup(iGFace);
- const dgGroupOfElements *groupLeft = &group->getGroupLeft();
- const dgGroupOfElements *groupRight = &group->getGroupRight();
+ const dgGroupOfConnections &left = group->getGroupOfConnections(0);
+ const dgGroupOfConnections &right = group->getGroupOfConnections(1);
+ const dgGroupOfElements *groupLeft = &left.getGroupOfElements();
+ const dgGroupOfElements *groupRight = &right.getGroupOfElements();
fullMatrix<double> &solleft = solution->getGroupProxy(groupLeft);
fullMatrix<double> &solright = solution->getGroupProxy(groupRight);
@@ -34,8 +36,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
for(int iFace=0 ; iFace<group->getNbElements();iFace++) {
- iElementL = group->getElementLeftId(iFace);
- iElementR = group->getElementRightId(iFace);
+ iElementL = left.getElementId(iFace);
+ iElementR = right.getElementId(iFace);
TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 07d94b7d9fa7f0cb23bbee469b6c25d445c036b9..c8f34f56478c338ff6d994dfffbc9be741e099fa 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -146,12 +146,14 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
- int nbNodesLeft = group.getGroupLeft().getNbNodes();
- int nbNodesRight = group.getGroupRight().getNbNodes();
+ const dgGroupOfConnections &left = group.getGroupOfConnections(0);
+ const dgGroupOfConnections &right = group.getGroupOfConnections(1);
+ const dgGroupOfElements &groupLeft = left.getGroupOfElements();
+ const dgGroupOfElements &groupRight = right.getGroupOfElements();
int nbFaces = group.getNbElements();
//get matrices needed to compute the gradient on both sides
- fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
- fullMatrix<double> &DPsiRightDx = group.getDPsiRightDxMatrix();
+ const fullMatrix<double> &DPsiLeftDx = left.getDPsiDx();
+ const fullMatrix<double> &DPsiRightDx = right.getDPsiDx();
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),nbFaces * _nbFields*2);
group.getCollocationMatrix().mult(solution, solutionQP);
@@ -161,26 +163,26 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
_cacheMapLeft->setNbEvaluationPoints(group.getNbIntegrationPoints());
_cacheMapRight->setNbEvaluationPoints(group.getNbIntegrationPoints());
fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
- int p = group.getGroupLeft().getOrder();
- int dim = group.getGroupLeft().getElement(0)->getDim();
+ int p = groupLeft.getOrder();
+ int dim = left.getElement(0)->getDim();
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
// B1 ) adjust the proxies for this element
// NB : the B1 section will almost completely disapear
// if we write a function (from the function class) for moving proxies across big matrices
// give the current elements to the dataCacheMap
- _cacheElementLeft.set(group.getElementLeft(iFace));
- _cacheElementRight.set(group.getElementRight(iFace));
+ _cacheElementLeft.set(left.getElement(iFace));
+ _cacheElementRight.set(left.getElement(iFace));
// proxies for the solution
_solutionQPLeft.setAsProxy(solutionQP, iFace*2*_nbFields, _nbFields);
_solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*_nbFields, _nbFields);
- _normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
+ _normals.setAsProxy(left.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
- dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
- dPsiRightDx.setAsProxy(DPsiRightDx,iFace*nbNodesRight,nbNodesRight);
- dofsLeft.setAsProxy(solutionLeft, _nbFields*group.getElementLeftId(iFace), _nbFields);
- dofsRight.setAsProxy(solutionRight, _nbFields*group.getElementRightId(iFace), _nbFields);
- _uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
- _uvwRight.setAsProxy(group.getIntegrationOnElementRight(iFace));
+ dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*groupLeft.getNbNodes(),groupLeft.getNbNodes());
+ dPsiRightDx.setAsProxy(DPsiRightDx,iFace*groupRight.getNbNodes(),groupRight.getNbNodes());
+ dofsLeft.setAsProxy(solutionLeft, _nbFields*left.getElementId(iFace), _nbFields);
+ dofsRight.setAsProxy(solutionRight, _nbFields*right.getElementId(iFace), _nbFields);
+ _uvwLeft.setAsProxy(left.getIntegrationPointsOnElement(iFace));
+ _uvwRight.setAsProxy(right.getIntegrationPointsOnElement(iFace));
// proxies for the flux
normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*2, _nbFields*2);
// B2 ) compute the gradient of the solution
@@ -200,8 +202,8 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
((dfl(iPt,k+_nbFields*0 )+dfr(iPt,k+_nbFields*0)) * _normals(0,iPt)
+(dfl(iPt,k+_nbFields*1 )+dfr(iPt,k+_nbFields*1)) * _normals(1,iPt)
+(dfl(iPt,k+_nbFields*2 )+dfr(iPt,k+_nbFields*2)) * _normals(2,iPt))/2;
- double minl = std::min(group.getElementVolumeLeft(iFace),
- group.getElementVolumeRight(iFace)
+ double minl = std::min(groupLeft.getElementVolume(left.getElementId(iFace)),
+ groupRight.getElementVolume(right.getElementId(iFace))
)/group.getInterfaceSurface(iFace);
double nu = std::max((*_maximumDiffusivityRight)(iPt,0),(*_maximumDiffusivityLeft)(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
@@ -229,9 +231,11 @@ void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofConta
{
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);
- compute1Group(faces,solInterface,solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupRight()),residuInterface);
- faces.mapFromInterface(_nbFields, residuInterface,residual.getGroupProxy(&faces.getGroupLeft()), residual.getGroupProxy(&faces.getGroupRight()));
+ const dgGroupOfElements *groupLeft = &faces.getGroupOfConnections(0).getGroupOfElements();
+ const dgGroupOfElements *groupRight = &faces.getGroupOfConnections(1).getGroupOfElements();
+ faces.mapToInterface(_nbFields, solution.getGroupProxy(groupLeft), solution.getGroupProxy(groupRight), solInterface);
+ compute1Group(faces,solInterface,solution.getGroupProxy(groupLeft), solution.getGroupProxy(groupRight),residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface,residual.getGroupProxy(groupLeft), residual.getGroupProxy(groupRight));
}
dgResidualInterface::~dgResidualInterface ()
@@ -253,10 +257,10 @@ void dgResidualBoundary::compute1Group(
{
//should be splitted like dgResidualInterface and dgResidualVolume
//but i do not do it know because dgResidualBoundary will probably disapear when we will have list of elements on interfaces
-
+ const dgGroupOfConnections &left = group.getGroupOfConnections(0);
+ const dgGroupOfElements &groupLeft = left.getGroupOfElements();
int nbFields = _claw.getNbFields();
- int nbNodesLeft = group.getGroupLeft().getNbNodes();
- const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getBoundaryTag());
+ const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getPhysicalTag());
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
group.getCollocationMatrix().mult(solution, solutionQP);
@@ -265,7 +269,7 @@ void dgResidualBoundary::compute1Group(
dataCacheMap cacheMapLeft;
cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
- fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
+ const fullMatrix<double> &DPsiLeftDx = left.getDPsiDx();
// provided dataCache
dataCacheDouble &uvw=cacheMapLeft.provideData("UVW",1,3);
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution",1,nbFields);
@@ -284,20 +288,20 @@ void dgResidualBoundary::compute1Group(
fullMatrix<double> normalFluxQP,dPsiLeftDx,dofsLeft;
- int p = group.getGroupLeft().getOrder();
- int dim = group.getGroupLeft().getElement(0)->getDim();
+ int p = groupLeft.getOrder();
+ int dim = left.getElement(0)->getDim();
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields, nbFields);
// ----- 2.3.1 --- provide the data to the cacheMap
solutionQPLeft.setAsProxy(solutionQP, iFace*nbFields, nbFields);
- normals.setAsProxy(group.getNormals(),iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
+ normals.setAsProxy(left.getNormals(),iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
- dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
- dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
+ dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*groupLeft.getNbNodes(), groupLeft.getNbNodes());
+ dofsLeft.setAsProxy(solutionLeft, nbFields*left.getElementId(iFace), nbFields);
- uvw.setAsProxy(group.getIntegrationOnElementLeft(iFace));
- cacheElementLeft.set(group.getElementLeft(iFace));
+ uvw.setAsProxy(left.getIntegrationPointsOnElement(iFace));
+ cacheElementLeft.set(left.getElement(iFace));
// compute the gradient of the solution
if(gradientSolutionLeft.somethingDependOnMe()){
@@ -317,7 +321,7 @@ void dgResidualBoundary::compute1Group(
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility :
for (int k=0;k<nbFields;k++) {
- double minl = group.getElementVolumeLeft(iFace)/group.getInterfaceSurface(iFace);
+ double minl = groupLeft.getElementVolume(left.getElementId(iFace))/group.getInterfaceSurface(iFace);
double nu = (*maximumDiffusivityLeft)(iPt,0);
if(maximumDiffusivityOut)
nu = std::max(nu,(*maximumDiffusivityOut)(iPt,0));
@@ -337,9 +341,10 @@ void dgResidualBoundary::computeAndMap1Group(dgGroupOfFaces &faces, dgDofContain
int _nbFields = _claw.getNbFields();
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
- faces.mapToInterface(_nbFields, solution.getGroupProxy(&faces.getGroupLeft()), solution.getGroupProxy(&faces.getGroupLeft()), solInterface);
- compute1Group(faces,solInterface,solution.getGroupProxy(&faces.getGroupLeft()),residuInterface);
- faces.mapFromInterface(_nbFields, residuInterface, residual.getGroupProxy(&faces.getGroupLeft()), residual.getGroupProxy(&faces.getGroupLeft()));
+ const dgGroupOfElements &groupLeft = faces.getGroupOfConnections(0).getGroupOfElements();
+ faces.mapToInterface(_nbFields, solution.getGroupProxy(&groupLeft), solution.getGroupProxy(&groupLeft), solInterface);
+ compute1Group(faces,solInterface,solution.getGroupProxy(&groupLeft),residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface, residual.getGroupProxy(&groupLeft), residual.getGroupProxy(&groupLeft));
}
void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual)
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index ffcf59692b78a40d116f8ec05d0262a056704622..75a12708b1fa4890d5141a61ba3bb3591e2f8472 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -151,7 +151,7 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
double dt,
dgDofContainer *solution,
int nStages,
- fullMatrix<double>&A, // c | A
+ fullMatrix<double> &A, // c | A
double *b, // Standard Butcher tableau notation :___|__
double *c // |b^T
)
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index b5686bdecb34ad90310fc1c43843c986a171c357..ac99bd9d2777daf808ea6ae546d2b90ec969e508 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -48,8 +48,8 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
const dgGroupOfElements *ge[2];
- ge[0]=&gf->getGroupLeft();
- ge[1]=&gf->getGroupRight();
+ ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
+ ge[1]=&gf->getGroupOfConnections(1).getGroupOfElements();
for(int i=0;i<2;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
@@ -64,7 +64,7 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
for(int iGroup=0;iGroup<gc->getNbBoundaryGroups();iGroup++){
dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
const dgGroupOfElements *ge[1];
- ge[0]=&gf->getGroupLeft();
+ ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
//ge[1]=&gf->getGroupRight();
for(int i=0;i<1;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
@@ -121,13 +121,13 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
gOFSet.insert(vfi.begin(),vfi.end());
for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
dgGroupOfFaces *faces=*it;
- if(faces->isBoundary()){
+ if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
//Msg::Info("Buffer face group %p is boundary in multirate",faces);
}
else{
- const dgGroupOfElements *gL=&faces->getGroupLeft();
- const dgGroupOfElements *gR=&faces->getGroupRight();
+ const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
+ const dgGroupOfElements *gR = &faces->getGroupOfConnections(1).getGroupOfElements();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
@@ -170,12 +170,12 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
std::vector<dgGroupOfFaces *> &vf=_bulkGroupsOfElements[exponent].second;
for(int i=0;i<vf.size();i++){
dgGroupOfFaces *faces=vf[i];
- if(faces->isBoundary()){
+ if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
}
else{
- const dgGroupOfElements *gL=&faces->getGroupLeft();
- const dgGroupOfElements *gR=&faces->getGroupRight();
+ const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
+ const dgGroupOfElements *gR = &faces->getGroupOfConnections(1).getGroupOfElements();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
diff --git a/Solver/function.h b/Solver/function.h
index bdfd9c64da7a4447e6402492021ad302549fa5cd..dbb42ef77ea8ea7a7ce9b0b8d6883d4701d12a82 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -80,14 +80,14 @@ class dataCacheDouble : public dataCache {
}
// take care if you use this you must ensure that the value pointed to are not modified
// without further call to setAsProxy because the dependencies won't be invalidate
- inline void setAsProxy(fullMatrix<double> &mat, int cShift, int c) {
+ inline void setAsProxy(const fullMatrix<double> &mat, int cShift, int c) {
_invalidateDependencies();
_value.setAsProxy(mat,cShift,c);
_valid=true;
}
// take care if you use this you must ensure that the value pointed to are not modified
// without further call to setAsProxy because the dependencies won't be invalidate
- inline void setAsProxy(fullMatrix<double> &mat) {
+ inline void setAsProxy(const fullMatrix<double> &mat) {
_invalidateDependencies();
_value.setAsProxy(mat,0,mat.size2());
_valid=true;