diff --git a/Solver/_ b/Solver/_
new file mode 100644
index 0000000000000000000000000000000000000000..0f2c8223facb9a348b52c3edb0df2d2cb909cfed
--- /dev/null
+++ b/Solver/_
@@ -0,0 +1,330 @@
+#include "dgGroupOfElements.h"
+#include "MElement.h"
+#include "functionSpace.h"
+#include "Numeric.h"
+#include "MTriangle.h"
+#include "MLine.h"
+
+static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
+ int npts;
+ IntPt *pts;
+ e->getIntegrationPoints(2*p+1, &npts, &pts);
+ fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
+ for (int i=0;i<npts;i++){
+ (*m)(i,0) = pts[i].pt[0];
+ (*m)(i,1) = pts[i].pt[1];
+ (*m)(i,2) = pts[i].pt[2];
+ (*m)(i,3) = pts[i].weight;
+ }
+ return m;
+}
+
+static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
+ const functionSpace *fsFace,
+ const fullMatrix<double> &intgFace,
+ const functionSpace *fsElement,
+ const std::vector <int> *closure) {
+ int npts=intgFace.size1();
+ fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
+ double f[256];
+ for (int i=0;i<npts;i++){
+ fsFace->f(intgFace(i,0),intgFace(i,1),intgFace(i,2),f);
+ for(size_t j=0; j<closure->size();j++){
+ int jNod=(*closure)[j];
+ (*m)(i,0) += f[j] * fsElement->points(jNod,0);
+ (*m)(i,1) += f[j] * fsElement->points(jNod,1);
+ (*m)(i,2) += f[j] * fsElement->points(jNod,2);
+ (*m)(i,3) = intgFace(i,3);
+ }
+ }
+ return m;
+}
+
+
+dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
+ : _elements(e),
+ _fs(*_elements[0]->getFunctionSpace(polyOrder)),
+ _integration(dgGetIntegrationRule (_elements[0], polyOrder)
+ )
+{
+ _dimUVW = _dimXYZ = e[0]->getDim();
+ // this is the biggest piece of data ... the mappings
+ 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);
+ 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];
+ (*_mapping)(i,10*j + 2) = ijac[0][2];
+ (*_mapping)(i,10*j + 3) = ijac[1][0];
+ (*_mapping)(i,10*j + 4) = ijac[1][1];
+ (*_mapping)(i,10*j + 5) = ijac[1][2];
+ (*_mapping)(i,10*j + 6) = ijac[2][0];
+ (*_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
+
+ for (int j=0;j<_integration->size1();j++) {
+ _fs.df((*_integration)(j,0),
+ (*_integration)(j,1),
+ (*_integration)(j,2), g);
+ _fs.f((*_integration)(j,0),
+ (*_integration)(j,1),
+ (*_integration)(j,2), f);
+ const double weight = (*_integration)(j,3);
+ for (int k=0;k<_fs.coefficients.size1();k++){
+ (*_redistributionFluxes[0])(k,j) = g[k][0] * weight;
+ (*_redistributionFluxes[1])(k,j) = g[k][1] * weight;
+ (*_redistributionFluxes[2])(k,j) = g[k][2] * weight;
+ (*_redistributionSource)(k,j) = f[k] * weight;
+ (*_collocation)(j,k) = f[k];
+ }
+ }
+}
+
+dgGroupOfElements::~dgGroupOfElements(){
+ delete _integration;
+ delete _redistributionFluxes[0];
+ delete _redistributionFluxes[1];
+ delete _redistributionFluxes[2];
+ delete _redistributionSource;
+ delete _mapping;
+ delete _collocation;
+ delete _imass;
+}
+
+
+void dgGroupOfFaces::computeFaceNormals () {
+ double g[256][3];
+ _normals = new fullMatrix<double> (3,_fsFace->points.size1()*_faces.size());
+ int index = 0;
+ for (size_t i=0; i<_faces.size();i++){
+ const std::vector<int> &closure=*_closuresLeft[i];
+ fullMatrix<double> *intLeft=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,&closure);
+ for (int j=0; j<intLeft->size1(); j++){
+ _fsLeft->df((*intLeft)(j,0),(*intLeft)(j,1),(*intLeft)(j,2),g);
+ double &nx=(*_normals)(0,index);
+ double &ny=(*_normals)(1,index);
+ double &nz=(*_normals)(2,index);
+ for (size_t k=0; k<closure.size(); k++){
+ nx += g[closure[k]][0];
+ ny += g[closure[k]][1];
+ nz += g[closure[k]][2];
+ }
+ double norm = sqrt(nx*nx+ny*ny+nz*nz);
+ nx/=norm;
+ ny/=norm;
+ nz/=norm;
+ index++;
+ }
+ delete intLeft;
+ }
+}
+
+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)));
+ // 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]) );
+ // triangular face
+ if (topoFace.getNumVertices() == 3){
+ switch(vertices.size()){
+ case 3 : _faces.push_back(new MTriangle (vertices) ); break;
+ case 6 : _faces.push_back(new MTriangle6 (vertices) ); break;
+ case 10 : _faces.push_back(new MTriangleN (vertices, 3) ); break;
+ case 15 : _faces.push_back(new MTriangleN (vertices, 4) ); break;
+ case 21 : _faces.push_back(new MTriangleN (vertices, 5) ); break;
+ default : throw;
+ }
+ }
+ // quad face 2 do
+ else throw;
+}
+static std::vector<int> *fakeClosure2d(const functionSpace *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);
+ std::vector<MVertex*> vertices;
+ for(int j=0;j<topoEdge.getNumVertices();j++)
+ vertices.push_back(topoEdge.getVertex(j));
+ for (int j=0; j<geomClosure.size() ; j++)
+ vertices.push_back( elRight.getVertex(geomClosure[j]) );
+ switch(vertices.size()){
+ case 2 : _faces.push_back(new MLine (vertices) ); break;
+ case 3 : _faces.push_back(new MLine3 (vertices) ); break;
+ default : _faces.push_back(new MLineN (vertices) ); break;
+ }
+}
+
+void dgGroupOfFaces::init(int pOrder) {
+ _fsFace = _faces[0]->getFunctionSpace (pOrder);
+ _integration=dgGetIntegrationRule (_faces[0],pOrder);
+ _redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
+ _collocation = new fullMatrix<double> (_fsFace->coefficients.size1(), _integration->size1());
+ _detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
+ double f[256];
+ for (int j=0;j<_integration->size1();j++) {
+ _fsFace->f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
+ const double weight = (*_integration)(j,3);
+ for (int k=0;k<_fsFace->coefficients.size1();k++){
+ (*_redistribution)(k,j) = f[j] * weight;
+ (*_collocation)(k,j) = f[k];
+ }
+ }
+ for (int i=0;i<_faces.size();i++){
+ MElement *f = _faces[i];
+ for (int j=0;j< _integration->size1() ; j++ ){
+ double jac[3][3],ijac[3][3],detjac;
+ f->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
+ const double weight = (*_integration)(j,3);
+ (*_detJac)(j,i) = inv3x3(jac,ijac);
+ }
+ }
+ computeFaceNormals();
+}
+
+dgGroupOfFaces::~dgGroupOfFaces()
+{
+
+}
+
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
+ _groupLeft(elGroup),_groupRight(elGroup)
+{
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
+ switch (_groupLeft.getElement(0)->getDim()) {
+ case 2 : {
+ std::map<MEdge,int,Less_Edge> edgeMap;
+ 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(edgeMap.find(edge) == edgeMap.end()){
+ edgeMap[edge] = i;
+ }else{
+ addEdge(edge,edgeMap[edge],i);
+ }
+ }
+ }
+ break;
+ }
+ case 3 : {
+ std::map<MFace,int,Less_Face> faceMap;
+ 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(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = i;
+ }else{
+ addFace(face,faceMap[face],i);
+ }
+ }
+ }
+ break;
+ }
+ default : throw;
+ }
+ init(pOrder);
+}
+
+void dgGroupOfFaces::mapToInterface ( int nFields,
+ const fullMatrix<double> &vLeft,
+ 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++){
+ printf("closure size=%i\n",closureLeft.size());
+ for(size_t j =0; j < closureLeft.size(); j++){
+ v(j, i*2*nFields + iField) = vLeft(closureLeft[j], _left[i]*nFields + iField);
+ printf("vv=%e\n",v(j, i*2*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);
+ }
+ }
+ v.print();
+}
+
+void dgGroupOfFaces::mapFromInterface ( int nFields,
+ const fullMatrix<double> &v,
+ fullMatrix<double> &vLeft,
+ 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);
+ }
+ }
+}
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 085b5d20dcf58427181afd46d52bbb98ac140624..8025824daed907c7fb0ebbe45d2daa5fd3ac1854 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -37,9 +37,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
// ----- 2.2 --- allocate parametric fluxes (computed in UVW coordinates) for all elements at all integration points
- fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
- fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
- fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints())};
+ fullMatrix<double> Fuvw[3] = {fullMatrix<double> ( group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
+ fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
+ fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields)};
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
@@ -48,11 +48,10 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );
else gradSolutionQPe = new fullMatrix<double>;
dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
- if(claw.convectiveFlux() || claw.diffusiveFlux()){
+ if(claw.convectiveFlux() || claw.diffusiveFlux()) {
// ----- 2.3.2 --- compute fluxes in XYZ coordinates
if (claw.convectiveFlux()) (*claw.convectiveFlux())(DGE,fConv);
if (claw.diffusiveFlux()) (*claw.diffusiveFlux())(DGE,fDiff);
-
// ----- 2.3.3 --- compute fluxes in UVW coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
// ----- 2.3.3.1 --- get a proxy on the big local flux matrix
@@ -61,14 +60,13 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
for (int iXYZ=0;iXYZ<group.getDimXYZ();iXYZ++) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
// get the right inv jacobian matrix dU/dX element
- const double invJ = group.getInvJ (iElement, iPt, iXYZ, iUVW);
+ const double invJ = group.getInvJ (iElement, iPt, iUVW, iXYZ);
// get the detJ at this point
const double detJ = group.getDetJ (iElement, iPt);
const double factor = invJ * detJ;
// compute fluxes in the reference coordinate system at this integration point
- for (int k=0;k<nbFields;k++) {
+ for (int k=0;k<nbFields;k++)
fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
- }
}
}
}
@@ -86,8 +84,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
if(claw.convectiveFlux() || claw.diffusiveFlux()){
- for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
+ for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
+ }
}
if(claw.sourceTerm()){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
@@ -103,7 +102,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
{
int nbFields = claw.nbFields();
// ----- 1 ---- get the solution at quadrature points
- fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields*2);
group.getCollocationMatrix().mult(solution, solutionQP);
// ----- 2 ---- compute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields*2);
@@ -113,17 +112,27 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
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());
+ solutionQPLeft, solutionQPRight, group.getIntegrationPointsMatrix(),group.getNormals(iFace));
// ----- 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++) {
+ 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);
}
+void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol)
+{
+ for(int i=0;i<group.getNbElements();i++) {
+ fullMatrix<double> residuEl(residu,i,1);
+ fullMatrix<double> solEl(sol,i,1);
+ solEl.gemm(group.getInverseMassMatrix(i),residuEl);
+ }
+}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 555b745c4246fd817aca13bbfd30d1e087dc1f1b..2af608c6df2d32098b24556d7dde4194e05a508c 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -6,6 +6,7 @@
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
+class dgTerm;
class dgAlgorithm {
public :
@@ -22,6 +23,10 @@ class dgAlgorithm {
void residualBoundaryConditions ( /*dofManager &dof,*/
const dgConservationLaw &law,
const dgGroupOfElements & group);
+ void multAddInverseMassMatrix ( /*dofManager &dof,*/
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol);
};
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index 3c1d9eb7d4cd8151d12402f2f220b7ebc1fb1ec8..269d43538a685a3740f2dc32c5b62909546815ab 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -3,24 +3,68 @@
#include "dgGroupOfElements.h"
#include "SPoint3.h"
#include "MElement.h"
-class testSourceTerm : public dgTerm {
- void operator () (const dgElement &el, fullMatrix<double> fcx[]) const{
+void getV(SPoint3 &p, double (&v)[3]){
+ 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;
+}
+
+class dgAdvectionFluxTerm : public dgTerm {
+ void operator () (const dgElement &el, fullMatrix<double> fcx[]) const
+ {
const fullMatrix<double> &sol = el.solution();
- const fullMatrix<double> &qp = el.integration();
- SPoint3 p;
+ const fullMatrix<double> &QP = el.integration();
for(int i=0; i< sol.size1(); i++) {
- el.element()->pnt(qp(i,0),qp(i,1),qp(i,2),p);
- //printf("%e - %e (%i)\n",p.x(),sol(i,0),sol.size2());
- fcx[0](i,0) = (p.x()*p.x()+p.y()*p.y()) - sol(i,0);
+ SPoint3 p;
+ el.element()->pnt(QP(i,0),QP(i,1),QP(i,2),p);
+ double v[3];
+ getV(p,v);
+ fcx[0](i,0) = sol(i,0)*v[0];
+ fcx[1](i,0) = sol(i,0)*v[1];
+ fcx[2](i,0) = sol(i,0)*v[2];
+ }
+ }
+};
+
+class dgAdvectionRiemanTerm : public dgFaceTerm{
+ public:
+ void operator () (const dgFace &face, fullMatrix<double> fcx[]) const
+ {
+ const fullMatrix<double> &solLeft = face.solutionLeft();
+ const fullMatrix<double> &solRight = face.solutionRight();
+ const fullMatrix<double> &normals = face.normals();
+ const fullMatrix<double> &QP = face.integration();
+ for(int i=0; i< solLeft.size1(); i++) {
+ SPoint3 p;
+ face.face()->pnt(QP(i,0),QP(i,1),QP(i,2),p);
+ double v[3];
+ getV(p,v);
+ 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;
+ }else{
+ fcx[0](i,0) = -solRight(i,0)*un;
+ fcx[0](i,1) = solRight(i,0)*un;
+ }
}
}
};
class dgConservationLawAdvection : public dgConservationLaw {
public:
- dgConservationLawAdvection() {
+ dgConservationLawAdvection()
+ {
_nbf = 1;
- _source = new testSourceTerm;
+ _convective = new dgAdvectionFluxTerm;
+ _riemannSolver = new dgAdvectionRiemanTerm;
+ }
+ ~dgConservationLawAdvection()
+ {
+ delete _convective;
+ delete _riemannSolver;
}
};
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index de8bfe3ddb0ed0699672a47a8bea0153fa8e0424..6d10a34e4ebbc1a8a716ca4af83cb84355afaea4 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -31,9 +31,8 @@ static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
fsFace->f(intgFace(i,0),intgFace(i,1),intgFace(i,2),f);
for(size_t j=0; j<closure->size();j++){
int jNod=(*closure)[j];
- (*m)(i,0) += f[j] * fsElement->points(jNod,0);
- (*m)(i,1) += f[j] * fsElement->points(jNod,1);
- (*m)(i,2) += f[j] * fsElement->points(jNod,2);
+ for(int k=0;k<fsElement->points.size2();k++)
+ (*m)(i,k) += f[j] * fsElement->points(jNod,k);
(*m)(i,3) = intgFace(i,3);
}
}
@@ -47,6 +46,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
_integration(dgGetIntegrationRule (_elements[0], polyOrder)
)
{
+ _dimUVW = _dimXYZ = e[0]->getDim();
// this is the biggest piece of data ... the mappings
int nbNodes = _fs.coefficients.size1();
_redistributionFluxes[0] = new fullMatrix<double> (nbNodes,_integration->size1());
@@ -120,21 +120,28 @@ dgGroupOfElements::~dgGroupOfElements(){
void dgGroupOfFaces::computeFaceNormals () {
double g[256][3];
- _normals = new fullMatrix<double> (_fsFace->points.size1()*_faces.size(),3);
+ _normals = new fullMatrix<double> (3,_fsFace->points.size1()*_faces.size());
int index = 0;
for (size_t i=0; i<_faces.size();i++){
const std::vector<int> &closure=*_closuresLeft[i];
fullMatrix<double> *intLeft=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,&closure);
+ double jac[3][3],ijac[3][3];
for (int j=0; j<intLeft->size1(); j++){
_fsLeft->df((*intLeft)(j,0),(*intLeft)(j,1),(*intLeft)(j,2),g);
- double &nx=(*_normals)(index,0);
- double &ny=(*_normals)(index,1);
- double &nz=(*_normals)(index,2);
+ getElementLeft(i)->getJacobian ((*intLeft)(j,0), (*intLeft)(j,1), (*intLeft)(j,2), jac);
+ inv3x3(jac,ijac);
+ double &nx=(*_normals)(0,index);
+ double &ny=(*_normals)(1,index);
+ double &nz=(*_normals)(2,index);
+ double nu=0,nv=0,nw=0;
for (size_t k=0; k<closure.size(); k++){
- nx += g[closure[k]][0];
- ny += g[closure[k]][1];
- nz += g[closure[k]][2];
+ nu += g[closure[k]][0];
+ nv += g[closure[k]][1];
+ nw += g[closure[k]][2];
}
+ nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
+ ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];
+ nz = nu*ijac[2][0]+nv*ijac[2][1]+nw*ijac[2][2];
double norm = sqrt(nx*nx+ny*ny+nz*nz);
nx/=norm;
ny/=norm;
@@ -179,6 +186,19 @@ 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);
@@ -187,13 +207,17 @@ void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight){
MElement &elLeft = *_groupLeft.getElement(iElLeft);
int ithEdge, sign;
elLeft.getEdgeInfo (topoEdge, ithEdge, sign);
- _closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
+ _closuresLeft.push_back(fakeClosure2d(_fsLeft, ithEdge, sign));
elRight.getEdgeInfo (topoEdge, ithEdge, sign);
- _closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));
+ _closuresRight.push_back(fakeClosure2d(_fsRight, ithEdge, sign));
const std::vector<int> &geomClosure = elRight.getFunctionSpace()->getEdgeClosure(ithEdge, sign);
std::vector<MVertex*> vertices;
- for(int j=0;j<topoEdge.getNumVertices();j++)
- vertices.push_back(topoEdge.getVertex(j));
+ if(sign==1)
+ for(int j=0;j<topoEdge.getNumVertices();j++)
+ vertices.push_back(topoEdge.getVertex(j));
+ else
+ 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]) );
switch(vertices.size()){
@@ -207,21 +231,32 @@ void dgGroupOfFaces::init(int pOrder) {
_fsFace = _faces[0]->getFunctionSpace (pOrder);
_integration=dgGetIntegrationRule (_faces[0],pOrder);
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
- _collocation = new fullMatrix<double> (_fsFace->coefficients.size1(), _integration->size1());
+ _collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
+ _detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
double f[256];
for (int j=0;j<_integration->size1();j++) {
_fsFace->f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
const double weight = (*_integration)(j,3);
for (int k=0;k<_fsFace->coefficients.size1();k++){
- (*_redistribution)(k,j) = f[j] * weight;
- (*_collocation)(k,j) = f[k];
+ (*_redistribution)(k,j) = f[k] * weight;
+ (*_collocation)(j,k) = f[k];
}
}
+ for (int i=0;i<_faces.size();i++){
+ MElement *f = _faces[i];
+ for (int j=0;j< _integration->size1() ; j++ ){
+ double jac[3][3];
+ (*_detJac)(j,i)=f->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
+ }
+ }
+ computeFaceNormals();
}
dgGroupOfFaces::~dgGroupOfFaces()
{
-
+ delete _redistribution;
+ delete _collocation;
+ delete _detJac;
}
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
@@ -274,8 +309,9 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
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++)
+ 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);
}
@@ -293,9 +329,9 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
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);
+ 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);
+ 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 f00b80a7b4802b3d9c0f365e8d14456b526cf6d2..a2f5035d5a72d4eb7074ad7fac40f8d547a9e708 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -86,18 +86,19 @@ class dgFace {
int nbGaussPoints;
MElement *_left, *_right;
MElement *_face;
- double *_normals;
- const fullMatrix<double> *_solutionRight, *_solutionLeft, *_integration;
+ const fullMatrix<double> *_solutionRight, *_solutionLeft, *_integration,*_normals;
public:
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)
+ const fullMatrix<double> &solRight,
+ const fullMatrix<double> &integration,
+ const fullMatrix<double> &normals
+ ) : _left(left), _right(right), _face(face),_solutionRight(&solRight),_solutionLeft(&solLeft),_integration(&integration),_normals(&normals)
{}
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 const fullMatrix<double> &normals() const { return *_normals; }
inline MElement *left() const { return _left;}
inline MElement *right() const { return _right;}
inline MElement *face() const { return _face;}
@@ -140,6 +141,7 @@ public:
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];}
+ inline fullMatrix<double> getNormals (int iFace) const {return fullMatrix<double>(*_normals,iFace*getNbIntegrationPoints(),getNbIntegrationPoints());}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
virtual ~dgGroupOfFaces ();
//this part is common with dgGroupOfElements, we should try polymorphism
@@ -149,7 +151,7 @@ public:
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);}
+ inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iGaussPoint,iElement);}
//keep this outside the Algorithm because this is the only place where data overlap
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);
diff --git a/Solver/dgMainTest.cc b/Solver/dgMainTest.cc
index a0ca1d73709dea69ec2b6f8142d5a15d90c56712..451bfde1b889cf96f8364e8e9a297c4eae78a8f6 100644
--- a/Solver/dgMainTest.cc
+++ b/Solver/dgMainTest.cc
@@ -11,29 +11,70 @@
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{
+ const fullMatrix<double> &sol = el.solution();
+ const fullMatrix<double> &qp = el.integration();
+ 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);
+ }
+ }
+};
+
+class dgConservationLawInitialCondition : public dgConservationLaw {
+ public:
+ dgConservationLawInitialCondition() {
+ _nbf = 1;
+ _source = new testSourceTerm;
+ }
+ ~dgConservationLawInitialCondition() {
+ delete _source;
+ }
+};
+
int main(int argc, char **argv){
GmshMergeFile("input/mesh1.msh");
+ //GmshMergeFile("square2.msh");
std::vector<MElement *> allTri=getAllTri(GModel::current());
int order=1;
dgGroupOfElements elements(allTri,order);
dgGroupOfFaces faces(elements,order);
- int nbNodes=elements.getNbNodes();
+ //dgGroupOfFaces faces(elements0,elements1,order);
+ int nbNodes = elements.getNbNodes();
fullMatrix<double> sol(nbNodes,elements.getNbElements());
- fullMatrix<double> solInterface(nbNodes,faces.getNbElements()*2);
fullMatrix<double> residu(nbNodes,elements.getNbElements());
- fullMatrix<double> residuInterface(nbNodes,faces.getNbElements()*2);
dgAlgorithm algo;
+ // initial condition
+ {
+ dgConservationLawInitialCondition initLaw;
+ algo.residualVolume(initLaw,elements,sol,residu);
+ algo.multAddInverseMassMatrix(elements,residu,sol);
+ }
+ print("init.pos",elements,&sol(0,0));
+ //advection
dgConservationLaw *law = dgNewConservationLawAdvection();
- algo.residualVolume(*law,elements,sol,residu);
- faces.mapToInterface(1, sol, sol, solInterface);
- algo.residualInterface(*law,faces,solInterface,residuInterface);
- faces.mapFromInterface(1, residuInterface, residu, residu);
- for(int i=0;i<elements.getNbElements();i++) {
- fullMatrix<double> residuEl(residu,i,1);
- fullMatrix<double> solEl(sol,i,1);
- solEl.gemm(elements.getInverseMassMatrix(i),residuEl);
+ 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);
+ }
+ residu.scale(0.01);
+ algo.multAddInverseMassMatrix(elements,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("test.pos",elements,&sol(0,0));
+ delete law;
}
std::vector<MElement *> getAllTri(GModel *model){