Inverse mass matrix directly multiply the sparse projection matrix

Solver/TESTCASES/Projection2D3D.lua

 model = GModel()
 -- model:load('tetrahedra.geo')
 -- model:load('tetrahedra.msh')
--- model:load('prisms.geo')
--- model:load('prisms.msh')
-model:load('mixed.geo')
-model:load('mixed.msh')
+model:load('prisms.geo')
+model:load('prisms.msh')
+-- model:load('mixed.geo')
+-- model:load('mixed.msh')
 order=1
 
 -- initial condition
@@ -13,11 +13,13 @@ function initial_condition( xyz , f )
       x = xyz:get(i,0)
       y = xyz:get(i,1)
       z = xyz:get(i,2)
-      f:set (i, 0, math.exp(-50*((x-0.4)^2+(y-0.3)^2)))
+      f:set (i, 0, math.exp(-50*((x-0.5)^2+(y-0.5)^2)))
 --       f:set (i, 0, math.exp(-50*((x-0.4)^2+(z+0.3)^2)))
    end
 end
 
+proj = linearSystemCSRGmmdouble()
+
 init_cond = functionLua(1, 'initial_condition', {'XYZ'}):getName()
 gc=dgGroupCollection(model,3,order)
 solTmp=dgDofContainer(gc,1)
@@ -25,10 +27,13 @@ solTmp:L2Projection(init_cond)
 gc:buildGroupsOfInterfaces(model,2,order)
 solution=dgDofContainer(gc,1)
 solution:L2Projection(init_cond)
-solution:exportMsh('sol3d.msh');
+solution:exportMsh('sol3d.msh')
+
+gc2D = gc:newByTag(model,2,order,{"top"})
+solution2d=dgDofContainer(gc2D,1)
+--solution2d:L2Projection(init_cond)
 
-gc2D = gc:newByTag(model,2,order,{"top","side1"});
-solution2d=dgDofContainer(gc2D,1);
-solution2d:L2Projection(init_cond);
-solution2d:exportMsh('sol2d.msh');
+solution2d:Mesh2Mesh_BuildL2Projection(proj,solution)
+solution2d:Mesh2Mesh_ApplyL2Projection(proj,solution)
+solution2d:exportMsh('sol2d.msh')
 

Solver/TESTCASES/prisms.geo

 //*********** prisms.geo *************//
 C = 1;
-Lup = 1;
 L = 1.;
-lc = .1;
+Lup = L;
+lc = 0.1;
 
 Point(1) = {0.0, 0.0, -Lup, lc};
 Point(2) = {C  , 0.0, -Lup, lc};
@@ -18,7 +18,7 @@ Line(4) = {4,1};
 Line Loop(5) = {1,2,3,4};
 Plane Surface(6) = {5};
 
-outpri[]= Extrude {0,0,1.2*L}{ Surface{6}; Layers{Ceil(1.2*L/lc)};Recombine;};
+outpri[]= Extrude {0,0,L}{ Surface{6}; Layers{Ceil(L/lc)};Recombine;};
 // outv[]= Extrude {0,0,0.5*L}{ Surface{7}; Layers{5};Recombine;};
 // Printf("top surface = %g", outpri[0]);
 // Printf("volume = %g", outpri[1]);

Solver/dgDofContainer.cpp

@@ -304,9 +304,39 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
       jGlobal += dGroup.getNbNodes();
     }
   }
+  multAddInverseMassMatrixL2Projection(projector);
 //   buffer.print();
 }
 
+void dgDofContainer::multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector){
+  dgGroupCollection* rGroups = this->getGroups();
+  int *startIndex;
+  int *columns;
+  double *values;
+  projector.getMatrix(startIndex,columns,values);
+
+  int iGlobal = 0;
+  for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {// for 2d.groups
+    const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
+    for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {// for elements
+      fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iGlobal+1]-startIndex[iGlobal]);
+      fullMatrix<double> buffer2 = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iGlobal+1]-startIndex[iGlobal]);
+      fullMatrix<double> iMassEl;
+      for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+        for (int i = startIndex[iGlobal+iNode]; i < startIndex[iGlobal+iNode+1] ; i++)
+          buffer(iNode,i-startIndex[iGlobal+iNode]) = values[i];
+      }
+      iMassEl.setAsProxy(rGroup.getInverseMassMatrix(),iElement*rGroup.getNbNodes(),rGroup.getNbNodes());
+      buffer2.gemm(iMassEl,buffer);
+      for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+        for (int i = startIndex[iGlobal+iNode]; i < startIndex[iGlobal+iNode+1] ; i++)
+          values[i] = buffer2(iNode,i-startIndex[iGlobal+iNode]);
+      }
+      iGlobal+=rGroup.getNbNodes();
+    }
+  }
+}
+
 void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor){
   scale(0.);
 
@@ -327,7 +357,7 @@ int nbFields = 1; // TO DEFINE !!!
 
   for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {// for 2d.groups
     const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
-    fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),rGroup.getNbElements() * _nbFields);
+//    fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),rGroup.getNbElements() * _nbFields);
     this->getGroupProxy(iGroup).scale(0);
     for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {// for elements
       for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
@@ -340,14 +370,14 @@ int nbFields = 1; // TO DEFINE !!!
           int jElement = (jGlobal-dGroupsStartIGlobal[jGroup])/dGroups->getElementGroup(jGroup)->getNbNodes();
           int jNode = jGlobal-dGroupsStartIGlobal[jGroup]-jElement*dGroups->getElementGroup(jGroup)->getNbNodes();
           for (int m = 0; m < nbFields; m++){
-            buffer(iNode,iElement*nbFields+m) += values[i]*(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m);
+            this->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += values[i]*(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m);
 //             printf("%d %d %d %f %f %f\n",iGlobal,jGlobal,i, buffer(iNode,iElement*nbFields+m),values[i],(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m));
           }
         }
       }
     }//*/
 //     buffer.print();
-    dgAlgorithm::multAddInverseMassMatrix(rGroup, buffer, this->getGroupProxy(iGroup));
+//    dgAlgorithm::multAddInverseMassMatrix(rGroup, buffer, this->getGroupProxy(iGroup));
 //     this->getGroupProxy(iGroup).print();
   }
 }

Solver/dgDofContainer.h

@@ -44,6 +44,7 @@ public:
   void setAll(double v);
   void L2Projection(std::string functionName);
   void Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
+  void multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector); // this method should be private
   void Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
   void exportMsh(const std::string name);
   void exportGroupIdMsh();