diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index d1523dc3deb62845e0cc6e8d7003890407b42c9a..82d4d65b00cf07a3a79d6cf3f06982ea0e687d93 100644
--- a/Common/LuaBindings.cpp
+++ b/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
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index e88ccd212ef704586aa5b5495ecce99a85fb8921..52d251eee9d922af2ffb32d32fb00452d49d28ac 100644
--- a/Solver/CMakeLists.txt
+++ b/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)
diff --git a/Solver/TESTCASES/Projection2D3D.lua b/Solver/TESTCASES/Projection2D3D.lua
index 19b3954c30920f429390e87c37387a2e6bd90c6d..d33e35719536730ff2a212542e177869234f6de0 100644
--- a/Solver/TESTCASES/Projection2D3D.lua
+++ b/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)
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 390f7b9f51e9ab0b2320ac316632c08c055e5ee9..ca0647e6ed94d2f6ba3a2d24323ce7f8feb44d2c 100644
--- a/Solver/dgDofContainer.cpp
+++ b/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);
}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 97569432a0465ac03718cae7b29296475a89c3a5..595d4dc1770172c23a0dfd2a9c78e4d5f315866a 100644
--- a/Solver/dgDofContainer.h
+++ b/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();
diff --git a/Solver/dgMesh2MeshProjection.cpp b/Solver/dgMesh2MeshProjection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..97146a6ecb8d50223b1f4d651b3de12543cec9c9
--- /dev/null
+++ b/Solver/dgMesh2MeshProjection.cpp
@@ -0,0 +1,289 @@
+#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);
+}
diff --git a/Solver/dgMesh2MeshProjection.h b/Solver/dgMesh2MeshProjection.h
new file mode 100644
index 0000000000000000000000000000000000000000..925ec53c219134777c1eeed3331c1a78495559aa
--- /dev/null
+++ b/Solver/dgMesh2MeshProjection.h
@@ -0,0 +1,60 @@
+#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