first version of dgMesh2MeshProjection class

Common/LuaBindings.cpp

@@ -31,6 +31,7 @@
 #include "drawContext.h"
 #include "GmshMessage.h"
 #include "linearSystemCSR.h"
+#include "dgMesh2MeshProjection.h"
 
 extern "C" {
   #include "lua.h"
@@ -369,6 +370,7 @@ binding::binding(){
   gmshOptions::registerBindings(this);
   Msg::registerBindings(this);
   linearSystemCSRGmm<double>::registerBindings(this);
+  dgMesh2MeshProjection::registerBindings(this);
 }
 binding *binding::_instance=NULL;
 #endif

Solver/CMakeLists.txt

@@ -32,6 +32,7 @@ set(SRC
   luaFunction.cpp
   functionSpace.cpp
   dgFunctionIntegrator.cpp
+  dgMesh2MeshProjection.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/TESTCASES/Projection2D3D.lua

@@ -5,7 +5,7 @@ model:load('prisms.geo')
 model:load('prisms.msh')
 -- model:load('mixed.geo')
 -- model:load('mixed.msh')
-order=1
+order=2
 
 -- initial condition
 function initial_condition( xyz , f )
@@ -13,27 +13,42 @@ 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.5)^2+(y-0.5)^2)))
---       f:set (i, 0, math.exp(-50*((x-0.4)^2+(z+0.3)^2)))
+--       f:set (i, 0, math.exp(x+3*y))
+      f:set (i, 0, math.exp(-20*((x-0.4)^2+(y-0.3)^2)))
    end
 end
 
-proj = linearSystemCSRGmmdouble()
-
-init_cond = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+init_cond = functionLua(1, 'initial_condition', {functionCoordinates.get()})
 gc=dgGroupCollection(model,3,order)
-solTmp=dgDofContainer(gc,1)
-solTmp:L2Projection(init_cond)
 gc:buildGroupsOfInterfaces(model,2,order)
 solution=dgDofContainer(gc,1)
 solution:L2Projection(init_cond)
+-- solution:setAll(3.0)
 solution:exportMsh('sol3d.msh')
 
 gc2D = gc:newByTag(model,2,order,{"top"})
 solution2d=dgDofContainer(gc2D,1)
---solution2d:L2Projection(init_cond)
-
-solution2d:Mesh2Mesh_BuildL2Projection(proj,solution)
-solution2d:Mesh2Mesh_ApplyL2Projection(proj,solution)
+-- solution2d:L2Projection(init_cond)
+solution2d:setAll(2.0)
 solution2d:exportMsh('sol2d.msh')
 
+print("** Build projection matrix **")
+pm = dgMesh2MeshProjection(solution,solution2d)
+print("** 3D->2D projection **")
+pm:projectFromTo(solution,solution2d)
+solution2d:exportMsh('proj3d2d.msh')
+print("** 3D->2D copy **")
+pm:copyFromTo(solution,solution2d)
+solution2d:exportMsh('copy3d2d.msh')
+
+print("** 2D->3D projection **")
+solution2d:setAll(2.0)
+pm:projectFromTo(solution2d,solution)
+solution:exportMsh('proj2d3d.msh')
+print("** 2D->3D copy **")
+pm:copyFromTo(solution2d,solution)
+solution:exportMsh('copy2d3d.msh')
+
+-- raises error
+-- pm:copyFromTo(solution,solution)
+-- pm:copyFromTo(solution2d,solution2d)

Solver/dgDofContainer.cpp

@@ -32,12 +32,14 @@ dgDofContainer::dgDofContainer (dgGroupCollection *groups, int nbFields):
   // create proxys for each group
   int offset = 0;
   _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
+  _groupFirstDofId.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
   for (int i=0;i<_groups.getNbElementGroups();i++){    
     dgGroupOfElements *group = _groups.getElementGroup(i);
     _groupId[group] = i;
     int nbNodes    = group->getNbNodes();
     int nbElements = group->getNbElements();
     _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, _nbFields*nbElements);
+    _groupFirstDofId[i] = offset;
     offset += nbNodes*_nbFields*nbElements;
   }  
 
@@ -51,16 +53,18 @@ dgDofContainer::dgDofContainer (dgGroupCollection *groups, int nbFields):
     _dataSizeGhost += nbNodes*_nbFields*nbElements;
   }
 
-  _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
   _ghostData = new fullVector<double>(_dataSizeGhost);
-  offset=0;
+  int ghostOffset=0;
   for (int i=0;i<_groups.getNbGhostGroups();i++){    
     dgGroupOfElements *group = _groups.getGhostGroup(i);
     int nbNodes    = group->getNbNodes();
     int nbElements = group->getNbElements();
-    _groupId[group] = i+_groups.getNbElementGroups();
-    _dataProxys[i+_groups.getNbElementGroups()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, _nbFields*nbElements);
+    int gid = i+_groups.getNbElementGroups();
+    _groupId[group] = gid;
+    _dataProxys[gid] = new fullMatrix<double> (&(*_ghostData)(ghostOffset),nbNodes, _nbFields*nbElements);
+    _groupFirstDofId[gid] = offset;
     offset += nbNodes*_nbFields*nbElements;
+    ghostOffset += nbNodes*_nbFields*nbElements;
   }  
 
 }
@@ -251,9 +255,10 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
 // evauate quad
 // add to matrix correct indices
 
-// For multigrid, the donor and receiver are inverted :S
   dgGroupCollection* dGroups = donor.getGroups();
   dgGroupCollection* rGroups = this->getGroups();
+  if (rGroups->getElementGroup(0)->getDimUVW() > dGroups->getElementGroup(0)->getDimUVW())
+    Msg::Fatal("Cannot build projection matrix from lower dim to higher dim. Use the transpose.");
   std::vector<int> rGroupsStartIGlobal(rGroups->getNbElementGroups() + 1);
   int jGlobal = 0; // indices in global container
   rGroupsStartIGlobal[0] = 0;
@@ -261,9 +266,9 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
     rGroupsStartIGlobal[i] = rGroupsStartIGlobal[i-1] + rGroups->getElementGroup(i-1)->getNbElements()*rGroups->getElementGroup(i-1)->getNbNodes();
   
   projector.allocate(rGroupsStartIGlobal[rGroupsStartIGlobal.size()-1]);
