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

dg non constant mass matrix and progress on interface term

parent 0c181d2b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/fullMatrix.h
+ 6
− 1
View file @ 3dea4b23
......@@ -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];
}
}
......
Solver/dgAlgorithm.cpp
+ 35
− 3
View file @ 3dea4b23
......@@ -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);
}
Solver/dgAlgorithm.h
+ 5
− 3
View file @ 3dea4b23
......@@ -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);
......
Solver/dgGroupOfElements.cpp
+ 37
− 33
View file @ 3dea4b23
......@@ -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);
......
Solver/dgGroupOfElements.h
+ 40
− 12
View file @ 3dea4b23
......@@ -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);}
};
......
Solver/dgMainTest.cc
+ 25
− 19
View file @ 3dea4b23
......@@ -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,"};");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment