diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 023995badfffa3c0cc3effc14801ab8fbc988096..427a9b532104bbd940b263e0dfee7e3dfc6f681a 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -91,6 +91,10 @@ class fullMatrix
_own_data = true;
scale(0.);
}
+ fullMatrix(int r, int c, double *data) : _r(r), _c(c), _data(data),_own_data(false)
+ {
+ scale(0.);
+ }
fullMatrix(const fullMatrix<scalar> &other) : _r(other._r), _c(other._c)
{
_data = new scalar[_r * _c];
@@ -106,12 +110,13 @@ class fullMatrix
if(this != &other){
_r = other._r;
_c = other._c;
- if (_data) delete[] _data;
+ if (_data && _own_data) delete[] _data;
if ((_r==0)||(_c==0))
_data=0;
else
{
_data = new scalar[_r * _c];
+ _own_data=true;
for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
}
}
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index b4e51694ff118e341a57d6a6e464a7dc4d278fdb..a0346fb8a74eacc80d33463be5c01eb7dfac23e7 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -76,12 +76,11 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
if (claw.sourceTerm()){
fullMatrix<double> source(Source, iElement*claw.nbFields(),claw.nbFields());
(*claw.sourceTerm())(DGE,&source);
- //we assume constant mass matrix and constant mapping for now
- /*for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
for (int k=0;k<nbFields;k++)
source(iPt,k) *= detJ;
- }*/
+ }
}
delete gradSolutionQPe;
}
@@ -94,3 +93,36 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
residual.gemm(group.getSourceRedistributionMatrix(),Source);
}
}
+
+void dgAlgorithm::residualInterface ( //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();
+ // ----- 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*2);
+ 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*2*nbFields, nbFields );
+ fullMatrix<double> solutionQPRight (solutionQP, (iFace*2+1)*nbFields, nbFields );
+ fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*2*nbFields, nbFields*2);
+ dgFace DGF( group.getFace(iFace), group.getElementLeft(iFace), group.getElementRight(iFace),
+ solutionQPLeft, solutionQPRight, group.getIntegrationPointsMatrix());
+ // ----- 2.3.2 --- compute fluxes
+ (*claw.riemannSolver())(DGF,&normalFluxQP);
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ for (int k=0;k<nbFields*2;k++) {
+ normalFluxQP(iPt,k) *= detJ;
+ }
+ }
+ }
+ // ----- 3 ---- do the redistribution at face nodes using BLAS3
+ residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 745cd42e5b980d367524d4f7e7ce9710a85b4ff0..555b745c4246fd817aca13bbfd30d1e087dc1f1b 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -14,9 +14,11 @@ class dgAlgorithm {
const dgGroupOfElements & group,
const fullMatrix<double> &solution,
fullMatrix<double> &residual);
- void residualInterface ( /*dofManager &dof,*/
- const dgConservationLaw &law,
- const dgGroupOfFaces & group);
+ void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ 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);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index e2f102bca491faf0c265524278a9c790ae389c84..420b6b667fcde6db3e579738ae4ea91815fda80c 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -48,16 +48,24 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
)
{
// this is the biggest piece of data ... the mappings
- _mapping = new fullMatrix<double> (_elements.size(), 10 * _integration->size1());
+ int nbNodes = _fs.coefficients.size1();
+ _redistributionFluxes[0] = new fullMatrix<double> (nbNodes,_integration->size1());
+ _redistributionFluxes[1] = new fullMatrix<double> (nbNodes,_integration->size1());
+ _redistributionFluxes[2] = new fullMatrix<double> (nbNodes,_integration->size1());
+ _redistributionSource = new fullMatrix<double> (nbNodes,_integration->size1());
+ _collocation = new fullMatrix<double> (_integration->size1(),nbNodes);
+ _mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
+ _imass = new fullMatrix<double> (nbNodes,nbNodes*e.size());
+ double g[256][3],f[256];
for (int i=0;i<_elements.size();i++){
MElement *e = _elements[i];
+ fullMatrix<double> imass(*_imass,nbNodes*i,nbNodes);
for (int j=0;j< _integration->size1() ; j++ ){
+ _fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
double jac[3][3],ijac[3][3],detjac;
(*_mapping)(i,10*j + 9) =
- e->getJacobian ((*_integration)(j,0),
- (*_integration)(j,1),
- (*_integration)(j,2),
- jac);
+ e->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
+ const double weight = (*_integration)(j,3);
detjac=inv3x3(jac,ijac);
(*_mapping)(i,10*j + 0) = ijac[0][0];
(*_mapping)(i,10*j + 1) = ijac[0][1];
@@ -69,19 +77,17 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
(*_mapping)(i,10*j + 7) = ijac[2][1];
(*_mapping)(i,10*j + 8) = ijac[2][2];
(*_mapping)(i,10*j + 9) = detjac;
+ for (int k=0;k<_fs.coefficients.size1();k++){
+ for (int l=0;l<_fs.coefficients.size1();l++) {
+ imass(k,l) += f[k]*f[l]*weight*detjac;
+ }
+ }
}
+ imass.invertInPlace();
}
// redistribution matrix
// quadrature weight x parametric gradients in quadrature points
- _redistributionFluxes[0] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
- _redistributionFluxes[1] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
- _redistributionFluxes[2] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
- _redistributionSource = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
- _collocation = new fullMatrix<double> (_integration->size1(),_fs.coefficients.size1());
- _imass = new fullMatrix<double> (_fs.coefficients.size1(),_fs.coefficients.size1());
-
- double g[256][3],f[256];
for (int j=0;j<_integration->size1();j++) {
_fs.df((*_integration)(j,0),
(*_integration)(j,1),
@@ -96,12 +102,8 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
(*_redistributionFluxes[2])(k,j) = g[k][2] * weight;
(*_redistributionSource)(k,j) = f[k] * weight;
(*_collocation)(j,k) = f[k];
- for (int l=0;l<_fs.coefficients.size1();l++) {
- (*_imass)(k,l) += f[k]*f[l]*weight;
- }
}
}
- _imass->invertInPlace();
}
dgGroupOfElements::~dgGroupOfElements(){
@@ -123,17 +125,17 @@ void dgGroupOfFaces::createFaceElements (const std::vector<MFace> &topo_faces){
for (int i=0;i<topo_faces.size();i++){
// compute closures for the interpolation
int ithFace, sign, rot;
- _left[i]->getFaceInfo (topo_faces[i], ithFace, sign, rot);
+ getElementLeft(i)->getFaceInfo (topo_faces[i], ithFace, sign, rot);
_closuresLeft.push_back(&(_fsLeft->getFaceClosure(ithFace, sign, rot)));
- _right[i]->getFaceInfo (topo_faces[i], 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 = _right[i]->getFunctionSpace()->getFaceClosure(ithFace, sign, rot);
+ 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 = _left[i]->getVertex(iNod);
+ MVertex *v = getElementLeft(i)->getVertex(iNod);
_vertices.push_back(v);
}
// triangular face
@@ -186,16 +188,16 @@ void dgGroupOfFaces::createEdgeElements (const std::vector<MEdge> &topo_edges){
// compute all closures
for (int i=0;i<topo_edges.size();i++){
int ithEdge, sign;
- _left[i]->getEdgeInfo (topo_edges[i], ithEdge, sign);
+ getElementLeft(i)->getEdgeInfo (topo_edges[i], ithEdge, sign);
_closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
- _right[i]->getEdgeInfo (topo_edges[i], ithEdge, sign);
+ getElementRight(i)->getEdgeInfo (topo_edges[i], ithEdge, sign);
_closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));
// get the vertices of the edge
- const std::vector<int> & geomClosure = _right[i]->getFunctionSpace()->getEdgeClosure(ithEdge, sign);
+ const std::vector<int> & geomClosure = getElementRight(i)->getFunctionSpace()->getEdgeClosure(ithEdge, sign);
std::vector<MVertex*> _vertices;
for (int j=0; j<geomClosure.size() ; j++){
int iNod = geomClosure[j];
- MVertex *v = _left[i]->getVertex(iNod);
+ MVertex *v = getElementLeft(i)->getVertex(iNod);
_vertices.push_back(v);
}
switch(_vertices.size()){
@@ -222,20 +224,22 @@ void dgGroupOfFaces::init(int pOrder) {
}
dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces,
- const std::vector<MElement*> &l,
- const std::vector<MElement*> &r,
- int pOrder) : _left(l), _right(r),
- _fsLeft(_left[0]->getFunctionSpace (pOrder)), _fsRight (_right[0]->getFunctionSpace (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) : _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,
- const std::vector<MElement*> &l,
- const std::vector<MElement*> &r,
- int pOrder) : _left(l), _right(r),
- _fsLeft(_left[0]->getFunctionSpace (pOrder)), _fsRight (_right[0]->getFunctionSpace (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) : _group_list(group_list),_left(l), _right(r),
+ _fsLeft(getElementLeft(0)->getFunctionSpace (pOrder)), _fsRight (getElementRight(0)->getFunctionSpace (pOrder))
{
createEdgeElements (topo_edges);
init(pOrder);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 61fc87e429054d8fbbffa542792e37f464c7a78a..2b34680d6c53ce28713439131345d9230f42942c 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -55,6 +55,7 @@ class dgGroupOfElements {
fullMatrix<double> *_redistributionFluxes[3];
// redistribution for the source term
fullMatrix<double> *_redistributionSource;
+ // inverse mass matrix of all elements
fullMatrix<double> *_imass;
// dimension of the parametric space and of the real space
// may be different if the domain is a surface in 3D (manifold)
@@ -75,24 +76,35 @@ public:
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
- inline const fullMatrix<double> & getInverseMassMatrix () const {return *_imass;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
inline double getInvJ (int iElement, int iGaussPoint, int i, int j) const {return (*_mapping)(iElement, 10*iGaussPoint + i + 3*j);}
- inline fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
+ inline fullMatrix<double> getInverseMassMatrix (int iElement) const {return fullMatrix<double>(*_imass,iElement*getNbNodes(),getNbNodes());}
+ inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
};
-/*class dgFace {
+class dgFace {
int nbGaussPoints;
MElement *_left, *_right;
MElement *_face;
double *_normals;
- double *_detJ;
- int *_closureLeft,*_closureRight;
+ const fullMatrix<double> *_solutionRight, *_solutionLeft, *_integration;
public:
- dgFace ();
+ dgFace (MElement *face,MElement *left, MElement *right,
+ const fullMatrix<double> &solRight,
+ const fullMatrix<double> &solLeft,
+ const fullMatrix<double> &integration
+ ) : _left(left), _right(right), _face(face),_solutionRight(&solRight),_solutionLeft(&solLeft),_integration(&integration)
+ {}
+ inline const fullMatrix<double> &solutionRight() const { return *_solutionRight; }
+ inline const fullMatrix<double> &solutionLeft() const { return *_solutionLeft; }
+ inline const fullMatrix<double> &integration() const { return *_integration; }
+ inline MElement *left() const { return _left;}
+ inline MElement *right() const { return _right;}
+ inline MElement *face() const { return _face;}
};
-*/
+
class dgGroupOfFaces {
+ std::vector<dgGroupOfElements*> _group_list;
void createFaceElements (const std::vector<MFace> &topo_faces);
void createEdgeElements (const std::vector<MEdge> &topo_faces);
void computeFaceNormals();
@@ -100,7 +112,8 @@ class dgGroupOfFaces {
// the group has always the same types for left and right
const functionSpace *_fsLeft,*_fsRight, *_fsFace;
// N elements in the group
- std::vector<MElement*> _left, _right, _faces;
+ std::vector<std::pair<int,int> >_left, _right;
+ 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
@@ -122,15 +135,30 @@ class dgGroupOfFaces {
//common part of the 3 constructors
void init(int pOrder);
public:
+ inline MElement* getElementLeft (int i) const {return _group_list[_left[i].first]->getElement(_left[i].second);}
+ inline MElement* getElementRight (int i) const {return _group_list[_right[i].first]->getElement(_right[i].second);}
+ inline MElement* getFace (int iElement) const {return _faces[iElement];}
+ const std::vector<int> * getClosureLeft(int iFace) const{ return _closuresLeft[iFace];}
+ const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
dgGroupOfFaces (const std::vector<MFace> &faces,
- const std::vector<MElement*> &l,
- const std::vector<MElement*> &r,
+ 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<MElement*> &l,
- const std::vector<MElement*> &r,
+ 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 ();
+ //this part is common with dgGroupOfElements, we should try polymorphism
+ inline int getNbElements() const {return _faces.size();}
+ inline int getNbNodes() const {return _collocation->size1();}
+ inline int getNbIntegrationPoints() const {return _collocation->size1();}
+ inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
+ inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
+ inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
+ inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iElement, iGaussPoint);}
};
diff --git a/Solver/dgMainTest.cc b/Solver/dgMainTest.cc
index f93e71a079742bba2304ae82afccf7368357f22b..cc0ea70b0f7c678f7c18e069dfc7577b70d902a6 100644
--- a/Solver/dgMainTest.cc
+++ b/Solver/dgMainTest.cc
@@ -8,26 +8,33 @@
#include "MElement.h"
-void print (const char *filename, std::vector<MElement *> &els, double *v);
-std::vector<MElement *> get_all_tri(GModel *model);
-
+void print (const char *filename,const dgGroupOfElements &els, double *v);
+void buildGroupAllTri(GModel *model, int order, //in
+ std::vector<dgGroupOfElements*> &elements, //out
+ std::vector<dgGroupOfFaces*> &faces); //out
int main(int argc, char **argv){
GmshMergeFile("input/mesh1.msh");
- std::vector<MElement *> all_tri=get_all_tri(GModel::current());
- dgGroupOfElements group(all_tri,1);
- fullMatrix<double> sol(3,all_tri.size());
- fullMatrix<double> residu(3,all_tri.size());
+ std::vector<dgGroupOfElements*> elements;
+ std::vector<dgGroupOfFaces*> faces;
+ buildGroupAllTri(GModel::current(),1,elements,faces);
+ int nbNodes=elements[0]->getNbNodes();
+ fullMatrix<double> sol(nbNodes,elements[0]->getNbElements());
+ fullMatrix<double> residu(nbNodes,elements[0]->getNbElements());
dgAlgorithm algo;
dgConservationLaw *law = dgNewConservationLawAdvection();
- algo.residualVolume(*law,group,sol,residu);
- sol.gemm(group.getInverseMassMatrix(),residu);
- print("test.pos",all_tri,&sol(0,0));
+ algo.residualVolume(*law,*elements[0],sol,residu);
+ for(int i=0;i<elements[0]->getNbElements();i++) {
+ fullMatrix<double> residuEl(residu,i,1);
+ fullMatrix<double> solEl(sol,i,1);
+ solEl.gemm(elements[0]->getInverseMassMatrix(i),residuEl);
+ }
+ print("test.pos",*elements[0],&sol(0,0));
}
-
-
-std::vector<MElement *> get_all_tri(GModel *model){
+void buildGroupAllTri(GModel *model, int order, //in
+ std::vector<dgGroupOfElements*> &elements, //out
+ std::vector<dgGroupOfFaces*> &faces){ //out
std::vector<GEntity*> entities;
model->getEntities(entities);
std::vector<MElement *> all_tri;
@@ -36,15 +43,14 @@ std::vector<MElement *> get_all_tri(GModel *model){
for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
all_tri.push_back((*itent)->getMeshElement(iel));
}
- return all_tri;
+ elements.push_back(new dgGroupOfElements(all_tri,order));
}
-void print (const char *filename, std::vector<MElement *> &els, double *v) {
+void print (const char *filename,const dgGroupOfElements &els, double *v) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);
- int i=0;
- for(std::vector<MElement *>::iterator itel= els.begin();itel!=els.end();itel++){
- MElement *el = *itel;
+ for(int iel=0;iel<els.getNbElements();iel++){
+ MElement *el = els.getElement(iel);
fprintf(file,"ST (");
for (int iv=0; iv<el->getNumVertices(); iv++) {
MVertex *vertex = el->getVertex(iv);
@@ -53,7 +59,7 @@ void print (const char *filename, std::vector<MElement *> &els, double *v) {
}
fprintf(file,"{");
for (int iv=0; iv<el->getNumVertices(); iv++)
- fprintf(file,"%e%c ",v[i++],iv==el->getNumVertices()-1?'}':',');
+ fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
fprintf(file,";\n");
}
fprintf(file,"};");