-//   fullMatrix<double> buffer = fullMatrix<double> (3,6);
-  for (int jGroup=0;jGroup<dGroups->getNbElementGroups();jGroup++) {// for 3d.groups
-    const dgGroupOfElements &dGroup = *dGroups->getElementGroup(jGroup);
+  // integration over donor function space
+  for (int jGroup=0;jGroup<dGroups->getNbElementGroups();jGroup++) {// for donor groups
+    const dgGroupOfElements &dGroup = *dGroups->getElementGroup(jGroup); 
     fullMatrix<double> iPtsMatrix = dGroup.getIntegrationPointsMatrix();
     for (int jElement=0 ; jElement<dGroup.getNbElements() ;++jElement) {// for elements
       double rShapeFun[256];
@@ -271,7 +276,7 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
       GModel *rModel = rGroups->getModel();
       GModel *dModel = dGroups->getModel();
       MElement* dElem = dGroup.getElement(jElement);
-      for (int iPt =0; iPt< dGroup.getNbIntegrationPoints(); iPt++) {
+      for (int iPt =0; iPt< iPtsMatrix.size1(); iPt++) {
         double x=0,y=0,z=0;
         dElem->getFunctionSpace()->f(iPtsMatrix(iPt,0),iPtsMatrix(iPt,1),iPtsMatrix(iPt,2),dShapeFun);      
         for (int iVer=0; iVer < dElem->getNumVertices(); iVer++) {
@@ -279,25 +284,27 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
           y += dElem->getVertex(iVer)->y()*dShapeFun[iVer];
           z += dElem->getVertex(iVer)->z()*dShapeFun[iVer];
         }
-        z = 0; // dummy projection to 2d mesh
-        SPoint3 p(x,y,z); // find p in 2D mesh
+        if (rGroups->getElementGroup(0)->getDimUVW() == 2) {
+          z = 0; // dummy projection to 2d mesh
+        }
+        SPoint3 p(x,y,z); // find p in receiver mesh
         MElement *rElem = rGroups->getModel()->getMeshElementByCoord(p);
         int iGroup,iElement;
         rGroups->find(rElem,iGroup,iElement);
         if (iElement == -1)
-          Msg::Error("Integration point not found in receiver mesh.");
+          Msg::Fatal("Integration point not found in receiver mesh.");
         const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);      
         double U[3],X[3]={x,y,z};
         rElem->xyz2uvw(X,U);
         dGroup.getFunctionSpace().f(iPtsMatrix(iPt,0),iPtsMatrix(iPt,1),iPtsMatrix(iPt,2),dShapeFun);
         rGroup.getFunctionSpace().f(U[0],U[1],U[2],rShapeFun);
         const double detJ = dGroup.getDetJ (jElement, iPt);
-        int iGlobal = rGroupsStartIGlobal[iGroup]+rGroup.getNbNodes()*iElement;
+//         int iGlobal = rGroupsStartIGlobal[iGroup]+rGroup.getNbNodes()*iElement;
+        int iGlobal = this->getDofId(iGroup,iElement,0,0); // works only for one field !
         for (int jNode=0;jNode<dGroup.getNbNodes();jNode++){
           for (int iNode=0;iNode<rGroup.getNbNodes();iNode++){
             double val = rShapeFun[iNode]*dShapeFun[jNode]*iPtsMatrix(iPt,3)*detJ;
             projector.addToMatrix(iGlobal+iNode,jGlobal+jNode,val);
-//             buffer(iGlobal+iNode,jGlobal+jNode) += val;
           }
         }
       }
@@ -305,7 +312,6 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
     }
   }
   multAddInverseMassMatrixL2Projection(projector);
-//   buffer.print();
 }
 
 void dgDofContainer::multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector){
@@ -337,11 +343,17 @@ void dgDofContainer::multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<dou
   }
 }
 
-void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor){
-  scale(0.);
+void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor, int transpose, int copy){
+  dgDofContainer* dDof = &donor;
+  dgDofContainer* rDof = this;
+  if (transpose) {
+    rDof = &donor;
+    dDof = this;
+  }
+  setAll(0);
 
-  dgGroupCollection* dGroups = donor.getGroups();
-  dgGroupCollection* rGroups = this->getGroups();
+  dgGroupCollection* dGroups = dDof->getGroups();
+  dgGroupCollection* rGroups = rDof->getGroups();
   std::vector<int> dGroupsStartIGlobal(dGroups->getNbElementGroups() + 1);
   int iGlobal = 0; // indices in global container
   dGroupsStartIGlobal[0] = 0;
@@ -353,12 +365,11 @@ void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &pro
   double *values;
   projector.getMatrix(startIndex,columns,values);
 
-int nbFields = 1; // TO DEFINE !!!
+  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);
-    this->getGroupProxy(iGroup).scale(0);
     for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {// for elements
       for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
         int jGroup = 0;
@@ -370,7 +381,12 @@ 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++){
-            this->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += values[i]*(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m);
+            double val = values[i];
+            if (copy) val = (fabs(values[i]) > 1e-8) ? 1 : 0;
+            if (transpose)
+              this->getGroupProxy(jGroup)(jNode,jElement*nbFields+m) += val*donor.getGroupProxy(iGroup)(iNode,iElement*nbFields+m);
+            else
+              this->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += val*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));
           }
         }
@@ -529,7 +545,7 @@ void dgDofContainer::exportMsh(const std::string name)
       name_oss<<"_"<<Msg::GetCommRank();
     FILE *f = fopen (name_oss.str().c_str(),"w");
     if(!f){
-      Msg::Error("Unable to open export file '%s'", name.c_str());
+      Msg::Fatal("Unable to open export file '%s'", name.c_str());
     }
 
     int COUNT = 0;
