work on dg solver

Numeric/fullMatrix.cpp

@@ -39,7 +39,7 @@ extern "C" {
 }
 
 template<> 
-void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c)
+void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c)const
 {
   int M = c.size1(), N = c.size2(), K = _c;
   int LDA = _r, LDB = b.size1(), LDC = c.size1();
@@ -50,7 +50,7 @@ void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c
 
 template<> 
 void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<double> > &b, 
-                                             fullMatrix<std::complex<double> > &c)
+                                             fullMatrix<std::complex<double> > &c)const
 {
   int M = c.size1(), N = c.size2(), K = _c;
   int LDA = _r, LDB = b.size1(), LDC = c.size1();

Numeric/fullMatrix.h

@@ -79,11 +79,11 @@ class fullMatrix
   int _r, _c;
   scalar *_data;
  public:
-  fullMatrix(fullMatrix<scalar> &original, int r_start, int r){
-    _r = r;
-    _c = original._c;
+  fullMatrix(fullMatrix<scalar> &original, int c_start, int c){
+    _c = c;
+    _r = original._r;
     _own_data = false;
-    _data = original._data + r_start * _c;
+    _data = original._data + c_start * _r;
   }
   fullMatrix(int r, int c) : _r(r), _c(c)
   {
@@ -132,7 +132,7 @@ class fullMatrix
       for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
         (*this)(desti, destj) = a(i, j);
   }
-  void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)
+  void mult(const fullMatrix<scalar> &b, fullMatrix<scalar> &c)const
 #if !defined(HAVE_BLAS)
   {
     c.scale(0.);

Solver/CMakeLists.txt

@@ -12,6 +12,7 @@ set(SRC
   eigenSolver.cpp
   dgGroupOfElements.cpp
   dgAlgorithm.cpp
+  dgConservationLawAdvection.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/dgAlgorithm.cpp

@@ -9,22 +9,24 @@
 
 void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
 				   const dgConservationLaw &claw,   // the conservation law
-				   const dgGroupOfElements & group ) // the residual
+				   const dgGroupOfElements & group, 
+           const fullMatrix<double> &solution,
+           fullMatrix<double> &residual // the residual
+           )
 { 
   // ----- 1 ----  get the solution at quadrature points
   // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints) 
-  int nbFields = group.getNbFields();
-  fullMatrix<double> solutionQP (group.getNbElements() * nbFields, group.getNbIntegrationPoints());
+  int nbFields = claw.nbFields();
+  fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
   // ----- 1.2 --- multiply the solution by the collocation matrix
-  group.getSolution().mult ( group.getCollocationMatrix() , solutionQP); 
+  group.getCollocationMatrix().mult(solution  , solutionQP); 
   // ----- 1.3 --- if the conservation law is diffusive, compute the gradients too
   fullMatrix<double> *gradientSolutionQP= 0;
   if (claw.diffusiveFlux()){
-     gradientSolutionQP = new  fullMatrix<double> (group.getNbElements() * nbFields * 3, group.getNbIntegrationPoints());
-     group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP); 
+     gradientSolutionQP = new  fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
+     //group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP); 
   }
   
-
   // ----- 2 ----  compute all fluxes (diffusive and convective) at integration points
   // ----- 2.1 --- allocate elementary fluxes (compued in XYZ coordinates)
   fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
@@ -33,43 +35,62 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
   fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
 				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
 				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
-  // ----- 2.2 --- allocate parametric fluxes (compued in UVW coordinates) for all elements at all integration points
+  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())};
   // ----- 2.3 --- iterate on elements
-  for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {        
+  for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
     // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
-    fullMatrix<double> solutionQPe (solutionQP, iElement*claw.nbFields(),claw.nbFields() );
+    fullMatrix<double> solutionQPe (solutionQP, iElement*nbFields, nbFields );
     fullMatrix<double> *gradSolutionQPe;
-    if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );      
+    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());
-    // ----- 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
-      fullMatrix<double> fuvwe(Fuvw[iUVW], iElement*claw.nbFields(),claw.nbFields());
-      // ----- 2.3.3.1 --- compute \vec{f}^{UVW} =\vec{f}^{XYZ} J^{-1} and multiply bt detJ
-      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);
-	  // 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++) {
-	    fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
-	  }
-	}
-      } 
+    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
+        fullMatrix<double> fuvwe(Fuvw[iUVW], iElement*nbFields,nbFields);
+        // ----- 2.3.3.1 --- compute \vec{f}^{UVW} =\vec{f}^{XYZ} J^{-1} and multiply bt detJ
+        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);
+            // 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++) {
+              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); 
+      //we assume constant mass matrix and constant mapping for now
+      /*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;
   }
   // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
-  for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
-    group.getResidual().gemm(group.getFluxRedistributionMatrix(iUVW), Fuvw[iUVW]);
+  if(claw.convectiveFlux() || claw.diffusiveFlux()){
+    for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
+      residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
+  }
+  if(claw.sourceTerm()){
+    residual.gemm(group.getSourceRedistributionMatrix(),Source);
+  }
 }

Solver/dgAlgorithm.h

@@ -9,9 +9,11 @@ class dgConservationLaw;
 
 class dgAlgorithm {
  public :
-  void residualVolume ( /*dofManager &dof,*/
-			const dgConservationLaw &law,
-			const dgGroupOfElements & group);
+   void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
+				   const dgConservationLaw &claw,   // the conservation law
+				   const dgGroupOfElements & group, 
+           const fullMatrix<double> &solution,
+           fullMatrix<double> &residual);
   void residualInterface ( /*dofManager &dof,*/
 			   const dgConservationLaw &law,
 			   const dgGroupOfFaces & group);

