dg advection 2d p1 p1 ok

Solver/_

+#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);
+    }
+  }
+}

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
@@ -52,7 +52,6 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
       // ----- 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,18 +60,17 @@ 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;
           }
         } 
       }
     }
-    }
     if (claw.sourceTerm()){
       fullMatrix<double> source(Source, iElement*claw.nbFields(),claw.nbFields());
       (*claw.sourceTerm())(DGE,&source); 
@@ -86,9 +84,10 @@ 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);
+  }
+}

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);
 };
 
 

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;
   }
 };
 

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;
+  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);
     }
   }
 }

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);

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);
+  //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();
+  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);
-  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);
     }
-  print("test.pos",elements,&sol(0,0));
+    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));
+    }
+  }
+  delete law;
 }
 
 std::vector<MElement *> getAllTri(GModel *model){