diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index c2db76fa759e4024f3f65cf22d314406fe0fe221..3b12907cd20ebffe709244c831cd892b94201c67 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -30,6 +30,7 @@
#include "dgResidual.h"
#include "drawContext.h"
#include "GmshMessage.h"
+#include "linearSystemCSR.h"
extern "C" {
#include "lua.h"
@@ -368,6 +369,7 @@ binding::binding(){
functionLua::registerBindings(this);
gmshOptions::registerBindings(this);
Msg::registerBindings(this);
+ linearSystemCSRGmm<double>::registerBindings(this);
}
binding *binding::_instance=NULL;
#endif
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 32df24e05fe232ab9d5e454ef4508f4cf8b2cc0b..3911f6d59039cc80d2c9b02aa57314d2887e14ad 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -4,6 +4,7 @@
#include "dgGroupOfElements.h"
#include "dgConservationLaw.h"
#include "dgAlgorithm.h"
+#include "GModel.h"
#ifdef HAVE_MPI
#include "mpi.h"
#else
@@ -242,10 +243,77 @@ void dgDofContainer::L2Projection(std::string functionName){
}
}
-void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor){
-// projector.allocate();
-// projector.addToMatrix(,,);
-/*nt nbCoarseNodes = 0; // OLD CODE
+void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor){
+// for basis function i
+// for integration pts
+// search for element in 2d mesh
+// for basis function j
+// evauate quad
+// add to matrix correct indices
+
+ dgGroupCollection* dGroups = donor.getGroups();
+ dgGroupCollection* rGroups = this->getGroups();
+ std::vector<int> rGroupsStartIGlobal(rGroups->getNbElementGroups() + 1);
+ int jGlobal = 0; // indices in global container
+ rGroupsStartIGlobal[0] = 0;
+ for (int i = 1; i<rGroupsStartIGlobal.size(); i++)
+ rGroupsStartIGlobal[i] = rGroupsStartIGlobal[i-1] + rGroups->getElementGroup(i-1)->getNbElements()*rGroups->getElementGroup(i-1)->getNbNodes();
+
+ projector.allocate(rGroupsStartIGlobal[rGroupsStartIGlobal.size()-1]);
+ for (int iGroup=0;iGroup<dGroups->getNbElementGroups();iGroup++) {// for 3d.groups
+ const dgGroupOfElements &dGroup = *dGroups->getElementGroup(iGroup);
+ fullMatrix<double> iPtsMatrix = dGroup.getIntegrationPointsMatrix();
+ for (int iElement=0 ; iElement<dGroup.getNbElements() ;++iElement) {// for elements
+ double rShapeFun[256];
+ double dShapeFun[256];
+ GModel *rModel = rGroups->getModel();
+ GModel *dModel = dGroups->getModel();
+// dataCacheDouble xyz = cacheMap.get("XYZ");
+ MElement* dElem = dGroup.getElement(iElement);
+ for (int iPt =0; iPt< dGroup.getNbIntegrationPoints(); 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];
+ }
+ z = 0;
+ SPoint3 p(x,y,z); // find p in 2D mesh
+ MElement *rElem = rGroups->getModel()->getMeshElementByCoord(p);
+ int iGroupR,iElementR;
+ rGroups->find(rElem,iGroupR,iElementR);
+ if (iElementR == -1)
+ Msg::Error("Integration point not found in receiver mesh.");
+ const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroupR);
+ 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 (iElement, iPt);
+ int iGlobal = rGroupsStartIGlobal[iGroupR]+rGroup.getNbNodes()*iElementR;
+ for (int jPsi=0;jPsi<dGroup.getNbNodes();jPsi++){
+ for (int iPsi=0;iPsi<rGroup.getNbNodes();iPsi++){
+ // get global index of shape functions
+ // sum of nElems in previous groups
+ // iElem + i|j
+ double val = rShapeFun[iPsi]*dShapeFun[jPsi];//*iPtsMatrix(iPt,3);//*detJ;
+ projector.addToMatrix(iGlobal+iPsi,jGlobal+jPsi,val);
+ }
+ }
+ }
+ jGlobal += dGroup.getNbNodes();
+ }
+ }
+}
+
+void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor){
+ scale(0.);
+ int *startIndex;
+ int *columns;
+ double *values;
+ projector.getMatrix(startIndex,columns,values);
+/*nt nbCoarseNodes = 0; // OLD CODE
int nbCoarseElements = 0;
int r = 0;
for (int kFine = 0;kFine < fineGroups.getNbElementGroups();kFine++) {
@@ -272,14 +340,6 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSR<double> &projec
}*/
}
-void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor){
- scale(0.);
- int *startIndex;
- int *columns;
- double *values;
- projector.getMatrix(startIndex,columns,values);
-}
-
void dgDofContainer::exportGroupIdMsh(){
// the elementnodedata::export does not work !!
@@ -517,4 +577,7 @@ void dgDofContainer::registerBindings(binding *b){
cm = cb->addMethod("axpy",(void (dgDofContainer::*)(dgDofContainer &,double)) &dgDofContainer::axpy);
cm->setArgNames("dofContainer","a",NULL);
cm->setDescription("this = this+a*dofContainer");
+ cm = cb->addMethod("Mesh2Mesh_BuildL2Projection",&dgDofContainer::Mesh2Mesh_BuildL2Projection);
+ cm->setDescription("Build L2 projection matrix between donor and this dofContainer.");
+ cm->setArgNames("projector","donor",NULL);
}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 29318b62f1973f7761159ff062f3b0ba313edaae..2565dd5cbf26f24faa35b7451540685380b50aa6 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -43,8 +43,8 @@ public:
void load(const std::string name);
void setAll(double v);
void L2Projection(std::string functionName);
- void Mesh2Mesh_BuildL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor);
- void Mesh2Mesh_ApplyL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor);
+ void Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
+ void Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
void exportMsh(const std::string name);
void exportGroupIdMsh();
void exportMultirateGroupMsh();
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 570657b5c450bc6e952f7027f9b458f6f8c81b86..3f1f221edac2ecb18ef648bfe4910d02c142e4ff 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -301,10 +301,13 @@ template<>
#include "Bindings.h"
template<>
- void linearSystemCSR<double>::registerBindings(binding *b)
+ void linearSystemCSRGmm<double>::registerBindings(binding *b)
{
- classBinding *cb = b->addClass< linearSystemCSR<double> >("linearSystemCSRdouble");
-// cb->setDescription("A GModel contains a geometrycal and it's mesh.");
+ classBinding *cb = b->addClass< linearSystemCSRGmm<double> >("linearSystemCSRGmmdouble");
+ cb->setDescription("Sparse matrix representation.");
+ methodBinding *cm;
+ cm = cb->setConstructor<linearSystemCSRGmm<double> >();
+ cm->setDescription("Build an empty container");
}
#if defined(HAVE_GMM)
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 58e343741af4e59b276d249d98f5fb3c815caa37..9e75a2d3f378ed2a456759e8a72d190baffa6eb9 100644
--- a/Solver/linearSystemCSR.h
+++ b/Solver/linearSystemCSR.h
@@ -114,7 +114,6 @@ class linearSystemCSR : public linearSystem<scalar> {
{
for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
}
- static void registerBindings(binding *b);
};
template <class scalar>
@@ -136,6 +135,7 @@ class linearSystemCSRGmm : public linearSystemCSR<scalar> {
}
#endif
;
+ static void registerBindings(binding *b);
};
template <class scalar>