From 825d207be2c2eb3f99cfe41e381d883374657a2e Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Mon, 17 Jun 2013 13:38:19 +0000
Subject: [PATCH] Switch to new algo for generation of Points/Monomials.
Monomials order changed ! Points of Prism p>2 changed ! 3D Serendipity
elements have now nodes only on edges !
---
FunctionSpace/BasisLagrange.h | 4 +-
Numeric/fullMatrix.h | 13 +-
Numeric/nodalBasis.cpp | 928 +---------------------------------
Numeric/nodalBasis.h | 6 +-
Numeric/pointsGenerators.cpp | 9 +
Numeric/pointsGenerators.h | 2 +-
Numeric/polynomialBasis.cpp | 473 +----------------
Numeric/polynomialBasis.h | 12 +-
Numeric/pyramidalBasis.cpp | 83 +--
Numeric/pyramidalBasis.h | 17 -
Post/PViewData.cpp | 64 +++
Post/PViewData.h | 13 +
12 files changed, 128 insertions(+), 1496 deletions(-)
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 37d85d31e3..7e340a9e96 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -58,7 +58,7 @@ class BasisLagrange: public BasisLocal{
getPreEvaluatedDerivatives(unsigned int orientation) const;
const fullMatrix<double>& getCoefficient(void) const;
- const fullMatrix<double>& getMonomial(void) const;
+ const fullMatrix<int>& getMonomial(void) const;
std::vector<double>
project(const MElement& element,
@@ -126,7 +126,7 @@ getCoefficient(void) const{
return lBasis->coefficients;
}
-inline const fullMatrix<double>& BasisLagrange::
+inline const fullMatrix<int>& BasisLagrange::
getMonomial(void) const{
return lBasis->monomials;
}
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index f278fe9836..420e6e4df5 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -475,6 +475,12 @@ class fullMatrix
setAll(scalar(0.));
return false; // no reallocation
}
+ void reshape(int nbRows, int nbColumns){
+ if (nbRows*nbColumns != size1()*size2())
+ Msg::Error("Invalid reshape, total number of entries must be equal");
+ _r = nbRows;
+ _c = nbColumns;
+ }
void setAsProxy(const fullMatrix<scalar> &original)
{
if(_data && _own_data)
@@ -788,12 +794,5 @@ class fullMatrix
#endif
;
- void reshape(int nbRows, int nbColumns){
- if (nbRows*nbColumns != size1()*size2())
- Msg::Error("Invalid reshape, total number of entries must be equal");
- _r = nbRows;
- _c = nbColumns;
- }
-
};
#endif
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index f6375623b8..874d8d1e13 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -7,786 +7,7 @@
#include <cmath>
#include "nodalBasis.h"
#include "BasisFactory.h"
-//#include "pointsGenerators.h"
-
-
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-
-//KH : caveat : node coordinates are not yet coherent with node numbering associated
-// to numbering of principal vertices of face !!!!
-
-// uv surface - orientation v0-v2-v1
-
-
-
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= order; v[2]= 0.; w[2] = 0.;
-
- // edges
- for (int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k] = k + 1;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = order - 1 - k;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
-
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3* (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-
-
-// uw surface - orientation v0-v1-v3
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0] = 0.; v[0]= 0.; w[0] = 0.;
- u[1] = order; v[1]= 0.; w[1] = 0.;
- u[2] = 0.; v[2]= 0.; w[2] = order;
- // edges
- for (int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = 0.;
- w[3 + order - 1+ k ] = k + 1;
-
- u[3 + 2 * (order - 1) + k] = 0. ;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-
-
-// vw surface - orientation v0-v3-v2
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= 0.; w[1] = order;
- u[2]= 0.; v[2]= order; w[2] = 0.;
- // edges
- for (int k = 0; k < (order - 1); k++){
-
- u[3 + k] = 0.;
- v[3 + k] = 0.;
- w[3 + k] = k + 1;
-
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = order - 1 - k;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = 0.;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-
-
-// uvw surface - orientation v3-v1-v2
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- u[0]= 0.; v[0]= 0.; w[0] = order;
- u[1]= order; v[1]= 0.; w[1] = 0.;
- u[2]= 0.; v[2]= order; w[2] = 0.;
- // edges
- for (int k = 0; k < (order - 1); k++){
-
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = order - 1 - k;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = k + 1;
- w[3 + order - 1 + k] = 0.;
-
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
- w[3 + 2 * (order - 1) + k] = k + 1;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
- &w[3 + 3 * (order - 1)]);
- for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-
-
-static fullMatrix<double> generate1DPoints(int order)
-{
- fullMatrix<double> line(order + 1, 1);
- line(0,0) = 0;
- if (order > 0) {
- line(0, 0) = -1.;
- line(1, 0) = 1.;
- double dd = 2. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
- }
- return line;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
-{
- int nbPoints =
- (serendip ?
- 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
- (order + 1) * (order + 2) * (order + 3) / 6);
-
- fullMatrix<double> point(nbPoints, 3);
-
- double overOrder = (order == 0 ? 1. : 1. / order);
-
- point(0, 0) = 0.;
- point(0, 1) = 0.;
- point(0, 2) = 0.;
-
- if (order > 0) {
- point(1, 0) = order;
- point(1, 1) = 0;
- point(1, 2) = 0;
-
- point(2, 0) = 0.;
- point(2, 1) = order;
- point(2, 2) = 0.;
-
- point(3, 0) = 0.;
- point(3, 1) = 0.;
- point(3, 2) = order;
-
- // edges e5 and e6 switched in original version, opposite direction
- // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
- if (order > 1) {
- for (int k = 0; k < (order - 1); k++) {
- point(4 + k, 0) = k + 1;
- point(4 + order - 1 + k, 0) = order - 1 - k;
- point(4 + 2 * (order - 1) + k, 0) = 0.;
- point(4 + 3 * (order - 1) + k, 0) = 0.;
- // point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
- // point(4 + 5 * (order - 1) + k, 0) = 0.;
- point(4 + 4 * (order - 1) + k, 0) = 0.;
- point(4 + 5 * (order - 1) + k, 0) = k+1;
-
- point(4 + k, 1) = 0.;
- point(4 + order - 1 + k, 1) = k + 1;
- point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 3 * (order - 1) + k, 1) = 0.;
- // point(4 + 4 * (order - 1) + k, 1) = 0.;
- // point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 1) = k+1;
- point(4 + 5 * (order - 1) + k, 1) = 0.;
-
- point(4 + k, 2) = 0.;
- point(4 + order - 1 + k, 2) = 0.;
- point(4 + 2 * (order - 1) + k, 2) = 0.;
- point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if (order > 2) {
- int ns = 4 + 6 * (order - 1);
- int nbdofface = nbdoftriangle(order - 3);
-
- double *u = new double[nbdofface];
- double *v = new double[nbdofface];
- double *w = new double[nbdofface];
-
- nodepositionface0(order - 3, u, v, w);
-
- // u-v plane
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3);
- }
-
- ns = ns + nbdofface;
-
- // u-w plane
-
- nodepositionface1(order - 3, u, v, w);
-
- for (int i=0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) ;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface2(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3);
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- // u-v-w plane
-
- ns = ns + nbdofface;
-
- nodepositionface3(order - 3, u, v, w);
-
- for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
- point(ns + i, 1) = v[i] * (order - 3) + 1.;
- point(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (!serendip && order > 3) {
-
- fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
- for (int k = 0; k < interior.size1(); k++) {
- point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
- point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
- point(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
- }
- }
- }
- }
- }
-
- point.scale(overOrder);
- return point;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
-{
- int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
- fullMatrix<double> point(nbPoints, 2);
-
- point(0, 0) = 0;
- point(0, 1) = 0;
-
- if (order > 0) {
- double dd = 1. / order;
-
- point(1, 0) = 1;
- point(1, 1) = 0;
- point(2, 0) = 0;
- point(2, 1) = 1;
-
- int index = 3;
-
- if (order > 1) {
-
- double ksi = 0;
- double eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- ksi = 1.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi -= dd;
- eta += dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- eta = 1.;
- ksi = 0.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta -= dd;
- point(index, 0) = ksi;
- point(index, 1) = eta;
- }
-
- if (order > 2 && !serendip) {
- fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
- inner.scale(1. - 3. * dd);
- inner.add(dd);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- return point;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
-{
- const double prism18Pts[18][3] = {
- {0, 0, -1}, // 0
- {1, 0, -1}, // 1
- {0, 1, -1}, // 2
- {0, 0, 1}, // 3
- {1, 0, 1}, // 4
- {0, 1, 1}, // 5
- {0.5, 0, -1}, // 6
- {0, 0.5, -1}, // 7
- {0, 0, 0}, // 8
- {0.5, 0.5, -1}, // 9
- {1, 0, 0}, // 10
- {0, 1, 0}, // 11
- {0.5, 0, 1}, // 12
- {0, 0.5, 1}, // 13
- {0.5, 0.5, 1}, // 14
- {0.5, 0, 0}, // 15
- {0, 0.5, 0}, // 16
- {0.5, 0.5, 0}, // 17
- };
-
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> point(nbPoints, 3);
-
- int index = 0;
- fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
- fullMatrix<double> linePoints = generate1DPoints(order);
-
- if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- point(i,j) = prism18Pts[i][j];
- else
- for (int j = 0; j <linePoints.size1() ; j++) {
- for (int i = 0; i < triPoints.size1(); i++) {
- point(index,0) = triPoints(i,0);
- point(index,1) = triPoints(i,1);
- point(index,2) = linePoints(j,0);
- index ++;
- }
- }
-
- return point;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
-{
- int nbPoints = serendip ? order*4 : (order+1)*(order+1);
- fullMatrix<double> point(nbPoints, 2);
-
- if (order > 0) {
- point(0, 0) = -1;
- point(0, 1) = -1;
- point(1, 0) = 1;
- point(1, 1) = -1;
- point(2, 0) = 1;
- point(2, 1) = 1;
- point(3, 0) = -1;
- point(3, 1) = 1;
-
- if (order > 1) {
- int index = 4;
- const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
- for (int iedge=0; iedge<4; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
- point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
- }
- }
- // FIXIT it was > 2 and I beleive it is >= 2 !!
- if (order >= 2 && !serendip) {
- fullMatrix<double> inner = gmshGeneratePointsQuad(order - 2, false);
- inner.scale(1. - 2./order);
- point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- else {
- point(0, 0) = 0;
- point(0, 1) = 0;
- }
- return point;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsHex(int order, bool serendip)
-{
- int nbPoints = (order+1)*(order+1)*(order+1);
- if (serendip) nbPoints -= (order-1)*(order-1)*(order-1);
- fullMatrix<double> point(nbPoints, 3);
-
- // should be a public member of MHexahedron, just copied
- static const int edges[12][2] = {
- {0, 1},
- {0, 3},
- {0, 4},
- {1, 2},
- {1, 5},
- {2, 3},
- {2, 6},
- {3, 7},
- {4, 5},
- {4, 7},
- {5, 6},
- {6, 7}
- };
- static const int faces[6][4] = {
- {0, 3, 2, 1},
- {0, 1, 5, 4},
- {0, 4, 7, 3},
- {1, 2, 6, 5},
- {2, 3, 7, 6},
- {4, 5, 6, 7}
- };
-
- if (order > 0) {
- point(0, 0) = -1;
- point(0, 1) = -1;
- point(0, 2) = -1;
- point(1, 0) = 1;
- point(1, 1) = -1;
- point(1, 2) = -1;
- point(2, 0) = 1;
- point(2, 1) = 1;
- point(2, 2) = -1;
- point(3, 0) = -1;
- point(3, 1) = 1;
- point(3, 2) = -1;
-
- point(4, 0) = -1;
- point(4, 1) = -1;
- point(4, 2) = 1;
- point(5, 0) = 1;
- point(5, 1) = -1;
- point(5, 2) = 1;
- point(6, 0) = 1;
- point(6, 1) = 1;
- point(6, 2) = 1;
- point(7, 0) = -1;
- point(7, 1) = 1;
- point(7, 2) = 1;
-
- if (order > 1) {
- int index = 8;
- for (int iedge=0; iedge<12; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
- point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
- point(index, 2) = point(p0, 2) + i*(point(p1,2)-point(p0,2))/order;
- }
- }
-
- fullMatrix<double> fp = gmshGeneratePointsQuad(order - 2, false);
- fp.scale(1. - 2./order);
- for (int iface=0; iface<6; iface++) {
- int p0 = faces[iface][0];
- int p1 = faces[iface][1];
- int p2 = faces[iface][2];
- int p3 = faces[iface][3];
- for (int i=0;i<fp.size1();i++, index++){
- const double u = fp(i,0);
- const double v = fp(i,1);
- const double newU =
- 0.25*(1.-u)*(1.-v)*point(p0,0) +
- 0.25*(1.+u)*(1.-v)*point(p1,0) +
- 0.25*(1.+u)*(1.+v)*point(p2,0) +
- 0.25*(1.-u)*(1.+v)*point(p3,0) ;
- const double newV =
- 0.25*(1.-u)*(1.-v)*point(p0,1) +
- 0.25*(1.+u)*(1.-v)*point(p1,1) +
- 0.25*(1.+u)*(1.+v)*point(p2,1) +
- 0.25*(1.-u)*(1.+v)*point(p3,1) ;
- const double newW =
- 0.25*(1.-u)*(1.-v)*point(p0,2) +
- 0.25*(1.+u)*(1.-v)*point(p1,2) +
- 0.25*(1.+u)*(1.+v)*point(p2,2) +
- 0.25*(1.-u)*(1.+v)*point(p3,2) ;
- point(index, 0) = newU;
- point(index, 1) = newV;
- point(index, 2) = newW;
- }
- }
- if (!serendip) {
- fullMatrix<double> inner = gmshGeneratePointsHex(order - 2, false);
- inner.scale(1. - 2./order);
- point.copy(inner, 0, nbPoints - index, 0, 3, index, 0);
- }
- }
- }
- else if (order == 0){
- point(0, 0) = 0;
- point(0, 1) = 0;
- point(0, 2) = 0;
- }
- return point;
-}
-
-
-
-static fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
-{
- int nbPoints = (order+1)*((order+1)+1)*(2*(order+1)+1)/6;
-
- // Remove volume points if incomplete basis
- if (serendip && order > 2) nbPoints -= (order-2)*((order-2)+1)*(2*(order-2)+1)/6;
-
- fullMatrix<double> points(nbPoints, 3);
-
- static const int edges[8][2] = {
- {0, 1},
- {0, 3},
- {0, 4},
- {1, 2},
- {1, 4},
- {2, 3},
- {2, 4},
- {3, 4},
- };
- static const int faces[5][4] = {
- {0, 1, 4, -1},
- {3, 0, 4, -1},
- {1, 2, 4, -1},
- {2, 3, 4, -1},
- {0, 3, 2, 1}
- };
-
- if (order == 0) {
- points(0,0) = 0.0;
- points(0,1) = 0.0;
- points(0,2) = 0.0;
- return points;
- }
-
- // Principal vertices of the pyramid
- points(0,0) = -1.0; points(0,1) = -1.0; points(0,2) = 0.0;
- points(1,0) = 1.0; points(1,1) = -1.0; points(1,2) = 0.0;
- points(2,0) = 1.0; points(2,1) = 1.0; points(2,2) = 0.0;
- points(3,0) = -1.0; points(3,1) = 1.0; points(3,2) = 0.0;
- points(4,0) = 0.0; points(4,1) = 0.0; points(4,2) = 1.0;
-
- if (order > 1) {
- int p = 5;
-
- // Edges
- for (int e = 0; e < 8; e++) {
- double vec[3];
- vec[0] = points(edges[e][1],0) - points(edges[e][0],0);
- vec[1] = points(edges[e][1],1) - points(edges[e][0],1);
- vec[2] = points(edges[e][1],2) - points(edges[e][0],2);
- for (int i = 0; i < order - 1; i++) {
- points(p,0) = points(edges[e][0],0) + vec[0]/order * (i+1);
- points(p,1) = points(edges[e][0],1) + vec[1]/order * (i+1);
- points(p,2) = points(edges[e][0],2) + vec[2]/order * (i+1);
- p++;
- }
- }
-
- // Faces
- for (int f = 0; f < 4; f++) {
- fullMatrix<double> fp = gmshGeneratePointsTriangle(order, false);
-
- fullVector<double> u(4);
- fullVector<double> v(4);
- fullVector<double> w(4);
- fullVector<double> y(4);
-
- fullVector<double> up(4);
- fullVector<double> vp(4);
- fullVector<double> wp(4);
- fullVector<double> yp(4);
-
- for (int c = 0; c < 3; c++) {
- u(c) = fp(0,c);
- v(c) = fp(1,c);
- w(c) = fp(2,c);
- up(c) = points(faces[f][0],c);
- vp(c) = points(faces[f][1],c);
- wp(c) = points(faces[f][2],c);
- }
-
- u(2) = 0.0;
- v(2) = 0.0;
- w(2) = 0.0;
- u(3) = 1.0;
- v(3) = 1.0;
- w(3) = 1.0;
-
- y(0) = u(0) + ((v(1)-u(1))*(w(2)-u(2)) - (v(2)-u(2))*(w(1)-u(1)));
- y(1) = u(1) + ((v(2)-u(2))*(w(0)-u(0)) - (v(0)-u(0))*(w(2)-u(2)));
- y(2) = u(2) + ((v(0)-u(0))*(w(1)-u(1)) - (v(1)-u(1))*(w(0)-u(0)));
- y(3) = 1.0;
-
- up(3) = 1.0;
- vp(3) = 1.0;
- wp(3) = 1.0;
-
- yp(0) = up(0) + ((vp(1)-up(1))*(wp(2)-up(2)) - (vp(2)-up(2))*(wp(1)-up(1)));
- yp(1) = up(1) + ((vp(2)-up(2))*(wp(0)-up(0)) - (vp(0)-up(0))*(wp(2)-up(2)));
- yp(2) = up(2) + ((vp(0)-up(0))*(wp(1)-up(1)) - (vp(1)-up(1))*(wp(0)-up(0)));
- yp(3) = 1.0;
-
- fullMatrix<double> M(4,4);
- fullMatrix<double> Mp(4,4);
-
- M(0,0) = u(0); M(0,1) = v(0); M(0,2) = w(0); M(0,3) = y(0);
- M(1,0) = u(1); M(1,1) = v(1); M(1,2) = w(1); M(1,3) = y(1);
- M(2,0) = u(2); M(2,1) = v(2); M(2,2) = w(2); M(2,3) = y(2);
- M(3,0) = u(3); M(3,1) = v(3); M(3,2) = w(3); M(3,3) = y(3);
-
- Mp(0,0) = up(0); Mp(0,1) = vp(0); Mp(0,2) = wp(0); Mp(0,3) = yp(0);
- Mp(1,0) = up(1); Mp(1,1) = vp(1); Mp(1,2) = wp(1); Mp(1,3) = yp(1);
- Mp(2,0) = up(2); Mp(2,1) = vp(2); Mp(2,2) = wp(2); Mp(2,3) = yp(2);
- Mp(3,0) = up(3); Mp(3,1) = vp(3); Mp(3,2) = wp(3); Mp(3,3) = yp(3);
-
-
- M.invertInPlace();
- fullMatrix<double> T(4,4);
- Mp.mult(M, T);
-
- for(int k = 3*order; k < fp.size1(); k++) {
- fullVector<double> pnt(4);
- pnt(0) = fp(k,0);
- pnt(1) = fp(k,1);
- pnt(2) = 0.0;
- pnt(3) = 1.0;
-
- fullVector<double> res(4);
- T.mult(pnt, res);
-
- points(p, 0) = res(0);
- points(p, 1) = res(1);
- points(p, 2) = res(2);
-
- p++;
-
- }
- }
-
- fullMatrix<double> fpq = gmshGeneratePointsQuad(order-2, false);
- fullMatrix<double> transform(2,2);
- fullMatrix<double> scaled(fpq.size1(), fpq.size2());
- transform(0,0) = (order-2)/(double)order; transform(0,1) = 0.0;
- transform(1,0) = 0.0; transform(1,1) = (order-2)/(double)order;
- fpq.mult(transform, scaled);
-
- for (int i = 0; i < scaled.size1(); i++) {
- points(p, 0) = scaled(i, 1);
- points(p, 1) = scaled(i, 0);
- points(p, 2) = 0.0;
- p++;
- }
-
- // Volume
- if ((!serendip) && (order > 2)) {
- fullMatrix<double> volume_points = gmshGeneratePointsPyramid(order - 3, false);
- // scale to order-3/order
- fullMatrix<double> T(3,3);
- T(0,0) = (order - 3.0) / order; T(0,1) = 0.0; T(0,2) = 0.0;
- T(1,0) = 0.0; T(1,1) = (order - 3.0) / order; T(1,2) = 0.0;
- T(2,0) = 0.0; T(2,1) = 0.0; T(2,2) = (order - 3.0) / order;
- fullMatrix<double> scaled(volume_points.size1(), volume_points.size2());
- volume_points.mult(T, scaled);
-
- // translate 1/order
- for (int i = 0; i < scaled.size1(); i++) {
- points(p, 0) = scaled(i,0);
- points(p, 1) = scaled(i,1);
- points(p, 2) = scaled(i,2) + 1.0/order;
- p++;
- }
- }
- }
- return points;
-}
-
+#include "pointsGenerators.h"
static void generate1dVertexClosure(nodalBasis::clCont &closure, int order)
@@ -1378,153 +599,21 @@ nodalBasis::nodalBasis(int tag)
parentType = MElement::ParentTypeFromTag(tag);
order = MElement::OrderFromTag(tag);
serendip = MElement::SerendipityFromTag(tag) > 1;
- /*
- switch (tag) {
- case MSH_PNT : parentType = TYPE_PNT; order = 0; serendip = false; break;
- case MSH_LIN_1 : parentType = TYPE_LIN; order = 0; serendip = false; break;
- case MSH_LIN_2 : parentType = TYPE_LIN; order = 1; serendip = false; break;
- case MSH_LIN_3 : parentType = TYPE_LIN; order = 2; serendip = false; break;
- case MSH_LIN_4 : parentType = TYPE_LIN; order = 3; serendip = false; break;
- case MSH_LIN_5 : parentType = TYPE_LIN; order = 4; serendip = false; break;
- case MSH_LIN_6 : parentType = TYPE_LIN; order = 5; serendip = false; break;
- case MSH_LIN_7 : parentType = TYPE_LIN; order = 6; serendip = false; break;
- case MSH_LIN_8 : parentType = TYPE_LIN; order = 7; serendip = false; break;
- case MSH_LIN_9 : parentType = TYPE_LIN; order = 8; serendip = false; break;
- case MSH_LIN_10 : parentType = TYPE_LIN; order = 9; serendip = false; break;
- case MSH_LIN_11 : parentType = TYPE_LIN; order = 10;serendip = false; break;
- case MSH_TRI_1 : parentType = TYPE_TRI; order = 0; serendip = false; break;
- case MSH_TRI_3 : parentType = TYPE_TRI; order = 1; serendip = false; break;
- case MSH_TRI_6 : parentType = TYPE_TRI; order = 2; serendip = false; break;
- case MSH_TRI_10 : parentType = TYPE_TRI; order = 3; serendip = false; break;
- case MSH_TRI_15 : parentType = TYPE_TRI; order = 4; serendip = false; break;
- case MSH_TRI_21 : parentType = TYPE_TRI; order = 5; serendip = false; break;
- case MSH_TRI_28 : parentType = TYPE_TRI; order = 6; serendip = false; break;
- case MSH_TRI_36 : parentType = TYPE_TRI; order = 7; serendip = false; break;
- case MSH_TRI_45 : parentType = TYPE_TRI; order = 8; serendip = false; break;
- case MSH_TRI_55 : parentType = TYPE_TRI; order = 9; serendip = false; break;
- case MSH_TRI_66 : parentType = TYPE_TRI; order = 10;serendip = false; break;
- case MSH_TRI_9 : parentType = TYPE_TRI; order = 3; serendip = true; break;
- case MSH_TRI_12 : parentType = TYPE_TRI; order = 4; serendip = true; break;
- case MSH_TRI_15I : parentType = TYPE_TRI; order = 5; serendip = true; break;
- case MSH_TRI_18 : parentType = TYPE_TRI; order = 6; serendip = true; break;
- case MSH_TRI_21I : parentType = TYPE_TRI; order = 7; serendip = true; break;
- case MSH_TRI_24 : parentType = TYPE_TRI; order = 8; serendip = true; break;
- case MSH_TRI_27 : parentType = TYPE_TRI; order = 9; serendip = true; break;
- case MSH_TRI_30 : parentType = TYPE_TRI; order = 10;serendip = true; break;
- case MSH_TET_1 : parentType = TYPE_TET; order = 0; serendip = false; break;
- case MSH_TET_4 : parentType = TYPE_TET; order = 1; serendip = false; break;
- case MSH_TET_10 : parentType = TYPE_TET; order = 2; serendip = false; break;
- case MSH_TET_20 : parentType = TYPE_TET; order = 3; serendip = false; break;
- case MSH_TET_35 : parentType = TYPE_TET; order = 4; serendip = false; break;
- case MSH_TET_56 : parentType = TYPE_TET; order = 5; serendip = false; break;
- case MSH_TET_84 : parentType = TYPE_TET; order = 6; serendip = false; break;
- case MSH_TET_120 : parentType = TYPE_TET; order = 7; serendip = false; break;
- case MSH_TET_165 : parentType = TYPE_TET; order = 8; serendip = false; break;
- case MSH_TET_220 : parentType = TYPE_TET; order = 9; serendip = false; break;
- case MSH_TET_286 : parentType = TYPE_TET; order = 10;serendip = false; break;
- case MSH_TET_34 : parentType = TYPE_TET; order = 4; serendip = true; break;
- case MSH_TET_52 : parentType = TYPE_TET; order = 5; serendip = true; break;
- case MSH_TET_74 : parentType = TYPE_TET; order = 6; serendip = true; break;
- case MSH_TET_100 : parentType = TYPE_TET; order = 7; serendip = true; break;
- case MSH_TET_130 : parentType = TYPE_TET; order = 8; serendip = true; break;
- case MSH_TET_164 : parentType = TYPE_TET; order = 9; serendip = true; break;
- case MSH_TET_202 : parentType = TYPE_TET; order = 10;serendip = true; break;
- case MSH_QUA_1 : parentType = TYPE_QUA; order = 0; serendip = false; break;
- case MSH_QUA_4 : parentType = TYPE_QUA; order = 1; serendip = false; break;
- case MSH_QUA_9 : parentType = TYPE_QUA; order = 2; serendip = false; break;
- case MSH_QUA_16 : parentType = TYPE_QUA; order = 3; serendip = false; break;
- case MSH_QUA_25 : parentType = TYPE_QUA; order = 4; serendip = false; break;
- case MSH_QUA_36 : parentType = TYPE_QUA; order = 5; serendip = false; break;
- case MSH_QUA_49 : parentType = TYPE_QUA; order = 6; serendip = false; break;
- case MSH_QUA_64 : parentType = TYPE_QUA; order = 7; serendip = false; break;
- case MSH_QUA_81 : parentType = TYPE_QUA; order = 8; serendip = false; break;
- case MSH_QUA_100 : parentType = TYPE_QUA; order = 9; serendip = false; break;
- case MSH_QUA_121 : parentType = TYPE_QUA; order = 10;serendip = false; break;
- case MSH_QUA_8 : parentType = TYPE_QUA; order = 2; serendip = true; break;
- case MSH_QUA_12 : parentType = TYPE_QUA; order = 3; serendip = true; break;
- case MSH_QUA_16I : parentType = TYPE_QUA; order = 4; serendip = true; break;
- case MSH_QUA_20 : parentType = TYPE_QUA; order = 5; serendip = true; break;
- case MSH_QUA_24 : parentType = TYPE_QUA; order = 6; serendip = true; break;
- case MSH_QUA_28 : parentType = TYPE_QUA; order = 7; serendip = true; break;
- case MSH_QUA_32 : parentType = TYPE_QUA; order = 8; serendip = true; break;
- case MSH_QUA_36I : parentType = TYPE_QUA; order = 9; serendip = true; break;
- case MSH_QUA_40 : parentType = TYPE_QUA; order = 10;serendip = true; break;
- case MSH_PRI_1 : parentType = TYPE_PRI; order = 0; serendip = false; break;
- case MSH_PRI_6 : parentType = TYPE_PRI; order = 1; serendip = false; break;
- case MSH_PRI_18 : parentType = TYPE_PRI; order = 2; serendip = false; break;
- case MSH_PRI_40 : parentType = TYPE_PRI; order = 3; serendip = false; break;
- case MSH_PRI_75 : parentType = TYPE_PRI; order = 4; serendip = false; break;
- case MSH_PRI_126 : parentType = TYPE_PRI; order = 5; serendip = false; break;
- case MSH_PRI_196 : parentType = TYPE_PRI; order = 6; serendip = false; break;
- case MSH_PRI_288 : parentType = TYPE_PRI; order = 7; serendip = false; break;
- case MSH_PRI_405 : parentType = TYPE_PRI; order = 8; serendip = false; break;
- case MSH_PRI_550 : parentType = TYPE_PRI; order = 9; serendip = false; break;
- case MSH_PRI_15 : parentType = TYPE_PRI; order = 2; serendip = true; break;
- case MSH_PRI_38 : parentType = TYPE_PRI; order = 3; serendip = true; break;
- case MSH_PRI_66 : parentType = TYPE_PRI; order = 4; serendip = true; break;
- case MSH_PRI_102 : parentType = TYPE_PRI; order = 5; serendip = true; break;
- case MSH_PRI_146 : parentType = TYPE_PRI; order = 6; serendip = true; break;
- case MSH_PRI_198 : parentType = TYPE_PRI; order = 7; serendip = true; break;
- case MSH_PRI_258 : parentType = TYPE_PRI; order = 8; serendip = true; break;
- case MSH_PRI_326 : parentType = TYPE_PRI; order = 9; serendip = true; break;
- case MSH_HEX_1 : parentType = TYPE_HEX; order = 0; serendip = false; break;
- case MSH_HEX_8 : parentType = TYPE_HEX; order = 1; serendip = false; break;
- case MSH_HEX_27 : parentType = TYPE_HEX; order = 2; serendip = false; break;
- case MSH_HEX_64 : parentType = TYPE_HEX; order = 3; serendip = false; break;
- case MSH_HEX_125 : parentType = TYPE_HEX; order = 4; serendip = false; break;
- case MSH_HEX_216 : parentType = TYPE_HEX; order = 5; serendip = false; break;
- case MSH_HEX_343 : parentType = TYPE_HEX; order = 6; serendip = false; break;
- case MSH_HEX_512 : parentType = TYPE_HEX; order = 7; serendip = false; break;
- case MSH_HEX_729 : parentType = TYPE_HEX; order = 8; serendip = false; break;
- case MSH_HEX_1000: parentType = TYPE_HEX; order = 9; serendip = false; break;
- case MSH_HEX_20 : parentType = TYPE_HEX; order = 2; serendip = true; break;
- case MSH_HEX_56 : parentType = TYPE_HEX; order = 3; serendip = true; break;
- case MSH_HEX_98 : parentType = TYPE_HEX; order = 4; serendip = true; break;
- case MSH_HEX_152 : parentType = TYPE_HEX; order = 5; serendip = true; break;
- case MSH_HEX_218 : parentType = TYPE_HEX; order = 6; serendip = true; break;
- case MSH_HEX_296 : parentType = TYPE_HEX; order = 7; serendip = true; break;
- case MSH_HEX_386 : parentType = TYPE_HEX; order = 8; serendip = true; break;
- case MSH_HEX_488 : parentType = TYPE_HEX; order = 9; serendip = true; break;
- case MSH_PYR_1 : parentType = TYPE_PYR; order = 0; serendip = false; break;
- case MSH_PYR_5 : parentType = TYPE_PYR; order = 1; serendip = false; break;
- case MSH_PYR_14 : parentType = TYPE_PYR; order = 2; serendip = false; break;
- case MSH_PYR_30 : parentType = TYPE_PYR; order = 3; serendip = false; break;
- case MSH_PYR_55 : parentType = TYPE_PYR; order = 4; serendip = false; break;
- case MSH_PYR_91 : parentType = TYPE_PYR; order = 5; serendip = false; break;
- case MSH_PYR_140 : parentType = TYPE_PYR; order = 6; serendip = false; break;
- case MSH_PYR_204 : parentType = TYPE_PYR; order = 7; serendip = false; break;
- case MSH_PYR_285 : parentType = TYPE_PYR; order = 8; serendip = false; break;
- case MSH_PYR_385 : parentType = TYPE_PYR; order = 9; serendip = false; break;
- case MSH_PYR_13 : parentType = TYPE_PYR; order = 2; serendip = false; break;
- case MSH_PYR_29 : parentType = TYPE_PYR; order = 3; serendip = true; break;
- case MSH_PYR_50 : parentType = TYPE_PYR; order = 4; serendip = true; break;
- case MSH_PYR_77 : parentType = TYPE_PYR; order = 5; serendip = true; break;
- case MSH_PYR_110 : parentType = TYPE_PYR; order = 6; serendip = true; break;
- case MSH_PYR_149 : parentType = TYPE_PYR; order = 7; serendip = true; break;
- case MSH_PYR_194 : parentType = TYPE_PYR; order = 8; serendip = true; break;
- case MSH_PYR_245 : parentType = TYPE_PYR; order = 9; serendip = true; break;
- default :
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- parentType = TYPE_TET; order = 1; serendip = false; break;
- }
- */
+ dimension = MElement::DimensionFromTag(tag);
switch (parentType) {
case TYPE_PNT :
numFaces = 1;
- dimension = 0;
- points = generate1DPoints(0);
+ points = gmshGeneratePointsLine(0);
break;
case TYPE_LIN :
numFaces = 2;
- dimension = 1;
- points = generate1DPoints(order);
+ points = gmshGeneratePointsLine(order);
generate1dVertexClosure(closures, order);
generate1dVertexClosureFull(fullClosures, closureRef, order);
break;
case TYPE_TRI :
numFaces = 3;
- dimension = 2;
points = gmshGeneratePointsTriangle(order, serendip);
if (order == 0) {
generateClosureOrder0(closures, 6);
@@ -1538,8 +627,7 @@ nodalBasis::nodalBasis(int tag)
break;
case TYPE_QUA :
numFaces = 4;
- dimension = 2;
- points = gmshGeneratePointsQuad(order, serendip);
+ points = gmshGeneratePointsQuadrangle(order, serendip);
if (order == 0) {
generateClosureOrder0(closures, 8);
generateClosureOrder0(fullClosures, 8);
@@ -1552,7 +640,6 @@ nodalBasis::nodalBasis(int tag)
break;
case TYPE_TET :
numFaces = 4;
- dimension = 3;
points = gmshGeneratePointsTetrahedron(order, serendip);
if (order == 0) {
generateClosureOrder0(closures,24);
@@ -1566,7 +653,6 @@ nodalBasis::nodalBasis(int tag)
break;
case TYPE_PRI :
numFaces = 5;
- dimension = 3;
points = gmshGeneratePointsPrism(order, serendip);
if (order == 0) {
generateClosureOrder0(closures,48);
@@ -1580,14 +666,12 @@ nodalBasis::nodalBasis(int tag)
break;
case TYPE_HEX :
numFaces = 6;
- dimension = 3;
- points = gmshGeneratePointsHex(order, serendip);
+ points = gmshGeneratePointsHexahedron(order, serendip);
generateFaceClosureHex(closures, order, serendip, points);
generateFaceClosureHexFull(fullClosures, closureRef, order, serendip, points);
break;
case TYPE_PYR :
numFaces = 5;
- dimension = 3;
points = gmshGeneratePointsPyramid(order, serendip);
break;
}
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index 789783f948..8e5cf93ac8 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -14,9 +14,8 @@ class nodalBasis {
int type, parentType, order, dimension, numFaces;
bool serendip;
fullMatrix<double> points;
- fullMatrix<double> points_newAlgo;
- fullMatrix<int> monomials_newAlgo;
- fullMatrix<double> coefficients_newAlgo;
+ fullMatrix<int> monomials;
+ fullMatrix<double> coefficients;
nodalBasis(int tag);
virtual ~nodalBasis() {}
@@ -25,7 +24,6 @@ class nodalBasis {
// Basis functions & gradients evaluation
virtual void f(double u, double v, double w, double *sf) const = 0;
- virtual void fnew(double u, double v, double w, double *sf) const = 0;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
virtual void df(double u, double v, double w, double grads[][3]) const = 0;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index b56dea5312..d2c2b50571 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -218,6 +218,15 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoint
return monomials;
}
+/*
+00 10 20 30 40 â..
+01 11 21 31 41 â..
+02 12
+03 13
+04 14
+â. â.
+*/
+
fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order)
{
int nbMonomials = order ? order*4 : 1;
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index dc47ba7189..15a4501d97 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -45,7 +45,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = f
fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order);
fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
-//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //FIXME
+//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //TODO
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 24826d3845..9c1c98fb1b 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -45,361 +45,8 @@ static void printClosure(polynomialBasis::clCont &fullClosure,
*/
-
-fullMatrix<double> generate1DMonomials(int order)
-{
- fullMatrix<double> monomials(order + 1, 1);
- for (int i = 0; i < order + 1; i++) monomials(i, 0) = i;
- return monomials;
-}
-
-
-
-static fullMatrix<double> generatePascalTriangle(int order)
-{
- fullMatrix<double> monomials((order + 1) * (order + 2) / 2, 2);
- int index = 0;
- for (int i = 0; i <= order; i++) {
- for (int j = 0; j <= i; j++) {
- monomials(index, 0) = i - j;
- monomials(index, 1) = j;
- index++;
- }
- }
- return monomials;
-}
-
-
-
-// generate the exterior hull of the Pascal triangle for Serendipity element
-
-static fullMatrix<double> generatePascalSerendipityTriangle(int order)
-{
- fullMatrix<double> monomials(3 * order, 2);
-
- monomials(0, 0) = 0;
- monomials(0, 1) = 0;
-
- int index = 1;
- for (int i = 1; i <= order; i++) {
- if (i == order) {
- for (int j = 0; j <= i; j++) {
- monomials(index, 0) = i - j;
- monomials(index, 1) = j;
- index++;
- }
- }
- else {
- monomials(index, 0) = i;
- monomials(index, 1) = 0;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = i;
- index++;
- }
- }
- return monomials;
-}
-
-
-
-// generate all monomials xi^m * eta^n with n and m <= order
-
-static fullMatrix<double> generatePascalQuad(int order)
-{
-
- fullMatrix<double> monomials( (order+1)*(order+1), 2);
- int index = 0;
- for (int p = 0; p <= order; p++) {
- for(int i = 0; i < p; i++, index++) {
- monomials(index, 0) = p;
- monomials(index, 1) = i;
- }
- for(int i = 0; i <= p; i++, index++) {
- monomials(index, 0) = p-i;
- monomials(index, 1) = p;
- }
- }
- return monomials;
-}
-
-/*
-00 10 20 30 40 ⋯
-01 11 21 31 41 ⋯
-02 12
-03 13
-04 14
-â‹® â‹®
-*/
-
-
-
-// generate all monomials xi^m * eta^n with n and m <= order
-
-static fullMatrix<double> generatePascalHex(int order, bool serendip)
-{
- if (false && serendip) {
- fullMatrix<double> monomials( 8+(order-1)*12, 3);
- monomials(0,0)=0.;
- monomials(0,1)=0.;
- monomials(0,2)=0.;
- monomials(1,0)=1.;
- monomials(1,1)=0.;
- monomials(1,2)=0.;
- monomials(2,0)=0.;
- monomials(2,1)=1.;
- monomials(2,2)=0.;
- monomials(3,0)=1.;
- monomials(3,1)=1.;
- monomials(3,2)=0.;
-
- monomials(4,0)=0.;
- monomials(4,1)=0.;
- monomials(4,2)=1.;
- monomials(5,0)=1.;
- monomials(5,1)=0.;
- monomials(5,2)=1.;
- monomials(6,0)=0.;
- monomials(6,1)=1.;
- monomials(6,2)=1.;
- monomials(7,0)=1.;
- monomials(7,1)=1.;
- monomials(7,2)=1.;
- int index = 8;
- for (int p = 2; p <= order; p++) {
- monomials(index, 0) = p;
- monomials(index, 1) = 0;
- monomials(index, 2) = 0;
- index++;
- monomials(index, 0) = p;
- monomials(index, 1) = 1;
- monomials(index, 2) = 0;
- index++;
- monomials(index, 0) = p;
- monomials(index, 1) = 1;
- monomials(index, 2) = 1;
- index++;
- monomials(index, 0) = p;
- monomials(index, 1) = 0;
- monomials(index, 2) = 1;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = p;
- monomials(index, 2) = 0;
- index++;
- monomials(index, 0) = 1;
- monomials(index, 1) = p;
- monomials(index, 2) = 0;
- index++;
- monomials(index, 0) = 1;
- monomials(index, 1) = p;
- monomials(index, 2) = 1;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = p;
- monomials(index, 2) = 1;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = 0;
- monomials(index, 2) = p;
- index++;
- monomials(index, 0) = 1;
- monomials(index, 1) = 0;
- monomials(index, 2) = p;
- index++;
- monomials(index, 0) = 1;
- monomials(index, 1) = 1;
- monomials(index, 2) = p;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = 1;
- monomials(index, 2) = p;
- index++;
- }
- return monomials;
- }
-
- int siz = (order+1)*(order+1)*(order+1);
- if (serendip) siz -= (order-1)*(order-1)*(order-1);
- fullMatrix<double> monomials( siz, 3);
- int index = 0;
- for (int i = 0; i <= order; i++) {
- for (int j = 0; j <= order; j++) {
- for (int k = 0; k <= order; k++) {
- if (!serendip || i<2 || j<2 || k<2){
- monomials(index,0) = i;
- monomials(index,1) = j;
- monomials(index,2) = k;
- index++;
- }
- }
- }
- }
- return monomials;
-}
-
-
-
-static fullMatrix<double> generatePascalQuadSerendip(int order)
-{
- fullMatrix<double> monomials( (order)*4, 2);
- monomials(0,0)=0.;
- monomials(0,1)=0.;
- monomials(1,0)=1.;
- monomials(1,1)=0.;
- monomials(2,0)=0.;
- monomials(2,1)=1.;
- monomials(3,0)=1.;
- monomials(3,1)=1.;
- int index = 4;
- for (int p = 2; p <= order; p++) {
- monomials(index, 0) = p;
- monomials(index, 1) = 0;
- index++;
- monomials(index, 0) = 0;
- monomials(index, 1) = p;
- index++;
- monomials(index, 0) = p;
- monomials(index, 1) = 1;
- index++;
- monomials(index, 0) = 1;
- monomials(index, 1) = p;
- index++;
- }
- return monomials;
-}
-
-
-
-// generate the monomials subspace of all monomials of order exactly == p
-
-static fullMatrix<double> generateMonomialSubspace(int dim, int p)
-{
- fullMatrix<double> monomials;
-
- switch (dim) {
- case 1:
- monomials = fullMatrix<double>(1, 1);
- monomials(0, 0) = p;
- break;
- case 2:
- monomials = fullMatrix<double>(p + 1, 2);
- for (int k = 0; k <= p; k++) {
- monomials(k, 0) = p - k;
- monomials(k, 1) = k;
- }
- break;
- case 3:
- monomials = fullMatrix<double>((p + 1) * (p + 2) / 2, 3);
- int index = 0;
- for (int i = 0; i <= p; i++) {
- for (int k = 0; k <= p - i; k++) {
- monomials(index, 0) = p - i - k;
- monomials(index, 1) = k;
- monomials(index, 2) = i;
- index++;
- }
- }
- break;
- }
- return monomials;
-}
-
-
-
-// generate external hull of the Pascal tetrahedron
-
-static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
-{
- int nbMonomials = 4 + 6 * std::max(0, order - 1) +
- 4 * std::max(0, (order - 2) * (order - 1) / 2);
- fullMatrix<double> monomials(nbMonomials, 3);
-
- monomials.setAll(0);
- int index = 1;
- for (int p = 1; p < order; p++) {
- for (int i = 0; i < 3; i++) {
- int j = (i + 1) % 3;
- int k = (i + 2) % 3;
- for (int ii = 0; ii < p; ii++) {
- monomials(index, i) = p - ii;
- monomials(index, j) = ii;
- monomials(index, k) = 0;
- index++;
- }
- }
- }
- fullMatrix<double> monomialsMaxOrder = generateMonomialSubspace(3, order);
- int nbMaxOrder = monomialsMaxOrder.size1();
- monomials.copy(monomialsMaxOrder, 0, nbMaxOrder, 0, 3, index, 0);
- return monomials;
-}
-
-
-
-// generate Pascal tetrahedron
-
-static fullMatrix<double> generatePascalTetrahedron(int order)
-{
- int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
-
- fullMatrix<double> monomials(nbMonomials, 3);
-
- int index = 0;
- for (int p = 0; p <= order; p++) {
- fullMatrix<double> monOrder = generateMonomialSubspace(3, p);
- int nb = monOrder.size1();
- monomials.copy(monOrder, 0, nb, 0, 3, index, 0);
- index += nb;
- }
-
- return monomials;
-}
-
-
-
-// generate Pascal prism
-
-static fullMatrix<double> generatePascalPrism(int order)
-{
- int nbMonomials = (order + 1) * (order + 1) * (order + 2) / 2;
-
- fullMatrix<double> monomials(nbMonomials, 3);
- int index = 0;
- fullMatrix<double> lineMonoms = generate1DMonomials(order);
- fullMatrix<double> triMonoms = generatePascalTriangle(order);
- // store monomials in right order
- for (int currentOrder = 0; currentOrder <= order; currentOrder++) {
- int orderT = currentOrder, orderL = currentOrder;
- for (orderL = 0; orderL < currentOrder; orderL++) {
- // do all permutations of monoms for orderL, orderT
- int iL = orderL;
- for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
- monomials(index,0) = triMonoms(iT,0);
- monomials(index,1) = triMonoms(iT,1);
- monomials(index,2) = lineMonoms(iL,0);
- index ++;
- }
- }
- orderL = currentOrder;
- for (orderT = 0; orderT <= currentOrder; orderT++) {
- int iL = orderL;
- for (int iT = (orderT)*(orderT+1)/2; iT < (orderT+1)*(orderT+2)/2 ;iT++) {
- monomials(index,0) = triMonoms(iT,0);
- monomials(index,1) = triMonoms(iT,1);
- monomials(index,2) = lineMonoms(iL,0);
- index ++;
- }
- }
- }
-
- return monomials;
-}
-
-
-
static fullMatrix<double> generateLagrangeMonomialCoefficients
- (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
+ (const fullMatrix<int>& monomial, const fullMatrix<double>& point)
{
if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
@@ -415,7 +62,7 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
for (int i = 0; i < ndofs; i++) {
for (int j = 0; j < ndofs; j++) {
double dd = 1.;
- for (int k = 0; k < dim; k++) dd *= pow(point(j, k), monomial(i, k));
+ for (int k = 0; k < dim; k++) dd *= pow_int(point(j, k), monomial(i, k));
Vandermonde(i, j) = dd;
}
}
@@ -426,108 +73,37 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
return coefficient;
}
-static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
-{
- b.resize(orig.size1(), orig.size2());
- for (int i = 0; i < orig.size1(); ++i) {
- for (int j = 0; j < orig.size2(); ++j) {
- b(i, j) = static_cast<double>(orig(i, j));
- }
- }
-}
-
polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
{
switch (parentType) {
case TYPE_PNT :
- monomials = generate1DMonomials(0);
+ monomials = gmshGenerateMonomialsLine(0);
break;
case TYPE_LIN :
- monomials = generate1DMonomials(order);
+ monomials = gmshGenerateMonomialsLine(order);
break;
case TYPE_TRI :
- monomials = serendip ? generatePascalSerendipityTriangle(order) :
- generatePascalTriangle(order);
+ monomials = gmshGenerateMonomialsTriangle(order, serendip);
break;
case TYPE_QUA :
- monomials = serendip ? generatePascalQuadSerendip(order) :
- generatePascalQuad(order);
+ monomials = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
+ gmshGenerateMonomialsQuadrangle(order);
break;
case TYPE_TET :
- monomials = serendip ? generatePascalSerendipityTetrahedron(order) :
- generatePascalTetrahedron(order);
+ monomials = gmshGenerateMonomialsTetrahedron(order, serendip);
break;
case TYPE_PRI :
- monomials = generatePascalPrism(order);
+ monomials = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
+ gmshGenerateMonomialsPrism(order);
break;
case TYPE_HEX :
- monomials = generatePascalHex(order, serendip);
+ monomials = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
+ gmshGenerateMonomialsHexahedron(order);
break;
}
coefficients = generateLagrangeMonomialCoefficients(monomials, points);
-
-
- // NEW ALGO POINTS / MONOMIALS
-
- switch (parentType) {
- case TYPE_PNT :
- points_newAlgo = gmshGeneratePointsLine(0);
- monomials_newAlgo = gmshGenerateMonomialsLine(0);
- break;
- case TYPE_LIN :
- points_newAlgo = gmshGeneratePointsLine(order);
- monomials_newAlgo = gmshGenerateMonomialsLine(order);
- break;
- case TYPE_TRI :
- points_newAlgo = gmshGeneratePointsTriangle(order, serendip);
- monomials_newAlgo = gmshGenerateMonomialsTriangle(order, serendip);
- break;
- case TYPE_QUA :
- points_newAlgo = gmshGeneratePointsQuadrangle(order, serendip);
- if (!serendip)
- monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order);
- else
- monomials_newAlgo = gmshGenerateMonomialsQuadSerendipity(order);
- break;
- case TYPE_TET :
- points_newAlgo = gmshGeneratePointsTetrahedron(order, serendip);
- monomials_newAlgo = gmshGenerateMonomialsTetrahedron(order, serendip);
- break;
- case TYPE_HEX :
- points_newAlgo = gmshGeneratePointsHexahedron(order, serendip);
- if (!serendip)
- monomials_newAlgo = gmshGenerateMonomialsHexahedron(order);
- else
- monomials_newAlgo = gmshGenerateMonomialsHexaSerendipity(order);
- break;
- case TYPE_PRI :
- points_newAlgo = gmshGeneratePointsPrism(order, serendip);
- if (!serendip)
- monomials_newAlgo = gmshGenerateMonomialsPrism(order);
- else
- monomials_newAlgo = gmshGenerateMonomialsPrismSerendipity(order);
- break;
- }
-
- fullMatrix<double> monDouble;
- switch (parentType) {
- case TYPE_PNT :
- case TYPE_LIN :
- case TYPE_TRI :
- case TYPE_TET :
- copy(monomials_newAlgo, monDouble);
- break;
- case TYPE_QUA :
- case TYPE_HEX :
- case TYPE_PRI :
- //if (serendip) monDouble = monomials;
- /*else*/ copy(monomials_newAlgo, monDouble);
- break;
- }
-
- coefficients_newAlgo = generateLagrangeMonomialCoefficients(monDouble, points_newAlgo);
}
@@ -539,17 +115,14 @@ polynomialBasis::~polynomialBasis()
-void polynomialBasis::evaluateMonomialsNew(double u, double v, double w, double p[]) const
+void polynomialBasis::evaluateMonomials(double u, double v, double w, double p[]) const
{
- for (int j = 0; j < monomials_newAlgo.size1(); ++j) {
+ for (int j = 0; j < monomials.size1(); ++j) {
p[j] = 1.;
switch (dimension) {
- case 3 :
- p[j] *= pow_int(w, monomials_newAlgo(j, 2));
- case 2 :
- p[j] *= pow_int(v, monomials_newAlgo(j, 1));
- case 1 :
- p[j] *= pow_int(u, monomials_newAlgo(j, 0));
+ case 3 : p[j] *= pow_int(w, monomials(j, 2));
+ case 2 : p[j] *= pow_int(v, monomials(j, 1));
+ case 1 : p[j] *= pow_int(u, monomials(j, 0));
default :
break;
}
@@ -569,18 +142,6 @@ void polynomialBasis::f(double u, double v, double w, double *sf) const
}
}
}
-void polynomialBasis::fnew(double u, double v, double w, double *sf) const
-{
- double p[1256];
- evaluateMonomialsNew(u, v, w, p);
- for (int i = 0; i < coefficients_newAlgo.size1(); i++) {
- sf[i] = 0.0;
- for (int j = 0; j < coefficients_newAlgo.size2(); j++) {
- sf[i] += coefficients_newAlgo(i, j) * p[j];
- }
- }
-}
-
void polynomialBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index aa91e6594f..5e247e3f41 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -69,8 +69,6 @@ class polynomialBasis : public nodalBasis
public:
// for now the only implemented polynomial basis are nodal poly
// basis, we use the type of the corresponding gmsh element as type
- fullMatrix<double> monomials;
- fullMatrix<double> coefficients;
polynomialBasis(int tag);
~polynomialBasis();
@@ -78,21 +76,13 @@ class polynomialBasis : public nodalBasis
virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
virtual void f(double u, double v, double w, double *sf) const;
- virtual void fnew(double u, double v, double w, double *sf) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void ddf(double u, double v, double w, double hess[][3][3]) const;
virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
- inline void evaluateMonomials(double u, double v, double w, double p[]) const {
- for (int j = 0; j < monomials.size1(); j++) {
- p[j] = pow_int(u, (int)monomials(j, 0));
- if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
- if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
- }
- }
- void evaluateMonomialsNew(double u, double v, double w, double p[]) const;
+ inline void evaluateMonomials(double u, double v, double w, double p[]) const;
};
#endif
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 2a3efeea5e..069d027599 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -7,54 +7,6 @@
#include "pointsGenerators.h"
#include <cmath>
-
-static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
-{
- b.resize(orig.size1(), orig.size2());
- for (int i = 0; i < orig.size1(); ++i) {
- for (int j = 0; j < orig.size2(); ++j) {
- b(i, j) = static_cast<double>(orig(i, j));
- }
- }
-}
-
-static fullMatrix<double> generateLagrangeMonomialCoefficientsPyr
- (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
-{
- if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
- monomial.size1(),point.size1(),
- monomial.size2(),point.size2() );
- return fullMatrix<double>(1, 1);
- }
-
- int ndofs = monomial.size1();
- int dim = monomial.size2();
- fullMatrix<double> ppoint(point.size1(), point.size2());
-
- for (int i = 0; i < ndofs; ++i) {
- ppoint(i, 2) = 1. - point(i, 2);
- if (ppoint(i, 2) != .0) {
- ppoint(i, 0) = point(i, 0) / ppoint(i, 2);
- ppoint(i, 1) = point(i, 1) / ppoint(i, 2);
- }
- }
-
- fullMatrix<double> Vandermonde(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
- for (int k = 0; k < dim; k++) dd *= pow(ppoint(j, k), monomial(i, k));
- Vandermonde(i, j) = dd;
- }
- }
-
- fullMatrix<double> coefficient(ndofs, ndofs);
- Vandermonde.invert(coefficient);
-
- return coefficient;
-}
-
pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
{
@@ -62,7 +14,7 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
int num_points = points.size1();
- VDMinv.resize(num_points, num_points);
+ coefficients.resize(num_points, num_points);
double *fval = new double[num_points];
// Invert the Vandermonde matrix
@@ -71,19 +23,9 @@ 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(VDMinv);
+ VDM.invert(coefficients);
delete[] fval;
-
-
- // TEST NEW ALGO POINTS / MONOMIAL
-
- monomials_newAlgo = gmshGenerateMonomialsPyramid(order, serendip);
- points_newAlgo = gmshGeneratePointsPyramid(order, serendip);
-
- fullMatrix<double> monDouble;
- copy(monomials_newAlgo, monDouble);
- coefficients_newAlgo = generateLagrangeMonomialCoefficientsPyr(monDouble, points_newAlgo);
}
@@ -109,23 +51,12 @@ 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] += VDMinv(i,j)*fval[j];
+ for (int j = 0; j < N; j++) val[i] += coefficients(i,j)*fval[j];
}
delete[] fval;
}
-void pyramidalBasis::fnew(double u, double v, double w, double *sf) const
-{
- double p[1256];
- evaluateMonomialsNew(u, v, w, p);
- for (int i = 0; i < coefficients_newAlgo.size1(); i++) {
- sf[i] = 0.0;
- for (int j = 0; j < coefficients_newAlgo.size2(); j++) {
- sf[i] += coefficients_newAlgo(i, j) * p[j];
- }
- }
-}
@@ -141,7 +72,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) += VDMinv(i,j)*fval[j];
+ for (int j = 0; j < N; j++) sf(iPt,i) += coefficients(i,j)*fval[j];
}
}
@@ -162,9 +93,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] += VDMinv(i,j)*dfval[j][0];
- grads[i][1] += VDMinv(i,j)*dfval[j][1];
- grads[i][2] += VDMinv(i,j)*dfval[j][2];
+ 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];
}
}
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index 998f49d4d7..8b0982b02f 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -16,9 +16,6 @@
class pyramidalBasis: public nodalBasis
{
private:
- // Inverse of the Vandermonde matrix
- fullMatrix<double> VDMinv;
-
// Orthogonal basis for the pyramid
BergotBasis *bergot;
@@ -27,25 +24,11 @@ class pyramidalBasis: public nodalBasis
~pyramidalBasis();
virtual void f(double u, double v, double w, double *val) const;
- virtual void fnew(double u, double v, double w, double *val) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual int getNumShapeFunctions() const;
-
- inline void evaluateMonomialsNew(double u, double v, double w, double p[]) const {
- w = 1-w;
- if (w != .0) {
- u = u/w;
- v = v/w;
- }
- for (int j = 0; j < monomials_newAlgo.size1(); j++) {
- p[j] = pow_int(u, monomials_newAlgo(j, 0));
- if (dimension > 1) p[j] *= pow_int(v, monomials_newAlgo(j, 1));
- if (dimension > 2) p[j] *= pow_int(w, monomials_newAlgo(j, 2));
- }
- }
};
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 6bcc7c11af..1ce16372cd 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -122,6 +122,23 @@ void PViewData::setInterpolationMatrices(int type,
_interpolation[type].push_back(new fullMatrix<double>(expVal));
}
+void PViewData::setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<int> &expVal)
+{
+ if(!type || _interpolation[type].size()) return;
+
+ fullMatrix<double> *dExpVal;
+ dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
+ for (int i = 0; i < expVal.size1(); ++i) {
+ for (int j = 0; i < expVal.size2(); ++j) {
+ (*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
+ }
+ }
+ _interpolation[type].push_back(new fullMatrix<double>(coefVal));
+ _interpolation[type].push_back(dExpVal);
+}
+
void PViewData::setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal,
@@ -134,6 +151,53 @@ void PViewData::setInterpolationMatrices(int type,
_interpolation[type].push_back(new fullMatrix<double>(coefGeo));
_interpolation[type].push_back(new fullMatrix<double>(expGeo));
}
+void PViewData::setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<int> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<int> &expGeo)
+{
+ if(!type || _interpolation[type].size()) return;
+
+ fullMatrix<double> *dExpVal, *dExpGeo;
+ dExpVal = new fullMatrix<double>(expVal.size1(), expVal.size2());
+ dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
+
+ for (int i = 0; i < expVal.size1(); ++i) {
+ for (int j = 0; i < expVal.size2(); ++j) {
+ (*dExpVal)(i, j) = static_cast<double>(expVal(i, j));
+ }
+ }
+ for (int i = 0; i < expGeo.size1(); ++i) {
+ for (int j = 0; i < expGeo.size2(); ++j) {
+ (*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
+ }
+ }
+ _interpolation[type].push_back(new fullMatrix<double>(coefVal));
+ _interpolation[type].push_back(dExpVal);
+ _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
+ _interpolation[type].push_back(dExpGeo);
+}
+void PViewData::setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<int> &expGeo)
+{
+ if(!type || _interpolation[type].size()) return;
+
+ fullMatrix<double> *dExpGeo;
+ dExpGeo = new fullMatrix<double>(expGeo.size1(), expGeo.size2());
+ for (int i = 0; i < expGeo.size1(); ++i) {
+ for (int j = 0; i < expGeo.size2(); ++j) {
+ (*dExpGeo)(i, j) = static_cast<double>(expGeo(i, j));
+ }
+ }
+ _interpolation[type].push_back(new fullMatrix<double>(coefVal));
+ _interpolation[type].push_back(new fullMatrix<double>(expVal));
+ _interpolation[type].push_back(new fullMatrix<double>(coefGeo));
+ _interpolation[type].push_back(dExpGeo);
+}
int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p)
{
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 5136bf93d4..998863cd89 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -204,11 +204,24 @@ class PViewData {
void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal);
+ void setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<int> &expVal);
void setInterpolationMatrices(int type,
const fullMatrix<double> &coefVal,
const fullMatrix<double> &expVal,
const fullMatrix<double> &coefGeo,
const fullMatrix<double> &expGeo);
+ void setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<int> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<int> &expGeo);
+ void setInterpolationMatrices(int type,
+ const fullMatrix<double> &coefVal,
+ const fullMatrix<double> &expVal,
+ const fullMatrix<double> &coefGeo,
+ const fullMatrix<int> &expGeo);
int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
bool haveInterpolationMatrices(int type=0);
--
GitLab