Solver/dgConservationLaw.h

@@ -6,7 +6,7 @@
                     + \nabla \cdot (\vec{g}(u,\nabla u,forcings)  -> diffusive flux g
 		    + r(u,forcings)                       -> source term r
 */
-
+#include "fullMatrix.h"
 class dgElement; 
 class dgFace;
 
@@ -23,6 +23,7 @@ class dgFaceTerm{
 };
 
 class dgConservationLaw {
+  protected :
   int _nbf;
   dgTerm *_diffusive, *_convective, *_source, *_maxConvectiveSpeed;
   dgFaceTerm *_riemannSolver;
@@ -42,5 +43,6 @@ public:
 };
 
 
+dgConservationLaw *dgNewConservationLawAdvection();
 
 #endif

Solver/dgConservationLawAdvection.cpp

+#include "dgConservationLaw.h"
+#include "fullMatrix.h"
+#include "dgGroupOfElements.h"
+#include "SPoint3.h"
+#include "MElement.h"
+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);
+      //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);
+    }
+  }
+};
+
+class dgConservationLawAdvection : public dgConservationLaw {
+  public:
+  dgConservationLawAdvection() {
+    _nbf = 1;
+    _source = new testSourceTerm;
+  }
+};
+
+dgConservationLaw *dgNewConservationLawAdvection() {
+  return new dgConservationLawAdvection;
+}

Solver/dgGroupOfElements.cpp

@@ -44,7 +44,8 @@ static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
 dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
   : _elements(e), 
     _fs(*_elements[0]->getFunctionSpace(polyOrder)),
-    _integration(dgGetIntegrationRule (_elements[0], polyOrder))
+    _integration(dgGetIntegrationRule (_elements[0], polyOrder)
+    )
 {
   // this is the biggest piece of data ... the mappings
   _mapping = new fullMatrix<double> (_elements.size(), 10 * _integration->size1());
@@ -56,7 +57,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
 	e->getJacobian ((*_integration)(j,0),
 			(*_integration)(j,1),
 			(*_integration)(j,2),
-			ijac);
+			jac);
       detjac=inv3x3(jac,ijac);
       (*_mapping)(i,10*j + 0) = ijac[0][0]; 
       (*_mapping)(i,10*j + 1) = ijac[0][1]; 
@@ -77,10 +78,10 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
   _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> (_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),
@@ -90,13 +91,17 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
 	   (*_integration)(j,2), f);
     const double weight = (*_integration)(j,3);
     for (int k=0;k<_fs.coefficients.size1();k++){ 
-      (*_redistributionFluxes[0])(k,j) = g[j][0] * weight;
-      (*_redistributionFluxes[1])(k,j) = g[j][1] * weight;
-      (*_redistributionFluxes[2])(k,j) = g[j][2] * weight;
-      (*_redistributionSource)(k,j) = f[j] * weight;
-      (*_collocation)(k,j) = f[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];
+      for (int l=0;l<_fs.coefficients.size1();l++) { 
+        (*_imass)(k,l) += f[k]*f[l]*weight;
+      }
     }
   }
