From 4192ff7553b09b4c6fa7af93b3fc5fd27d16795d Mon Sep 17 00:00:00 2001
From: Koen Hillewaert <koen.hillewaert@cenaero.be>
Date: Mon, 4 Aug 2014 16:26:49 +0000
Subject: [PATCH] Integrated adaptive views for the pyramids
---
Numeric/pyramidalBasis.cpp | 35 ++++-
Numeric/pyramidalBasis.h | 5 +
Plugin/Levelset.cpp | 48 +++++++
Post/PViewDataGModel.cpp | 53 ++++++--
Post/PViewDataList.cpp | 2 +-
Post/adaptiveData.cpp | 262 +++++++++++++++++++++++++++++++++++++
Post/adaptiveData.h | 51 ++++++++
7 files changed, 435 insertions(+), 21 deletions(-)
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 634ab598fe..aa9c8b2c93 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -18,7 +18,7 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
int num_points = points.size1();
- coefficients.resize(num_points, num_points);
+ bergotCoefficients.resize(num_points, num_points);
double *fval = new double[num_points];
// Invert the Vandermonde matrix
@@ -27,6 +27,29 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
bergot->f(points(j,0), points(j,1), points(j, 2), fval);
for (int i = 0; i < num_points; i++) VDM(i,j) = fval[i];
}
+ VDM.invert(bergotCoefficients);
+
+ coefficients.resize(num_points,num_points);
+ monomials.resize(num_points,3);
+
+ int idx = 0;
+ for (int i=0;i<=order;i++) {
+ for (int j=0;j<=order;j++) {
+ for (int k=0;k<=order-std::max(i,j);k++,idx++) {
+ monomials(idx,0) = i;
+ monomials(idx,1) = j;
+ monomials(idx,2) = k;
+
+ for (int l=0;l<num_points;l++) {
+ double oneMinW = std::max(1e-14,1-points(l,2));
+ VDM(idx,l) = std::pow(points(l,0),i);
+ VDM(idx,l) *= std::pow(points(l,1),j);
+ VDM(idx,l) *= std::pow(points(l,2),k);
+ VDM(idx,l) *= std::pow(oneMinW,std::max(i,j)-i-j);
+ }
+ }
+ }
+ }
VDM.invert(coefficients);
delete[] fval;
@@ -55,7 +78,7 @@ void pyramidalBasis::f(double u, double v, double w, double *val) const
for (int i = 0; i < N; i++) {
val[i] = 0.;
- for (int j = 0; j < N; j++) val[i] += coefficients(i,j)*fval[j];
+ for (int j = 0; j < N; j++) val[i] += bergotCoefficients(i,j)*fval[j];
}
delete[] fval;
@@ -76,7 +99,7 @@ void pyramidalBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf)
bergot->f(coord(iPt,0), coord(iPt,1), coord(iPt,2), fval);
for (int i = 0; i < N; i++) {
sf(iPt,i) = 0.;
- for (int j = 0; j < N; j++) sf(iPt,i) += coefficients(i,j)*fval[j];
+ for (int j = 0; j < N; j++) sf(iPt,i) += bergotCoefficients(i,j)*fval[j];
}
}
@@ -97,9 +120,9 @@ void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
for (int i = 0; i < N; i++) {
grads[i][0] = 0.; grads[i][1] = 0.; grads[i][2] = 0.;
for (int j = 0; j < N; j++) {
- grads[i][0] += coefficients(i,j)*dfval[j][0];
- grads[i][1] += coefficients(i,j)*dfval[j][1];
- grads[i][2] += coefficients(i,j)*dfval[j][2];
+ grads[i][0] += bergotCoefficients(i,j)*dfval[j][0];
+ grads[i][1] += bergotCoefficients(i,j)*dfval[j][1];
+ grads[i][2] += bergotCoefficients(i,j)*dfval[j][2];
}
}
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index 610f48005d..98157b1135 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -18,7 +18,12 @@ class pyramidalBasis: public nodalBasis
private:
// Orthogonal basis for the pyramid
BergotBasis *bergot;
+ fullMatrix<double> bergotCoefficients;
+
+ public:
+
fullMatrix<double> coefficients;
+ fullMatrix<double> monomials;
public:
pyramidalBasis(int tag);
diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp
index 8036dace51..56c1e17e06 100644
--- a/Plugin/Levelset.cpp
+++ b/Plugin/Levelset.cpp
@@ -702,6 +702,50 @@ static bool recur_sign_change(adaptivePrism *t,
}
}
+static bool recur_sign_change(adaptivePyramid *t,
+ const GMSH_LevelsetPlugin *plug)
+{
+ if (!t->e[0] || t->visible){
+ double v1 = plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
+ double v2 = plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
+ double v3 = plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
+ double v4 = plug->levelset(t->p[3]->X, t->p[3]->Y, t->p[3]->Z, t->p[3]->val);
+ double v5 = plug->levelset(t->p[4]->X, t->p[4]->Y, t->p[4]->Z, t->p[4]->val);
+ if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0 && v1 * v5 > 0)
+ t->visible = false;
+ else
+ t->visible = true;
+ return t->visible;
+ }
+ else{
+ bool sc1 = recur_sign_change(t->e[0], plug);
+ bool sc2 = recur_sign_change(t->e[1], plug);
+ bool sc3 = recur_sign_change(t->e[2], plug);
+ bool sc4 = recur_sign_change(t->e[3], plug);
+ bool sc5 = recur_sign_change(t->e[4], plug);
+ bool sc6 = recur_sign_change(t->e[5], plug);
+ bool sc7 = recur_sign_change(t->e[6], plug);
+ bool sc8 = recur_sign_change(t->e[7], plug);
+ bool sc9 = recur_sign_change(t->e[8], plug);
+ bool sc10 = recur_sign_change(t->e[9], plug);
+ if(sc1 || sc2 || sc3 || sc4 || sc5 || sc6 || sc7 || sc8 || sc9 || sc10){
+ if (!sc1) t->e[0]->visible = true;
+ if (!sc2) t->e[1]->visible = true;
+ if (!sc3) t->e[2]->visible = true;
+ if (!sc4) t->e[3]->visible = true;
+ if (!sc5) t->e[4]->visible = true;
+ if (!sc6) t->e[5]->visible = true;
+ if (!sc7) t->e[6]->visible = true;
+ if (!sc8) t->e[7]->visible = true;
+ if (!sc9) t->e[8]->visible = true;
+ if (!sc10) t->e[9]->visible = true;
+ return true;
+ }
+ t->visible = false;
+ return false;
+ }
+}
+
void GMSH_LevelsetPlugin::assignSpecificVisibility() const
{
if(adaptiveTriangle::all.size()){
@@ -724,4 +768,8 @@ void GMSH_LevelsetPlugin::assignSpecificVisibility() const
adaptivePrism *p = *adaptivePrism::all.begin();
if(!p->visible) p->visible = !recur_sign_change(p, this);
}
+ if(adaptivePyramid::all.size()){
+ adaptivePyramid *p = *adaptivePyramid::all.begin();
+ if(!p->visible) p->visible = !recur_sign_change(p, this);
+ }
}
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 42a878e2c2..5eaace718e 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -15,6 +15,7 @@
#include "MElementCut.h"
#include "Numeric.h"
#include "GmshMessage.h"
+#include "pyramidalBasis.h"
PViewDataGModel::PViewDataGModel(DataType type)
: PViewData(), _min(VAL_INF), _max(-VAL_INF), _type(type)
@@ -168,38 +169,62 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
// type, assume isoparametric elements (except for ElementData,
// for which we know the interpolation: it's constant)
int types[] = {TYPE_PNT, TYPE_LIN, TYPE_TRI, TYPE_QUA, TYPE_TET, TYPE_HEX,
- TYPE_PRI, /*TYPE_PYR - Not polynomial!,*/ TYPE_POLYG, TYPE_POLYH};
+ TYPE_PRI,TYPE_PYR, TYPE_POLYG, TYPE_POLYH};
for(unsigned 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 = (polynomialBasis*) e->getFunctionSpace();
+
+ const polynomialBasis *fs = dynamic_cast<const polynomialBasis*> (e->getFunctionSpace());
if(fs){
if(e->getPolynomialOrder() > 1){
if(_type == ElementData){
// data is constant per element: force the interpolation matrix
fullMatrix<double> coef(1, 1);
- coef(0, 0) = 1.;
- fullMatrix<double> mono(1, 3);
- mono(0, 0) = 0.;
- mono(0, 1) = 0.;
- mono(0, 2) = 0.;
- setInterpolationMatrices(types[i], coef, mono,
- fs->coefficients, fs->monomials);
- }
+ coef(0, 0) = 1.;
+ fullMatrix<double> mono(1, 3);
+ mono(0, 0) = 0.;
+ mono(0, 1) = 0.;
+ mono(0, 2) = 0.;
+ setInterpolationMatrices(types[i], coef, mono,
+ fs->coefficients, fs->monomials);
+ }
else
- setInterpolationMatrices(types[i], fs->coefficients, fs->monomials,
+ setInterpolationMatrices(types[i],
+ fs->coefficients, fs->monomials,
fs->coefficients, fs->monomials);
- }
+ }
else
setInterpolationMatrices(types[i], fs->coefficients, fs->monomials);
}
+ else {
+ const pyramidalBasis *fs = dynamic_cast<const pyramidalBasis*> (e->getFunctionSpace());
+ if(fs){
+ if(e->getPolynomialOrder() > 1){
+ if(_type == ElementData){
+ // data is constant per element: force the interpolation matrix
+ fullMatrix<double> coef(1, 1);
+ coef(0, 0) = 1.;
+ fullMatrix<double> mono(1, 3);
+ mono(0, 0) = 0.;
+ mono(0, 1) = 0.;
+ mono(0, 2) = 0.;
+ setInterpolationMatrices(types[i], coef, mono,
+ fs->coefficients, fs->monomials);
+ }
+ else
+ setInterpolationMatrices(types[i],
+ fs->coefficients, fs->monomials,
+ fs->coefficients, fs->monomials);
+ }
+ else
+ setInterpolationMatrices(types[i], fs->coefficients, fs->monomials);
+ }
+ }
}
}
}
-
}
-
return PViewData::finalize();
}
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 1886fee2ff..de18c1ca85 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -905,7 +905,7 @@ void PViewDataList::setOrder2(int type)
case TYPE_TET: typeMSH = MSH_TET_10; break;
case TYPE_HEX: typeMSH = MSH_HEX_27; break;
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
- // case TYPE_PYR: typeMSH = MSH_PYR_14; break;
+ case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
const polynomialBasis *fs = (polynomialBasis*)BasisFactory::getNodalBasis(typeMSH);
if(!fs){
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index 2d91ede6ed..5871cb6ed4 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -20,6 +20,7 @@ std::set<adaptiveVertex> adaptiveQuadrangle::allVertices;
std::set<adaptiveVertex> adaptiveTetrahedron::allVertices;
std::set<adaptiveVertex> adaptiveHexahedron::allVertices;
std::set<adaptiveVertex> adaptivePrism::allVertices;
+std::set<adaptiveVertex> adaptivePyramid::allVertices;
std::list<adaptivePoint*> adaptivePoint::all;
std::list<adaptiveLine*> adaptiveLine::all;
@@ -28,6 +29,7 @@ std::list<adaptiveQuadrangle*> adaptiveQuadrangle::all;
std::list<adaptiveTetrahedron*> adaptiveTetrahedron::all;
std::list<adaptiveHexahedron*> adaptiveHexahedron::all;
std::list<adaptivePrism*> adaptivePrism::all;
+std::list<adaptivePyramid*> adaptivePyramid::all;
int adaptivePoint::numNodes = 1;
int adaptiveLine::numNodes = 2;
@@ -36,6 +38,7 @@ int adaptiveQuadrangle::numNodes = 4;
int adaptivePrism::numNodes = 6;
int adaptiveTetrahedron::numNodes = 4;
int adaptiveHexahedron::numNodes = 8;
+int adaptivePyramid::numNodes = 5;
int adaptivePoint::numEdges = 0;
int adaptiveLine::numEdges = 1;
@@ -44,6 +47,7 @@ int adaptiveQuadrangle::numEdges = 4;
int adaptivePrism::numEdges = 9;
int adaptiveTetrahedron::numEdges = 6;
int adaptiveHexahedron::numEdges = 12;
+int adaptivePyramid::numEdges = 8;
template <class T>
static void cleanElement()
@@ -66,6 +70,40 @@ static void computeShapeFunctions(fullMatrix<double> *coeffs, fullMatrix<double>
coeffs->mult(*tmp, *sf);
}
+/*! Bergot space is characterised by polynomials
+ \f$ \mathcal B_{ijk} =
+ \mathcal P_i \left(\frac{\xi }{1-\zeta}\right)
+ \mathcal P_j \left(\frac{\eta}{1-\zeta}\right)
+ \left(1-\zeta\right)^{max(i,j)}
+ \mathcal P^{2 max(i,j),0}_k \left(2 \zeta -1\right)~|~i,j \leq p, k \leq p - max(i,j) \f$
+ and hence by the "monomials"
+ \f$ \mu_{ijk} =
+ \left(\frac{\xi }{\1-\zeta}\right)^i
+ \left(\frac{\eta}{\1-\zeta}\right)^j
+ \left(1-\zeta\right)^{max(i,j)} \zeta^k~|~i,j \leq p~,~k \leq p-max(i,j)
+ \f$
+*/
+static void computeShapeFunctionsPyramid(fullMatrix<double> *coeffs,
+ fullMatrix<double> *eexps,
+ double u, double v, double w,
+ fullVector<double> *sf,
+ fullVector<double> *tmp)
+{
+
+ double oneMinW = (w==1) ? 1e-12:1-w;
+ for(int l = 0; l < eexps->size1(); l++) {
+ int i = (*eexps)(l,0);
+ int j = (*eexps)(l,1);
+ int k = (*eexps)(l,2);
+ int m = std::max(i,j);
+ (*tmp)(l) = pow(u,i);
+ (*tmp)(l) *= pow(v,j);
+ (*tmp)(l) *= pow(w,k);
+ (*tmp)(l) *= pow(oneMinW,m-i-j);
+ }
+ coeffs->mult(*tmp, *sf);
+}
+
adaptiveVertex *adaptiveVertex::add(double x, double y, double z,
std::set<adaptiveVertex> &allVertices)
{
@@ -883,6 +921,163 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
}
}
+void adaptivePyramid::create(int maxlevel)
+{
+ cleanElement<adaptivePyramid>();
+ adaptiveVertex *p1 = adaptiveVertex::add(-1, -1, 0, allVertices);
+ adaptiveVertex *p2 = adaptiveVertex::add(1, -1, 0, allVertices);
+ adaptiveVertex *p3 = adaptiveVertex::add(1, 1, 0, allVertices);
+ adaptiveVertex *p4 = adaptiveVertex::add(-1, 1, 0, allVertices);
+ adaptiveVertex *p5 = adaptiveVertex::add(0, 0, 1, allVertices);
+ adaptivePyramid *p = new adaptivePyramid(p1, p2, p3, p4, p5);
+ recurCreate(p, maxlevel, 0);
+}
+
+void adaptivePyramid::recurCreate(adaptivePyramid *p, int maxlevel, int level)
+{
+ all.push_back(p);
+ if(level++ >= maxlevel) return;
+
+ // quad points
+ adaptiveVertex *p1 = p->p[0];
+ adaptiveVertex *p2 = p->p[1];
+ adaptiveVertex *p3 = p->p[2];
+ adaptiveVertex *p4 = p->p[3];
+
+ // apex
+ adaptiveVertex *p5 = p->p[4];
+
+ // center of the quad
+
+ adaptiveVertex *p1234 = adaptiveVertex::add
+ ((p1->x + p2->x + p3->x + p4->x)*0.25,
+ (p1->y + p2->y + p3->y + p4->y)*0.25,
+ (p1->z + p2->z + p3->z + p4->z)*0.25,allVertices);
+
+ // quad edge points
+
+ adaptiveVertex *p12 = adaptiveVertex::add
+ ((p1->x + p2->x)*0.5,
+ (p1->y + p2->y)*0.5,
+ (p1->z + p2->z)*0.5,allVertices);
+
+ adaptiveVertex *p23 = adaptiveVertex::add
+ ((p2->x + p3->x)*0.5,
+ (p2->y + p3->y)*0.5,
+ (p2->z + p3->z)*0.5,allVertices);
+
+ adaptiveVertex *p34 = adaptiveVertex::add
+ ((p3->x + p4->x)*0.5,
+ (p3->y + p4->y)*0.5,
+ (p3->z + p4->z)*0.5,allVertices);
+
+ adaptiveVertex *p41 = adaptiveVertex::add
+ ((p4->x + p1->x)*0.5,
+ (p4->y + p1->y)*0.5,
+ (p4->z + p1->z)*0.5,allVertices);
+
+ // quad vertex to apex edge points
+
+ adaptiveVertex *p15 = adaptiveVertex::add
+ ((p1->x + p5->x)*0.5,
+ (p1->y + p5->y)*0.5,
+ (p1->z + p5->z)*0.5,allVertices);
+
+ adaptiveVertex *p25 = adaptiveVertex::add
+ ((p2->x + p5->x)*0.5,
+ (p2->y + p5->y)*0.5,
+ (p2->z + p5->z)*0.5,allVertices);
+
+ adaptiveVertex *p35 = adaptiveVertex::add
+ ((p3->x + p5->x)*0.5,
+ (p3->y + p5->y)*0.5,
+ (p3->z + p5->z)*0.5,allVertices);
+
+ adaptiveVertex *p45 = adaptiveVertex::add
+ ((p4->x + p5->x)*0.5,
+ (p4->y + p5->y)*0.5,
+ (p4->z + p5->z)*0.5,allVertices);
+
+ // four base pyramids on the quad base
+
+ p->e[0] = new adaptivePyramid(p1, p12, p1234, p41, p15);
+ recurCreate(p->e[0], maxlevel, level);
+ p->e[1] = new adaptivePyramid(p2, p23, p1234, p12, p25);
+ recurCreate(p->e[1], maxlevel, level);
+ p->e[2] = new adaptivePyramid(p3, p34, p1234, p23, p35);
+ recurCreate(p->e[2], maxlevel, level);
+ p->e[3] = new adaptivePyramid(p4, p41, p1234, p34, p45);
+ recurCreate(p->e[3], maxlevel, level);
+
+ // top pyramids
+
+ p->e[4] = new adaptivePyramid(p15,p25,p35,p45,p5);
+ recurCreate(p->e[4], maxlevel, level);
+ p->e[5] = new adaptivePyramid(p15,p45,p35,p25,p1234);
+ recurCreate(p->e[5], maxlevel, level);
+
+ // degenerated pyramids to replace the remaining tetrahedral holes
+ // degenerated quad in the interior of the element, apices on the quad edges
+
+ p->e[6] = new adaptivePyramid(p1234,p25,p15,p1234,p12);
+ recurCreate(p->e[6], maxlevel, level);
+ p->e[7] = new adaptivePyramid(p1234,p35,p25,p1234,p23);
+ recurCreate(p->e[7], maxlevel, level);
+ p->e[8] = new adaptivePyramid(p1234,p45,p35,p1234,p34);
+ recurCreate(p->e[8], maxlevel, level);
+ p->e[9] = new adaptivePyramid(p1234,p15,p45,p1234,p41);
+ recurCreate(p->e[9], maxlevel, level);
+}
+
+void adaptivePyramid::error(double AVG, double tol)
+{
+ adaptivePyramid *p = *all.begin();
+ recurError(p, AVG, tol);
+}
+
+void adaptivePyramid::recurError(adaptivePyramid *p, double AVG, double tol)
+{
+ if(!p->e[0])
+ p->visible = true;
+ else {
+ double vi[10];
+ for (int i = 0; i < 10; i++) vi[i] = p->e[i]->V();
+ double vr = 0;
+ for (int i = 0; i < 6 ; i++) vr += vi[i]; // pyramids have volume V/8
+ for (int i = 6; i < 10; i++) vr += vi[i]*0.5; // tetrahedra have volume V/16
+ vr /= 8.;
+ const double v = p->V();
+ if(!p->e[0]->e[0]) {
+ if(fabs(v - vr) > AVG * tol){
+ p->visible = false;
+ for (int i = 0; i < 10; i++) recurError(p->e[i],AVG,tol);
+ }
+ else
+ p->visible = true;
+ }
+ else {
+ bool err = false;
+ for(int i = 0; i < 10; i++){
+ double vj[10];
+ for (int j = 0; j < 10; j++) vj[j] = p->e[i]->e[j]->V();
+ double vri = 0;
+ for (int j = 0; j < 6 ; j++) vri += vj[j];
+ for (int j = 6; j < 10; j++) vri += vj[j]*0.5;
+ vri /= 8.;
+ err |= (fabs((vi[i] - vri)) > AVG * tol);
+ }
+ err |= (fabs((v - vr)) > AVG * tol);
+ if(err) {
+ p->visible = false;
+ for(int i = 0; i < 10; i++)
+ recurError(p->e[i], AVG, tol);
+ }
+ else
+ p->visible = true;
+ }
+ }
+}
+
template <class T>
adaptiveElements<T>::adaptiveElements(std::vector<fullMatrix<double>*> &p)
: _coeffsVal(0), _eexpsVal(0), _interpolVal(0),
@@ -960,6 +1155,61 @@ void adaptiveElements<T>::init(int level)
#endif
}
+template <>
+void adaptiveElements<adaptivePyramid>::init(int level)
+{
+#ifdef TIMER
+ double t1 = GetTimeInSeconds();
+#endif
+
+ adaptivePyramid::create(level);
+ int numVals = _coeffsVal ? _coeffsVal->size1() : adaptivePyramid::numNodes;
+ int numNodes = _coeffsGeom ? _coeffsGeom->size1() : adaptivePyramid::numNodes;
+
+ if(_interpolVal) delete _interpolVal;
+ _interpolVal = new fullMatrix<double>(adaptivePyramid::allVertices.size(), numVals);
+
+ if(_interpolGeom) delete _interpolGeom;
+ _interpolGeom = new fullMatrix<double>(adaptivePyramid::allVertices.size(), numNodes);
+
+ fullVector<double> sfv(numVals), *tmpv = 0;
+ fullVector<double> sfg(numNodes), *tmpg = 0;
+ if(_eexpsVal) tmpv = new fullVector<double>(_eexpsVal->size1());
+ if(_eexpsGeom) tmpg = new fullVector<double>(_eexpsGeom->size1());
+
+ int i = 0;
+ for(std::set<adaptiveVertex>::iterator it = adaptivePyramid::allVertices.begin();
+ it != adaptivePyramid::allVertices.end(); ++it) {
+
+ if(_coeffsVal && _eexpsVal)
+ computeShapeFunctionsPyramid(_coeffsVal, _eexpsVal,
+ it->x, it->y, it->z, &sfv, tmpv);
+ else
+ adaptivePyramid::GSF(it->x, it->y, it->z, sfv);
+
+ for(int j = 0; j < numVals; j++)
+ (*_interpolVal)(i, j) = sfv(j);
+
+ if(_coeffsGeom && _eexpsGeom)
+ computeShapeFunctionsPyramid(_coeffsGeom, _eexpsGeom,
+ it->x, it->y, it->z, &sfg, tmpg);
+ else
+ adaptivePyramid::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;
+
+#ifdef TIMER
+ adaptiveData::timerInit += GetTimeInSeconds() - t1;
+ return;
+#endif
+}
+
template <class T>
void adaptiveElements<T>::adapt(double tol, int numComp,
std::vector<PCoords> &coords,
@@ -1134,6 +1384,11 @@ void adaptiveElements<T>::addInView(double tol, int step,
outNb = (numComp == 1) ? &out->NbSI : &out->NbVI;
outList = (numComp == 1) ? &out->SI : &out->VI;
break;
+ case 8:
+ numEle = in->getNumPyramids();
+ outNb = (numComp == 1) ? &out->NbSY : &out->NbVY;
+ outList = (numComp == 1) ? &out->SY : &out->VY;
+ break;
case 12:
numEle = in->getNumHexahedra();
outNb = (numComp == 1) ? &out->NbSH : &out->NbVH;
@@ -1227,6 +1482,10 @@ adaptiveData::adaptiveData(PViewData *data)
_inData->getInterpolationMatrices(TYPE_HEX, p);
_hexahedra = new adaptiveElements<adaptiveHexahedron>(p);
}
+ if(_inData->getNumPyramids()){
+ _inData->getInterpolationMatrices(TYPE_PYR, p);
+ _pyramids = new adaptiveElements<adaptivePyramid>(p);
+ }
}
adaptiveData::~adaptiveData()
@@ -1238,6 +1497,7 @@ adaptiveData::~adaptiveData()
if(_tetrahedra) delete _tetrahedra;
if(_prisms) delete _prisms;
if(_hexahedra) delete _hexahedra;
+ if(_pyramids) delete _pyramids;
delete _outData;
}
@@ -1257,6 +1517,7 @@ void adaptiveData::changeResolution(int step, int level, double tol,
if(_tetrahedra) _tetrahedra->init(level);
if(_prisms) _prisms->init(level);
if(_hexahedra) _hexahedra->init(level);
+ if(_pyramids) _pyramids->init(level);
}
if(plug || _step != step || _level != level || _tol != tol){
_outData->setDirty(true);
@@ -1267,6 +1528,7 @@ void adaptiveData::changeResolution(int step, int level, double tol,
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);
+ if(_pyramids) _pyramids->addInView(tol, step, _inData, _outData, plug);
_outData->finalize();
}
_step = step;
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index ff8a28e97b..15f5758a0f 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -9,6 +9,7 @@
#include <list>
#include <set>
#include <vector>
+#include <cstdlib>
#include "fullMatrix.h"
class PViewData;
@@ -287,6 +288,55 @@ class adaptiveHexahedron {
static void recurError(adaptiveHexahedron *h, double AVG, double tol);
};
+// modif koen.hillewaert@cenaero.be, 31/07/2014
+
+class adaptivePyramid {
+ public:
+ bool visible;
+ adaptiveVertex *p[5];
+ adaptivePyramid *e[10];
+ static std::list<adaptivePyramid*> all;
+ static std::set<adaptiveVertex> allVertices;
+ static int numNodes, numEdges;
+ public:
+ adaptivePyramid(adaptiveVertex *p1,
+ adaptiveVertex *p2,
+ adaptiveVertex *p3,
+ adaptiveVertex *p4,
+ adaptiveVertex *p5)
+ : visible(false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ p[3] = p4;
+ p[4] = p5;
+ for (int i=0;i<10;i++) e[i] = NULL;
+ }
+ inline double V() const
+ {
+ return (p[0]->val +
+ p[1]->val +
+ p[2]->val +
+ p[3]->val +
+ p[4]->val) / 5.;
+ }
+ // barycentric coordinates ?
+ inline static void GSF(double u, double v, double w, fullVector<double> &sf)
+ {
+ double ww = 0.25 / std::max(1e-14,1.-w);
+ sf(0) = (1 - u - w) * (1 - v - w) * ww;
+ sf(1) = (1 + u - w) * (1 - v - w) * ww;
+ sf(2) = (1 + u - w) * (1 + v - w) * ww;
+ sf(3) = (1 - u - w) * (1 + v - w) * ww;
+ sf(4) = w;
+ }
+ static void create(int maxlevel);
+ static void recurCreate(adaptivePyramid *h, int maxlevel, int level);
+ static void error(double AVG, double tol);
+ static void recurError(adaptivePyramid *h, double AVG, double tol);
+};
+
class PCoords {
public:
double c[3];
@@ -346,6 +396,7 @@ class adaptiveData {
adaptiveElements<adaptiveTetrahedron> *_tetrahedra;
adaptiveElements<adaptiveHexahedron> *_hexahedra;
adaptiveElements<adaptivePrism> *_prisms;
+ adaptiveElements<adaptivePyramid> *_pyramids;
public:
static double timerInit, timerAdapt;
adaptiveData(PViewData *data);
--
GitLab