towards initial dg code ...

Numeric/functionSpace.h

@@ -23,10 +23,10 @@ struct functionSpace
   fullMatrix<double> coefficients;
   // for a given face/edge, with both a sign and a rotation,
   // give an ordered list of nodes on this face/edge
-  std::vector<int> & getFaceClosure (int iFace, int iSign, int iRot){
+  inline const std::vector<int> & getFaceClosure (int iFace, int iSign, int iRot)const{
     return faceClosure[iFace+4*(iSign==1?0:1)+8*iRot];
   }
-  inline std::vector<int> & getEdgeClosure (int iEdge, int iSign){
+  inline const std::vector<int> & getEdgeClosure (int iEdge, int iSign) const{
     return edgeClosure[iSign == 1 ? iEdge : 3+iEdge];
   }
   inline void evaluateMonomials(double u, double v, double w, double p[]) const 

Solver/dgAlgorithm.cpp

@@ -18,7 +18,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
   fullMatrix<double> *gradientSolutionQP= 0;
   if (claw.diffusiveFlux()){
      gradientSolutionQP = new  fullMatrix<double> (group.getNbElements() * group.getNbFields() * 3, group.getNbIntegrationPoints());
-     group.getGradientOfSolution().mult ( group.getCollocationMatrix() , solutionQP); 
+     group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP); 
   }
   
 
@@ -40,7 +40,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
     fullMatrix<double> solutionQPe (solutionQP, iElement*claw.nbFields(),claw.nbFields() );
     fullMatrix<double> *gradSolutionQPe;
     if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );      
-    else gradSolutionQPe = new fullMatrix<double>
+    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);

Solver/dgGroupOfElements.cpp

@@ -199,6 +199,21 @@ void dgGroupOfFaces::createEdgeElements (const std::vector<MEdge> &topo_edges){
     }
   }
 }
+void dgGroupOfFace::init() {
+  _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());
+  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<_fs.coefficients.size1();k++){ 
+      (*_redistribution)(k,j) = f[j] * weight;
+      (*_collocation)(k,j) = f[k];
+    }
+  }
+}
 
 dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces, 		  
 				const std::vector<MElement*> &l, 
@@ -207,8 +222,7 @@ dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces,
         _fsLeft(_left[0]->getFunctionSpace (pOrder)), _fsRight (_right[0]->getFunctionSpace (pOrder))
 {
   createFaceElements (topo_faces);
-  _fsFace = _faces[0]->getFunctionSpace (pOrder);
-  dgGetIntegrationRule (_faces[0],pOrder);
+  init();
 }
 
 dgGroupOfFaces::dgGroupOfFaces (const std::vector<MEdge> &topo_edges, 		  
@@ -218,6 +232,5 @@ dgGroupOfFaces::dgGroupOfFaces (const std::vector<MEdge> &topo_edges,
         _fsLeft(_left[0]->getFunctionSpace (pOrder)), _fsRight (_right[0]->getFunctionSpace (pOrder))
 {
   createEdgeElements (topo_edges);
-  _fsFace = _faces[0]->getFunctionSpace (pOrder);
-  dgGetIntegrationRule (_faces[0],pOrder);
+  init();
 }

Solver/dgGroupOfElements.h

@@ -121,10 +121,12 @@ class dgGroupOfFaces {
   // detJac at integration points (N*Ni) x 1
   fullMatrix<double> *_detJac;
   // collocation matrices \psi_i (GP_j) 
-  fullMatrix<double> *_collocationLeft, *_collocationRight;
+  fullMatrix<double> *_collocation;
+  //fullMatrix<double> *_collocationLeft, *_collocationRight;
   // redistribution matrices \psi_i (GP_j) * weight_j
-  fullMatrix<double> *_redistributionLeft, *_redistributionRight;
-  
+  fullMatrix<double> *_redistribution;
+  //common part of the 3 constructors
+  void init();
 public:
   dgGroupOfFaces (const std::vector<MFace> &faces, 		  
 		  const std::vector<MElement*> &l,