Skip to content
Snippets Groups Projects
Commit ced88a13 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

JFR has added a function that enables to d mesh to mesh interpolation

an example (mesh2mesh.lua) has been added
parent bdf82447
No related branches found
No related tags found
No related merge requests found
......@@ -20,8 +20,6 @@ void lloydAlgorithm::operator () ( GFace * gf) {
}
}
// assume every
// Create a triangulator
DocRecord triangulator (all.size());
......
......@@ -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)
......
-- 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');
......@@ -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);
......
......@@ -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");
......
......@@ -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);
......
......@@ -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);
......
......@@ -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>();
}
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment