Skip to content
Snippets Groups Projects
Commit 6216cadd authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

refined algorithm to set interpolation matrices for model-based views

parent 0a8e8561
No related branches found
No related tags found
No related merge requests found
...@@ -118,6 +118,14 @@ int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double> ...@@ -118,6 +118,14 @@ int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>
return 0; return 0;
} }
bool PViewData::haveInterpolationMatrices(int type)
{
if(!type)
return !_interpolation.empty();
else
return _interpolation.count(type) ? true : false;
}
void PViewData::removeInterpolationScheme(const std::string &name) void PViewData::removeInterpolationScheme(const std::string &name)
{ {
std::map<std::string, interpolationMatrices>::iterator it = _interpolationSchemes.find(name); std::map<std::string, interpolationMatrices>::iterator it = _interpolationSchemes.find(name);
......
...@@ -200,7 +200,7 @@ class PViewData { ...@@ -200,7 +200,7 @@ class PViewData {
const fullMatrix<double> &coefGeo, const fullMatrix<double> &coefGeo,
const fullMatrix<double> &expGeo); const fullMatrix<double> &expGeo);
int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p); int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
inline bool haveInterpolationMatrices(){ return !_interpolation.empty(); } bool haveInterpolationMatrices(int type=0);
// access to global interpolation schemes // access to global interpolation schemes
static void removeInterpolationScheme(const std::string &name); static void removeInterpolationScheme(const std::string &name);
......
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
// bugs and problems to <gmsh@geuz.org>. // bugs and problems to <gmsh@geuz.org>.
#include "PViewDataGModel.h" #include "PViewDataGModel.h"
#include "MPoint.h"
#include "MLine.h" #include "MLine.h"
#include "MTriangle.h" #include "MTriangle.h"
#include "MQuadrangle.h" #include "MQuadrangle.h"
...@@ -25,6 +26,63 @@ PViewDataGModel::~PViewDataGModel() ...@@ -25,6 +26,63 @@ PViewDataGModel::~PViewDataGModel()
for(unsigned int i = 0; i < _steps.size(); i++) delete _steps[i]; for(unsigned int i = 0; i < _steps.size(); i++) delete _steps[i];
} }
static MElement *_getOneElementOfGivenType(GModel *m, int type)
{
switch(type){
case TYPE_PNT:
for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); it++){
if((*it)->points.size()) return (*it)->points[0];
}
break;
case TYPE_LIN:
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
if((*it)->lines.size()) return (*it)->lines[0];
}
break;
case TYPE_TRI:
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
if((*it)->triangles.size()) return (*it)->triangles[0];
}
break;
case TYPE_QUA:
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
if((*it)->quadrangles.size()) return (*it)->quadrangles[0];
}
break;
case TYPE_POLYG:
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
if((*it)->polygons.size()) return (*it)->polygons[0];
}
break;
case TYPE_TET:
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->tetrahedra.size()) return (*it)->tetrahedra[0];
}
break;
case TYPE_HEX:
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->hexahedra.size()) return (*it)->hexahedra[0];
}
break;
case TYPE_PRI:
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->prisms.size()) return (*it)->prisms[0];
}
break;
case TYPE_PYR:
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->pyramids.size()) return (*it)->pyramids[0];
}
break;
case TYPE_POLYH:
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->polyhedra.size()) return (*it)->polyhedra[0];
}
break;
}
return 0;
}
bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolationScheme) bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolationScheme)
{ {
if(computeMinMax){ if(computeMinMax){
...@@ -64,6 +122,12 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat ...@@ -64,6 +122,12 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
} }
} }
// set up interpolation matrices
if(!haveInterpolationMatrices()){
GModel *model = _steps[0]->getModel();
// if an interpolation scheme is explicitly provided, use it
if(interpolationScheme.size()){ if(interpolationScheme.size()){
interpolationMatrices m = _interpolationSchemes[interpolationScheme]; interpolationMatrices m = _interpolationSchemes[interpolationScheme];
if(m.size()) if(m.size())
...@@ -74,66 +138,53 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat ...@@ -74,66 +138,53 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
interpolationScheme.c_str()); interpolationScheme.c_str());
for(interpolationMatrices::iterator it = m.begin(); it != m.end(); it++){ for(interpolationMatrices::iterator it = m.begin(); it != m.end(); it++){
if(it->second.size() == 2){ if(it->second.size() == 2){
// use provided interpolation matrices for field interpolation
// and use geometrical interpolation matrices from the mesh if
// the mesh is curved
MElement *e = _getOneElementOfGivenType(model, it->first);
if(e && e->getPolynomialOrder() > 1 && e->getFunctionSpace()){
const polynomialBasis *fs = e->getFunctionSpace();
setInterpolationMatrices(it->first, *(it->second[0]), *(it->second[1]),
fs->coefficients, fs->monomials);
}
else
setInterpolationMatrices(it->first, *(it->second[0]), *(it->second[1])); setInterpolationMatrices(it->first, *(it->second[0]), *(it->second[1]));
// we should maybe automatically add the geometrical
// interpolation matrices here (if geometrically the elements
// are of order > 1)
} }
else if(it->second.size() == 4) else if(it->second.size() == 4){
// use provided matrices for field and geometry
setInterpolationMatrices(it->first, *it->second[0], *it->second[1], setInterpolationMatrices(it->first, *it->second[0], *it->second[1],
*it->second[2], *it->second[3]); *it->second[2], *it->second[3]);
}
else else
Msg::Error("Wrong number of interpolation matrices (%d) for scheme '%s'", Msg::Error("Wrong number of interpolation matrices (%d) for scheme '%s'",
(int)it->second.size(), interpolationScheme.c_str()); (int)it->second.size(), interpolationScheme.c_str());
} }
} }
else if(!haveInterpolationMatrices()){
// add default interpolation matrices for known element types
for(int step = 0; step < getNumTimeSteps(); step++){
if(!_steps[step]->getNumData()) continue;
GModel *m = _steps[step]->getModel();
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
if((*it)->lines.size())
_addInterpolationMatricesForElement((*it)->lines[0]);
}
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
if((*it)->triangles.size())
_addInterpolationMatricesForElement((*it)->triangles[0]);
if((*it)->quadrangles.size())
_addInterpolationMatricesForElement((*it)->quadrangles[0]);
if((*it)->polygons.size())
_addInterpolationMatricesForElement((*it)->polygons[0]);
}
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
if((*it)->tetrahedra.size())
_addInterpolationMatricesForElement((*it)->tetrahedra[0]);
if((*it)->hexahedra.size())
_addInterpolationMatricesForElement((*it)->hexahedra[0]);
if((*it)->prisms.size())
_addInterpolationMatricesForElement((*it)->prisms[0]);
if((*it)->pyramids.size())
_addInterpolationMatricesForElement((*it)->pyramids[0]);
if((*it)->polyhedra.size())
_addInterpolationMatricesForElement((*it)->polyhedra[0]);
}
}
}
return PViewData::finalize(); // if we don't have interpolation matrices for a given element type,
} // assume isoparametric elements
int types[] = {TYPE_PNT, TYPE_LIN, TYPE_TRI, TYPE_QUA, TYPE_TET, TYPE_HEX,
void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e) TYPE_PRI, TYPE_PYR, TYPE_POLYG, TYPE_POLYH};
{ for(int i = 0; i < sizeof(types) / sizeof(types[0]); i++){
int type = e->getType(); if(!haveInterpolationMatrices(types[i])){
MElement *e = _getOneElementOfGivenType(model, types[i]);
if(e){
const polynomialBasis *fs = e->getFunctionSpace(); const polynomialBasis *fs = e->getFunctionSpace();
if(fs){ if(fs){
if(e->getPolynomialOrder() > 1) if(e->getPolynomialOrder() > 1)
setInterpolationMatrices(type, fs->coefficients, fs->monomials, setInterpolationMatrices(types[i], fs->coefficients, fs->monomials,
fs->coefficients, fs->monomials); fs->coefficients, fs->monomials);
else else
setInterpolationMatrices(type, fs->coefficients, fs->monomials); setInterpolationMatrices(types[i], fs->coefficients, fs->monomials);
} }
} }
}
}
}
return PViewData::finalize();
}
MElement *PViewDataGModel::_getElement(int step, int ent, int ele) MElement *PViewDataGModel::_getElement(int step, int ent, int ele)
{ {
......
...@@ -166,8 +166,6 @@ class PViewDataGModel : public PViewData { ...@@ -166,8 +166,6 @@ class PViewDataGModel : public PViewData {
// cache last element to speed up loops // cache last element to speed up loops
MElement *_getElement(int step, int ent, int ele); MElement *_getElement(int step, int ent, int ele);
MVertex *_getNode(MElement *e, int nod); MVertex *_getNode(MElement *e, int nod);
// helper function to populate the interpolation matrix list
void _addInterpolationMatricesForElement(MElement *e);
public: public:
PViewDataGModel(DataType type=NodeData); PViewDataGModel(DataType type=NodeData);
~PViewDataGModel(); ~PViewDataGModel();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment