diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index a15aef284afed00c97f0a3513951424d9751dff4..b72bb8b93db1221b897de28001776b43b6c369ca 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -993,6 +993,8 @@ StringXNumber MeshOptions_Number[] = {
"Labels display frequency?" },
{ F|O, "LabelType" , opt_mesh_label_type , 0. ,
"Type of element label (0=element number, 1=elementary entity number, 2=physical entity number, 3=partition number, 4=coordinates)" },
+ { F|O, "LcIntegrationPrecision" , opt_mesh_lc_integration_precision, 1.e-9 ,
+ "Accuracy of evaluation of the LC field for 1D mesh generation" },
{ F|O, "Light" , opt_mesh_light , 1. ,
"Enable lighting for the mesh" },
{ F|O, "LightLines" , opt_mesh_light_lines , 1. ,
@@ -1098,8 +1100,6 @@ StringXNumber MeshOptions_Number[] = {
"Create incomplete second order elements? (8-node quads, 20-node hexas, etc.)" },
{ F|O, "SecondOrderLinear" , opt_mesh_second_order_linear , 0. ,
"Should second order vertices simply be created by linear interpolation?" },
- { F|O, "LcIntegrationPrecision" , opt_mesh_lc_integration_precision, 1.e-9 ,
- "Accuracy of evaluation of the LC field for 1D mesh generation" },
{ F|O, "Smoothing" , opt_mesh_nb_smoothing , 1. ,
"Number of smoothing steps applied to the final mesh" },
{ F|O, "SmoothInternalEdges" , opt_mesh_smooth_internal_edges , 0 ,
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 4415d2c8de4fc53f3bd24b12061c3dd20713dd27..35adb58388ee30e041ff9f3ddaf5cfaeb88f6521 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -31,7 +31,7 @@ bool PViewData::finalize()
void PViewData::initAdaptiveData(int step, int level, double tol)
{
- if(!_adaptive && _interpolation.size()){
+ if(!_adaptive){
Msg::Info("Initializing adaptive data %p interp size=%d", this, _interpolation.size());
_adaptive = new adaptiveData(this);
_adaptive->changeResolution(step, level, tol);
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index f7546c15f7237a98eb858f3f521c3c5c17482a04..e653811e66b278b0cc9e7f21124bfa42417018f3 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -94,23 +94,9 @@ void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
_interpolation[edg].push_back(new Double_Matrix(fs->monomials));
}
}
- else if(edg == 4 && e->getPolynomialOrder() == 1){
- // quad4
- Double_Matrix *c = new Double_Matrix(4, 4);
- (*c)(0, 0) = .25; (*c)(0, 1) = -.25; (*c)(0, 2) = .25; (*c)(0, 3) = -.25;
- (*c)(1, 0) = .25; (*c)(1, 1) = .25; (*c)(1, 2) = -.25; (*c)(1, 3) = -.25;
- (*c)(2, 0) = .25; (*c)(2, 1) = .25; (*c)(2, 2) = .25; (*c)(2, 3) = .25;
- (*c)(3, 0) = .25; (*c)(3, 1) = -.25; (*c)(3, 2) = -.25; (*c)(3, 3) = .25;
- Double_Matrix *m = new Double_Matrix(4, 2);
- (*m)(0, 0) = 0; (*m)(0, 1) = 0;
- (*m)(1, 0) = 1; (*m)(1, 1) = 0;
- (*m)(2, 0) = 1; (*m)(2, 1) = 1;
- (*m)(3, 0) = 0; (*m)(3, 1) = 1;
- _interpolation[edg].push_back(c);
- _interpolation[edg].push_back(m);
- }
else{
- Msg::Warning("need to add interpol matrices for ele type %d", e->getTypeForMSH());
+ Msg::Info("No interpolation matrix for element type %d: using first order",
+ e->getTypeForMSH());
}
}
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index 74ed6623d307abae9ca73f1f401e6e5fa8fc255a..f367cb821f24bc5176bec0afa07aa623ac9b0007 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -47,20 +47,15 @@ static void cleanElement()
}
static void computeShapeFunctions(Double_Matrix *coeffs, Double_Matrix *eexps,
- double u, double v, double w, double *sf)
+ double u, double v, double w, Double_Vector *sf,
+ Double_Vector *tmp)
{
- static double powsuvw[256];
- for(int j = 0; j < eexps->size1(); ++j) {
- powsuvw[j] = pow(u, (*eexps)(j, 0));
- if(eexps->size2() > 1) powsuvw[j] *= pow(v, (*eexps)(j, 1));
- if(eexps->size2() > 2) powsuvw[j] *= pow(w, (*eexps)(j, 2));
- }
- for(int i = 0; i < coeffs->size1(); ++i) {
- sf[i] = 0.;
- for(int j = 0; j < coeffs->size2(); ++j) {
- sf[i] += (*coeffs)(i, j) * powsuvw[j];
- }
+ for(int i = 0; i < eexps->size1(); i++) {
+ (*tmp)(i) = pow(u, (*eexps)(i, 0));
+ if(eexps->size2() > 1) (*tmp)(i) *= pow(v, (*eexps)(i, 1));
+ if(eexps->size2() > 2) (*tmp)(i) *= pow(w, (*eexps)(i, 2));
}
+ coeffs->mult(*tmp, *sf);
}
adaptivePoint *adaptivePoint::add(double x, double y, double z,
@@ -801,9 +796,9 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
p->visible = true;
}
else {
- bool err=false;
+ bool err = false;
double ve[8];
- for(int i=0; i<8; i++){
+ for(int i = 0; i < 8; i++){
double v1 = p->e[i]->e[0]->V();
double v2 = p->e[i]->e[1]->V();
double v3 = p->e[i]->e[2]->V();
@@ -821,7 +816,7 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
err |= (fabs((p->V() - vr))>AVG*tol);
if(err) {
p->visible = false;
- for(int i=0;i<8;i++)
+ for(int i = 0; i < 8; i++)
recurError(p->e[i], AVG, tol);
}
else
@@ -830,6 +825,21 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
}
}
+template <class T>
+adaptiveElements<T>::adaptiveElements(std::vector<Double_Matrix*> &p)
+ : _coeffsVal(0), _eexpsVal(0), _interpolVal(0),
+ _coeffsGeom(0), _eexpsGeom(0), _interpolGeom(0)
+{
+ if(p.size() >= 2){
+ _coeffsVal = p[0];
+ _eexpsVal = p[1];
+ }
+ if(p.size() == 4){
+ _coeffsGeom = p[2];
+ _eexpsGeom = p[3];
+ }
+}
+
template <class T>
adaptiveElements<T>::~adaptiveElements()
{
@@ -842,7 +852,7 @@ template <class T>
void adaptiveElements<T>::init(int level)
{
T::create(level);
- int numVals = _coeffsVal->size1();
+ int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
if(_interpolVal) delete _interpolVal;
@@ -851,25 +861,34 @@ void adaptiveElements<T>::init(int level)
if(_interpolGeom) delete _interpolGeom;
_interpolGeom = new Double_Matrix(T::allPoints.size(), numNodes);
- double sf[100];
- int kk = 0;
+ Double_Vector sfv(numVals), *tmpv = 0;
+ Double_Vector sfg(numNodes), *tmpg = 0;
+ if(_eexpsVal) tmpv = new Double_Vector(_eexpsVal->size1());
+ if(_eexpsGeom) tmpg = new Double_Vector(_eexpsGeom->size1());
+
+ int i = 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);
+
+ if(_coeffsVal && _eexpsVal)
+ computeShapeFunctions(_coeffsVal, _eexpsVal, it->x, it->y, it->z, &sfv, tmpv);
else
- T::GSF(p->x, p->y, p->z, sf);
- for(int k = 0; k < numNodes; k++)
- (*_interpolGeom)(kk, k) = sf[k];
-
- kk++;
+ T::GSF(it->x, it->y, it->z, sfv);
+ for(int j = 0; j < numVals; j++)
+ (*_interpolVal)(i, j) = sfv(j);
+
+ if(_coeffsGeom && _eexpsGeom)
+ computeShapeFunctions(_coeffsGeom, _eexpsGeom, it->x, it->y, it->z, &sfg, tmpg);
+ else
+ T::GSF(it->x, it->y, it->z, sfg);
+ for(int j = 0; j < numNodes; j++)
+ (*_interpolGeom)(i, j) = sfg(j);
+
+ i++;
}
+
+ if(tmpv) delete tmpv;
+ if(tmpg) delete tmpg;
}
template <class T>
@@ -891,7 +910,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
return;
}
- int numVals = _coeffsVal->size1();
+ int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
if(numVals != values.size()){
Msg::Error("Wrong number of values in adaptation %d != %i",
numVals, values.size());
@@ -900,16 +919,18 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
Double_Vector val(numVals), res(numPoints);
if(numComp == 1){
- for(int k = 0; k < numVals; ++k)
- val(k) = values[k].v[0];
+ for(int i = 0; i < numVals; i++)
+ val(i) = values[i].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];
+ for(int i = 0; i < numVals; i++)
+ val(i) = values[i].v[0] * values[i].v[0] + values[i].v[1] * values[i].v[1] +
+ values[i].v[2] * values[i].v[2];
}
_interpolVal->mult(val, res);
-
+
+ //minVal = VAL_INF;
+ //maxVal = -VAL_INF;
for(int i = 0; i < numPoints; i++){
minVal = std::min(minVal, res(i));
maxVal = std::max(maxVal, res(i));
@@ -922,10 +943,10 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
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];
+ for(int i = 0; i < numVals; i++){
+ valx(i) = values[i].v[0];
+ valy(i) = values[i].v[1];
+ valz(i) = values[i].v[2];
}
_interpolVal->mult(valx, *resx);
_interpolVal->mult(valy, *resy);
@@ -940,27 +961,28 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
}
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];
+ for(int i = 0; i < numNodes; i++){
+ xyz(i, 0) = coords[i].c[0];
+ xyz(i, 1) = coords[i].c[1];
+ xyz(i, 2) = coords[i].c[2];
}
_interpolGeom->mult(xyz, XYZ);
- int k = 0;
+ int i = 0;
for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
it != T::allPoints.end(); ++it){
+ // ok because we know this will not change the set ordering
adaptivePoint *p = (adaptivePoint*)&(*it);
- p->val = res(k);
+ p->val = res(i);
if(resx){
- p->valx = (*resx)(k);
- p->valy = (*resy)(k);
- p->valz = (*resz)(k);
+ p->valx = (*resx)(i);
+ p->valy = (*resy)(i);
+ p->valz = (*resz)(i);
}
- p->X = XYZ(k, 0);
- p->Y = XYZ(k, 1);
- p->Z = XYZ(k, 2);
- k++;
+ p->X = XYZ(i, 0);
+ p->Y = XYZ(i, 1);
+ p->Z = XYZ(i, 2);
+ i++;
}
if(resx){
@@ -972,7 +994,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
(*it)->visible = false;
if(!plug || tol != 0.)
- T::error(maxVal - minVal, tol);
+ T::error(fabs(maxVal - minVal), tol);
if(plug)
plug->assignSpecificVisibility();
@@ -983,12 +1005,12 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
it != T::all.end(); it++){
if((*it)->visible){
adaptivePoint **p = (*it)->p;
- for(int k = 0; k < T::numNodes; ++k) {
- coords.push_back(PCoords(p[k]->X, p[k]->Y, p[k]->Z));
+ for(int i = 0; i < T::numNodes; i++) {
+ coords.push_back(PCoords(p[i]->X, p[i]->Y, p[i]->Z));
if(numComp == 1)
- values.push_back(PValues(p[k]->val));
+ values.push_back(PValues(p[i]->val));
else
- values.push_back(PValues(p[k]->valx, p[k]->valy, p[k]->valz));
+ values.push_back(PValues(p[i]->valx, p[i]->valy, p[i]->valz));
}
}
}
@@ -1094,24 +1116,30 @@ adaptiveData::adaptiveData(PViewData *data)
{
_outData = new PViewDataList(true);
std::vector<Double_Matrix*> p;
- if(_inData->getNumLines() && _inData->getInterpolationMatrices(1, p) >= 2)
- _lines = new adaptiveElements<adaptiveLine>
- (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>
- (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>
- (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>
- (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>
- (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>
- (p[0], p[1], (p.size() == 4) ? p[2] : 0, (p.size() == 4) ? p[3] : 0);
+ if(_inData->getNumLines()){
+ _inData->getInterpolationMatrices(1, p);
+ _lines = new adaptiveElements<adaptiveLine>(p);
+ }
+ if(_inData->getNumTriangles()){
+ _inData->getInterpolationMatrices(3, p);
+ _triangles = new adaptiveElements<adaptiveTriangle>(p);
+ }
+ if(_inData->getNumQuadrangles()){
+ _inData->getInterpolationMatrices(4, p);
+ _quadrangles = new adaptiveElements<adaptiveQuadrangle>(p);
+ }
+ if(_inData->getNumTetrahedra()){
+ _inData->getInterpolationMatrices(6, p);
+ _tetrahedra = new adaptiveElements<adaptiveTetrahedron>(p);
+ }
+ if(_inData->getNumPrisms()){
+ _inData->getInterpolationMatrices(9, p);
+ _prisms = new adaptiveElements<adaptivePrism>(p);
+ }
+ if(_inData->getNumHexahedra()){
+ _inData->getInterpolationMatrices(12, p);
+ _hexahedra = new adaptiveElements<adaptiveHexahedron>(p);
+ }
}
adaptiveData::~adaptiveData()
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 087116e9af39f11396bf0f39d15bc50e18fb345b..500fceb7ae266ebd2c7f4e408daeac9774f8756d 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -53,10 +53,10 @@ class adaptiveLine {
{
return (p[0]->val + p[1]->val) / 2.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = (1 - u) / 2.;
- sf[1] = (1 + u) / 2.;
+ sf(0) = (1 - u) / 2.;
+ sf(1) = (1 + u) / 2.;
}
static void create(int maxlevel);
static void recurCreate(adaptiveLine *e, int maxlevel, int level);
@@ -85,11 +85,11 @@ class adaptiveTriangle {
{
return (p[0]->val + p[1]->val + p[2]->val) / 3.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = 1. - u - v;
- sf[1] = u;
- sf[2] = v;
+ sf(0) = 1. - u - v;
+ sf(1) = u;
+ sf(2) = v;
}
static void create(int maxlevel);
static void recurCreate(adaptiveTriangle *t, int maxlevel, int level);
@@ -120,12 +120,12 @@ class adaptiveQuadrangle {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = 0.25 * (1. - u) * (1. - v);
- sf[1] = 0.25 * (1. + u) * (1. - v);
- sf[2] = 0.25 * (1. + u) * (1. + v);
- sf[3] = 0.25 * (1. - u) * (1. + v);
+ sf(0) = 0.25 * (1. - u) * (1. - v);
+ sf(1) = 0.25 * (1. + u) * (1. - v);
+ sf(2) = 0.25 * (1. + u) * (1. + v);
+ sf(3) = 0.25 * (1. - u) * (1. + v);
}
static void create(int maxlevel);
static void recurCreate(adaptiveQuadrangle *q, int maxlevel, int level);
@@ -161,14 +161,14 @@ class adaptivePrism {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val + p[4]->val + p[5]->val) / 6.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = (1. - u - v) * (1 - w) / 2;
- sf[1] = u * (1-w)/2;
- sf[2] = v*(1-w)/2;
- sf[3] = (1. - u - v)*(1+w)/2;
- sf[4] = u*(1+w)/2;
- sf[5] = v*(1+w)/2;
+ sf(0) = (1. - u - v) * (1 - w) / 2;
+ sf(1) = u * (1-w)/2;
+ sf(2) = v*(1-w)/2;
+ sf(3) = (1. - u - v)*(1+w)/2;
+ sf(4) = u*(1+w)/2;
+ sf(5) = v*(1+w)/2;
}
static void create(int maxlevel);
static void recurCreate(adaptivePrism *p, int maxlevel, int level);
@@ -200,12 +200,12 @@ class adaptiveTetrahedron {
{
return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = 1. - u - v - w;
- sf[1] = u;
- sf[2] = v;
- sf[3] = w;
+ sf(0) = 1. - u - v - w;
+ sf(1) = u;
+ sf(2) = v;
+ sf(3) = w;
}
static void create(int maxlevel);
static void recurCreate(adaptiveTetrahedron *t, int maxlevel, int level);
@@ -243,16 +243,16 @@ class adaptiveHexahedron {
return (p[0]->val + p[1]->val + p[2]->val+ p[3]->val +
p[4]->val + p[5]->val + p[6]->val+ p[7]->val) / 8.;
}
- inline static void GSF(double u, double v, double w, double sf[])
+ inline static void GSF(double u, double v, double w, Double_Vector &sf)
{
- sf[0] = 0.125 * (1 - u) * (1 - v) * (1 - w);
- sf[1] = 0.125 * (1 + u) * (1 - v) * (1 - w);
- sf[2] = 0.125 * (1 + u) * (1 + v) * (1 - w);
- sf[3] = 0.125 * (1 - u) * (1 + v) * (1 - w);
- sf[4] = 0.125 * (1 - u) * (1 - v) * (1 + w);
- sf[5] = 0.125 * (1 + u) * (1 - v) * (1 + w);
- sf[6] = 0.125 * (1 + u) * (1 + v) * (1 + w);
- sf[7] = 0.125 * (1 - u) * (1 + v) * (1 + w);
+ sf(0) = 0.125 * (1 - u) * (1 - v) * (1 - w);
+ sf(1) = 0.125 * (1 + u) * (1 - v) * (1 - w);
+ sf(2) = 0.125 * (1 + u) * (1 + v) * (1 - w);
+ sf(3) = 0.125 * (1 - u) * (1 + v) * (1 - w);
+ sf(4) = 0.125 * (1 - u) * (1 - v) * (1 + w);
+ sf(5) = 0.125 * (1 + u) * (1 - v) * (1 + w);
+ sf(6) = 0.125 * (1 + u) * (1 + v) * (1 + w);
+ sf(7) = 0.125 * (1 - u) * (1 + v) * (1 + w);
}
static void create(int maxlevel);
static void recurCreate(adaptiveHexahedron *h, int maxlevel, int level);
@@ -288,10 +288,7 @@ class adaptiveElements {
Double_Matrix *_coeffsVal, *_eexpsVal, *_interpolVal;
Double_Matrix *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
public:
- 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(std::vector<Double_Matrix*> &interpolationMatrices);
~adaptiveElements();
// create the _interpolVal and _interpolGeom matrices at the given
// refinement level