From 6216caddfac70cf02b29e51f0446cc48b4938631 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 10 May 2011 18:15:00 +0000
Subject: [PATCH] refined algorithm to set interpolation matrices for
model-based views
---
Post/PViewData.cpp | 8 ++
Post/PViewData.h | 2 +-
Post/PViewDataGModel.cpp | 177 +++++++++++++++++++++++++--------------
Post/PViewDataGModel.h | 2 -
4 files changed, 123 insertions(+), 66 deletions(-)
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 9263c9c1b3..dd56dc81ca 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -118,6 +118,14 @@ int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>
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)
{
std::map<std::string, interpolationMatrices>::iterator it = _interpolationSchemes.find(name);
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 2dca280807..bbf47d872c 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -200,7 +200,7 @@ class PViewData {
const fullMatrix<double> &coefGeo,
const fullMatrix<double> &expGeo);
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
static void removeInterpolationScheme(const std::string &name);
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index ecba42296c..96fadd77b9 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -4,6 +4,7 @@
// bugs and problems to <gmsh@geuz.org>.
#include "PViewDataGModel.h"
+#include "MPoint.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
@@ -25,6 +26,63 @@ PViewDataGModel::~PViewDataGModel()
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)
{
if(computeMinMax){
@@ -64,75 +122,68 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
}
}
- if(interpolationScheme.size()){
- interpolationMatrices m = _interpolationSchemes[interpolationScheme];
- if(m.size())
- Msg::Info("Setting interpolation matrices from scheme '%s'",
- interpolationScheme.c_str());
- else
- Msg::Error("Could not find interpolation scheme '%s'",
- interpolationScheme.c_str());
- for(interpolationMatrices::iterator it = m.begin(); it != m.end(); it++){
- if(it->second.size() == 2){
- 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)
- setInterpolationMatrices(it->first, *it->second[0], *it->second[1],
- *it->second[2], *it->second[3]);
+ // set up interpolation matrices
+ if(!haveInterpolationMatrices()){
+
+ GModel *model = _steps[0]->getModel();
+
+ // if an interpolation scheme is explicitly provided, use it
+ if(interpolationScheme.size()){
+ interpolationMatrices m = _interpolationSchemes[interpolationScheme];
+ if(m.size())
+ Msg::Info("Setting interpolation matrices from scheme '%s'",
+ interpolationScheme.c_str());
else
- Msg::Error("Wrong number of interpolation matrices (%d) for scheme '%s'",
- (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]);
+ Msg::Error("Could not find interpolation scheme '%s'",
+ interpolationScheme.c_str());
+ for(interpolationMatrices::iterator it = m.begin(); it != m.end(); it++){
+ 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]));
+ }
+ else if(it->second.size() == 4){
+ // use provided matrices for field and geometry
+ setInterpolationMatrices(it->first, *it->second[0], *it->second[1],
+ *it->second[2], *it->second[3]);
+ }
+ else
+ Msg::Error("Wrong number of interpolation matrices (%d) for scheme '%s'",
+ (int)it->second.size(), interpolationScheme.c_str());
}
- 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]);
+ }
+
+ // 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,
+ TYPE_PRI, TYPE_PYR, TYPE_POLYG, TYPE_POLYH};
+ for(int i = 0; i < sizeof(types) / sizeof(types[0]); i++){
+ if(!haveInterpolationMatrices(types[i])){
+ MElement *e = _getOneElementOfGivenType(model, types[i]);
+ if(e){
+ const polynomialBasis *fs = e->getFunctionSpace();
+ if(fs){
+ if(e->getPolynomialOrder() > 1)
+ setInterpolationMatrices(types[i], fs->coefficients, fs->monomials,
+ fs->coefficients, fs->monomials);
+ else
+ setInterpolationMatrices(types[i], fs->coefficients, fs->monomials);
+ }
+ }
}
}
- }
- return PViewData::finalize();
-}
-
-void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
-{
- int type = e->getType();
- const polynomialBasis *fs = e->getFunctionSpace();
- if(fs){
- if(e->getPolynomialOrder() > 1)
- setInterpolationMatrices(type, fs->coefficients, fs->monomials,
- fs->coefficients, fs->monomials);
- else
- setInterpolationMatrices(type, fs->coefficients, fs->monomials);
}
+
+ return PViewData::finalize();
}
MElement *PViewDataGModel::_getElement(int step, int ent, int ele)
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 5acdad0557..783c7cad6b 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -166,8 +166,6 @@ class PViewDataGModel : public PViewData {
// cache last element to speed up loops
MElement *_getElement(int step, int ent, int ele);
MVertex *_getNode(MElement *e, int nod);
- // helper function to populate the interpolation matrix list
- void _addInterpolationMatricesForElement(MElement *e);
public:
PViewDataGModel(DataType type=NodeData);
~PViewDataGModel();
--
GitLab