diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 71ebfbc725bf32c27bdd12b74c90679f16603382..78167989a72703e3f11e2bea41aa2b2ce32bddd3 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -269,7 +269,7 @@ void statisticsWindow::compute(bool elementQuality)
Field *f = fields->get(fields->background_field);
int nbEdges = edges.size();
printf("nb edges =%d \n", nbEdges);
- system("rm qualEdges.txt");
+ if(system("rm qualEdges.txt"));
FILE *fp = fopen("qualEdges.txt", "w");
std::vector<int> qualE;
int nbS = 50;
@@ -347,7 +347,7 @@ void statisticsWindow::compute(bool elementQuality)
FieldManager *fields = GModel::current()->getFields();
Field *f = fields->get(fields.getBackgroundField());
int nbEdges = edges.size();
- system("rm qualEdges.txt");
+ if(system("rm qualEdges.txt"));
FILE *fp = fopen("qualEdges.txt", "w");
std::vector<int> qualE;
int nbS = 50;
diff --git a/Plugin/Curl.cpp b/Plugin/Curl.cpp
index bec30666b9dbd5736d9613eb4fbd2279e4db3553..e10c5e3fda8510ec45535950f8aa02da0e7e017a 100644
--- a/Plugin/Curl.cpp
+++ b/Plugin/Curl.cpp
@@ -60,9 +60,9 @@ PView *GMSH_CurlPlugin::execute(PView *v)
int numComp = data1->getNumComponents(firstNonEmptyStep, ent, ele);
if(numComp != 3) continue;
int type = data1->getType(firstNonEmptyStep, ent, ele);
- std::vector<double> *out = data2->incrementList(3, type);
- if(!out) continue;
int numNodes = data1->getNumNodes(firstNonEmptyStep, ent, ele);
+ std::vector<double> *out = data2->incrementList(3, type, numNodes);
+ if(!out) continue;
double x[8], y[8], z[8], val[8 * 3];
for(int nod = 0; nod < numNodes; nod++)
data1->getNode(firstNonEmptyStep, ent, ele, nod, x[nod], y[nod], z[nod]);
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index fdf557d39085538e7a5a2dc7978331894287e38a..f0b79894884fae3bdc3349e74e614293cf94946a 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -114,19 +114,28 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
if(_entities[ii]->dim() == _maxDim) {
for(unsigned int i = 0; i < _entities[ii]->getNumMeshElements(); i++){
MElement *e = _entities[ii]->getMeshElement(i);
-
+
int numNodes = e->getNumVertices();
+ if(e->getNumChildren()) numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
- std::vector<double> *out = _data->incrementList(1, e->getType());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
- for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
-
+ std::vector<double> *out = _data->incrementList(1, e->getType(), numNodes);
+ std::vector<MVertex*> nods;
+ if(!e->getNumChildren())
+ for(int i = 0; i < numNodes; i++)
+ nods.push_back(e->getVertex(i));
+ else
+ for(int i = 0; i < e->getNumChildren(); i++)
+ for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
+ nods.push_back(e->getChild(i)->getVertex(j));
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->x());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->y());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->z());
+
if (_maxDim == 3) fprintf(fName,"SS(");
else if (_maxDim == 2) fprintf(fName,"ST(");
std::vector<double> dist;
for(int j = 0; j < numNodes; j++) {
- MVertex *v = e->getVertex(j);
+ MVertex *v = nods[j];
if(j) fprintf(fName,",%.16g,%.16g,%.16g",v->x(),v->y(), v->z());
else fprintf(fName,"%.16g,%.16g,%.16g", v->x(),v->y(), v->z());
std::map<MVertex*, double>::iterator it = _distance_map.find(v);
@@ -157,7 +166,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int id_face = (int) DistanceOptions_Number[2].def;
double type = (double) DistanceOptions_Number[3].def;
int ortho = (int) DistanceOptions_Number[6].def;
-
+
PView *view = new PView();
_data = getDataList(view);
#if defined(HAVE_SOLVER)
@@ -190,7 +199,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
Msg::Error("This plugin needs a mesh !");
return view;
}
-
+
std::vector<SPoint3> pts;
std::vector<double> distances;
std::vector<MVertex* > pt2Vertex;
@@ -285,13 +294,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
}
printView(_entities, _distance_map);
}
-
+
// Compute PDE for distance function
//-----------------------------------
else if (type > 0.0){
-
+
#if defined(HAVE_SOLVER)
-
+
bool existEntity = false;
SBoundingBox3d bbox;
for(unsigned int i = 0; i < _entities.size(); i++){
@@ -339,24 +348,24 @@ PView *GMSH_DistancePlugin::execute(PView *v)
}
}
}
-
+
double L = norm(SVector3(bbox.max(), bbox.min()));
double mu = type*L;
-
+
simpleFunction<double> DIFF(mu*mu), ONE(1.0);
distanceTerm distance(GModel::current(), 1, &DIFF, &ONE);
-
+
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
SElement se((*it));
distance.addToMatrix(*dofView, &se);
}
groupOfElements gr(allElems);
distance.addToRightHandSide(*dofView, gr);
-
+
Msg::Info("Distance Computation: Assembly done");
lsys->systemSolve();
Msg::Info("Distance Computation: System solved");
-
+
for(std::map<MVertex*,double >::iterator itv = _distance_map.begin();
itv != _distance_map.end() ; ++itv){
MVertex *v = itv->first;
@@ -366,9 +375,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
double dist = -mu * log(1. - value);
itv->second = dist;
}
-
+
printView(_entities, _distance_map);
-
+
#endif
}
@@ -376,14 +385,14 @@ PView *GMSH_DistancePlugin::execute(PView *v)
_data->Time.push_back(0);
_data->setFileName(_fileName.c_str());
_data->finalize();
-
-
+
+
//compute also orthogonal vector to distance field
// A Uortho = -C DIST
//------------------------------------------------
if (ortho > 0){
#if defined(HAVE_SOLVER)
-
+
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys2 = new linearSystemCSRTaucs<double>;
#else
@@ -394,9 +403,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
#endif
dofManager<double> myAssembler(lsys2);
simpleFunction<double> ONE(1.0);
-
+
double dMax = 1.0; //EMI TO CHANGE
-
+
std::vector<MElement *> allElems;
for(unsigned int ii = 0; ii < _entities.size(); ii++){
if(_entities[ii]->dim() == _maxDim) {
@@ -418,53 +427,62 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MElement *e = allElems[mid];
MVertex *vFIX = e->getVertex(0);
myAssembler.fixVertex(vFIX, 0, 1, 0.0);
-
+
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
MElement *t = *it;
for(int k = 0; k < t->getNumVertices(); k++)
myAssembler.numberVertex(t->getVertex(k), 0, 1);
}
-
+
orthogonalTerm *ortho;
ortho = new orthogonalTerm(GModel::current(), 1, &ONE, &_distance_map);
// if (type < 0)
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, view);
// else
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, dofView);
-
+
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
SElement se((*it));
ortho->addToMatrix(myAssembler, &se);
}
groupOfElements gr(allElems);
ortho->addToRightHandSide(myAssembler, gr);
-
+
Msg::Info("Orthogonal Computation: Assembly done");
lsys2->systemSolve();
Msg::Info("Orthogonal Computation: System solved");
-
+
PView *view2 = new PView();
PViewDataList *data2 = getDataList(view2);
data2->setName("ortogonal field");
-
+
Msg::Info("Writing orthogonal.pos");
FILE * f5 = fopen("orthogonal.pos","w");
fprintf(f5,"View \"orthogonal\"{\n");
for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
MElement *e = *it;
-
+
int numNodes = e->getNumVertices();
+ if(e->getType() == TYPE_POLYG) numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
- std::vector<double> *out2 = data2->incrementList(1, e->getType());
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->x());
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->y());
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((e->getVertex(nod))->z());
-
+ std::vector<double> *out2 = data2->incrementList(1, e->getType(), numNodes);
+ std::vector<MVertex*> nods;
+ if(!e->getNumChildren())
+ for(int i = 0; i < numNodes; i++)
+ nods.push_back(e->getVertex(i));
+ else
+ for(int i = 0; i < e->getNumChildren(); i++)
+ for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
+ nods.push_back(e->getChild(i)->getVertex(j));
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->x());
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->y());
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->z());
+
if (_maxDim == 3) fprintf(f5,"SS(");
else if (_maxDim == 2) fprintf(f5,"ST(");
std::vector<double> orth;
for(int j = 0; j < numNodes; j++) {
- MVertex *v = e->getVertex(j);
+ MVertex *v = nods[j];
if(j) fprintf(f5,",%g,%g,%g",v->x(),v->y(), v->z());
else fprintf(f5,"%g,%g,%g", v->x(),v->y(), v->z());
double value;
@@ -488,9 +506,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
data2->Time.push_back(0);
data2->setFileName("orthogonal.pos");
data2->finalize();
-
+
#endif
}
-
+
return view;
}
diff --git a/Plugin/Divergence.cpp b/Plugin/Divergence.cpp
index ffabf5a367b97227a25f7d7f6307665f4a7dc131..77719944946854e8bc58732f42c7b13296ec0d3c 100644
--- a/Plugin/Divergence.cpp
+++ b/Plugin/Divergence.cpp
@@ -60,9 +60,9 @@ PView *GMSH_DivergencePlugin::execute(PView *v)
int numComp = data1->getNumComponents(firstNonEmptyStep, ent, ele);
if(numComp != 3) continue;
int type = data1->getType(firstNonEmptyStep, ent, ele);
- std::vector<double> *out = data2->incrementList(1, type);
- if(!out) continue;
int numNodes = data1->getNumNodes(firstNonEmptyStep, ent, ele);
+ std::vector<double> *out = data2->incrementList(1, type, numNodes);
+ if(!out) continue;
double x[8], y[8], z[8], val[8 * 3];
for(int nod = 0; nod < numNodes; nod++)
data1->getNode(firstNonEmptyStep, ent, ele, nod, x[nod], y[nod], z[nod]);
diff --git a/Plugin/Eigenvalues.cpp b/Plugin/Eigenvalues.cpp
index 559c53128e6fd4bb5ede78ae118bea8bfa531c79..f65f9859b7eb71b6b994d9fc14a09e151f2690f0 100644
--- a/Plugin/Eigenvalues.cpp
+++ b/Plugin/Eigenvalues.cpp
@@ -64,11 +64,11 @@ PView *GMSH_EigenvaluesPlugin::execute(PView *v)
int numComp = data1->getNumComponents(0, ent, ele);
if(numComp != 9) continue;
int type = data1->getType(0, ent, ele);
- std::vector<double> *outmin = dmin->incrementList(1, type);
- std::vector<double> *outmid = dmid->incrementList(1, type);
- std::vector<double> *outmax = dmax->incrementList(1, type);
- if(!outmin || !outmid || !outmax) continue;
int numNodes = data1->getNumNodes(0, ent, ele);
+ std::vector<double> *outmin = dmin->incrementList(1, type, numNodes);
+ std::vector<double> *outmid = dmid->incrementList(1, type, numNodes);
+ std::vector<double> *outmax = dmax->incrementList(1, type, numNodes);
+ if(!outmin || !outmid || !outmax) continue;
double xyz[3][8];
for(int nod = 0; nod < numNodes; nod++)
data1->getNode(0, ent, ele, nod, xyz[0][nod], xyz[1][nod], xyz[2][nod]);
diff --git a/Plugin/Eigenvectors.cpp b/Plugin/Eigenvectors.cpp
index 70ce1b22f26e4692f0653f479925673b00ec9ac2..4172e427e1bcb18128fcf205c1afcec5fb3afe3e 100644
--- a/Plugin/Eigenvectors.cpp
+++ b/Plugin/Eigenvectors.cpp
@@ -75,11 +75,11 @@ PView *GMSH_EigenvectorsPlugin::execute(PView *v)
int numComp = data1->getNumComponents(0, ent, ele);
if(numComp != 9) continue;
int type = data1->getType(0, ent, ele);
- std::vector<double> *outmin = dmin->incrementList(3, type);
- std::vector<double> *outmid = dmid->incrementList(3, type);
- std::vector<double> *outmax = dmax->incrementList(3, type);
- if(!outmin || !outmid || !outmax) continue;
int numNodes = data1->getNumNodes(0, ent, ele);
+ std::vector<double> *outmin = dmin->incrementList(3, type, numNodes);
+ std::vector<double> *outmid = dmid->incrementList(3, type, numNodes);
+ std::vector<double> *outmax = dmax->incrementList(3, type, numNodes);
+ if(!outmin || !outmid || !outmax) continue;
double xyz[3][8];
for(int nod = 0; nod < numNodes; nod++)
data1->getNode(0, ent, ele, nod, xyz[0][nod], xyz[1][nod], xyz[2][nod]);
diff --git a/Plugin/ExtractElements.cpp b/Plugin/ExtractElements.cpp
index a957962534be0f52f7a8fd04f168af677201cba8..44537ae6d9782f864380354c001db63f2f0cc7e5 100644
--- a/Plugin/ExtractElements.cpp
+++ b/Plugin/ExtractElements.cpp
@@ -87,7 +87,7 @@ PView *GMSH_ExtractElementsPlugin::execute(PView *v)
}
int type = data1->getType(step, ent, ele);
int numComp = data1->getNumComponents(step, ent, ele);
- std::vector<double> *out = data2->incrementList(numComp, type);
+ std::vector<double> *out = data2->incrementList(numComp, type, numNodes);
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
std::vector<double> v(numNodes * numComp);
for(int nod = 0; nod < numNodes; nod++){
diff --git a/Plugin/Gradient.cpp b/Plugin/Gradient.cpp
index 98fe67a2985aca278c06d95099032f29be3428b7..a9fe9077a574b58492d73ba3f428f1dedb54d4de 100644
--- a/Plugin/Gradient.cpp
+++ b/Plugin/Gradient.cpp
@@ -60,9 +60,9 @@ PView *GMSH_GradientPlugin::execute(PView *v)
int numComp = data1->getNumComponents(firstNonEmptyStep, ent, ele);
if(numComp != 1 && numComp != 3) continue;
int type = data1->getType(firstNonEmptyStep, ent, ele);
- std::vector<double> *out = data2->incrementList((numComp == 1) ? 3 : 9, type);
- if(!out) continue;
int numNodes = data1->getNumNodes(firstNonEmptyStep, ent, ele);
+ std::vector<double> *out = data2->incrementList((numComp == 1) ? 3 : 9, type, numNodes);
+ if(!out) continue;
double x[8], y[8], z[8], val[8 * 3];
for(int nod = 0; nod < numNodes; nod++)
data1->getNode(firstNonEmptyStep, ent, ele, nod, x[nod], y[nod], z[nod]);
diff --git a/Plugin/HarmonicToTime.cpp b/Plugin/HarmonicToTime.cpp
index 2cfcc6270f3b56d4ecdef3cf7b750c4c5acaf11e..fa5cfb8f736ef7f28e7640c6285d7df378c563c5 100644
--- a/Plugin/HarmonicToTime.cpp
+++ b/Plugin/HarmonicToTime.cpp
@@ -79,7 +79,7 @@ PView *GMSH_HarmonicToTimePlugin::execute(PView * v)
int numNodes = data1->getNumNodes(0, ent, ele);
int type = data1->getType(0, ent, ele);
int numComp = data1->getNumComponents(0, ent, ele);
- std::vector<double> *out = data2->incrementList(numComp, type);
+ std::vector<double> *out = data2->incrementList(numComp, type, numNodes);
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
std::vector<double> vr(numNodes * numComp), vi(numNodes * numComp);
for(int nod = 0; nod < numNodes; nod++){
diff --git a/Plugin/MathEval.cpp b/Plugin/MathEval.cpp
index 492392b58b136bfe1486268713e1a7aa2e3144d9..548c06542a8a97a21984614cb9209fa5c4693050 100644
--- a/Plugin/MathEval.cpp
+++ b/Plugin/MathEval.cpp
@@ -204,7 +204,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
int numComp = data1->getNumComponents(timeBeg, ent, ele);
int otherNumComp = (!otherData || octree) ? 9 :
otherData->getNumComponents(timeBeg, ent, ele);
- std::vector<double> *out = data2->incrementList(numComp2, type);
+ std::vector<double> *out = data2->incrementList(numComp2, type, numNodes);
std::vector<double> w(std::max(9, otherNumComp), 0.);
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
for(int nod = 0; nod < numNodes; nod++)
diff --git a/Plugin/Scal2Vec.cpp b/Plugin/Scal2Vec.cpp
index e7ecf449bbffb8e4b20e57a5deadb659c0014a60..5d68c68d561ca13bdb57e9ec077942051b93750b 100644
--- a/Plugin/Scal2Vec.cpp
+++ b/Plugin/Scal2Vec.cpp
@@ -116,9 +116,9 @@ PView *GMSH_Scal2VecPlugin::execute(PView *v)
if(dataRef->skipElement(0, ent, ele)) continue;
int numComp = 3; // The 3 components of the new view: x,y,z
int type = dataRef->getType(0, ent, ele);
- std::vector<double> *out = dataNew->incrementList(numComp, type); // Pointer in data of the new view
- if(!out) continue;
int numNodes = dataRef->getNumNodes(0, ent, ele);
+ std::vector<double> *out = dataNew->incrementList(numComp, type, numNodes); // Pointer in data of the new view
+ if(!out) continue;
double x[8], y[8], z[8], valX, valY, valZ;
for(int nod = 0; nod < numNodes; nod++)
dataRef->getNode(0, ent, ele, nod, x[nod], y[nod], z[nod]);
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 80e3360757cdbc18f7019aa63bbce11d3b80d1ef..5516811652de4050992e7414e04556bf3c9dd86b 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -15,14 +15,17 @@ PViewDataList::PViewDataList(bool isAdapted)
: PViewData(), NbTimeStep(0), Min(VAL_INF), Max(-VAL_INF),
NbSP(0), NbVP(0), NbTP(0), NbSL(0), NbVL(0), NbTL(0),
NbST(0), NbVT(0), NbTT(0), NbSQ(0), NbVQ(0), NbTQ(0),
+ NbSG(0), NbVG(0), NbTG(0),
NbSS(0), NbVS(0), NbTS(0), NbSH(0), NbVH(0), NbTH(0),
NbSI(0), NbVI(0), NbTI(0), NbSY(0), NbVY(0), NbTY(0),
- NbT2(0), NbT3(0), _lastElement(-1), _lastDimension(-1),
- _lastNumNodes(-1), _lastNumComponents(-1), _lastNumValues(-1),
- _lastNumEdges(-1), _lastType(-1), _lastXYZ(0), _lastVal(0),
- _isAdapted(isAdapted)
+ NbSD(0), NbVD(0), NbTD(0), NbT2(0), NbT3(0),
+ _lastElement(-1), _lastDimension(-1), _lastNumNodes(-1),
+ _lastNumComponents(-1), _lastNumValues(-1), _lastNumEdges(-1),
+ _lastType(-1), _lastXYZ(0), _lastVal(0), _isAdapted(isAdapted)
{
- for(int i = 0; i < 24; i++) _index[i] = 0;
+ for(int i = 0; i < 30; i++) _index[i] = 0;
+ polyTotNumNodes[0] = 0.; polyTotNumNodes[1] = 0.;
+ polyAgNumNodes[0].push_back(0.); polyAgNumNodes[1].push_back(0.);
}
void PViewDataList::setXY(std::vector<double> &x, std::vector<double> &y)
@@ -69,6 +72,10 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
_stat(TI, 9, NbTI, 6, TYPE_PRI);
_stat(SY, 1, NbSY, 5, TYPE_PYR); _stat(VY, 3, NbVY, 5, TYPE_PYR);
_stat(TY, 9, NbTY, 5, TYPE_PYR);
+ _stat(SG, 1, NbSG, 3, TYPE_POLYG); _stat(VG, 3, NbVG, 3, TYPE_POLYG);
+ _stat(TG, 9, NbTG, 3, TYPE_POLYG);
+ _stat(SD, 1, NbSD, 4, TYPE_POLYH); _stat(VD, 3, NbVD, 4, TYPE_POLYH);
+ _stat(TD, 9, NbTD, 4, TYPE_POLYH);
// add dummy time values if none (or too few) time values are
// provided (e.g. using the old parsed format)
@@ -78,10 +85,11 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
}
// compute starting element indices
- int nb[24] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
+ int nb[30] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT,
NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH,
- NbSI, NbVI, NbTI, NbSY, NbVY, NbTY};
- for(int i = 0; i < 24; i++){
+ NbSI, NbVI, NbTI, NbSY, NbVY, NbTY, NbSG, NbVG, NbTG,
+ NbSD, NbVD, NbTD};
+ for(int i = 0; i < 30; i++){
_index[i] = 0;
for(int j = 0; j <= i; j++)
_index[i] += nb[j];
@@ -93,18 +101,18 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
}
int PViewDataList::getNumScalars(int step)
-{
- return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY;
+{
+ return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY + NbSG + NbSD;
}
int PViewDataList::getNumVectors(int step)
{
- return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY;
+ return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVG + NbVD;
}
int PViewDataList::getNumTensors(int step)
{
- return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY;
+ return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTG + NbTD;
}
int PViewDataList::getNumElements(int step, int ent)
@@ -207,7 +215,15 @@ void PViewDataList::_stat(std::vector<double> &list, int nbcomp, int nbelm,
}
int nb = list.size() / nbelm;
- for(unsigned int i = 0; i < list.size(); i += nb){
+ for(unsigned int ele = 0; ele < nbelm; ele ++){
+ int i = ele * nb;
+ if(type == TYPE_POLYG || type == TYPE_POLYH){
+ int t = (type == TYPE_POLYG) ? 0 : 1;
+ nbnod = polyNumNodes[t][ele];
+ nb = list.size() / polyTotNumNodes[t] * nbnod;
+ i = polyAgNumNodes[t][ele] * nb / nbnod;
+ nbval = nbcomp * nbnod;
+ }
int N = nb - 3 * nbnod;
double *X = &list[i];
double *Y = &list[i + 1 * nbnod];
@@ -261,10 +277,16 @@ void PViewDataList::_setLast(int ele, int dim, int nbnod, int nbcomp, int nbedg,
_lastNumComponents = nbcomp;
_lastNumEdges = nbedg;
_lastType = type;
- int nb = list.size() / nblist;
- _lastXYZ = &list[ele * nb];
- _lastVal = &list[ele * nb + 3 * _lastNumNodes];
+ int nb = list.size() / nblist; // number of coords and values for the element
+ int nbAg = ele * nb; // number of coords and values before the ones of the element
+ if(type == TYPE_POLYG || type == TYPE_POLYH){
+ int t = (type == TYPE_POLYG) ? 0 : 1;
+ nb = list.size() / polyTotNumNodes[t] * nbnod;
+ nbAg = polyAgNumNodes[t][ele] * nb / nbnod;
+ }
_lastNumValues = (nb - 3 * nbnod) / NbTimeStep;
+ _lastXYZ = &list[nbAg];
+ _lastVal = &list[nbAg + 3 * _lastNumNodes];
}
void PViewDataList::_setLast(int ele)
@@ -305,11 +327,23 @@ void PViewDataList::_setLast(int ele)
else if(ele < _index[19]) _setLast(ele - _index[18], 3, 6, 3, 9, TYPE_PRI, VI, NbVI);
else _setLast(ele - _index[19], 3, 6, 9, 9, TYPE_PRI, TI, NbTI);
}
- else{ // pyramids
+ else if(ele < _index[23]){ // pyramids
if(ele < _index[21]) _setLast(ele - _index[20], 3, 5, 1, 8, TYPE_PYR, SY, NbSY);
else if(ele < _index[22]) _setLast(ele - _index[21], 3, 5, 3, 8, TYPE_PYR, VY, NbVY);
else _setLast(ele - _index[22], 3, 5, 9, 8, TYPE_PYR, TY, NbTY);
}
+ else if(ele < _index[26]){ // polygons
+ int nN = polyNumNodes[0][ele - _index[23]];
+ if(ele < _index[24]) _setLast(ele - _index[23], 2, nN, 1, nN, TYPE_POLYG, SG, NbSG);
+ else if(ele < _index[25]) _setLast(ele - _index[24], 2, nN, 3, nN, TYPE_POLYG, VG, NbVG);
+ else _setLast(ele - _index[25], 2, nN, 9, nN, TYPE_POLYG, TG, NbTG);
+ }
+ else if(ele < _index[29]){ // polyhedra
+ int nN = polyNumNodes[1][ele - _index[26]];
+ if(ele < _index[27]) _setLast(ele - _index[26], 3, nN, 1, nN*1.5, TYPE_POLYH, SD, NbSD);
+ else if(ele < _index[28]) _setLast(ele - _index[27], 3, nN, 3, nN*1.5, TYPE_POLYH, VD, NbVD);
+ else _setLast(ele - _index[28], 3, nN, 9, nN*1.5, TYPE_POLYH, TD, NbTD);
+ }
}
int PViewDataList::getDimension(int step, int ent, int ele)
@@ -794,6 +828,9 @@ int PViewDataList::_getRawData(int idxtype, std::vector<double> **l, int **ne,
int *nc, int *nn)
{
int type = 0;
+ // No constant nn for polygons!
+ if(idxtype > 23 && idxtype < 30)
+ Msg::Warning("No constant number of nodes for polygons and polyhedra.");
switch(idxtype){
case 0 : *l = &SP; *ne = &NbSP; *nc = 1; *nn = 1; type = TYPE_PNT; break;
case 1 : *l = &VP; *ne = &NbVP; *nc = 3; *nn = 1; type = TYPE_PNT; break;
@@ -819,6 +856,12 @@ int PViewDataList::_getRawData(int idxtype, std::vector<double> **l, int **ne,
case 21: *l = &SY; *ne = &NbSY; *nc = 1; *nn = 5; type = TYPE_PYR; break;
case 22: *l = &VY; *ne = &NbVY; *nc = 3; *nn = 5; type = TYPE_PYR; break;
case 23: *l = &TY; *ne = &NbTY; *nc = 9; *nn = 5; type = TYPE_PYR; break;
+ case 24: *l = &SG; *ne = &NbSG; *nc = 1; *nn = 3; type = TYPE_POLYG; break;
+ case 25: *l = &VG; *ne = &NbVG; *nc = 3; *nn = 3; type = TYPE_POLYG; break;
+ case 26: *l = &TG; *ne = &NbTG; *nc = 9; *nn = 3; type = TYPE_POLYG; break;
+ case 27: *l = &SD; *ne = &NbSD; *nc = 1; *nn = 4; type = TYPE_POLYH; break;
+ case 28: *l = &VD; *ne = &NbVD; *nc = 3; *nn = 4; type = TYPE_POLYH; break;
+ case 29: *l = &TD; *ne = &NbTD; *nc = 9; *nn = 4; type = TYPE_POLYH; break;
default: Msg::Error("Wrong type in PViewDataList"); break;
}
@@ -851,8 +894,9 @@ void PViewDataList::setOrder2(int type)
fs->coefficients, fs->monomials);
}
-std::vector<double> *PViewDataList::incrementList(int numComp, int type)
+std::vector<double> *PViewDataList::incrementList(int numComp, int type, int numNodes)
{
+ int nb;
switch(type){
case TYPE_PNT:
if (numComp == 1){ NbSP++; return &SP; }
@@ -894,6 +938,24 @@ std::vector<double> *PViewDataList::incrementList(int numComp, int type)
else if(numComp == 3){ NbVY++; return &VY; }
else if(numComp == 9){ NbTY++; return &TY; }
break;
+ case TYPE_POLYG:
+ polyNumNodes[0].push_back(numNodes);
+ nb = (polyAgNumNodes[0].size()) ? polyAgNumNodes[0].back() : 0;
+ polyAgNumNodes[0].push_back(numNodes + nb);
+ polyTotNumNodes[0] += numNodes;
+ if (numComp == 1){ NbSG++; return &SG; }
+ else if(numComp == 3){ NbVG++; return &VG; }
+ else if(numComp == 9){ NbTG++; return &TG; }
+ break;
+ case TYPE_POLYH:
+ polyNumNodes[1].push_back(numNodes);
+ nb = (polyAgNumNodes[1].size()) ? polyAgNumNodes[1].back() : 0;
+ polyAgNumNodes[1].push_back(numNodes + nb);
+ polyTotNumNodes[1] += numNodes;
+ if (numComp == 1){ NbSD++; return &SD; }
+ else if(numComp == 3){ NbVD++; return &VD; }
+ else if(numComp == 9){ NbTD++; return &TD; }
+ break;
}
return 0;
}
diff --git a/Post/PViewDataList.h b/Post/PViewDataList.h
index 593285819f62152bcef901aadcabe7df7727c225..a15b4dfe266d9796923e39d33d49beef3541e520 100644
--- a/Post/PViewDataList.h
+++ b/Post/PViewDataList.h
@@ -31,6 +31,8 @@ class PViewDataList : public PViewData {
std::vector<double> ST, VT, TT; // triangles
int NbSQ, NbVQ, NbTQ;
std::vector<double> SQ, VQ, TQ; // quadrangles
+ int NbSG, NbVG, NbTG;
+ std::vector<double> SG, VG, TG; // polygons
int NbSS, NbVS, NbTS;
std::vector<double> SS, VS, TS; // tetrahedra
int NbSH, NbVH, NbTH;
@@ -39,11 +41,16 @@ class PViewDataList : public PViewData {
std::vector<double> SI, VI, TI; // prisms
int NbSY, NbVY, NbTY;
std::vector<double> SY, VY, TY; // pyramids
+ int NbSD, NbVD, NbTD;
+ std::vector<double> SD, VD, TD; // polyhedra
int NbT2, NbT3;
std::vector<double> T2D, T3D; // 2D and 3D text strings
std::vector<char> T2C, T3C;
+ std::vector<int> polyNumNodes[2];
+ std::vector<int> polyAgNumNodes[2];
+ int polyTotNumNodes[2];
private:
- int _index[24];
+ int _index[30];
int _lastElement, _lastDimension;
int _lastNumNodes, _lastNumComponents, _lastNumValues, _lastNumEdges, _lastType;
double *_lastXYZ, *_lastVal;
@@ -78,10 +85,12 @@ class PViewDataList : public PViewData {
int getNumLines(int step=-1){ return NbSL + NbVL + NbTL; }
int getNumTriangles(int step=-1){ return NbST + NbVT + NbTT; }
int getNumQuadrangles(int step=-1){ return NbSQ + NbVQ + NbTQ; }
+ int getNumPolygons(int step=-1){ return NbSG + NbVG + NbTG; }
int getNumTetrahedra(int step=-1){ return NbSS + NbVS + NbTS; }
int getNumHexahedra(int step=-1){ return NbSH + NbVH + NbTH; }
int getNumPrisms(int step=-1){ return NbSI + NbVI + NbTI; }
int getNumPyramids(int step=-1){ return NbSY + NbVY + NbTY; }
+ int getNumPolyhedra(int step=-1){ return NbSD + NbVD + NbTD; }
int getNumEntities(int step=-1){ return 1; }
int getNumElements(int step=-1, int ent=-1);
int getDimension(int step, int ent, int ele);
@@ -109,7 +118,7 @@ class PViewDataList : public PViewData {
// specific to list-based data sets
void setOrder2(int type);
- std::vector<double> *incrementList(int numComp, int type);
+ std::vector<double> *incrementList(int numComp, int type, int numNodes = 0);
// I/O routines
bool readPOS(FILE *fp, double version, bool binary);