diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index fbe10d714a617eb499bc3584f2062030f1dff22c..8978c5d8780e705356a1d223d821f0a03b07eeae 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -20,8 +20,6 @@ void lloydAlgorithm::operator () ( GFace * gf) {
}
}
- // assume every
-
// Create a triangulator
DocRecord triangulator (all.size());
diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index e11f1cd7600d6b4e86790911f47a4b1c4e78df8a..286b6d85e282720cc36db1d4e5b1271b086716eb 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -8,6 +8,7 @@ groups = dgGroupCollection(model, dimension, order)
-- groups:split...
groups:buildGroupsOfInterfaces()
solution = dgDofContainer(groups, claw:getNbFields())
+--[[
a = fullMatrix(3,1);
a:set(0,0, 1.);
a:set(1,0, 2.);
@@ -15,8 +16,9 @@ a:set(2,0, 3.);
f = functionConstant(a):getName()
solution:L2Projection(f)
solution:exportMsh('output/init')
+]]--
-for i=1,100000 do
+for i=1,00000 do
norm = dg:RK44(150*(3/(2.*order+1)))
if ( i%100 ==0 ) then
print ('iter ', i, norm)
diff --git a/Solver/TESTCASES/mesh2mesh.lua b/Solver/TESTCASES/mesh2mesh.lua
new file mode 100644
index 0000000000000000000000000000000000000000..40cd190f1f661ee7d9419f26960b78b01d36e699
--- /dev/null
+++ b/Solver/TESTCASES/mesh2mesh.lua
@@ -0,0 +1,47 @@
+-- initial condition
+function initial_condition( xyz , f )
+ for i=0,xyz:size1()-1 do
+ x = xyz:get(i,0)
+ y = xyz:get(i,1)
+ z = xyz:get(i,2)
+ f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
+ f:set (i, 1, math.exp(-100*((x-0.4)^2 +(y-0.3)^2)))
+ f:set (i, 2, math.exp(-100*((x+0.2)^2 +(y-0.3)^2)))
+ end
+end
+
+dimension = 2
+order1 = 1
+order2 = 2
+
+model1 = GModel()
+model2 = GModel()
+-- load a mesh
+model1:load('square1.msh')
+-- load another mesh
+model2:load('square2.msh');
+
+-- create 2 groups related to the 2 models with different orders
+groups1 = dgGroupCollection (model1 , dimension, order1)
+groups1:buildGroupsOfInterfaces()
+groups2 = dgGroupCollection (model2 , dimension, order2)
+groups2:buildGroupsOfInterfaces()
+
+-- create 2 solutions
+solution1 = dgDofContainer(groups1, 3);
+solution2 = dgDofContainer(groups2, 3);
+
+-- apply initial solution to solution1
+solution1:L2Projection(functionLua(3,'initial_condition',{'XYZ'}):getName())
+-- save the stuff
+solution1:exportMsh('solution1');
+
+-- this is the function that interpolates solution1
+F = functionMesh2Mesh(solution1)
+
+-- project solution 1 onto mesh 2
+solution2:L2Projection(F:getName())
+-- save the stuff
+solution2:exportMsh('solution2');
+
+
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index bcbdb732ea72709d17dd4fd3a3c75e64ead4cfe3..c9659752ffd174eee80934a3a50c1c34ca6db940 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -32,11 +32,13 @@ public:
dgDofContainer (dgGroupCollection *groups, int nbFields);
~dgDofContainer ();
int getNbElements() {return _totalNbElements;}
+ int getNbFields() {return _nbFields;}
void scatter();
void save(const std::string name);
void load(const std::string name);
void setAll(double v);
void L2Projection(std::string functionName);
+ void Mesh2Mesh_L2Projection(dgDofContainer &other);
void exportMsh(const std::string name);
static void registerBindings(binding *b);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c1cfdab539c5420ac78400f46a3bffebe418f94f..551752deec2e9efce10c6fb581497e7ca05eb1ec 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -53,7 +53,9 @@ static fullMatrix<double> dgGetFaceIntegrationRuleOnElement (
}
-dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder,int ghostPartition)
+dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
+ int polyOrder,
+ int ghostPartition)
: _elements(e),
_fs(*_elements[0]->getFunctionSpace(polyOrder)),
_integration(dgGetIntegrationRule (_elements[0], polyOrder)),
@@ -76,6 +78,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
for (int i=0;i<_elements.size();i++){
MElement *e = _elements[i];
+ element_to_index[e] = i;
fullMatrix<double> imass(*_imass,nbPsi*i,nbPsi);
for (int j=0;j< _integration->size1() ; j++ ){
_fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
@@ -921,6 +924,13 @@ dgGroupCollection::~dgGroupCollection() {
delete _ghostGroups[i];
}
+void dgGroupCollection::find (MElement*e, int &ig, int &index){
+ for (ig=0;ig<_elementGroups.size();ig++){
+ index = _elementGroups[ig]->getIndexOfElement(e);
+ if (index != -1)return;
+ }
+}
+
#include "LuaBindings.h"
void dgGroupCollection::registerBindings(binding *b){
classBinding *cb = b->addClass<dgGroupCollection>("dgGroupCollection");
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 063b0afd6ee3356ed14905f3176da40160c778c9..41c6f46450d61d9d12ed016a5ae609f363453bf8 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -46,6 +46,8 @@ class dgGroupOfElements {
int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition
// N elements in the group
std::vector<MElement*> _elements;
+ // inverse map that gives the index of an element in the group
+ std::map<MElement*,int> element_to_index;
// the ONLY function space that is used to
// inerpolated the fields (may be different to the
// one of the elements)
@@ -99,6 +101,11 @@ public:
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
inline int getOrder() const {return _order;}
inline int getGhostPartition() const {return _ghostPartition;}
+ inline int getIndexOfElement (MElement *e) const {
+ std::map<MElement*,int>::const_iterator it = element_to_index.find(e);
+ if (it == element_to_index.end())return -1;
+ return it->second;
+ }
};
class dgGroupOfFaces {
@@ -199,6 +206,7 @@ class dgGroupCollection {
bool _groupsOfElementsBuilt;
bool _groupsOfInterfacesBuilt;
public:
+ inline GModel* getModel() {return _model;}
inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();}
inline int getNbBoundaryGroups() const {return _boundaryGroups.size();}
@@ -217,6 +225,8 @@ class dgGroupCollection {
void splitGroupsForMultirate(double refDt,dgConservationLaw *claw);
+ void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
+
dgGroupCollection(GModel *model, int dimension, int order);
dgGroupCollection();
static void registerBindings(binding *b);
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 0cafe0ced7b0c185175c5047268f1239a879456a..391402b9a5014cc7c1b64a696ef6b0b84a95468e 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -34,7 +34,6 @@ class dgSystemOfEquations {
public:
const dgConservationLaw * getLaw() const {return _claw;}
const GModel * getModel() const {return _gm;}
- std::pair<dgGroupOfElements*,int> getElementPosition (MElement *);
void setOrder (int order); // set the polynomial order
void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
dgSystemOfEquations(GModel *_gm);
diff --git a/Solver/function.cpp b/Solver/function.cpp
index d0f5fbb8e558ab40e86f32bc418575908294b979..15afbfa243e6ff81182ba4dc73754a9e7e161621 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -3,6 +3,9 @@
#include "function.h"
#include "SPoint3.h"
#include "MElement.h"
+#include "dgDofContainer.h"
+#include "dgGroupOfElements.h"
+#include "GModel.h"
// dataCache members
dataCache::dataCache(dataCacheMap *cacheMap) : _valid(false) {
@@ -214,6 +217,7 @@ dataCacheDouble *functionConstant::newDataCache(dataCacheMap *m)
{
return new data(this,m);
}
+
functionConstant::functionConstant(const fullMatrix<double> *source){
_source = *source;
static int c=0;
@@ -225,40 +229,66 @@ functionConstant::functionConstant(const fullMatrix<double> *source){
// function that enables to interpolate a DG solution using
// geometrical search in a mesh
-/*
-class functionSystemOfEquations::data : public dataCacheDouble {
- dgSystemOfEquations *_sys;
+
+functionMesh2Mesh::functionMesh2Mesh(dgDofContainer *dofc)
+ : _dofContainer(dofc)
+{
+ static int c=0;
+ std::ostringstream oss;
+ oss<<"FunctionMesh2Mesh_"<<c++;
+ _name = oss.str();
+ function::add(_name,this);
+}
+
+class functionMesh2Mesh::data : public dataCacheDouble {
+ dgDofContainer *_dofContainer;
dataCacheDouble &xyz;
public:
- data(dataCacheMap &m, dgSystemOfEquations *sys) :
- _sys(sys),
- xyz(cacheMap.get("Solution",this)),
- dataCacheDouble(m.getNbEvaluationPoints(), sys->getLaw()->getNbFields())
+ data(dataCacheMap &m, dgDofContainer *sys) :
+ _dofContainer(sys),
+ xyz(m.get("XYZ",this)),
+ dataCacheDouble(m,m.getNbEvaluationPoints(), sys->getNbFields())
{
}
void _eval() {
int nP =xyz().size1();
if(_value.size1() != nP)
- _value = fullMatrix<double>(nP,sys->getLaw()->getNbFields());
+ _value = fullMatrix<double>(nP,_dofContainer->getNbFields());
_value.setAll(0.0);
double fs[256];
+ fullMatrix<double> solEl;
+ GModel *m = _dofContainer->getGroups()->getModel();
for (int i=0;i<_value.size1();i++){
const double x = xyz(i,0);
const double y = xyz(i,1);
const double z = xyz(i,2);
- MElement *e = _sys->getModel()->getMeshElementByCoord(SPoint3(x,y,z));
- std::pair<dgGroupOfElements*,int> location = _sys->getElementPosition(e);
- double U[3],X[3]={xyz(i,0),xyz(i,1),xyz(i,2)};
+ SPoint3 p(x,y,z);
+ MElement *e = m->getMeshElementByCoord(p);
+ int ig,index;
+ _dofContainer->getGroups()->find (e,ig,index);
+ dgGroupOfElements *group = _dofContainer->getGroups()->getElementGroup(ig);
+ double U[3],X[3]={x,y,z};
e->xyz2uvw (X,U);
- location.first->getFunctionSpace().f(U[0],U[1],U[2],fs);
+ group->getFunctionSpace().f(U[0],U[1],U[2],fs);
+ fullMatrix<double> &sol = _dofContainer->getGroupProxy(ig);
+ solEl.setAsProxy(sol,index*_dofContainer->getNbFields(),_dofContainer->getNbFields());
+ int fSize = group->getNbNodes();
+ for (int k=0;k<_dofContainer->getNbFields();k++){
+ _value(i,k) = 0.0;
+ for (int j=0;j<fSize;j++){
+ _value(i,k) += solEl(j,k)*fs[j];
+ }
+ }
}
}
- dataCacheDouble *newDataCache(dataCacheMap *m)
- {
- return new data(this,_sys);
- }
};
-*/
+
+dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
+{
+ printf("coussdo %d %d\n",m->getNbEvaluationPoints(),_dofContainer->getNbFields());
+ return new data(*m,_dofContainer);
+}
+
#include "Bindings.h"
@@ -287,5 +317,13 @@ void function::registerBindings(binding *b){
mb->setArgNames("fileName","coordinateFunction",NULL);
mb->setDescription("Tri-linearly interpolate through data in file 'fileName' at coordinate given by 'coordinateFunction'.\nThe file format is :\nx0 y0 z0\ndx dy dz\nnx ny nz\nv(0,0,0) v(0,0,1) v(0 0 2) ...");
+ cb = b->addClass<functionMesh2Mesh>("functionMesh2Mesh");
+ cb->setDescription("A function that can be used to interpolate into a given mesh");
+ mb = cb->setConstructor<functionMesh2Mesh,dgDofContainer*>();
+ mb->setArgNames("solution",NULL);
+ mb->setDescription("A solution.");
+ cb->setParentClass<function>();
+
+
}
diff --git a/Solver/function.h b/Solver/function.h
index 9a81cbfa7d2311d1224fad62bcb0b600a29e2bd4..333dca348a904abb485900467a66ec52ea71939e 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -8,7 +8,7 @@
class dataCacheMap;
class MElement;
class binding;
-class dgSystemOfEquations;
+class dgDofContainer;
// those classes manage complex function dependencies and keep their values in cache so that they are not recomputed when it is not necessary. To do this, we use three classes : function, dataCache and dataCacheMap. The workflow is :
//
@@ -195,12 +195,12 @@ class functionConstant : public function {
functionConstant(const fullMatrix<double> *source);
};
-class functionSystemOfEquations : public function {
- dgSystemOfEquations *_sys;
+class functionMesh2Mesh : public function {
+ dgDofContainer *_dofContainer;
public:
class data ;
dataCacheDouble *newDataCache(dataCacheMap *m);
- functionSystemOfEquations(dgSystemOfEquations *sys) : _sys(sys) {}
+ functionMesh2Mesh(dgDofContainer *dofc) ;
};
#endif