@@ -592,7 +608,6 @@ void dgDofContainer::exportMsh(const std::string name)
   */
 }
 
-
 #include "LuaBindings.h"
 void dgDofContainer::registerBindings(binding *b){
   classBinding *cb = b->addClass<dgDofContainer>("dgDofContainer");
@@ -627,5 +642,5 @@ void dgDofContainer::registerBindings(binding *b){
   cm->setArgNames("projector","donor",NULL);
   cm = cb->addMethod("Mesh2Mesh_ApplyL2Projection",&dgDofContainer::Mesh2Mesh_ApplyL2Projection);
   cm->setDescription("Apply L2 projection matrix from donor to this dofContainer.");
-  cm->setArgNames("projector","donor",NULL);
+  cm->setArgNames("projector","donor","transpose","copy",NULL);
 }

Solver/dgDofContainer.h

@@ -25,15 +25,25 @@ private:
   double *sendBuf, *recvBuf;
   std::vector<fullMatrix<double> *> _dataProxys; // proxys 
   std::map<const dgGroupOfElements*,int> _groupId;
+  std::vector<int> _groupFirstDofId;
   int _mshStep;
 public:
   void scale(double f);
+  inline int getDofId (int groupId, int elementId, int fieldId, int nodeId) const {
+    const fullMatrix<double> &proxy = getGroupProxy(groupId);
+    return _groupFirstDofId[groupId]+(elementId*_nbFields+fieldId)*proxy.size1()+nodeId;
+  }
+  inline int getGroupFirstDofId(int groupId) {
+    return _groupFirstDofId[groupId];
+  }
+  inline fullVector<double> &getVector() {return *_data;}
   void scale(std::vector<dgGroupOfElements*>groups, double f);
   double norm();
   double norm(std::vector<dgGroupOfElements*>groups);
   void axpy(dgDofContainer &x, double a=1.);
   void axpy(std::vector<dgGroupOfElements*>groups,dgDofContainer &x, double a=1.);
   inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
+  inline const fullMatrix<double> &getGroupProxy(int gId) const { return *(_dataProxys[gId]); }
   inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); }
   dgDofContainer (dgGroupCollection *groups, int nbFields);
   ~dgDofContainer ();  
@@ -46,7 +56,7 @@ public:
   void L2Projection(const function *f);
   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 Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor, int transpose, int copy);
   void exportMsh(const std::string name);
   void exportGroupIdMsh();
   void exportMultirateGroupMsh();

Solver/dgMesh2MeshProjection.cpp

