From aa98966ceb7e1785d49188a5d6fabd8a23853129 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 14 Dec 2008 12:43:08 +0000
Subject: [PATCH] rewrote adaptive post in terms of local adapt() operator
---
Mesh/Makefile | 59 +++--
Plugin/Makefile | 8 +-
Post/adaptiveData.cpp | 513 ++++++++++++++++++++----------------------
Post/adaptiveData.h | 63 +++---
4 files changed, 301 insertions(+), 342 deletions(-)
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 15a3dec75e..8d1e113c4f 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -182,8 +182,8 @@ meshGFace${OBJEXT}: meshGFace.cpp meshGFace.h meshGFaceBDS.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h ../Numeric/Numeric.h \
- ../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h qualityMeasures.h \
- Field.h ../Common/OS.h HighOrder.h
+ ../Numeric/NumericEmbedded.h BDS.h qualityMeasures.h Field.h \
+ ../Post/PView.h ../Common/OS.h HighOrder.h
meshGFaceTransfinite${OBJEXT}: meshGFaceTransfinite.cpp meshGFace.h \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -227,8 +227,8 @@ meshGFaceBDS${OBJEXT}: meshGFaceBDS.cpp meshGFace.h meshGFaceOptimize.h \
../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h ../Geo/GModel.h \
../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h \
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h \
- qualityMeasures.h Field.h ../Common/OS.h
+ ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h \
+ qualityMeasures.h Field.h ../Post/PView.h ../Common/OS.h
meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
BackgroundMesh.h meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
@@ -256,18 +256,18 @@ meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
../Geo/Pair.h BackgroundMesh.h Generator.h
meshGFaceQuadrilateralize${OBJEXT}: meshGFaceQuadrilateralize.cpp \
meshGFaceQuadrilateralize.h ../Common/GmshMessage.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
- meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
- ../Numeric/GmshMatrix.h meshGFaceOptimize.h meshGFaceBDS.h BDS.h \
- ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Geo/GFace.h \
+ ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/Pair.h ../Post/PView.h
+ ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
+ ../Geo/SVector3.h ../Geo/Pair.h meshGFaceDelaunayInsertion.h \
+ ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h meshGFaceOptimize.h \
+ meshGFaceBDS.h BDS.h
meshGRegion${OBJEXT}: meshGRegion.cpp meshGRegion.h \
meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
@@ -288,8 +288,8 @@ meshGRegion${OBJEXT}: meshGRegion.cpp meshGRegion.h \
../Geo/SPoint3.h ../Geo/SVector3.h ../Geo/SBoundingBox3d.h \
../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
../Common/ListUtils.h ../Geo/SPoint2.h ../Geo/ExtrudeParams.h \
- ../Common/SmoothData.h ../Geo/GRegion.h BDS.h ../Post/PView.h \
- ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
+ ../Common/SmoothData.h ../Geo/GRegion.h BDS.h ../Common/Context.h \
+ ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
meshGRegionDelaunayInsertion${OBJEXT}: meshGRegionDelaunayInsertion.cpp \
../Common/OS.h BackgroundMesh.h meshGRegion.h meshGRegionLocalMeshMod.h \
meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
@@ -374,17 +374,12 @@ BackgroundMesh${OBJEXT}: BackgroundMesh.cpp ../Common/GmshMessage.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h Field.h ../Post/PView.h
qualityMeasures${OBJEXT}: qualityMeasures.cpp qualityMeasures.h BDS.h \
- ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
- ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
- ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
- ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
- ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h
+ ../Common/GmshMessage.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MElement.h ../Common/GmshDefines.h \
+ ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h
BoundaryLayers${OBJEXT}: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -400,12 +395,12 @@ BoundaryLayers${OBJEXT}: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Numeric/GmshMatrix.h BoundaryLayers.h ../Geo/ExtrudeParams.h \
../Common/SmoothData.h meshGEdge.h meshGFace.h
BDS${OBJEXT}: BDS.cpp ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h \
- ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
- ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
- ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
- ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
+ ../Common/GmshMessage.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+ ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h \
+ ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+ ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
diff --git a/Plugin/Makefile b/Plugin/Makefile
index b540066178..7c6019118b 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -216,13 +216,7 @@ ExtractEdges${OBJEXT}: ExtractEdges.cpp ExtractEdges.h Plugin.h \
../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
- ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Mesh/BDS.h \
- ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
- ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/GEntity.h \
- ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/SPoint3.h \
- ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/Pair.h
+ ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Mesh/BDS.h
MakeSimplex${OBJEXT}: MakeSimplex.cpp MakeSimplex.h Plugin.h ../Common/Options.h \
../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index 8df5037f3d..78d39946f7 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -8,7 +8,6 @@
#include <set>
#include "Plugin.h"
#include "adaptiveData.h"
-#include "OS.h"
std::set<adaptivePoint> adaptiveLine::allPoints;
std::set<adaptivePoint> adaptiveTriangle::allPoints;
@@ -832,323 +831,288 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
}
template <class T>
-int adaptiveElements<T>::_zoomElement(int ielem, int level, GMSH_Post_Plugin *plug)
+adaptiveElements<T>::~adaptiveElements()
{
- const int N = _coeffs->size1();
- Double_Vector val(N), res(T::allPoints.size());
- Double_Vector valx(N), resx(T::allPoints.size());
- Double_Vector valy(N), resy(T::allPoints.size());
- Double_Vector valz(N), resz(T::allPoints.size());
- Double_Matrix xyz(_posX->size2(), 3);
- Double_Matrix XYZ(T::allPoints.size(), 3);
+ if(_interpolVal) delete _interpolVal;
+ if(_interpolGeom) delete _interpolGeom;
+ cleanElement<T>();
+}
- for(int k = 0; k < _posX->size2(); ++k){
- xyz(k, 0) = (*_posX)(ielem, k);
- xyz(k, 1) = (*_posY)(ielem, k);
- xyz(k, 2) = (*_posZ)(ielem, k);
+template <class T>
+void adaptiveElements<T>::init(int level)
+{
+ T::create(level);
+ int numVals = _coeffsVal->size1();
+ int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
+
+ if(_interpolVal) delete _interpolVal;
+ _interpolVal = new Double_Matrix(T::allPoints.size(), numVals);
+
+ if(_interpolGeom) delete _interpolGeom;
+ _interpolGeom = new Double_Matrix(T::allPoints.size(), numNodes);
+
+ double sf[100];
+ int kk = 0;
+ for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
+ it != T::allPoints.end(); ++it) {
+ adaptivePoint *p = (adaptivePoint*)&(*it);
+
+ computeShapeFunctions(_coeffsVal, _eexpsVal, p->x, p->y, p->z, sf);
+ for(int k = 0; k < numVals; k++)
+ (*_interpolVal)(kk, k) = sf[k];
+
+ if(_coeffsGeom)
+ computeShapeFunctions(_coeffsGeom, _eexpsGeom, p->x, p->y, p->z, sf);
+ else
+ T::GSF(p->x, p->y, p->z, sf);
+ for(int k = 0; k < numNodes; k++)
+ (*_interpolGeom)(kk, k) = sf[k];
+
+ kk++;
}
+}
- for(int k = 0; k < N; ++k)
- val(k) = (*_val)(ielem, k);
-
- _interpolate->mult(val, res);
-
- if(_valX){
- for(int k = 0; k < N; ++k){
- valx(k) = (*_valX)(ielem, k);
- valy(k) = (*_valY)(ielem, k);
- valz(k) = (*_valZ)(ielem, k);
- }
- _interpolate->mult(valx, resx);
- _interpolate->mult(valy, resy);
- _interpolate->mult(valz, resz);
+template <class T>
+void adaptiveElements<T>::adapt(double tol, int numComp,
+ std::vector<PCoords> &coords,
+ std::vector<PValues> &values,
+ double &minVal, double &maxVal,
+ GMSH_Post_Plugin *plug,
+ bool onlyComputeMinMax)
+{
+ if(numComp != 1 && numComp != 3){
+ Msg::Error("Can only adapt scalar or vector data");
+ return;
}
-
- _geometry->mult(xyz, XYZ);
-
+
+ int numPoints = T::allPoints.size();
+ if(!numPoints){
+ Msg::Error("No adapted points to interpolate");
+ return;
+ }
+
+ int numVals = _coeffsVal->size1();
+ if(numVals != values.size()){
+ Msg::Error("Wrong number of values in adaptation %d != %i",
+ numVals, values.size());
+ return;
+ }
+
+ Double_Vector val(numVals), res(numPoints);
+ if(numComp == 1){
+ for(int k = 0; k < numVals; ++k)
+ val(k) = values[k].v[0];
+ }
+ else{
+ for(int k = 0; k < numVals; ++k)
+ val(k) = values[k].v[0] * values[k].v[0] + values[k].v[1] * values[k].v[1] +
+ values[k].v[2] * values[k].v[2];
+ }
+ _interpolVal->mult(val, res);
+
+ for(int i = 0; i < numPoints; i++){
+ minVal = std::min(minVal, res(i));
+ maxVal = std::max(maxVal, res(i));
+ }
+ if(onlyComputeMinMax) return;
+
+ Double_Vector *resx = 0, *resy = 0, *resz = 0;
+ if(numComp == 3){
+ Double_Vector valx(numVals), valy(numVals), valz(numVals);
+ resx = new Double_Vector(numPoints);
+ resy = new Double_Vector(numPoints);
+ resz = new Double_Vector(numPoints);
+ for(int k = 0; k < numVals; ++k){
+ valx(k) = values[k].v[0];
+ valy(k) = values[k].v[1];
+ valz(k) = values[k].v[2];
+ }
+ _interpolVal->mult(valx, *resx);
+ _interpolVal->mult(valy, *resy);
+ _interpolVal->mult(valz, *resz);
+ }
+
+ int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
+ if(numNodes != coords.size()){
+ Msg::Error("Wrong number of nodes in adaptation %d != %i",
+ numNodes, coords.size());
+ return;
+ }
+
+ Double_Matrix xyz(numNodes, 3), XYZ(numPoints, 3);
+ for(int k = 0; k < numNodes; ++k){
+ xyz(k, 0) = coords[k].c[0];
+ xyz(k, 1) = coords[k].c[1];
+ xyz(k, 2) = coords[k].c[2];
+ }
+ _interpolGeom->mult(xyz, XYZ);
+
int k = 0;
for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
it != T::allPoints.end(); ++it){
adaptivePoint *p = (adaptivePoint*)&(*it);
p->val = res(k);
- if(_valX){
- p->valx = resx(k);
- p->valy = resy(k);
- p->valz = resz(k);
+ if(resx){
+ p->valx = (*resx)(k);
+ p->valy = (*resy)(k);
+ p->valz = (*resz)(k);
}
p->X = XYZ(k, 0);
p->Y = XYZ(k, 1);
p->Z = XYZ(k, 2);
- if(_minVal > p->val) _minVal = p->val;
- if(_maxVal < p->val) _maxVal = p->val;
k++;
}
-
+
+ if(resx){
+ delete resx; delete resy; delete resz;
+ }
+
for(typename std::list<T*>::iterator it = T::all.begin();
- it != T::all.end(); it++)
+ it != T::all.end(); it++)
(*it)->visible = false;
-
- if(!plug || _tol != 0.)
- T::error(_maxVal - _minVal, _tol);
-
+
+ if(!plug || tol != 0.)
+ T::error(maxVal - minVal, tol);
+
if(plug)
plug->assignSpecificVisibility();
-
- for(typename std::list<T*>::iterator it = T::all.begin();
- it != T::all.end(); it++){
- if((*it)->visible && !(*it)->e[0] && level != _level)
- return 0;
- }
- for(typename std::list<T*>::iterator it = T::all.begin();
+ coords.clear();
+ values.clear();
+ for(typename std::list<T*>::iterator it = T::all.begin();
it != T::all.end(); it++){
if((*it)->visible){
adaptivePoint **p = (*it)->p;
- for(int k = 0; k < T::numNodes; ++k) List_Add(_listEle, &p[k]->X);
- for(int k = 0; k < T::numNodes; ++k) List_Add(_listEle, &p[k]->Y);
- for(int k = 0; k < T::numNodes; ++k) List_Add(_listEle, &p[k]->Z);
- if(_valX){
- for(int k = 0; k < T::numNodes; ++k){
- List_Add(_listEle, &p[k]->valx);
- List_Add(_listEle, &p[k]->valy);
- List_Add(_listEle, &p[k]->valz);
- }
- }
- else{
- for (int k = 0; k < T::numNodes; ++k)
- List_Add(_listEle, &p[k]->val);
+ for(int k = 0; k < T::numNodes; ++k) {
+ coords.push_back(PCoords(p[k]->X, p[k]->Y, p[k]->Z));
+ if(numComp == 1)
+ values.push_back(PValues(p[k]->val));
+ else
+ values.push_back(PValues(p[k]->valx, p[k]->valy, p[k]->valz));
}
- (*_numEle)++;
}
}
- return 1;
-}
-
-template <class T>
-void adaptiveElements<T>::_changeResolution(int level, GMSH_Post_Plugin *plug, int *done)
-{
- T::create(level);
-
- if(!T::allPoints.size()) return;
-
- const int N = _coeffs->size1();
- if(_interpolate) delete _interpolate;
- _interpolate = new Double_Matrix(T::allPoints.size(), N);
-
- if(_geometry) delete _geometry;
- _geometry = new Double_Matrix(T::allPoints.size(), _posX->size2());
-
- double sf[100];
- int kk = 0;
- for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
- it != T::allPoints.end(); ++it) {
- adaptivePoint *p = (adaptivePoint*)&(*it);
-
- computeShapeFunctions(_coeffs, _eexps, p->x, p->y, p->z, sf);
- for(int k = 0; k < N; ++k)
- (*_interpolate)(kk, k) = sf[k];
-
- if(_coeffsGeom)
- computeShapeFunctions(_coeffsGeom, _eexpsGeom, p->x, p->y, p->z, sf);
- else
- T::GSF(p->x, p->y, p->z, sf);
- for(int k = 0; k < _posX->size2(); k++)
- (*_geometry)(kk, k) = sf[k];
-
- kk++;
- }
-
- const int nbelm = _posX->size1();
- for(int i = 0; i < nbelm; ++i)
- done[i] = _zoomElement(i, level, plug);
-}
-
-template <class T>
-adaptiveElements<T>::~adaptiveElements()
-{
- if(_posX) delete _posX;
- if(_posY) delete _posY;
- if(_posZ) delete _posZ;
- if(_val) delete _val;
- if(_valX) delete _valX;
- if(_valY) delete _valY;
- if(_valZ) delete _valZ;
- if(_interpolate) delete _interpolate;
- if(_geometry) delete _geometry;
- cleanElement<T>();
}
template <class T>
-void adaptiveElements<T>::initData(PViewData *data, int step)
+void adaptiveElements<T>::addInView(double tol, int step,
+ PViewData *in, PViewDataList *out,
+ GMSH_Post_Plugin *plug)
{
- int numComp = data->getNumComponents(0, 0, 0);
+ int numComp = in->getNumComponents(0, 0, 0);
if(numComp != 1 && numComp != 3) return;
-
- int numEle = 0;
+
+ int numEle = 0, *outNb = 0;
+ List_T *outList = 0;
switch(T::numEdges){
- case 1: numEle = data->getNumLines(); break;
- case 3: numEle = data->getNumTriangles(); break;
- case 4: numEle = data->getNumQuadrangles(); break;
- case 6: numEle = data->getNumTetrahedra(); break;
- case 9: numEle = data->getNumPrisms(); break;
- case 12: numEle = data->getNumHexahedra(); break;
+ case 1:
+ numEle = in->getNumLines();
+ outNb = (numComp == 1) ? &out->NbSL : &out->NbVL;
+ outList = (numComp == 1) ? out->SL : out->VL;
+ break;
+ case 3:
+ numEle = in->getNumTriangles();
+ outNb = (numComp == 1) ? &out->NbST : &out->NbVT;
+ outList = (numComp == 1) ? out->ST : out->VT;
+ break;
+ case 4:
+ numEle = in->getNumQuadrangles();
+ outNb = (numComp == 1) ? &out->NbSQ : &out->NbVQ;
+ outList = (numComp == 1) ? out->SQ : out->VQ;
+ break;
+ case 6:
+ numEle = in->getNumTetrahedra();
+ outNb = (numComp == 1) ? &out->NbSS : &out->NbVS;
+ outList = (numComp == 1) ? out->SS : out->VS;
+ break;
+ case 9:
+ numEle = in->getNumPrisms();
+ outNb = (numComp == 1) ? &out->NbSI : &out->NbVI;
+ outList = (numComp == 1) ? out->SI : out->VI;
+ break;
+ case 12:
+ numEle = in->getNumHexahedra();
+ outNb = (numComp == 1) ? &out->NbSH : &out->NbVH;
+ outList = (numComp == 1) ? out->SH : out->VH;
+ break;
}
if(!numEle) return;
-
- int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
- int numVal = _coeffs->size1() * numComp;
- if(!numNodes || !numVal) return;
-
- _minVal = VAL_INF;
- _maxVal = -VAL_INF;
-
- if(_posX) delete _posX;
- if(_posY) delete _posY;
- if(_posZ) delete _posZ;
- if(_val) delete _val;
- if(_valX) delete _valX;
- if(_valY) delete _valY;
- if(_valZ) delete _valZ;
- _posX = new Double_Matrix(numEle, numNodes);
- _posY = new Double_Matrix(numEle, numNodes);
- _posZ = new Double_Matrix(numEle, numNodes);
- _val = new Double_Matrix(numEle, numVal);
- if(numComp == 3){
- _valX = new Double_Matrix(numEle, numVal);
- _valY = new Double_Matrix(numEle, numVal);
- _valZ = new Double_Matrix(numEle, numVal);
- }
-
- // store non-interpolated data
+
+ List_Reset(outList);
+ *outNb = 0;
+
int k = 0;
- for(int ent = 0; ent < data->getNumEntities(step); ent++){
- for(int ele = 0; ele < data->getNumElements(step, ent); ele++){
- if(data->skipElement(step, ent, ele) ||
- data->getNumEdges(step, ent, ele) != T::numEdges) continue;
- if(numVal != data->getNumValues(step, ent, ele)){
- Msg::Error("Wrong number of values (%d) in element %d (expecting %d)",
- numVal, ele, data->getNumValues(step, ent, ele));
- continue;
- }
- if(numNodes != data->getNumNodes(step, ent, ele)){
- Msg::Error("Wrong number of nodes (%d) in element %d (expecting %d)",
- numNodes, ele, data->getNumNodes(step, ent, ele));
- continue;
- }
- for(int nod = 0; nod < numNodes; nod++){
- double x, y, z;
- data->getNode(step, ent, ele, nod, x, y, z);
- (*_posX)(k, nod) = x;
- (*_posY)(k, nod) = y;
- (*_posZ)(k, nod) = z;
+ for(int ent = 0; ent < in->getNumEntities(step); ent++){
+ for(int ele = 0; ele < in->getNumElements(step, ent); ele++){
+ if(in->skipElement(step, ent, ele) ||
+ in->getNumEdges(step, ent, ele) != T::numEdges) continue;
+ int numNodes = in->getNumNodes(step, ent, ele);
+ std::vector<PCoords> coords;
+ for(int i = 0; i < numNodes; i++){
+ double x, y, z;
+ in->getNode(step, ent, ele, i, x, y, z);
+ coords.push_back(PCoords(x, y, z));
}
+ int numVal = in->getNumValues(step, ent, ele);
+ std::vector<PValues> values;
if(numComp == 1){
- for(int i = 0; i < numVal; i++){
- double val;
- data->getValue(step, ent, ele, i, val);
- (*_val)(k, i) = val;
- }
+ for(int i = 0; i < numVal; i++){
+ double val;
+ in->getValue(step, ent, ele, i, val);
+ values.push_back(PValues(val));
+ }
}
else if(numComp == 3){
- for(int i = 0; i < numVal / 3; i++){
- double val[3];
- // adaptation of the visualization mesh is based on the norm
- // squared of the vector
- data->getValue(step, ent, ele, 3 * i, val[0]);
- data->getValue(step, ent, ele, 3 * i + 1, val[1]);
- data->getValue(step, ent, ele, 3 * i + 2, val[2]);
- (*_val)(k, i) = (val[0] * val[0] + val[1] * val[1] + val[2] * val[2]);
- (*_valX)(k, i) = val[0];
- (*_valY)(k, i) = val[1];
- (*_valZ)(k, i) = val[2];
- }
+ for(int i = 0; i < numVal / 3; i++){
+ double vx, vy, vz;
+ in->getValue(step, ent, ele, 3 * i, vx);
+ in->getValue(step, ent, ele, 3 * i + 1, vy);
+ in->getValue(step, ent, ele, 3 * i + 2, vz);
+ values.push_back(PValues(vx, vy, vz));
+ }
+ }
+ adapt(tol, numComp, coords, values, out->Min, out->Max, plug);
+ *outNb += coords.size() / T::numNodes;
+ for(unsigned int i = 0; i < coords.size() / T::numNodes; i++){
+ for(int k = 0; k < T::numNodes; ++k)
+ List_Add(outList, &coords[T::numNodes * i + k].c[0]);
+ for(int k = 0; k < T::numNodes; ++k)
+ List_Add(outList, &coords[T::numNodes * i + k].c[1]);
+ for(int k = 0; k < T::numNodes; ++k)
+ List_Add(outList, &coords[T::numNodes * i + k].c[2]);
+ for(int k = 0; k < T::numNodes; ++k)
+ for(int l = 0; l < numComp; ++l)
+ List_Add(outList, &values[T::numNodes * i + k].v[l]);
}
- k++;
}
}
}
-template <class T>
-void adaptiveElements<T>::changeResolution(int level, double tol, GMSH_Post_Plugin *plug)
-{
- if(!_val){
- Msg::Error("No data available to change adaptive resolution");
- return;
- }
-
- _level = level;
- _tol = tol;
-
- List_Reset(_listEle);
- *_numEle = 0;
-
- std::vector<int> done(_posX->size1(), 0);
- // We first do the adaptive stuff at level 2 and will only process
- // elements that have reached the maximal recursion level
- int level_act = (level > 2) ? 2 : level;
- while(1){
- _changeResolution(level_act, plug, &done[0]);
- int numDone = 0;
- for(int i = 0; i < _posX->size1(); ++i) numDone += done[i];
- if(numDone == _posX->size1()) break;
- if(level_act >= level) break;
- level_act++;
- }
- //_changeResolution(level, plug, &done[0]);
-}
-
adaptiveData::adaptiveData(PViewData *data)
: _step(-1), _level(-1), _tol(-1.), _inData(data),
_lines(0), _triangles(0), _quadrangles(0),
_tetrahedra(0), _hexahedra(0), _prisms(0)
{
- // We could do this, but it's a bit tricky (need to set a flag in
- // the view to say "don't use the adaptive stuff anymore!")
- /*
- if(dynamic_cast<PViewDataList*>(_inData) && _inData->getNumTimeSteps() == 1)
- _outData = (PViewDataList*)_inData;
- else
- */
_outData = new PViewDataList(true);
-
- int numComp = _inData->getNumComponents(0, 0, 0);
std::vector<Double_Matrix*> p;
- if(_inData->getNumLines() && _inData->getInterpolationMatrices(1, p) >= 2){
+ if(_inData->getNumLines() && _inData->getInterpolationMatrices(1, p) >= 2)
_lines = new adaptiveElements<adaptiveLine>
- ((numComp == 1) ? _outData->SL : _outData->VL,
- (numComp == 1) ? &_outData->NbSL : &_outData->NbVL,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
- if(_inData->getNumTriangles() && _inData->getInterpolationMatrices(3, p) >= 2){
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumTriangles() && _inData->getInterpolationMatrices(3, p) >= 2)
_triangles = new adaptiveElements<adaptiveTriangle>
- ((numComp == 1) ? _outData->ST : _outData->VT,
- (numComp == 1) ? &_outData->NbST : &_outData->NbVT,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
- if(_inData->getNumQuadrangles() && _inData->getInterpolationMatrices(4, p) >= 2){
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumQuadrangles() && _inData->getInterpolationMatrices(4, p) >= 2)
_quadrangles = new adaptiveElements<adaptiveQuadrangle>
- ((numComp == 1) ? _outData->SQ : _outData->VQ,
- (numComp == 1) ? &_outData->NbSQ : &_outData->NbVQ,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
- if(_inData->getNumTetrahedra() && _inData->getInterpolationMatrices(6, p) >= 2){
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumTetrahedra() && _inData->getInterpolationMatrices(6, p) >= 2)
_tetrahedra = new adaptiveElements<adaptiveTetrahedron>
- ((numComp == 1) ? _outData->SS : _outData->VS,
- (numComp == 1) ? &_outData->NbSS : &_outData->NbVS,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
- if(_inData->getNumPrisms() && _inData->getInterpolationMatrices(9, p) >= 2){
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumPrisms() && _inData->getInterpolationMatrices(9, p) >= 2)
_prisms = new adaptiveElements<adaptivePrism>
- ((numComp == 1) ? _outData->SI : _outData->VI,
- (numComp == 1) ? &_outData->NbSI : &_outData->NbVI,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
- if(_inData->getNumHexahedra() && _inData->getInterpolationMatrices(12, p) >= 2){
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumHexahedra() && _inData->getInterpolationMatrices(12, p) >= 2)
_hexahedra = new adaptiveElements<adaptiveHexahedron>
- ((numComp == 1) ? _outData->SH : _outData->VH,
- (numComp == 1) ? &_outData->NbSH : &_outData->NbVH,
- p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
- }
+ (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
}
adaptiveData::~adaptiveData()
@@ -1159,31 +1123,30 @@ adaptiveData::~adaptiveData()
if(_tetrahedra) delete _tetrahedra;
if(_prisms) delete _prisms;
if(_hexahedra) delete _hexahedra;
- if(_inData != _outData) delete _outData;
+ delete _outData;
}
void adaptiveData::changeResolution(int step, int level, double tol, GMSH_Post_Plugin *plug)
{
- if(_step != step){
- if(_lines) _lines->initData(_inData, step);
- if(_triangles) _triangles->initData(_inData, step);
- if(_quadrangles) _quadrangles->initData(_inData, step);
- if(_tetrahedra) _tetrahedra->initData(_inData, step);
- if(_prisms) _prisms->initData(_inData, step);
- if(_hexahedra) _hexahedra->initData(_inData, step);
+ if(_level != level){
+ if(_lines) _lines->init(level);
+ if(_triangles) _triangles->init(level);
+ if(_quadrangles) _quadrangles->init(level);
+ if(_tetrahedra) _tetrahedra->init(level);
+ if(_prisms) _prisms->init(level);
+ if(_hexahedra) _hexahedra->init(level);
}
if(plug || _step != step || _level != level || _tol != tol){
_outData->setDirty(true);
- if(_lines) _lines->changeResolution(level, tol, plug);
- if(_triangles) _triangles->changeResolution(level, tol, plug);
- if(_quadrangles) _quadrangles->changeResolution(level, tol, plug);
- if(_tetrahedra) _tetrahedra->changeResolution(level, tol, plug);
- if(_prisms) _prisms->changeResolution(level, tol, plug);
- if(_hexahedra) _hexahedra->changeResolution(level, tol, plug);
+ if(_lines) _lines->addInView(tol, step, _inData, _outData, plug);
+ if(_triangles) _triangles->addInView(tol, step, _inData, _outData, plug);
+ if(_quadrangles) _quadrangles->addInView(tol, step, _inData, _outData, plug);
+ if(_tetrahedra) _tetrahedra->addInView(tol, step, _inData, _outData, plug);
+ if(_prisms) _prisms->addInView(tol, step, _inData, _outData, plug);
+ if(_hexahedra) _hexahedra->addInView(tol, step, _inData, _outData, plug);
_outData->finalize();
}
_step = step;
_level = level;
_tol = tol;
}
-
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 200339a221..8477bdb8bb 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -260,39 +260,46 @@ class adaptiveHexahedron {
static void recurError(adaptiveHexahedron *h, double AVG, double tol);
};
+class PCoords {
+ public:
+ double c[3];
+ PCoords(double x, double y, double z)
+ {
+ c[0] = x; c[1] = y; c[2] = z;
+ }
+};
+
+class PValues{
+ public:
+ double v[3];
+ PValues(double vx)
+ {
+ v[0] = vx;
+ }
+ PValues(double vx, double vy, double vz)
+ {
+ v[0] = vx; v[1] = vy; v[2] = vz;
+ }
+};
+
template <class T>
class adaptiveElements {
private:
- int _level;
- double _tol, _minVal, _maxVal;
- List_T *_listEle;
- int *_numEle;
- Double_Matrix *_coeffs, *_eexps;
- Double_Matrix *_coeffsGeom, *_eexpsGeom;
- Double_Matrix *_posX, *_posY, *_posZ;
- Double_Matrix *_val, *_valX, *_valY, *_valZ;
- Double_Matrix *_interpolate, *_geometry;
- void _changeResolution(int level, GMSH_Post_Plugin *plug, int *done);
- int _zoomElement(int ielem, int level, GMSH_Post_Plugin *plug);
+ Double_Matrix *_coeffsVal, *_eexpsVal, *_interpolVal;
+ Double_Matrix *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
public:
- adaptiveElements(List_T *listEle, int *numEle,
- Double_Matrix *coeffs, Double_Matrix *eexps,
- Double_Matrix *coeffsGeom=0, Double_Matrix *eexpsGeom=0)
- : _level(-1), _tol(-1.), _minVal(0.), _maxVal(0.),
- _listEle(listEle), _numEle(numEle), _coeffs(coeffs), _eexps(eexps),
- _coeffsGeom(coeffsGeom), _eexpsGeom(eexpsGeom), _posX(0), _posY(0), _posZ(0),
- _val(0), _valX(0), _valY(0), _valZ(0), _interpolate(0), _geometry(0){}
+ adaptiveElements(Double_Matrix *coeffsVal, Double_Matrix *eexpsVal,
+ Double_Matrix *coeffsGeom=0, Double_Matrix *eexpsGeom=0)
+ : _coeffsVal(coeffsVal), _eexpsVal(eexpsVal), _interpolVal(0),
+ _coeffsGeom(coeffsGeom), _eexpsGeom(eexpsGeom), _interpolGeom(0) {}
~adaptiveElements();
- // store data's step-th timestep data into _val and _pos arrays
- // (This makes adaptive views even more memory hungry than what they
- // already are due to their discontinous nature. We should evaluate
- // the performance hit incurred if we loop directly in data in
- // _zoomElement instead)
- void initData(PViewData *data, int step);
- // compute the adaptive representation and store it in _listEle
- void changeResolution(int level, double tol, GMSH_Post_Plugin *plug=0);
- // The number of nodes is supposed to be fixed in an adaptive view
- inline int getNumNodes () const { return _coeffsGeom ? _coeffsGeom->size1() : T::numNodes; }
+ void init(int level);
+ void adapt(double tol, int numComp,
+ std::vector<PCoords> &coords, std::vector<PValues> &values,
+ double &minVal, double &maxVal, GMSH_Post_Plugin *plug=0,
+ bool onlyComputeMinMax=false);
+ void addInView(double tol, int step, PViewData *in, PViewDataList *out,
+ GMSH_Post_Plugin *plug=0);
};
class adaptiveData {
--
GitLab