+  _imass->invertInPlace();
 }
 
 dgGroupOfElements::~dgGroupOfElements(){
@@ -107,6 +112,7 @@ dgGroupOfElements::~dgGroupOfElements(){
   delete _redistributionSource;
   delete _mapping;
   delete _collocation;
+  delete _imass;
 }
 
 // dgGroupOfFaces

Solver/dgGroupOfElements.h

@@ -19,18 +19,23 @@ class MEdge;
 class functionSpace;
 
 class dgElement {
-  const MElement *_element;
+  MElement *_element;
   // solution at points
   const fullMatrix<double> &_solution, &_integration, &_gradients;
 public:
-  dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
+  dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
     : _element(e), _solution(sol), _integration(integ), _gradients(sol)
   {}
-  dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
+  dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
     : _element(e), _solution(sol), _integration(integ), _gradients(grads)
   {}
+  const fullMatrix<double> &solution() const { return _solution; }
+  const fullMatrix<double> &integration() const { return _integration; }
+  MElement *element() const { return _element;}
 };
 
+
+// store topological and geometrical data for 1 group for 1 discretisation
 class dgGroupOfElements {
   // N elements in the group
   std::vector<MElement*> _elements;
@@ -50,15 +55,10 @@ class dgGroupOfElements {
   fullMatrix<double> *_redistributionFluxes[3];
   // redistribution for the source term
   fullMatrix<double> *_redistributionSource;
-  // the "finite element" gradient of the solution for this group (not owned)
-  fullMatrix<double> *_gradSolution;
-  // the solution for this group (not owned)
-  fullMatrix<double> *_solution;
-  // the residual for this group (not owned)
-  fullMatrix<double> *_residual;
+  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)
-  int _dimUVW, _dimXYZ, _nbFields;
+  int _dimUVW, _dimXYZ;
   // forbid the copy 
   //  dgGroupOfElements (const dgGroupOfElements &e, int order) {}
   //  dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
@@ -66,23 +66,16 @@ public:
   dgGroupOfElements (const std::vector<MElement*> &e, int pOrder);
   virtual ~dgGroupOfElements ();
   inline int getNbElements() const {return _elements.size();}
-  inline int getNbFields() const {return _nbFields;}
-  inline int getNbNodes() const {return _collocation->size1();}
-  inline int getNbIntegrationPoints() const {return _collocation->size2();}
+  inline int getNbNodes() const {return _collocation->size2();}
+  inline int getNbIntegrationPoints() const {return _collocation->size1();}
   inline int getDimUVW () const {return _dimUVW;}
   inline int getDimXYZ () const {return _dimXYZ;}
-  inline const MElement* getElement (int iElement) const {return _elements[iElement];}  
+  inline MElement* getElement (int iElement) const {return _elements[iElement];}  
   inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
   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 fullMatrix<double> & getSolution () const {return *_solution;}
-  inline fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
-  // get a proxy on the solution for element iElement
-  inline fullMatrix<double> getSolution (int iElement) const {return fullMatrix<double>(*_solution, iElement*_nbFields, _nbFields);}
-  inline fullMatrix<double> & getResidual () const {return *_solution;}
-  // get a proxy on the residual for element iElement
-  inline fullMatrix<double> getResidual (int iElement) const {return fullMatrix<double>(*_residual, iElement*_nbFields, _nbFields);}
+  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);}
@@ -141,4 +134,5 @@ public:
 };
 
 
+
 #endif

Solver/dgMainTest.cc

+#include <stdio.h>
+#include <vector>
+#include "GModel.h"
+#include "dgGroupOfElements.h"
+#include "dgAlgorithm.h"
+#include "dgConservationLaw.h"
+#include "Gmsh.h"
+
+
+#include "MElement.h"
+void print (const char *filename, std::vector<MElement *> &els, double *v);
+std::vector<MElement *> get_all_tri(GModel *model);
+
+
+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());
+  dgAlgorithm algo;
+  dgConservationLaw *law = dgNewConservationLawAdvection();
+  algo.residualVolume(*law,group,sol,residu);
+  sol.gemm(group.getInverseMassMatrix(),residu);
+  print("test.pos",all_tri,&sol(0,0));
+}
+
+
+
+std::vector<MElement *> get_all_tri(GModel *model){
+  std::vector<GEntity*> entities;
+  model->getEntities(entities);
+  std::vector<MElement *> all_tri;
+  for(std::vector<GEntity *>::iterator itent=entities.begin(); itent!=entities.end(); itent++){
+    if ((*itent)->dim()!=2) continue;
+    for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
+      all_tri.push_back((*itent)->getMeshElement(iel));
+  }
+  return all_tri;
+}
+
+void print (const char *filename, std::vector<MElement *> &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;
+    fprintf(file,"ST (");
+    for (int iv=0; iv<el->getNumVertices(); iv++) {
+      MVertex *vertex = el->getVertex(iv);
+      SPoint3 coord = vertex->point();
+      fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
+    }
+    fprintf(file,"{");
+    for (int iv=0; iv<el->getNumVertices(); iv++)
+      fprintf(file,"%e%c ",v[i++],iv==el->getNumVertices()-1?'}':',');
+    fprintf(file,";\n");
+  }
+  fprintf(file,"};");
+  fclose(file);
+}