+#include "dgMesh2MeshProjection.h"
+#include "MElementOctree.h"
+#include "Octree.h"
+
+dgMesh2MeshProjection::dgMesh2MeshProjection()
+{
+  _projectionMatrixIsBuilt = false;
+  _projectionMatrix = new linearSystemCSRGmm<double>;
+  _numInDofs = _numOutDofs = _inMaxDim = _outMaxDim = 0;
+  _inDofPattern = _outDofPattern = NULL;
+  _useDofContainerOctree = true;
+  _octree = NULL;
+}
+
+dgMesh2MeshProjection::dgMesh2MeshProjection(dgDofContainer* donor, dgDofContainer* receiver)
+{
+  _projectionMatrixIsBuilt = false;
+  _projectionMatrix = new linearSystemCSRGmm<double>;
+  _numInDofs = _numOutDofs = 0;
+  _inDofPattern = _outDofPattern = NULL;
+  _useDofContainerOctree = true;
+  _octree = NULL;
+  _buildProjectionMatrix(donor,receiver);
+}
+
+dgMesh2MeshProjection::dofStoragePattern::dofStoragePattern(dgDofContainer* dofContainer)
+{
+  _nbGroups = dofContainer->getGroups()->getNbElementGroups();
+  _nbFields = dofContainer->getNbFields();
+  _nbElems.resize(_nbGroups);
+  _nbNodes.resize(_nbGroups);
+  _groupFirstProjectionDofId.resize(_nbGroups);
+
+  int offset = 0;
+  for (int i = 0; i < _nbGroups; i++) {
+    int nbElems = dofContainer->getGroups()->getElementGroup(i)->getNbElements();
+    int nbNodes = dofContainer->getGroups()->getElementGroup(i)->getNbNodes();
+    _nbElems[i] = nbElems;
+    _nbNodes[i] = nbNodes;
+    _groupFirstProjectionDofId[i] = offset;
+    offset += nbElems*nbNodes;
+  }
+  _nbProjDofs = offset;
+}
+
+bool dgMesh2MeshProjection::dofStoragePattern::operator==(const dgMesh2MeshProjection::dofStoragePattern& other) const
+{
+  if(this->_nbGroups != other._nbGroups)
+    return false;
+  for (int i =0; i<this->_nbGroups; i++) {
+    if(this->_nbElems[i] != other._nbElems[i])
+      return false;
+    if(this->_nbNodes[i] != other._nbNodes[i])
+      return false;
+  }
+  return true;
+}
+
+bool dgMesh2MeshProjection::_checkTranspose(dgDofContainer* donor, dgDofContainer* receiver)
+{
+  // check if transpose is necessary
+  dofStoragePattern* dDofPattern = new dofStoragePattern(donor);
+  dofStoragePattern* rDofPattern = new dofStoragePattern(receiver);
+  bool transpose;
+  if (*dDofPattern == *_inDofPattern && *rDofPattern == *_outDofPattern)
+    transpose = false;
+  else if (*dDofPattern == *_outDofPattern && *rDofPattern == *_inDofPattern)
+    transpose = true;
+  else
+    Msg::Fatal("the DOF storage patterns do not match the projection matrix.");
+  return transpose;
+}
+
+void dgMesh2MeshProjection::projectFromTo(dgDofContainer* donor, dgDofContainer* receiver)
+{
+  if(!_projectionMatrixIsBuilt)
+    _buildProjectionMatrix(donor,receiver);
+  _applyProjection(donor,receiver,_checkTranspose(donor,receiver),0);
+}
+
+void dgMesh2MeshProjection::copyFromTo(dgDofContainer* donor, dgDofContainer* receiver)
+{
+  if(!_projectionMatrixIsBuilt)
+    _buildProjectionMatrix(donor,receiver);
+  _applyProjection(donor,receiver,_checkTranspose(donor,receiver),1);
+}
+
+
+void dgMesh2MeshProjection::_buildReceiverOctree(dgDofContainer* receiver)
+{
+  std::vector<MElement*> allElements;
+  for (int i = 0; i < receiver->getGroups()->getNbElementGroups(); i++ ) {
+    dgGroupOfElements* group = receiver->getGroups()->getElementGroup(i);
+    for (int iElem = 0; iElem < group->getNbElements(); iElem++ ) {
+      allElements.push_back(group->getElement(iElem));
+    }
+  }
+  _octree = buildMElementOctree(allElements);
+}
+
+MElement* dgMesh2MeshProjection::_getReceiverElementByCoord(SPoint3 point)
+{
+  if (!_octree)
+    Msg::Fatal("MElement Octree is not built");
+  double P[3] = {point.x(), point.y(), point.z()};
+  return (MElement*)Octree_Search(P, _octree);
+}
+
+void dgMesh2MeshProjection::_buildProjectionMatrix(dgDofContainer* donor, dgDofContainer* receiver)
+{
+  _inDofPattern = new dofStoragePattern(donor);
+  _outDofPattern = new dofStoragePattern(receiver);
+  _inMaxDim = donor->getGroups()->getElementGroup(0)->getDimUVW(); // not necessary max...
+  _outMaxDim = receiver->getGroups()->getElementGroup(0)->getDimUVW();
+  dgGroupCollection* dGroups = donor->getGroups();
+  dgGroupCollection* rGroups = receiver->getGroups();
+  if (_outMaxDim == 2, receiver->getGroups()->getModel()->getDim() == 3) {
+    // need to find 2D elements in 3D geometry, build special octree
+    _useDofContainerOctree = false;
+    _buildReceiverOctree(receiver);
+  }
+  if (_outMaxDim > _inMaxDim)
+    Msg::Fatal("Cannot build projection matrix from lower dim to higher dim. Swap donor and receiver.");
+  
+  int jProj = 0; // indices in projection matrix
+  int totalNbOutDofs = 0;
+  for (int i = 0; i<rGroups->getNbElementGroups(); i++)
+    totalNbOutDofs += rGroups->getElementGroup(i)->getNbElements()*rGroups->getElementGroup(i)->getNbNodes();
+  _projectionMatrix->allocate(totalNbOutDofs);
+  
+  // integration over donor function space
+  for (int jGroup=0;jGroup<dGroups->getNbElementGroups();jGroup++) {// for donor groups
+    const dgGroupOfElements &dGroup = *dGroups->getElementGroup(jGroup); 
+    fullMatrix<double> iPtsMatrix = dGroup.getIntegrationPointsMatrix();
+    for (int jElement=0 ; jElement<dGroup.getNbElements() ;++jElement) {// for elements
+      double rShapeFun[256];
+      double dShapeFun[256];
+      GModel *rModel = rGroups->getModel();
+      GModel *dModel = dGroups->getModel();
+      MElement* dElem = dGroup.getElement(jElement);
+      for (int iPt =0; iPt< iPtsMatrix.size1(); iPt++) {
+        double x=0,y=0,z=0;
+        dElem->getFunctionSpace()->f(iPtsMatrix(iPt,0),iPtsMatrix(iPt,1),iPtsMatrix(iPt,2),dShapeFun);      
+        for (int iVer=0; iVer < dElem->getNumVertices(); iVer++) {
+          x += dElem->getVertex(iVer)->x()*dShapeFun[iVer];
+          y += dElem->getVertex(iVer)->y()*dShapeFun[iVer];
+          z += dElem->getVertex(iVer)->z()*dShapeFun[iVer];
+        }
+        if (_inMaxDim==3 && _outMaxDim == 2) {
+          z = 0; // dummy projection to 2d mesh
+        }
+        // find p in receiver mesh
+        SPoint3 p(x,y,z);
+        MElement *rElem;
+        if (_useDofContainerOctree)
+          rElem = rGroups->getModel()->getMeshElementByCoord(p);
+        else
+          rElem = _getReceiverElementByCoord(p);
+          
+        int iGroup,iElement;
+        rGroups->find(rElem,iGroup,iElement);
+        if (iElement == -1)
+          Msg::Fatal("Integration point (%g,%g,%g) not found in receiver mesh",p.x(),p.y(),p.z()); 
+        const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);      
+        double U[3],X[3]={x,y,z};
+        rElem->xyz2uvw(X,U);
+        dGroup.getFunctionSpace().f(iPtsMatrix(iPt,0),iPtsMatrix(iPt,1),iPtsMatrix(iPt,2),dShapeFun);
+        rGroup.getFunctionSpace().f(U[0],U[1],U[2],rShapeFun);
+        const double detJ = dGroup.getDetJ (jElement, iPt);
+        int iProj = _outDofPattern->getProjDofId(iGroup,iElement,0);
+        for (int jNode=0;jNode<dGroup.getNbNodes();jNode++){
+          for (int iNode=0;iNode<rGroup.getNbNodes();iNode++){
+            double val = rShapeFun[iNode]*dShapeFun[jNode]*iPtsMatrix(iPt,3)*detJ;
+            _projectionMatrix->addToMatrix(iProj+iNode,jProj+jNode,val);
+          }
+        }
+      }
+      jProj += dGroup.getNbNodes();
+    }
+  }
+  _multiplyByInverseMassMatrix(receiver);
+  _projectionMatrixIsBuilt = true;
+}
+
+// could be improved by using a better matrix representation and an efficient M*P routine
+void dgMesh2MeshProjection::_multiplyByInverseMassMatrix(dgDofContainer* donor)
+{
+  dgGroupCollection* rGroups = donor->getGroups();
+  int *startIndex;
+  int *columns;
+  double *values;
+  _projectionMatrix->getMatrix(startIndex,columns,values);
+
+  int iProj = 0;
+  for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {
+    const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
+    for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {
+      fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iProj+1]-startIndex[iProj]);
+      fullMatrix<double> buffer2 = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iProj+1]-startIndex[iProj]);
+      fullMatrix<double> iMassEl;
+      for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+        for (int i = startIndex[iProj+iNode]; i < startIndex[iProj+iNode+1] ; i++)
+          buffer(iNode,i-startIndex[iProj+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[iProj+iNode]; i < startIndex[iProj+iNode+1] ; i++)
+          values[i] = buffer2(iNode,i-startIndex[iProj+iNode]);
+      }
+      iProj+=rGroup.getNbNodes();
+    }
+  }
+}
+
+// could be improved by using a better matrix representation and an efficient M*V routine
+void dgMesh2MeshProjection::_applyProjection(dgDofContainer* donor, dgDofContainer* receiver, bool transpose, bool copy)
+{
+  dgDofContainer* dDof = donor;
+  dgDofContainer* rDof = receiver;
+  if (transpose) {
+    rDof = donor;
+    dDof = receiver;
+  }
+  receiver->setAll(0);
+
+  dgGroupCollection* dGroups = dDof->getGroups();
+  dgGroupCollection* rGroups = rDof->getGroups();
+  std::vector<int> dGroupsStartIGlobal(dGroups->getNbElementGroups() + 1);
+  int iProj = 0; // indices in projection matrix
+  dGroupsStartIGlobal[0] = 0;
+  for (int i = 1; i<dGroupsStartIGlobal.size(); i++)
+    dGroupsStartIGlobal[i] = dGroupsStartIGlobal[i-1] + dGroups->getElementGroup(i-1)->getNbElements()*dGroups->getElementGroup(i-1)->getNbNodes();
+
+  int *startIndex;
+  int *columns;
+  double *values;
+  _projectionMatrix->getMatrix(startIndex,columns,values);
+
+  int nbFields = receiver->getNbFields();
+  if (nbFields != donor->getNbFields())
+    Msg::Fatal("Number of fields in donor and receiver dofContainer must match (for now)");
+  // full dofContainer data
+//   fullVector<double>& rVector = receiver->getVector();
+//   fullVector<double>& dVector = donor->getVector();
+  
+  for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {
+    const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
+    for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {
+      for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+        int jGroup = 0;
+        for (int i = startIndex[iProj++]; i < startIndex[iProj] ; i++){
+          int jProj = columns[i];
+          while (jProj > dGroupsStartIGlobal[jGroup+1]) {
+            jGroup++;
+          }
+          int jElement = (jProj-dGroupsStartIGlobal[jGroup])/dGroups->getElementGroup(jGroup)->getNbNodes();
+          int jNode = jProj-dGroupsStartIGlobal[jGroup]-jElement*dGroups->getElementGroup(jGroup)->getNbNodes();
+          for (int m = 0; m < nbFields; m++){
+            double val = values[i];
+            if (copy) val = (fabs(values[i]) > 1e-8) ? 1 : 0;
+            if (transpose)
+              receiver->getGroupProxy(jGroup)(jNode,jElement*nbFields+m) += val*donor->getGroupProxy(iGroup)(iNode,iElement*nbFields+m);
+//               rVector(jProj) += val*donor->getGroupProxy(iGroup)(iNode,iElement*nbFields+m); // only for one field
+            else
+              receiver->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += val*donor->getGroupProxy(jGroup)(jNode,jElement*nbFields+m);
+//               receiver->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += val*dVector(jProj); // only for one field
+          }
+        }
+      }
+    }
+  }
+}
+
+#include "LuaBindings.h"
+void dgMesh2MeshProjection::registerBindings(binding *b){
+  classBinding *cb = b->addClass<dgMesh2MeshProjection>("dgMesh2MeshProjection");
+  cb->setDescription("The dgMesh2MeshProjection class provides methods for projecting fields from one mesh on another.");
+  methodBinding *cm;
+  cm = cb->setConstructor<dgMesh2MeshProjection,dgDofContainer*,dgDofContainer*>();
+  cm->setDescription("Build a projection instance between two dgDofContainers. Integration is performend in the donor Dof function space.");
+  cm->setArgNames("donorDof","receiverDof",NULL);
+  cm = cb->addMethod("projectFromTo",&dgMesh2MeshProjection::projectFromTo);
+  cm->setDescription("Projects all the fields in one container to the other");
+  cm->setArgNames("fromContainer","toContainer",NULL);
+  cm = cb->addMethod("copyFromTo",&dgMesh2MeshProjection::copyFromTo);
+  cm->setDescription("Copies all the fields on one container to the other (projection with binary projection matrix)");
+  cm->setArgNames("fromContainer","toContainer",NULL);
+}

Solver/dgMesh2MeshProjection.h

+#ifndef _DG_MESH2MESH_PROJECTION_H_
+#define _DG_MESH2MESH_PROJECTION_H_
+#include "GmshConfig.h"
+#include "dgDofContainer.h"
+// #include "function.h"
+#include "dgGroupOfElements.h"
+// #include "dgConservationLaw.h"
+// #include "dgAlgorithm.h"
+#include "GModel.h"
+#ifdef HAVE_MPI
+#include "mpi.h"
+#else
+#include "string.h"
+#endif
+// #include <sstream>
+#include "MElement.h"
+// #include <benchmarks/misc/stipple.pos>
+
+// TODO: better representation of matrices, faster matrix-vector matrix-matrix products
+// projecting 3D integration points to 2D plane, instead of assuming z=0
+// possibility to use spesific integration rules
+class dgMesh2MeshProjection {
+public: class dofStoragePattern;
+private:
+  bool _projectionMatrixIsBuilt, _useDofContainerOctree;
+  Octree* _octree;
+  int _numInDofs, _numOutDofs, _inMaxDim, _outMaxDim;
+  dofStoragePattern* _inDofPattern;
+  dofStoragePattern* _outDofPattern;
+  linearSystemCSRGmm<double>* _projectionMatrix;
+  void _buildProjectionMatrix(dgDofContainer* donor, dgDofContainer* receiver);
+  void _multiplyByInverseMassMatrix(dgDofContainer* donor);
+  void _applyProjection(dgDofContainer* donor, dgDofContainer* receiver, bool transpose, bool copy);
+  bool _checkTranspose(dgDofContainer* donor,dgDofContainer* receiver);
+  void _buildReceiverOctree(dgDofContainer* receiver);
+  MElement* _getReceiverElementByCoord(SPoint3 point);
+public:
+  // Contains the number of groups, elements and nodes to check if DOF containers are mutually compatible
+  class dofStoragePattern {
+    private:
+      int _nbGroups, _nbFields, _nbProjDofs;
+      std::vector<int> _nbElems, _nbNodes, _groupFirstProjectionDofId;
+    public:
+      dofStoragePattern(dgDofContainer* dofContainer);
+      bool operator == (const dofStoragePattern& other) const;
+      // indices needed when accessing the projection matrix
+      inline int getProjDofId (int groupId, int elementId, int nodeId) const {
+        // this is different from dgDofContainer::getDofId() because nbFields is ignored
+        return _groupFirstProjectionDofId[groupId]+elementId*_nbNodes[groupId]+nodeId;
+      }
+  };
+  dgMesh2MeshProjection();
+  dgMesh2MeshProjection(dgDofContainer* donor,dgDofContainer* receiver);
+  void setIntegrationRule();
+  void projectFromTo(dgDofContainer* donor,dgDofContainer* receiver);
+  void copyFromTo(dgDofContainer* donor,dgDofContainer* receiver);
+  static void registerBindings(binding *b);
+};
+
+#endif