diff --git a/CMakeLists.txt b/CMakeLists.txt
index 05fd38f7b36ab98f7a47898f0ad442d27f81b912..dd7ea8784e3b83e31e807ed8b956bb5ff64245c1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -685,6 +685,8 @@ if(ENABLE_MED OR ENABLE_CGNS)
endif(ZLIB_FOUND)
endif(NOT HAVE_LIBZ)
endif(MED_LIB OR CGNS_LIB)
+ else(HDF5_LIB)
+ message(STATUS "HDF5 not found")
endif(HDF5_LIB)
endif(ENABLE_MED OR ENABLE_CGNS)
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 93223fa539c2dc224e966a7111a5093078fb7835..8952556a28b2b5058de493e237ddcb8c93b30c7a 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -609,6 +609,7 @@ int GModel::readCGNS(const std::string &name)
double ielem = irmax[0] - 1;
double jelem = irmax[1] - 1;
double kelem = irmax[2] - 1;
+ printf("Elems %g %g %g\n", ielem, jelem, kelem);
int order = 1;
bool done = false;
while(fmod(ielem / 2.0, 1.0) == 0.0 && fmod(jelem / 2.0, 1.0) == 0.0 && fmod(kelem / 2.0, 1.0) == 0.0 and order < 5) {
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6be81f73551efbf2ef336a54881b9d88e986fcd3..9b7280cc1655405861c9e9a64299db117fdcaa1f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -147,7 +147,7 @@ void MElement::getNode(int num, double &u, double &v, double &w)
void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
{
- const polynomialBasis* fs = getFunctionSpace(o);
+ const nodalBasis* fs = getFunctionSpace(o);
if(fs) fs->f(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
@@ -155,7 +155,7 @@ void MElement::getShapeFunctions(double u, double v, double w, double s[], int o
void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
int o)
{
- const polynomialBasis* fs = getFunctionSpace(o);
+ const nodalBasis* fs = getFunctionSpace(o);
if(fs) fs->df(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
@@ -163,15 +163,15 @@ void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3]
void MElement::getHessShapeFunctions(double u, double v, double w, double s[][3][3],
int o)
{
- const polynomialBasis* fs = getFunctionSpace(o);
+ const nodalBasis* fs = getFunctionSpace(o);
if(fs) fs->ddf(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
-void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w,
+void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w,
double s[][3][3][3], int o)
{
- const polynomialBasis* fs = getFunctionSpace(o);
+ const nodalBasis* fs = getFunctionSpace(o);
if(fs) fs->dddf(u, v, w, s);
else Msg::Error("Function space not implemented for this type of element");
}
@@ -1230,6 +1230,13 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_5 : if(name) *name = "Pyramid 5"; return 5;
case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8;
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
+ case MSH_PYR_30 : if(name) *name = "Pyramid 30"; return 5 + 8*2 + 4*1 + 1*4 + 1;
+ case MSH_PYR_55 : if(name) *name = "Pyramid 55"; return 5 + 8*3 + 4*3 + 1*9 + 5;
+ case MSH_PYR_91 : if(name) *name = "Pyramid 91"; return 91;
+ case MSH_PYR_140 : if(name) *name = "Pyramid 140"; return 140;
+ case MSH_PYR_204 : if(name) *name = "Pyramid 204"; return 204;
+ case MSH_PYR_285 : if(name) *name = "Pyramid 285"; return 285;
+ case MSH_PYR_385 : if(name) *name = "Pyramid 385"; return 385;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
case MSH_PNT_SUB : if(name) *name = "Point Xfem"; return 1;
case MSH_LIN_SUB : if(name) *name = "Line Xfem"; return 2;
@@ -1320,6 +1327,96 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
return newEl;
}
+// Gives the parent type corresponding to
+// any element type.
+// Ex. : MSN_TRI_3 -> TYPE_TRI
+int MElement::ParentTypeFromTag(int tag)
+{
+ switch(tag) {
+ case(MSH_PNT):
+ return TYPE_PNT;
+ case(MSH_LIN_2): case(MSH_LIN_3):
+ case(MSH_LIN_4): case(MSH_LIN_5):
+ case(MSH_LIN_6): case(MSH_LIN_7):
+ case(MSH_LIN_8): case(MSH_LIN_9):
+ case(MSH_LIN_10): case(MSH_LIN_11):
+ case(MSH_LIN_B): case(MSH_LIN_C):
+ case(MSH_LIN_1):
+ return TYPE_LIN;
+ case(MSH_TRI_3): case(MSH_TRI_6):
+ case(MSH_TRI_9): case(MSH_TRI_10):
+ case(MSH_TRI_12): case(MSH_TRI_15):
+ case(MSH_TRI_15I): case(MSH_TRI_21):
+ case(MSH_TRI_28): case(MSH_TRI_36):
+ case(MSH_TRI_45): case(MSH_TRI_55):
+ case(MSH_TRI_66): case(MSH_TRI_18):
+ case(MSH_TRI_21I): case(MSH_TRI_24):
+ case(MSH_TRI_27): case(MSH_TRI_30):
+ case(MSH_TRI_B): case(MSH_TRI_1):
+ return TYPE_TRI;
+ case(MSH_QUA_4): case(MSH_QUA_9):
+ case(MSH_QUA_8): case(MSH_QUA_16):
+ case(MSH_QUA_25): case(MSH_QUA_36):
+ case(MSH_QUA_12): case(MSH_QUA_16I):
+ case(MSH_QUA_20): case(MSH_QUA_49):
+ case(MSH_QUA_64): case(MSH_QUA_81):
+ case(MSH_QUA_100): case(MSH_QUA_121):
+ case(MSH_QUA_24): case(MSH_QUA_28):
+ case(MSH_QUA_32): case(MSH_QUA_36I):
+ case(MSH_QUA_40): case(MSH_QUA_1):
+ return TYPE_QUA;
+ case(MSH_TET_4): case(MSH_TET_10):
+ case(MSH_TET_20): case(MSH_TET_35):
+ case(MSH_TET_56): case(MSH_TET_34):
+ case(MSH_TET_52): case(MSH_TET_84):
+ case(MSH_TET_120): case(MSH_TET_165):
+ case(MSH_TET_220): case(MSH_TET_286):
+ case(MSH_TET_74): case(MSH_TET_100):
+ case(MSH_TET_130): case(MSH_TET_164):
+ case(MSH_TET_202): case(MSH_TET_1):
+ return TYPE_TET;
+ case(MSH_PYR_5): case(MSH_PYR_14):
+ case(MSH_PYR_13): case(MSH_PYR_30):
+ case(MSH_PYR_55): case(MSH_PYR_91):
+ case(MSH_PYR_140): case(MSH_PYR_204):
+ case(MSH_PYR_285): case(MSH_PYR_385):
+ case(MSH_PYR_29): case(MSH_PYR_50):
+ case(MSH_PYR_77): case(MSH_PYR_110):
+ case(MSH_PYR_149): case(MSH_PYR_194):
+ case(MSH_PYR_245): case(MSH_PYR_1):
+ return TYPE_PYR;
+ case(MSH_PRI_6): case(MSH_PRI_18):
+ case(MSH_PRI_15): case(MSH_PRI_1):
+ case(MSH_PRI_40): case(MSH_PRI_75):
+ case(MSH_PRI_126): case(MSH_PRI_196):
+ case(MSH_PRI_288): case(MSH_PRI_405):
+ case(MSH_PRI_550): case(MSH_PRI_38):
+ case(MSH_PRI_66): case(MSH_PRI_102):
+ case(MSH_PRI_146): case(MSH_PRI_198):
+ case(MSH_PRI_258): case(MSH_PRI_326):
+ return TYPE_PRI;
+ case(MSH_HEX_8): case(MSH_HEX_27):
+ case(MSH_HEX_20): case(MSH_HEX_1):
+ case(MSH_HEX_64): case(MSH_HEX_125):
+ case(MSH_HEX_216): case(MSH_HEX_343):
+ case(MSH_HEX_512): case(MSH_HEX_729):
+ case(MSH_HEX_1000): case(MSH_HEX_56):
+ case(MSH_HEX_98): case(MSH_HEX_152):
+ case(MSH_HEX_222): case(MSH_HEX_296):
+ case(MSH_HEX_386): case(MSH_HEX_488):
+ return TYPE_HEX;
+ case(MSH_POLYG_): case(MSH_POLYG_B):
+ return TYPE_POLYG;
+ case(MSH_POLYH_):
+ return TYPE_POLYH;
+ default:
+ Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
+ return TYPE_TET;
+ }
+}
+
+
+
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part, bool owner, MElement *parent,
MElement *d1, MElement *d2)
@@ -1400,6 +1497,13 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_LIN_SUB: return new MSubLine(v, num, part, owner, parent);
case MSH_TRI_SUB: return new MSubTriangle(v, num, part, owner, parent);
case MSH_TET_SUB: return new MSubTetrahedron(v, num, part, owner, parent);
+ case MSH_PYR_30: return new MPyramidN(v, 3, num, part);
+ case MSH_PYR_55: return new MPyramidN(v, 4, num, part);
+ case MSH_PYR_91: return new MPyramidN(v, 5, num, part);
+ case MSH_PYR_140: return new MPyramidN(v, 6, num, part);
+ case MSH_PYR_204: return new MPyramidN(v, 7, num, part);
+ case MSH_PYR_285: return new MPyramidN(v, 8, num, part);
+ case MSH_PYR_385: return new MPyramidN(v, 9, num, part);
default: return 0;
}
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 67d48a6ca1d891c3923d1950a89f5a242d0c20b4..7abe28e7bbb911cc9e7b3ccb3c503fa9081f2830 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -14,6 +14,7 @@
#include "MVertex.h"
#include "MEdge.h"
#include "MFace.h"
+#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "JacobianBasis.h"
#include "GaussIntegration.h"
@@ -219,7 +220,7 @@ class MElement
virtual std::string getInfoString();
// get the function space for the element
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
+ virtual const nodalBasis* getFunctionSpace(int order=-1) const { return 0; }
// get the function space for the jacobian of the element
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; }
@@ -353,6 +354,8 @@ class MElement
std::map<MElement*, MElement*> &newParents,
std::map<MElement*, MElement*> &newDomains);
+ static int ParentTypeFromTag(int tag);
+
};
class MElementFactory{
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index ebfe8eab5ab99ee3d5b8622adf81fe70fe241752..fa143ddfb4d71e1a28abe70ad531ff8603a63211 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -117,7 +117,7 @@ class MPolyhedron : public MElement {
vol += _parts[i]->getVolume();
return vol;
}
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
+ virtual const nodalBasis* getFunctionSpace(int order=-1) const
{
return (_orig ? _orig->getFunctionSpace(order) : 0);
}
@@ -258,7 +258,7 @@ class MPolygon : public MElement {
virtual int getNumChildren() const { return _parts.size(); }
virtual MElement *getChild(int i) const { return _parts[i]; }
virtual bool ownsParent() const { return _owner; }
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
+ virtual const nodalBasis* getFunctionSpace(int order=-1) const
{
return (_orig ? _orig->getFunctionSpace(order) : 0);
}
@@ -323,7 +323,7 @@ class MLineChild : public MLine {
delete _orig;
}
virtual int getTypeForMSH() const { return MSH_LIN_C; }
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
+ virtual const nodalBasis* getFunctionSpace(int order=-1) const
{
if(_orig) return _orig->getFunctionSpace(order);
return 0;
@@ -354,6 +354,7 @@ class MLineChild : public MLine {
virtual bool ownsParent() const { return _owner; }
};
+
// -------------------- Border classes
class MTriangleBorder : public MTriangle {
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 070e21af38ddcd60ba398ac7b556dad4aa201885..38dbab818806c213a7b42583929e9425158e89c9 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -6,9 +6,12 @@
#include "MHexahedron.h"
#include "Numeric.h"
#include "Context.h"
+#include "BasisFactory.h"
+#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "MQuadrangle.h"
#include "qualityMeasures.h"
+
int MHexahedron::getVolumeSign()
{
double mat[3][3];
@@ -162,8 +165,8 @@ int MHexahedronN::getNumEdgesRep()
{
return 12 * CTX::instance()->mesh.numSubEdges;
}
-
-const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
+
+const nodalBasis* MHexahedron::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -171,31 +174,31 @@ const polynomialBasis* MHexahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return polynomialBases::find(MSH_HEX_1);
- case 1: return polynomialBases::find(MSH_HEX_8);
- case 2: return polynomialBases::find(MSH_HEX_20);
- case 3: return polynomialBases::find(MSH_HEX_56);
- case 4: return polynomialBases::find(MSH_HEX_98);
- case 5: return polynomialBases::find(MSH_HEX_152);
- case 6: return polynomialBases::find(MSH_HEX_222);
- case 7: return polynomialBases::find(MSH_HEX_296);
- case 8: return polynomialBases::find(MSH_HEX_386);
- case 9: return polynomialBases::find(MSH_HEX_488);
+ case 0: return BasisFactory::create(MSH_HEX_1);
+ case 1: return BasisFactory::create(MSH_HEX_8);
+ case 2: return BasisFactory::create(MSH_HEX_20);
+ case 3: return BasisFactory::create(MSH_HEX_56);
+ case 4: return BasisFactory::create(MSH_HEX_98);
+ case 5: return BasisFactory::create(MSH_HEX_152);
+ case 6: return BasisFactory::create(MSH_HEX_222);
+ case 7: return BasisFactory::create(MSH_HEX_296);
+ case 8: return BasisFactory::create(MSH_HEX_386);
+ case 9: return BasisFactory::create(MSH_HEX_488);
default: Msg::Error("Order %d hex function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return polynomialBases::find(MSH_HEX_1);
- case 1: return polynomialBases::find(MSH_HEX_8);
- case 2: return polynomialBases::find(MSH_HEX_27);
- case 3: return polynomialBases::find(MSH_HEX_64);
- case 4: return polynomialBases::find(MSH_HEX_125);
- case 5: return polynomialBases::find(MSH_HEX_216);
- case 6: return polynomialBases::find(MSH_HEX_343);
- case 7: return polynomialBases::find(MSH_HEX_512);
- case 8: return polynomialBases::find(MSH_HEX_729);
- case 9: return polynomialBases::find(MSH_HEX_1000);
+ case 0: return BasisFactory::create(MSH_HEX_1);
+ case 1: return BasisFactory::create(MSH_HEX_8);
+ case 2: return BasisFactory::create(MSH_HEX_27);
+ case 3: return BasisFactory::create(MSH_HEX_64);
+ case 4: return BasisFactory::create(MSH_HEX_125);
+ case 5: return BasisFactory::create(MSH_HEX_216);
+ case 6: return BasisFactory::create(MSH_HEX_343);
+ case 7: return BasisFactory::create(MSH_HEX_512);
+ case 8: return BasisFactory::create(MSH_HEX_729);
+ case 9: return BasisFactory::create(MSH_HEX_1000);
default: Msg::Error("Order %d hex function space not implemented", order);
}
}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 4cadf664f2197904838c92b82453a9291a436012..d089dbae6141ed459a5c29144bdfd73f26d98860 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -57,7 +57,7 @@ class MHexahedron : public MElement {
virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return _v[num]; }
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
virtual MVertex *getVertexDIFF(int num)
{
static const int map[8] = {2, 3, 7, 6, 0, 1, 5, 4};
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 3ffe7f1e106e9efdb02c53b42d3a0f6c0740cf57..52d949135342946cfddef890c1fa71e8aa915a86 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -5,26 +5,28 @@
#include "GmshDefines.h"
#include "MLine.h"
+#include "nodalBasis.h"
+#include "BasisFactory.h"
#include "GaussLegendre1D.h"
#include "Context.h"
#include "qualityMeasures.h"
-const polynomialBasis* MLine::getFunctionSpace(int o) const
+const nodalBasis* MLine::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 0: return polynomialBases::find(MSH_LIN_1);
- case 1: return polynomialBases::find(MSH_LIN_2);
- case 2: return polynomialBases::find(MSH_LIN_3);
- case 3: return polynomialBases::find(MSH_LIN_4);
- case 4: return polynomialBases::find(MSH_LIN_5);
- case 5: return polynomialBases::find(MSH_LIN_6);
- case 6: return polynomialBases::find(MSH_LIN_7);
- case 7: return polynomialBases::find(MSH_LIN_8);
- case 8: return polynomialBases::find(MSH_LIN_9);
- case 9: return polynomialBases::find(MSH_LIN_10);
- case 10: return polynomialBases::find(MSH_LIN_11);
+ case 0: return BasisFactory::create(MSH_LIN_1);
+ case 1: return BasisFactory::create(MSH_LIN_2);
+ case 2: return BasisFactory::create(MSH_LIN_3);
+ case 3: return BasisFactory::create(MSH_LIN_4);
+ case 4: return BasisFactory::create(MSH_LIN_5);
+ case 5: return BasisFactory::create(MSH_LIN_6);
+ case 6: return BasisFactory::create(MSH_LIN_7);
+ case 7: return BasisFactory::create(MSH_LIN_8);
+ case 8: return BasisFactory::create(MSH_LIN_9);
+ case 9: return BasisFactory::create(MSH_LIN_10);
+ case 10: return BasisFactory::create(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 68a22428c501e3b84fcea1e7aa42dcf27d7787b3..0737024bef66db8b72561a58dbdc95ac10ae4475 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -7,6 +7,7 @@
#define _MLINE_H_
#include "MElement.h"
+#include "nodalBasis.h"
/*
* MLine
@@ -71,7 +72,7 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual bool isInside(double u, double v, double w)
{
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 0ae5009b7a03c13457917f57a721a07b65f5d352..7cc20180db60102d2cb438809e1d1857a49832b9 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -5,6 +5,9 @@
#include "MPyramid.h"
#include "Numeric.h"
+#include "Context.h"
+#include "BergotBasis.h"
+#include "BasisFactory.h"
int MPyramid::getVolumeSign()
{
@@ -23,3 +26,337 @@ int MPyramid::getVolumeSign()
else if(d > 0.) return 1;
else return 0;
}
+
+const nodalBasis* MPyramid::getFunctionSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return BasisFactory::create(MSH_PYR_5);
+ case 2: return BasisFactory::create(MSH_PYR_14);
+ case 0: return BasisFactory::create(MSH_PYR_1);
+ case 3: return BasisFactory::create(MSH_PYR_29);
+ case 4: return BasisFactory::create(MSH_PYR_50);
+ case 5: return BasisFactory::create(MSH_PYR_77);
+ case 6: return BasisFactory::create(MSH_PYR_110);
+ case 7: return BasisFactory::create(MSH_PYR_149);
+ case 8: return BasisFactory::create(MSH_PYR_194);
+ case 9: return BasisFactory::create(MSH_PYR_245);
+ default: Msg::Error("Order %d pyramid function space not implemented", order);
+ }
+ }
+ else {
+ switch (order) {
+ case 0: return BasisFactory::create(MSH_PYR_1);
+ case 1: return BasisFactory::create(MSH_PYR_5);
+ case 2: return BasisFactory::create(MSH_PYR_14);
+ case 3: return BasisFactory::create(MSH_PYR_30);
+ case 4: return BasisFactory::create(MSH_PYR_55);
+ case 5: return BasisFactory::create(MSH_PYR_91);
+ case 6: return BasisFactory::create(MSH_PYR_140);
+ case 7: return BasisFactory::create(MSH_PYR_204);
+ case 8: return BasisFactory::create(MSH_PYR_285);
+ case 9: return BasisFactory::create(MSH_PYR_385);
+ default: Msg::Error("Order %d pyramid function space not implemented", order);
+ }
+ }
+ return 0;
+}
+
+
+
+
+MPyramidN::~MPyramidN() {}
+
+double MPyramidN::distoShapeMeasure()
+{
+#if defined(HAVE_MESH)
+ _disto = 0.0; //qmDistorsionOfMapping(this);
+#else
+ _disto = 0.;
+#endif
+ return _disto;
+}
+
+int MPyramidN::getNumEdgesRep(){ return 8 * CTX::instance()->mesh.numSubEdges; }
+
+void MPyramidN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ int numSubEdges = CTX::instance()->mesh.numSubEdges;
+ static double pp[5][3] = {{-1,-1,0},{1,-1,0},{1,1,0},{-1,1,0},{0,0,1}};
+ static int ed [8][2] = {{0,1},{0,3},{0,4},{1,2},{1,4},{2,3},{2,4},{3,4}};
+ int iEdge = num / numSubEdges;
+ int iSubEdge = num % numSubEdges;
+
+ int iVertex1 = ed [iEdge][0];
+ int iVertex2 = ed [iEdge][1];
+ double t1 = (double) iSubEdge / (double) numSubEdges;
+ double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+
+ double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+ double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+
+ SPoint3 pnt1, pnt2;
+ pnt(u1,v1,w1,pnt1);
+ pnt(u2,v2,w2,pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
+
+ // not great, but better than nothing
+ static const int f[8] = {0, 1, 0, 2, 2, 3, 2, 3};
+ n[0] = n[1] = getFace(f[iEdge]).normal();
+}
+
+
+int MPyramidN::getNumFacesRep(){ return 6 * SQU(CTX::instance()->mesh.numSubEdges); }
+
+static void _myGetFaceRep(MPyramid *pyr, int num, double *x, double *y, double *z,
+ SVector3 *n, int numSubEdges)
+{
+ static double pp[5][3] = {
+ {-1,-1,0},
+ {1,-1,0},
+ {1,1,0},
+ {-1,1,0},
+ {0,0,1}
+ };
+ int iFace = num / (numSubEdges * numSubEdges);
+ int iSubFace = num % (numSubEdges * numSubEdges);
+
+ if (iFace > 3) {
+ iFace = 4;
+ iSubFace = num % (2*numSubEdges * numSubEdges);
+ }
+
+ int iVertex1, iVertex2, iVertex3, iVertex4;
+
+ if (iFace < 4) {
+ iVertex1 = pyr->faces_pyramid(iFace,0);
+ iVertex2 = pyr->faces_pyramid(iFace,1);
+ iVertex3 = pyr->faces_pyramid(iFace,2);
+ } else {
+ iVertex1 = 0;
+ iVertex2 = 1;
+ iVertex3 = 2;
+ iVertex4 = 3;
+ }
+
+ if (iFace < 4) {
+
+ int ix = 0, iy = 0;
+ int nbt = 0;
+ for (int i = 0; i < numSubEdges; i++){
+ int nbl = (numSubEdges - i - 1) * 2 + 1;
+ nbt += nbl;
+ if (nbt > iSubFace){
+ iy = i;
+ ix = nbl - (nbt - iSubFace);
+ break;
+ }
+ }
+
+ const double d = 1. / numSubEdges;
+
+ SPoint3 pnt1, pnt2, pnt3;
+ double u1, v1, u2, v2, u3, v3;
+ if (ix % 2 == 0){
+ u1 = ix / 2 * d; v1= iy*d;
+ u2 = (ix / 2 + 1) * d ; v2 = iy * d;
+ u3 = ix / 2 * d ; v3 = (iy+1) * d;
+ }
+ else{
+ u1 = (ix / 2 + 1) * d; v1= iy * d;
+ u2 = (ix / 2 + 1) * d; v2= (iy + 1) * d;
+ u3 = ix / 2 * d ; v3 = (iy + 1) * d;
+ }
+
+ double U1 = pp[iVertex1][0] * (1.-u1-v1) + pp[iVertex2][0] * u1 + pp[iVertex3][0] * v1;
+ double U2 = pp[iVertex1][0] * (1.-u2-v2) + pp[iVertex2][0] * u2 + pp[iVertex3][0] * v2;
+ double U3 = pp[iVertex1][0] * (1.-u3-v3) + pp[iVertex2][0] * u3 + pp[iVertex3][0] * v3;
+
+ double V1 = pp[iVertex1][1] * (1.-u1-v1) + pp[iVertex2][1] * u1 + pp[iVertex3][1] * v1;
+ double V2 = pp[iVertex1][1] * (1.-u2-v2) + pp[iVertex2][1] * u2 + pp[iVertex3][1] * v2;
+ double V3 = pp[iVertex1][1] * (1.-u3-v3) + pp[iVertex2][1] * u3 + pp[iVertex3][1] * v3;
+
+ double W1 = pp[iVertex1][2] * (1.-u1-v1) + pp[iVertex2][2] * u1 + pp[iVertex3][2] * v1;
+ double W2 = pp[iVertex1][2] * (1.-u2-v2) + pp[iVertex2][2] * u2 + pp[iVertex3][2] * v2;
+ double W3 = pp[iVertex1][2] * (1.-u3-v3) + pp[iVertex2][2] * u3 + pp[iVertex3][2] * v3;
+
+ pyr->pnt(U1, V1, W1, pnt1);
+ pyr->pnt(U2, V2, W2, pnt2);
+ pyr->pnt(U3, V3, W3, pnt3);
+
+ x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
+
+ SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
+ SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
+ n[0] = crossprod(d1, d2);
+ n[0].normalize();
+ n[1] = n[0];
+ n[2] = n[0];
+
+ } else {
+
+
+ SPoint3 pnt1, pnt2, pnt3;
+ double J1[3][3], J2[3][3], J3[3][3];
+
+ /*
+ 0
+ 0 1
+ 0 1 2
+ 0 1 2 3
+ 0 1 2 3 4
+ 0 1 2 3 4 5
+ */
+
+ // on each layer, we have (numSubEdges) * 2 triangles
+ // ix and iy are the coordinates of the sub-quadrangle
+
+ int io = iSubFace%2;
+ int ix = (iSubFace/2)/numSubEdges;
+ int iy = (iSubFace/2)%numSubEdges;
+
+ const double d = 2. / numSubEdges;
+ double ox = -1. + d*ix;
+ double oy = -1. + d*iy;
+
+ if (io == 0){
+ double U1 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V1 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W1 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox += d;
+
+ double U2 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V2 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W2 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ oy += d;
+
+ double U3 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V3 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W3 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ pyr->pnt(U1,V1,W1, pnt1);
+ pyr->pnt(U2,V2,W2, pnt2);
+ pyr->pnt(U3,V3,W3, pnt3);
+ }
+ else{
+ double U1 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V1 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W1 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox += d;
+ oy += d;
+
+ double U2 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V2 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W2 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox -= d;
+
+ double U3 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V3 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W3 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ pyr->pnt(U1,V1,W1, pnt1);
+ pyr->pnt(U2,V2,W2, pnt2);
+ pyr->pnt(U3,V3,W3, pnt3);
+ }
+
+ n[0] = 1;
+ n[1] = 1;
+ n[2] = 1;
+ x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
+ }
+}
+
+void MPyramidN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 065d1186dc74602c09b4c6d561e01984d6ebb580..1277b1d75533e7867ba7605d3c01cc6964c196df 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -67,6 +67,7 @@ class MPyramid : public MElement {
virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 5; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
virtual int getNumEdges(){ return 8; }
virtual MEdge getEdge(int num)
{
@@ -122,7 +123,7 @@ class MPyramid : public MElement {
MVertex *tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
}
virtual int getVolumeSign();
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+ /*virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
double r = (w != 1.) ? (u * v * w / (1. - w)) : 0.;
s[0] = 0.25 * ((1. - u) * (1. - v) - w + r);
@@ -131,6 +132,7 @@ class MPyramid : public MElement {
s[3] = 0.25 * ((1. - u) * (1. + v) - w - r);
s[4] = w;
}
+ */
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
{
if(w == 1.) {
@@ -450,4 +452,105 @@ class MPyramid14 : public MPyramid {
}
};
+//------------------------------------------------------------------------------
+
+class MPyramidN : public MPyramid {
+
+ protected:
+ std::vector<MVertex*> _vs;
+ const char _order;
+ double _disto;
+ public:
+ MPyramidN(MVertex* v0, MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MPyramid(v0, v1, v2, v3, v4, num, part), _vs(v), _order(order), _disto(-1e22)
+ {
+ for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ getFunctionSpace(order);
+ }
+
+ MPyramidN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MPyramid(v[0], v[1], v[2], v[3], v[4], num, part), _order(order), _disto(-1e22)
+ {
+ for (unsigned int i = 5; i < v.size(); i++ ) _vs.push_back(v[i]);
+ for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ getFunctionSpace(order);
+ }
+
+ ~MPyramidN();
+
+ virtual double distoShapeMeasure();
+ virtual int getPolynomialOrder() const { return _order; }
+ virtual int getNumVertices() const { return 5 + _vs.size(); }
+ virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
+ virtual int getNumEdgeVertices() const { return 8 * (_order - 1); }
+ virtual int getNumFaceVertices() const
+ {
+ return (_order-1)*(_order-1) + 4 * ((_order - 1) * (_order - 2)) / 2;
+ }
+ virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ v.resize(_order+1);
+ MPyramid::_getEdgeVertices(num, v);
+ int j = 2;
+ const int ie = (num + 1) * (_order - 1);
+ for(int i = num * (_order -1); i != ie; ++i) v[j++] = _vs[i];
+ }
+ virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
+ {
+ int j = 3;
+ if (num == 4) {
+ j = 4;
+ v.resize(_order * _order);
+ }
+ else {
+ v.resize(3 + 3 * (_order - 1) + (_order-1) * (_order - 2) /2);
+ }
+ MPyramid::_getFaceVertices(num, v);
+ int nbVQ = (_order-1)*(_order-1);
+ int nbVT = (_order - 1) * (_order - 2) / 2;
+ const int ie = (num == 4) ? 4*nbVT + nbVQ : (num+1)*nbVT;
+ for (int i = num*nbVT; i != ie; ++i) v[j++] = _vs[i];
+ }
+ virtual int getNumVolumeVertices() const
+ {
+ switch(getTypeForMSH()){
+ case MSH_PYR_30 : return 1;
+ case MSH_PYR_55 : return 5;
+ case MSH_PYR_91 : return 14;
+ case MSH_PYR_140 : return 30;
+ case MSH_PYR_204 : return 55;
+ case MSH_PYR_285 : return 91;
+ case MSH_PYR_385 : return 140;
+ default : return 0;
+ }
+ }
+ virtual int getTypeForMSH() const
+ {
+ if(_order == 3 && _vs.size() + 5 == 29) return MSH_PYR_29;
+ if(_order == 3 && _vs.size() + 5 == 30) return MSH_PYR_30;
+ if(_order == 4 && _vs.size() + 5 == 50) return MSH_PYR_50;
+ if(_order == 4 && _vs.size() + 5 == 55) return MSH_PYR_55;
+ if(_order == 5 && _vs.size() + 5 == 77) return MSH_PYR_77;
+ if(_order == 5 && _vs.size() + 5 == 91) return MSH_PYR_91;
+ if(_order == 6 && _vs.size() + 5 == 110) return MSH_PYR_110;
+ if(_order == 6 && _vs.size() + 5 == 140) return MSH_PYR_140;
+ if(_order == 7 && _vs.size() + 5 == 149) return MSH_PYR_149;
+ if(_order == 7 && _vs.size() + 5 == 204) return MSH_PYR_204;
+ if(_order == 8 && _vs.size() + 5 == 194) return MSH_PYR_194;
+ if(_order == 8 && _vs.size() + 5 == 285) return MSH_PYR_285;
+ if(_order == 9 && _vs.size() + 5 == 245) return MSH_PYR_245;
+ if(_order == 9 && _vs.size() + 5 == 385) return MSH_PYR_385;
+ return 0;
+ }
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumEdgesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumFacesRep();
+ virtual void getNode(int num, double &u, double &v, double &w)
+ {
+ num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+ }
+};
+
#endif
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 7640a3eb068d80bb52e076eafe1390e9b805e0d9..60c4e82e15a4f5797dd2ea448b192dd1fc653c39 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -18,7 +18,7 @@ MSubTetrahedron::~MSubTetrahedron()
const polynomialBasis* MSubTetrahedron::getFunctionSpace(int order) const
{
- if(_orig) return _orig->getFunctionSpace(order);
+ if(_orig) return (polynomialBasis*)_orig->getFunctionSpace(order);
return 0;
}
@@ -116,7 +116,7 @@ MSubTriangle::~MSubTriangle()
const polynomialBasis* MSubTriangle::getFunctionSpace(int order) const
{
- if(_orig) return _orig->getFunctionSpace(order);
+ if(_orig) return (polynomialBasis*)_orig->getFunctionSpace(order);
return 0;
}
@@ -211,7 +211,7 @@ MSubLine::~MSubLine()
const polynomialBasis* MSubLine::getFunctionSpace(int order) const
{
- if(_orig) return _orig->getFunctionSpace(order);
+ if(_orig) return (polynomialBasis*)_orig->getFunctionSpace(order);
return 0;
}
@@ -301,7 +301,7 @@ MSubPoint::~MSubPoint()
const polynomialBasis* MSubPoint::getFunctionSpace(int order) const
{
- if(_orig) return _orig->getFunctionSpace(order);
+ if(_orig) return (polynomialBasis*)_orig->getFunctionSpace(order);
return 0;
}
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 7d0cd1509838483e6f1a12328628018598bd71de..1dde8a9e9a5cd385dfd67fde9dc7241e74948d00 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -7,6 +7,7 @@
#include "MTetrahedron.h"
#include "Numeric.h"
#include "Context.h"
+#include "BasisFactory.h"
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
@@ -184,7 +185,7 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
sys3x3(mat, b, uvw, &det);
}
-const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
+const nodalBasis* MTetrahedron::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
@@ -192,33 +193,33 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return polynomialBases::find(MSH_TET_1);
- case 1: return polynomialBases::find(MSH_TET_4);
- case 2: return polynomialBases::find(MSH_TET_10);
- case 3: return polynomialBases::find(MSH_TET_20);
- case 4: return polynomialBases::find(MSH_TET_34);
- case 5: return polynomialBases::find(MSH_TET_52);
- case 6: return polynomialBases::find(MSH_TET_74);
- case 7: return polynomialBases::find(MSH_TET_100);
- case 8: return polynomialBases::find(MSH_TET_130);
- case 9: return polynomialBases::find(MSH_TET_164);
- case 10: return polynomialBases::find(MSH_TET_202);
+ case 0: return BasisFactory::create(MSH_TET_1);
+ case 1: return BasisFactory::create(MSH_TET_4);
+ case 2: return BasisFactory::create(MSH_TET_10);
+ case 3: return BasisFactory::create(MSH_TET_20);
+ case 4: return BasisFactory::create(MSH_TET_34);
+ case 5: return BasisFactory::create(MSH_TET_52);
+ case 6: return BasisFactory::create(MSH_TET_74);
+ case 7: return BasisFactory::create(MSH_TET_100);
+ case 8: return BasisFactory::create(MSH_TET_130);
+ case 9: return BasisFactory::create(MSH_TET_164);
+ case 10: return BasisFactory::create(MSH_TET_202);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return polynomialBases::find(MSH_TET_1);
- case 1: return polynomialBases::find(MSH_TET_4);
- case 2: return polynomialBases::find(MSH_TET_10);
- case 3: return polynomialBases::find(MSH_TET_20);
- case 4: return polynomialBases::find(MSH_TET_35);
- case 5: return polynomialBases::find(MSH_TET_56);
- case 6: return polynomialBases::find(MSH_TET_84);
- case 7: return polynomialBases::find(MSH_TET_120);
- case 8: return polynomialBases::find(MSH_TET_165);
- case 9: return polynomialBases::find(MSH_TET_220);
- case 10: return polynomialBases::find(MSH_TET_286);
+ case 0: return BasisFactory::create(MSH_TET_1);
+ case 1: return BasisFactory::create(MSH_TET_4);
+ case 2: return BasisFactory::create(MSH_TET_10);
+ case 3: return BasisFactory::create(MSH_TET_20);
+ case 4: return BasisFactory::create(MSH_TET_35);
+ case 5: return BasisFactory::create(MSH_TET_56);
+ case 6: return BasisFactory::create(MSH_TET_84);
+ case 7: return BasisFactory::create(MSH_TET_120);
+ case 8: return BasisFactory::create(MSH_TET_165);
+ case 9: return BasisFactory::create(MSH_TET_220);
+ case 10: return BasisFactory::create(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 0be96dcf11df3c5ff7a8977a7087c48ac1e0e879..2df127b1640bc3c12d7ebd9bae6fd3b709705576 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -128,7 +128,7 @@ class MTetrahedron : public MElement {
virtual double getCircumRadius();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual void getNode(int num, double &u, double &v, double &w)
{
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index f37cb864725d1760500d0e7821d7be03fcedf98b..746bff9db943ae6c90193e463ccdada64a704174 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -22,6 +22,7 @@
#include "Context.h"
#include "fullMatrix.h"
#include "polynomialBasis.h"
+#include "BasisFactory.h"
#define SQU(a) ((a)*(a))
@@ -723,14 +724,16 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
switch (nPts){
case 0: return;
case 1: return;
- case 2: points = polynomialBases::find(MSH_TET_20)->points; break;
- case 3: points = polynomialBases::find(MSH_TET_35)->points; break;
- case 4: points = polynomialBases::find(MSH_TET_56)->points; break;
- case 5: points = polynomialBases::find(MSH_TET_84)->points; break;
- case 6: points = polynomialBases::find(MSH_TET_120)->points; break;
- case 7: points = polynomialBases::find(MSH_TET_165)->points; break;
- case 8: points = polynomialBases::find(MSH_TET_220)->points; break;
- case 9: points = polynomialBases::find(MSH_TET_286)->points; break;
+ case 2:
+ BasisFactory::create(MSH_TET_20)->points.print();
+ points = BasisFactory::create(MSH_TET_20)->points; break;
+ case 3: points = BasisFactory::create(MSH_TET_35)->points; break;
+ case 4: points = BasisFactory::create(MSH_TET_56)->points; break;
+ case 5: points = BasisFactory::create(MSH_TET_84)->points; break;
+ case 6: points = BasisFactory::create(MSH_TET_120)->points; break;
+ case 7: points = BasisFactory::create(MSH_TET_165)->points; break;
+ case 8: points = BasisFactory::create(MSH_TET_220)->points; break;
+ case 9: points = BasisFactory::create(MSH_TET_286)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
@@ -755,17 +758,16 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ;
break;
case TYPE_PYR:
- printf("aaaaaand...\n");
switch (nPts){
case 0:
case 1: return;
- case 2: points = polynomialBases::find(MSH_PYR_30)->points; break;
- case 3: points = polynomialBases::find(MSH_PYR_55)->points; break;
- case 4: points = polynomialBases::find(MSH_PYR_91)->points; break;
- case 5: points = polynomialBases::find(MSH_PYR_140)->points; break;
- case 6: points = polynomialBases::find(MSH_PYR_204)->points; break;
- case 7: points = polynomialBases::find(MSH_PYR_285)->points; break;
- case 8: points = polynomialBases::find(MSH_PYR_385)->points; break;
+ case 2: points = BasisFactory::create(MSH_PYR_30)->points; break;
+ case 3: points = BasisFactory::create(MSH_PYR_55)->points; break;
+ case 4: points = BasisFactory::create(MSH_PYR_91)->points; break;
+ case 5: points = BasisFactory::create(MSH_PYR_140)->points; break;
+ case 6: points = BasisFactory::create(MSH_PYR_204)->points; break;
+ case 7: points = BasisFactory::create(MSH_PYR_285)->points; break;
+ case 8: points = BasisFactory::create(MSH_PYR_385)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
@@ -781,9 +783,11 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
const double t1 = points(k, 0);
const double t2 = points(k, 1);
const double t3 = points(k, 2);
+ //printf("inside getRegionVertices, point is %g %g %g\n", t1, t2, t3);
SPoint3 pos;
incomplete->pnt(t1, t2, t3, pos);
v = new MVertex(pos.x(), pos.y(), pos.z(), gr);
+ //printf("inside getRegionVertices, vertex is %g %g %g\n", v->x(), v->y(), v->z());
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
@@ -1011,39 +1015,167 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
gr->prisms = prisms2;
std::vector<MPyramid*> pyramids2;
- for(unsigned int i = 0; i < gr->pyramids.size(); i++){
+ for(unsigned int i = 0; i < gr->pyramids.size(); i++) {
MPyramid *p = gr->pyramids[i];
std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
if (nPts == 1) {
- if(incomplete){
- pyramids2.push_back
- (new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7]));
+ if(incomplete) {
+ pyramids2.push_back(new MPyramid13(p->getVertex(0), p->getVertex(1),
+ p->getVertex(2), p->getVertex(3),
+ p->getVertex(4), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7]));
}
- else{
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- pyramids2.push_back
- (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
+ else {
+ getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ pyramids2.push_back
+ (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
}
}
else {
- Msg::Error("PyramidN generation not implemented");
- /*
getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
- MPyramidN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve, nPts + 1);
- getRegionVertices(gr, &incpl, p, vr, linear, nPts);
+
+ // Creating quads to get the internal vertices
+ /*
+ * P3 : q1 - 21 22 23 24
+ * P4 : q1 - 29 30 32 33 35 36 38 39
+ * q2 - 31 34 37 40
+ * P5 : q1 - 37 40 38 43 46 44 49 52 50 55 58 56
+ * q2 - 42 41 48 47 54 53 60 59
+ * q3 - 39 45 51 57
+ */
+
+ vr.reserve((nPts-1)*(nPts)*(2*(nPts-1)+1)/6);
+ /*for (int tmp = 0; tmp < (nPts-1)*(nPts)*(2*(nPts-1)+1)/6; tmp++) {
+ vr.push_back(0);
+ }*/
+
+ int verts_lvl3[12] = {37,40,38,43,46,44,49,52,50,55,58,56};
+
+ int verts_lvl2[8];
+ if (nPts == 4) {
+ verts_lvl2[0] = 42; verts_lvl2[1] = 41;
+ verts_lvl2[2] = 48; verts_lvl2[3] = 47;
+ verts_lvl2[4] = 54; verts_lvl2[5] = 53;
+ verts_lvl2[6] = 60; verts_lvl2[7] = 59;
+ } else {
+ verts_lvl2[0] = 29; verts_lvl2[1] = 30;
+ verts_lvl2[2] = 35; verts_lvl2[3] = 36;
+ verts_lvl2[4] = 38; verts_lvl2[5] = 39;
+ verts_lvl2[6] = 32; verts_lvl2[7] = 33;
+ }
+
+ int verts_lvl1[4];
+ switch(nPts) {
+ case(4):
+ verts_lvl1[0] = 39;
+ verts_lvl1[1] = 45;
+ verts_lvl1[2] = 51;
+ verts_lvl1[3] = 57;
+ break;
+ case(3):
+ verts_lvl1[0] = 31;
+ verts_lvl1[1] = 37;
+ verts_lvl1[2] = 40;
+ verts_lvl1[3] = 34;
+ break;
+ case(2):
+ verts_lvl1[0] = 21;
+ verts_lvl1[1] = 23;
+ verts_lvl1[2] = 24;
+ verts_lvl1[3] = 22;
+ break;
+ }
+
+ for (int q = 0; q < nPts - 1; q++) {
+ std::vector<MVertex*> vq, veq;
+ vq.push_back(ve[2*nPts + q]);
+ vq.push_back(ve[4*nPts + q]);
+ vq.push_back(ve[6*nPts + q]);
+ vq.push_back(ve[7*nPts + q]);
+
+ int triverts = nPts*(nPts-1)/2;
+
+ if (nPts-q == 4)
+ for (int f = 0; f < 12; f++)
+ veq.push_back(ve[verts_lvl3[f]-5]);
+ else if (nPts-q == 3)
+ for (int f = 0; f < 8; f++)
+ veq.push_back(ve[verts_lvl2[f]-5]);
+ else if (nPts-q == 2)
+ for (int f = 0; f < 4; f++)
+ veq.push_back(ve[verts_lvl1[f]-5]);
+
+ if (nPts-q == 2) {
+ MQuadrangle8 incpl2(vq[0], vq[1], vq[2], vq[3],
+ veq[0], veq[1], veq[2], veq[3]);
+ SPoint3 pointz;
+ incpl2.pnt(0,0,0,pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+
+ gr->mesh_vertices.push_back(v);
+ std::vector<MVertex*>::iterator cursor = vr.begin();
+ cursor += nPts == 2 ? 0 : 4;
+ vr.insert(cursor, v);
+ }
+ else if (nPts-q == 3) {
+
+ MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 3);
+ int offsets[4] = {nPts == 4 ? 7 : 0,
+ nPts == 4 ? 9 : 1,
+ nPts == 4 ? 11 : 2,
+ nPts == 4 ? 12 : 3};
+ double quad_v [4][2] = {{-1.0/3.0, -1.0/3.0},
+ { 1.0/3.0, -1.0/3.0},
+ { 1.0/3.0, 1.0/3.0},
+ {-1.0/3.0, 1.0/3.0}};
+ SPoint3 pointz;
+ for (int k = 0; k<4; k++) {
+ incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ std::vector<MVertex*>::iterator cursor = vr.begin();
+ cursor += offsets[k];
+ vr.insert(cursor, v);
+ }
+ }
+ else if (nPts-q == 4) {
+ MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 4);
+ int offsets[9] = {0, 1, 2, 3, 5, 8, 10, 6, 13};
+ double quad_v [9][2] = {
+ { -0.5, -0.5},
+ { 0.5, -0.5},
+ { 0.5, 0.5},
+ { -0.5, 0.5},
+ { 0.0, -0.5},
+ { 0.5, 0.0},
+ { 0.0, 0.5},
+ { -0.5, 0.0},
+ { 0.0, 0.0}
+ };
+ SPoint3 pointz;
+ for (int k = 0; k<9; k++) {
+ incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ std::vector<MVertex*>::iterator cursor = vr.begin();
+ cursor += offsets[k];
+ vr.insert(cursor, v);
+ }
+
+
+ }
+
+ }
ve.insert(ve.end(), vr.begin(), vr.end());
- printf("Hmmm, let's see : %d\n", nPts);
MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve, nPts + 1);
+ p->getVertex(3), p->getVertex(4), ve, nPts + 1);
pyramids2.push_back(n);
- */
+ SPoint3 test_pnt;
+ n->pnt(-1,-1,0, test_pnt);
}
delete p;
}
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index 2e966677a2b9f32b08b7cecb78ad6ee20cbad761..b2d2cb06af000d66646449d7e74eaf6bc9f78a66 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -469,7 +469,7 @@ void highOrderTools::computeStraightSidedPositions ()
for (unsigned int i=0;i<(*it)->triangles.size();i++){
MTriangle *t = (*it)->triangles[i];
MFace face = t->getFace(0);
- const polynomialBasis* fs = t->getFunctionSpace();
+ const nodalBasis* fs = t->getFunctionSpace();
for (int j=0;j<t->getNumVertices();j++){
if (t->getVertex(j)->onWhat() == *it){
const double t1 = fs->points(j, 0);
@@ -484,7 +484,7 @@ void highOrderTools::computeStraightSidedPositions ()
// printf("coucou quad %d\n",i);
MQuadrangle *q = (*it)->quadrangles[i];
MFace face = q->getFace(0);
- const polynomialBasis* fs = q->getFunctionSpace();
+ const nodalBasis* fs = q->getFunctionSpace();
for (int j=0;j<q->getNumVertices();j++){
if (q->getVertex(j)->onWhat() == *it){
const double t1 = fs->points(j, 0);
@@ -509,7 +509,7 @@ void highOrderTools::computeStraightSidedPositions ()
(*it)->tetrahedra[i]->getVertex(1),
(*it)->tetrahedra[i]->getVertex(2),
(*it)->tetrahedra[i]->getVertex(3));
- const polynomialBasis* fs = t->getFunctionSpace();
+ const nodalBasis* fs = t->getFunctionSpace();
for (int j=0;j<t->getNumVertices();j++){
if (t->getVertex(j)->onWhat() == *it){
double t1 = fs->points(j, 0);
@@ -532,7 +532,7 @@ void highOrderTools::computeStraightSidedPositions ()
(*it)->hexahedra[i]->getVertex(5),
(*it)->hexahedra[i]->getVertex(6),
(*it)->hexahedra[i]->getVertex(7));
- const polynomialBasis* fs = h->getFunctionSpace();
+ const nodalBasis* fs = h->getFunctionSpace();
for (int j=0;j<h->getNumVertices();j++){
if (h->getVertex(j)->onWhat() == *it){
double t1 = fs->points(j, 0);
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..56f38633cab15172ff58a50d2dc7f5730238b7dd
--- /dev/null
+++ b/Numeric/BasisFactory.cpp
@@ -0,0 +1,216 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-B-> Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include "polynomialBasis.h"
+#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
+#include "BasisFactory.h"
+
+std::map<int, nodalBasis*> BasisFactory::fs;
+
+const nodalBasis* BasisFactory::create(int elementType) {
+
+ // If the Basis has already been built, return it.
+ std::map<int, nodalBasis*>::const_iterator it = fs.find(elementType);
+ if (it != fs.end()) {
+ return it->second;
+ }
+ // Get the parent type to see which kind of basis
+ // we want to create
+ int parentType = MElement::ParentTypeFromTag(elementType);
+ nodalBasis* B = 0;
+
+ switch(parentType) {
+ case(TYPE_PNT):
+ case(TYPE_LIN):
+ case(TYPE_TRI):
+ case(TYPE_QUA):
+ case(TYPE_PRI):
+ case(TYPE_TET):
+ case(TYPE_HEX):
+ B = new polynomialBasis();
+ break;
+ case(TYPE_PYR):
+ B = new pyramidalBasis();
+ break;
+ default:
+ Msg::Error("Unknown type of element.");
+ return 0;
+ }
+ B->parentType = parentType;
+
+ // There is currently only 1 type of basis for any
+ // element type, hence we use the element type as the
+ // basis type.
+ B->type = elementType;
+ switch(elementType) {
+ case MSH_PNT : B->order = 0; B->serendip = false; break;
+ case MSH_LIN_1 : B->order = 0; B->serendip = false; break;
+ case MSH_LIN_2 : B->order = 1; B->serendip = false; break;
+ case MSH_LIN_3 : B->order = 2; B->serendip = false; break;
+ case MSH_LIN_4 : B->order = 3; B->serendip = false; break;
+ case MSH_LIN_5 : B->order = 4; B->serendip = false; break;
+ case MSH_LIN_6 : B->order = 5; B->serendip = false; break;
+ case MSH_LIN_7 : B->order = 6; B->serendip = false; break;
+ case MSH_LIN_8 : B->order = 7; B->serendip = false; break;
+ case MSH_LIN_9 : B->order = 8; B->serendip = false; break;
+ case MSH_LIN_10 : B->order = 9; B->serendip = false; break;
+ case MSH_LIN_11 : B->order = 10;B->serendip = false; break;
+ case MSH_TRI_1 : B->order = 0; B->serendip = false; break;
+ case MSH_TRI_3 : B->order = 1; B->serendip = false; break;
+ case MSH_TRI_6 : B->order = 2; B->serendip = false; break;
+ case MSH_TRI_10 : B->order = 3; B->serendip = false; break;
+ case MSH_TRI_15 : B->order = 4; B->serendip = false; break;
+ case MSH_TRI_21 : B->order = 5; B->serendip = false; break;
+ case MSH_TRI_28 : B->order = 6; B->serendip = false; break;
+ case MSH_TRI_36 : B->order = 7; B->serendip = false; break;
+ case MSH_TRI_45 : B->order = 8; B->serendip = false; break;
+ case MSH_TRI_55 : B->order = 9; B->serendip = false; break;
+ case MSH_TRI_66 : B->order = 10;B->serendip = false; break;
+ case MSH_TRI_9 : B->order = 3; B->serendip = true; break;
+ case MSH_TRI_12 : B->order = 4; B->serendip = true; break;
+ case MSH_TRI_15I : B->order = 5; B->serendip = true; break;
+ case MSH_TRI_18 : B->order = 6; B->serendip = true; break;
+ case MSH_TRI_21I : B->order = 7; B->serendip = true; break;
+ case MSH_TRI_24 : B->order = 8; B->serendip = true; break;
+ case MSH_TRI_27 : B->order = 9; B->serendip = true; break;
+ case MSH_TRI_30 : B->order = 10;B->serendip = true; break;
+ case MSH_TET_1 : B->order = 0; B->serendip = false; break;
+ case MSH_TET_4 : B->order = 1; B->serendip = false; break;
+ case MSH_TET_10 : B->order = 2; B->serendip = false; break;
+ case MSH_TET_20 : B->order = 3; B->serendip = false; break;
+ case MSH_TET_35 : B->order = 4; B->serendip = false; break;
+ case MSH_TET_56 : B->order = 5; B->serendip = false; break;
+ case MSH_TET_84 : B->order = 6; B->serendip = false; break;
+ case MSH_TET_120 : B->order = 7; B->serendip = false; break;
+ case MSH_TET_165 : B->order = 8; B->serendip = false; break;
+ case MSH_TET_220 : B->order = 9; B->serendip = false; break;
+ case MSH_TET_286 : B->order = 10;B->serendip = false; break;
+ case MSH_TET_34 : B->order = 4; B->serendip = true; break;
+ case MSH_TET_52 : B->order = 5; B->serendip = true; break;
+ case MSH_TET_74 : B->order = 6; B->serendip = true; break;
+ case MSH_TET_100 : B->order = 7; B->serendip = true; break;
+ case MSH_TET_130 : B->order = 8; B->serendip = true; break;
+ case MSH_TET_164 : B->order = 9; B->serendip = true; break;
+ case MSH_TET_202 : B->order = 10;B->serendip = true; break;
+ case MSH_QUA_1 : B->order = 0; B->serendip = false; break;
+ case MSH_QUA_4 : B->order = 1; B->serendip = false; break;
+ case MSH_QUA_9 : B->order = 2; B->serendip = false; break;
+ case MSH_QUA_16 : B->order = 3; B->serendip = false; break;
+ case MSH_QUA_25 : B->order = 4; B->serendip = false; break;
+ case MSH_QUA_36 : B->order = 5; B->serendip = false; break;
+ case MSH_QUA_49 : B->order = 6; B->serendip = false; break;
+ case MSH_QUA_64 : B->order = 7; B->serendip = false; break;
+ case MSH_QUA_81 : B->order = 8; B->serendip = false; break;
+ case MSH_QUA_100 : B->order = 9; B->serendip = false; break;
+ case MSH_QUA_121 : B->order = 10;B->serendip = false; break;
+ case MSH_QUA_8 : B->order = 2; B->serendip = true; break;
+ case MSH_QUA_12 : B->order = 3; B->serendip = true; break;
+ case MSH_QUA_16I : B->order = 4; B->serendip = true; break;
+ case MSH_QUA_20 : B->order = 5; B->serendip = true; break;
+ case MSH_QUA_24 : B->order = 6; B->serendip = true; break;
+ case MSH_QUA_28 : B->order = 7; B->serendip = true; break;
+ case MSH_QUA_32 : B->order = 8; B->serendip = true; break;
+ case MSH_QUA_36I : B->order = 9; B->serendip = true; break;
+ case MSH_QUA_40 : B->order = 10;B->serendip = true; break;
+ case MSH_PRI_1 : B->order = 0; B->serendip = false; break;
+ case MSH_PRI_6 : B->order = 1; B->serendip = false; break;
+ case MSH_PRI_18 : B->order = 2; B->serendip = false; break;
+ case MSH_HEX_8 : B->order = 1; B->serendip = false; break;
+ case MSH_HEX_27 : B->order = 2; B->serendip = false; break;
+ case MSH_HEX_64 : B->order = 3; B->serendip = false; break;
+ case MSH_HEX_125 : B->order = 4; B->serendip = false; break;
+ case MSH_HEX_216 : B->order = 5; B->serendip = false; break;
+ case MSH_HEX_343 : B->order = 6; B->serendip = false; break;
+ case MSH_HEX_512 : B->order = 7; B->serendip = false; break;
+ case MSH_HEX_729 : B->order = 8; B->serendip = false; break;
+ case MSH_HEX_1000: B->order = 9; B->serendip = false; break;
+ case MSH_HEX_20 : B->order = 2; B->serendip = false; break;
+ case MSH_HEX_56 : B->order = 3; B->serendip = true; break;
+ case MSH_HEX_98 : B->order = 4; B->serendip = true; break;
+ case MSH_HEX_152 : B->order = 5; B->serendip = true; break;
+ case MSH_HEX_222 : B->order = 6; B->serendip = true; break;
+ case MSH_HEX_296 : B->order = 7; B->serendip = true; break;
+ case MSH_HEX_386 : B->order = 8; B->serendip = true; break;
+ case MSH_HEX_488 : B->order = 9; B->serendip = true; break;
+ case MSH_PYR_1 : B->order = 0; B->serendip = false; break;
+ case MSH_PYR_5 : B->order = 1; B->serendip = false; break;
+ case MSH_PYR_14 : B->order = 2; B->serendip = false; break;
+ case MSH_PYR_30 : B->order = 3; B->serendip = false; break;
+ case MSH_PYR_55 : B->order = 4; B->serendip = false; break;
+ case MSH_PYR_91 : B->order = 5; B->serendip = false; break;
+ case MSH_PYR_140 : B->order = 6; B->serendip = false; break;
+ case MSH_PYR_204 : B->order = 7; B->serendip = false; break;
+ case MSH_PYR_285 : B->order = 8; B->serendip = false; break;
+ case MSH_PYR_385 : B->order = 9; B->serendip = false; break;
+ case MSH_PYR_29 : B->order = 3; B->serendip = true; break;
+ case MSH_PYR_50 : B->order = 4; B->serendip = true; break;
+ case MSH_PYR_77 : B->order = 5; B->serendip = true; break;
+ case MSH_PYR_110 : B->order = 6; B->serendip = true; break;
+ case MSH_PYR_149 : B->order = 7; B->serendip = true; break;
+ case MSH_PYR_194 : B->order = 8; B->serendip = true; break;
+ case MSH_PYR_245 : B->order = 9; B->serendip = true; break;
+
+ default :
+ Msg::Error("Unknown function space %d: reverting to TET_4", elementType);
+ B->parentType = TYPE_TET; B->order = 1; B->serendip = false; break;
+ }
+
+ // Here is where we do create the basis.
+ switch (B->parentType) {
+ case TYPE_PNT :
+ B->numFaces = 1;
+ B->dimension = 0;
+ B->points = gmshGeneratePointsLine(0);
+ break;
+ case TYPE_LIN :
+ B->numFaces = 2;
+ B->dimension = 1;
+ B->points = gmshGeneratePointsLine(B->order);
+ break;
+ case TYPE_TRI :
+ B->numFaces = 3;
+ B->dimension = 2;
+ B->points = gmshGeneratePointsTriangle(B->order, B->serendip);
+ break;
+ case TYPE_QUA :
+ B->numFaces = 4;
+ B->dimension = 2;
+ B->points = gmshGeneratePointsQuadrangle(B->order, B->serendip);
+ break;
+ case TYPE_TET :
+ B->numFaces = 4;
+ B->dimension = 3;
+ B->points = gmshGeneratePointsTetrahedron(B->order, B->serendip);
+ B->points.print();
+ break;
+ case TYPE_PRI :
+ B->numFaces = 5;
+ B->dimension = 3;
+ B->points = gmshGeneratePointsPrism(B->order, B->serendip);
+ break;
+ case TYPE_HEX :
+ B->numFaces = 6;
+ B->dimension = 3;
+ B->points = gmshGeneratePointsHexahedron(B->order, B->serendip);
+ break;
+ case TYPE_PYR :
+ B->numFaces = 5;
+ B->dimension = 3;
+ B->points = gmshGeneratePointsPyramid(B->order, B->serendip);
+ break;
+ default:
+ Msg::Error("Unknown element type!\n");
+ }
+
+ B->initialize();
+
+ fs.insert(std::make_pair(elementType, B));
+
+ return fs[elementType];
+
+}
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..691a8340d29088da2d715d2a632e642041ce0afa
--- /dev/null
+++ b/Numeric/BasisFactory.h
@@ -0,0 +1,24 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef BASISFACTORY_H
+#define BASISFACTORY_H
+
+#include "MElement.h"
+#include "MPyramid.h"
+#include "nodalBasis.h"
+
+class BasisFactory
+{
+ private:
+ static std::map<int, nodalBasis*> fs;
+
+ public:
+ // Caution: the returned pointer should be
+ // checked against 0.
+ static const nodalBasis* create(int elementType);
+};
+
+#endif
\ No newline at end of file
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index bc8ce13da41f92770cebeafd83f486dae3c5fcfc..e361960798f64e78940a03ef197e3f4b28ce98b0 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -5,6 +5,11 @@
#include "BergotBasis.h"
+/*BergotBasis::BergotBasis() {
+
+
+}*/
+
BergotBasis::BergotBasis(int p): order(p)
{
// allocate function information and fill
@@ -13,7 +18,8 @@ BergotBasis::BergotBasis(int p): order(p)
jOrder = new int[size()];
kOrder = new int[size()];
- legendre = new LegendrePolynomials(order);
+ LegendrePolynomials lp(p);
+ legendre = lp;
int index = 0;
for (int i=0;i<=order;i++) {
@@ -26,22 +32,21 @@ BergotBasis::BergotBasis(int p): order(p)
kOrder[index] = k;
if (jacobi.find(mIJ) == jacobi.end()) {
- jacobi[mIJ] = new JacobiPolynomials(2*mIJ+2,0,order);
+ JacobiPolynomials jp(2*mIJ+2,0,order);
+ jacobi[mIJ] = jp;
}
}
}
}
}
+
BergotBasis::~BergotBasis()
{
delete [] iOrder;
delete [] jOrder;
delete [] kOrder;
- delete legendre;
- std::map<int,JacobiPolynomials*>::iterator jIter = jacobi.begin();
- for (;jIter!=jacobi.end();++jIter) delete jIter->second;
}
int BergotBasis::size() const
@@ -54,21 +59,22 @@ void BergotBasis::f(double u, double v, double w, double* val) const
double uhat = (w == 1.) ? 1 : u/(1.-w);
std::vector<double> uFcts(order+1);
- legendre->f(uhat,&(uFcts[0]));
+ //double uFcts[order+1];
+ legendre.f(uhat,&(uFcts[0]));
double vhat = (w == 1.) ? 1 : v/(1.-w);
std::vector<double> vFcts(order+1);
- legendre->f(vhat,&(vFcts[0]));
+ legendre.f(vhat,&(vFcts[0]));
double what = 2.*w - 1.;
std::map<int,double* > wFcts;
- std::map<int,JacobiPolynomials*>::const_iterator jIter=jacobi.begin();
+ std::map<int,JacobiPolynomials>::const_iterator jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
wFcts[a] = new double[order + 1];
double* wf = wFcts[a];
- jIter->second->f(what,wf);
+ jIter->second.f(what,wf);
}
// recombine to find shape function values
@@ -93,33 +99,32 @@ void BergotBasis::f(double u, double v, double w, double* val) const
}
}
-void BergotBasis::df(double u, double v, double w, double grads[][3]) const
+inline void BergotBasis::df(double u, double v, double w, double grads[][3]) const
{
std::vector<double> uFcts(order+1);
- legendre->f(u,&(uFcts[0]));
+ legendre.f(u,&(uFcts[0]));
std::vector<double> uGrads(order+1);
- legendre->df(u,&(uGrads[0]));
+ legendre.df(u,&(uGrads[0]));
std::vector<double> vFcts(order+1);
- legendre->f(v,&(vFcts[0]));
-
- std::vector<double> vGrads(order+1);
- legendre->df(v,&(vGrads[0]));
+ legendre.f(v,&(vFcts[0]));
+ std::vector<double> vGrads(order+1);
+ legendre.df(v,&(vGrads[0]));
std::map<int,double* > wFcts;
std::map<int,double* > wGrads;
- std::map<int,JacobiPolynomials*>::const_iterator jIter=jacobi.begin();
+ std::map<int,JacobiPolynomials>::const_iterator jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
wFcts[a] = new double[order+1];
double* wf = wFcts[a];
- jIter->second->f(w,wf);
+ jIter->second.f(w,wf);
wGrads[a] = new double[order+1];
double* wg = wGrads[a];
- jIter->second->df(w,wg);
+ jIter->second.df(w,wg);
}
// now recombine to find the shape function gradients
diff --git a/Numeric/BergotBasis.h b/Numeric/BergotBasis.h
index cecf7b1bbe7e6f0c7a2bee7eaacf82ac8bfbfb4e..0a68e68ce46ae34cb3db150fd8799ae812c54c5f 100644
--- a/Numeric/BergotBasis.h
+++ b/Numeric/BergotBasis.h
@@ -6,18 +6,22 @@
#ifndef BERGOTBASIS_H
#define BERGOTBASIS_H
+#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "jacobiPolynomials.h"
#include "legendrePolynomials.h"
-class BergotBasis : public polynomialBasis {
+class BergotBasis : public nodalBasis {
public:
+ BergotBasis() {};
BergotBasis(int p);
~BergotBasis();
void f(double u, double v, double w, double *val) const;
- void df(double u, double v, double w, double grads[][3]) const;
+ inline void df(double u, double v, double w, double grads[][3]) const;
+
+ void initialize() {};
private:
int order; //!< maximal order of surrounding functional spaces (on triangle / quad)
@@ -27,10 +31,10 @@ class BergotBasis : public polynomialBasis {
int *kOrder; //!< order of \f$\hat \zeta \f$ polynomial
//! list of Legendre polynomials up to order p
- LegendrePolynomials* legendre;
+ LegendrePolynomials legendre;
//! list of Jacobi polynomials up to order p in function of index i (\f$ \alpha = 2*i + 2\f$)
- std::map<int,JacobiPolynomials*> jacobi;
+ std::map<int,JacobiPolynomials> jacobi;
int size() const;
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 62bab4da4f03f835e9f700770b972d44bad3d9f4..aa2bb2bb932f5a5e7a23a0dd2df366f98ef1419a 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -6,9 +6,11 @@
set(SRC
Numeric.cpp
fullMatrix.cpp
+ BasisFactory.cpp
polynomialBasis.cpp
JacobianBasis.cpp
BergotBasis.cpp
+ pyramidalBasis.cpp
jacobiPolynomials.cpp
legendrePolynomials.cpp
pointsGenerators.cpp
diff --git a/Numeric/jacobiPolynomials.h b/Numeric/jacobiPolynomials.h
index 7b88a950c1be2b35727ca225166e8ef7f3c5c1f5..0b9919cc5f5afd01b2ea1cfda53c93a8d1fe5bb2 100644
--- a/Numeric/jacobiPolynomials.h
+++ b/Numeric/jacobiPolynomials.h
@@ -11,6 +11,7 @@
class JacobiPolynomials {
public:
+ JacobiPolynomials() {};
JacobiPolynomials(double a, double b, int o);
~JacobiPolynomials();
diff --git a/Numeric/legendrePolynomials.cpp b/Numeric/legendrePolynomials.cpp
index 3968e53c4eb48dee71f19fbfac5235cf273030c6..b0cc7281c14b5b3ed8bbab47c1925d395bfee216 100644
--- a/Numeric/legendrePolynomials.cpp
+++ b/Numeric/legendrePolynomials.cpp
@@ -5,8 +5,10 @@
#include "legendrePolynomials.h"
-LegendrePolynomials::LegendrePolynomials(int o): n(o) {}
-LegendrePolynomials::~LegendrePolynomials() {;}
+LegendrePolynomials::LegendrePolynomials(int o): n(o) {
+
+}
+LegendrePolynomials::~LegendrePolynomials() {}
void LegendrePolynomials::f(double u, double *val) const
{
diff --git a/Numeric/legendrePolynomials.h b/Numeric/legendrePolynomials.h
index ce7dd998b99f72adec4af6ca31e93aae01f6caa6..a7d9a497c1fe7dad03e4e106138b06805c933e03 100644
--- a/Numeric/legendrePolynomials.h
+++ b/Numeric/legendrePolynomials.h
@@ -11,6 +11,7 @@
class LegendrePolynomials {
public:
+ LegendrePolynomials() {};
LegendrePolynomials(int o);
~LegendrePolynomials();
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..94791d6bc7ba7be9fdb515446d9f7997ba032395
--- /dev/null
+++ b/Numeric/nodalBasis.h
@@ -0,0 +1,32 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef BASIS_H
+#define BASIS_H
+
+#include "fullMatrix.h"
+
+class nodalBasis {
+
+ public:
+ int type, parentType, order, dimension, numFaces;
+ bool serendip;
+ fullMatrix<double> points;
+
+ virtual void initialize() {};
+
+ // Basis functions evaluation
+ inline virtual void f(double u, double v, double w, double *sf) const {};
+ inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const;
+
+ // Basis functions gradients evaluation
+ inline virtual void df(double u, double v, double w, double grads[][3]) const {};
+ inline void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
+
+ inline void ddf(double u, double v, double w, double grads[][3][3]) const {};
+ inline void dddf(double u, double v, double w, double grads[][3][3][3]) const {};
+
+};
+#endif
\ No newline at end of file
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index c340ca51d66bb7ec7a7e0801d78b2f2c810532cf..aa14897fe739f2f96beef6a1b861183f8753ca36 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -430,7 +430,6 @@ fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
}
points.scale(overOrder);
- points.print();
return points;
}
@@ -778,9 +777,5 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
}
}
}
- for (int i = 0; i < points.size1(); i++) {
- printf("Point(%d) %g, %g, %g\n", i, points(i,0),points(i,1),points(i,2));
- }
- printf("number of points : %d\n\n", nbPoints);
return points;
}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 1922abd238990e89adeeaf916d690ec09fe4b04e..ea67a22634c9f32d7c237a81fdb582dc1f785000 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -42,6 +42,8 @@ static void printClosure(polynomialBasis::clCont &fullClosure,
}
*/
+
+
int polynomialBasis::getTag(int parentTag, int order, bool serendip)
{
switch (parentTag) {
@@ -131,7 +133,7 @@ int polynomialBasis::getTag(int parentTag, int order, bool serendip)
}
}
-static fullMatrix<double> generate1DMonomials(int order)
+fullMatrix<double> generate1DMonomials(int order)
{
fullMatrix<double> monomials(order + 1, 1);
for (int i = 0; i < order + 1; i++) monomials(i, 0) = i;
@@ -964,7 +966,6 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
(const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{
if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Info("Here");
Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
monomial.size1(),point.size1(),
monomial.size2(),point.size2() );
@@ -1538,6 +1539,82 @@ static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
std::map<int, polynomialBasis> polynomialBases::fs;
+
+void polynomialBasis::initialize()
+{
+ switch(parentType) {
+ case TYPE_PNT:
+ monomials = generate1DMonomials(0);
+ break;
+ case TYPE_LIN:
+ monomials = generate1DMonomials(order);
+ generate1dVertexClosure(closures, order);
+ generate1dVertexClosureFull(fullClosures, closureRef, order);
+ break;
+ case TYPE_TRI:
+ monomials = serendip ? generatePascalSerendipityTriangle(order) :
+ generatePascalTriangle(order);
+ if (order == 0) {
+ generateClosureOrder0(closures, 6);
+ generateClosureOrder0(fullClosures, 6);
+ closureRef.resize(6,0);
+ }
+ else {
+ generate2dEdgeClosure(closures, order);
+ generate2dEdgeClosureFull(fullClosures, closureRef, order, 3, serendip);
+ }
+ break;
+ case TYPE_QUA:
+ monomials = serendip ? generatePascalQuadSerendip(order) :
+ generatePascalQuad(order);
+ if (order == 0) {
+ generateClosureOrder0(closures, 8);
+ generateClosureOrder0(fullClosures, 8);
+ closureRef.resize(8, 0);
+ }
+ else {
+ generate2dEdgeClosure(closures, order, 4);
+ generate2dEdgeClosureFull(fullClosures, closureRef, order, 4, serendip);
+ }
+ break;
+ case TYPE_TET:
+ monomials = serendip ? generatePascalSerendipityTetrahedron(order) :
+ generatePascalTetrahedron(order);
+ if (order == 0) {
+ generateClosureOrder0(closures,24);
+ generateClosureOrder0(fullClosures, 24);
+ closureRef.resize(24, 0);
+ }
+ else {
+ generateFaceClosureTet(closures, order);
+ generateFaceClosureTetFull(fullClosures, closureRef, order, serendip);
+ }
+ break;
+ case TYPE_PRI:
+ monomials = generatePascalPrism(order);
+ if (order == 0) {
+ generateClosureOrder0(closures,48);
+ generateClosureOrder0(fullClosures,48);
+ closureRef.resize(48, 0);
+ }
+ else {
+ generateFaceClosurePrism(closures, order);
+ generateFaceClosurePrismFull(fullClosures, closureRef, order);
+ }
+ break;
+ case TYPE_HEX:
+ monomials = generatePascalHex(order, serendip);
+ generateFaceClosureHex(closures, order, serendip, points);
+ generateFaceClosureHexFull(fullClosures, closureRef, order, serendip, points);
+ break;
+ }
+
+ this->points.print();
+ coefficients = generateLagrangeMonomialCoefficients(monomials, points);
+
+}
+
+
const polynomialBasis *polynomialBases::find(int tag)
{
std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 7b1988c4e8d8ee16e4fb25fb3a4509f2417520e0..72e5a01c2d34b0f2edccd27aa4455013a4ab90ed 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -10,6 +10,7 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
+#include "nodalBasis.h"
#include <iostream>
#define SQU(a) ((a)*(a))
@@ -65,15 +66,20 @@ inline double pow_int(const double &a, const int &n)
}
}
-class polynomialBasis
+fullMatrix<double> generate1DMonomials(int order);
+
+
+
+
+class polynomialBasis : public nodalBasis
{
// integrationOrder, closureId => df/dXi
mutable std::map<int,std::vector<fullMatrix<double> > > _dfAtFace;
public:
// for now the only implemented polynomial basis are nodal poly
// basis, we use the type of the corresponding gmsh element as type
- int type, parentType, order, dimension;
- bool serendip;
+ //int type, parentType, order, dimension;
+ //bool serendip;
class closure : public std::vector<int> {
public:
int type;
@@ -89,10 +95,13 @@ class polynomialBasis
// quad face)
clCont closures, fullClosures;
std::vector<int> closureRef;
- fullMatrix<double> points;
+ //fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
int numFaces;
+
+ void initialize();
+
inline int getNumShapeFunctions() const {return coefficients.size1();}
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf84f8832f88de579e8366bc4eb2cd4e7a840fb4
--- /dev/null
+++ b/Numeric/pyramidalBasis.cpp
@@ -0,0 +1,60 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "pyramidalBasis.h"
+
+pyramidalBasis::pyramidalBasis()
+{
+
+}
+
+void pyramidalBasis::initialize()
+{
+
+ BergotBasis bb(order);
+ bergot = bb;
+
+ int n = order+1;
+ int num_points = n*(n+1)*(2*n+1)/6;
+ if (serendip and order > 2) {
+ num_points -= (order-2)*((order-2)+1)*(2*(order-2)+1)/6;
+ }
+
+ VDMinv = new fullMatrix<double>(num_points, num_points);
+
+ // Invert the Vandermonde matrix
+ fullMatrix<double> VDM = fullMatrix<double>(num_points, num_points);
+ for (int j = 0; j < num_points; j++) {
+ double f[num_points];
+ bergot.f(points(j,0), points(j,1), points(j, 2), f);
+ for (int i = 0; i < num_points; i++) {
+ VDM(i,j) = f[i];
+ }
+ }
+ VDM.invert(*VDMinv);
+}
+
+inline void pyramidalBasis::f(double u, double v, double w, double *val) const
+{
+ int n = order+1;
+ int num_points = n*(n+1)*(2*n+1)/6;
+
+ double f[n*n*n];
+ bergot.f(u,v,w, f);
+
+ for (int i = 0; i < num_points; i++) {
+ val[i] = 0.0;
+ for (int j = 0; j < num_points; j++) {
+ val[i] += (*VDMinv)(i,j)*f[j];
+ }
+ }
+}
+
+inline void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
+{
+
+}
+
+
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..1055d039e8980cf18bdbe424a3468ae763dc031e
--- /dev/null
+++ b/Numeric/pyramidalBasis.h
@@ -0,0 +1,34 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef PYRAMIDALBASIS_H
+#define PYRAMIDALBASIS_H
+
+#include "fullMatrix.h"
+
+#include "nodalBasis.h"
+#include "BergotBasis.h"
+
+class pyramidalBasis: public nodalBasis
+{
+ private:
+ // Inverse of the Vandermonde matrix
+ fullMatrix<double>* VDMinv;
+
+ // Orthogonal basis for the pyramid
+ BergotBasis bergot;
+
+ public:
+ pyramidalBasis();
+
+ inline void f(double u, double v, double w, double *val) const;
+ inline void df(double u, double v, double w, double grads[][3]) const;
+
+ void initialize();
+
+};
+
+
+#endif
\ No newline at end of file
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 3935fd52356d70585ca9a220704471966cf3b87d..b54f7b204b307d9bc41d02a50f77c0d270dfc1ee 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -143,7 +143,7 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
// the mesh is curved
MElement *e = _getOneElementOfGivenType(model, it->first);
if(e && e->getPolynomialOrder() > 1 && e->getFunctionSpace()){
- const polynomialBasis *fs = e->getFunctionSpace();
+ const polynomialBasis *fs = (polynomialBasis*) e->getFunctionSpace();
setInterpolationMatrices(it->first, *(it->second[0]), *(it->second[1]),
fs->coefficients, fs->monomials);
}
@@ -173,7 +173,7 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
if(!haveInterpolationMatrices(types[i])){
MElement *e = _getOneElementOfGivenType(model, types[i]);
if(e){
- const polynomialBasis *fs = e->getFunctionSpace();
+ const polynomialBasis *fs = (polynomialBasis*) e->getFunctionSpace();
if(fs){
if(e->getPolynomialOrder() > 1){
if(_type == ElementData){
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 20f7035bc0555ea2e53bdb976e5a9a07c286bb6c..e940e8b3951b390ea1331087a49b14e757a3cdad 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -70,7 +70,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
for(std::set<MElement*>::const_iterator it = els.begin(); it != els.end(); ++it, ++iEl) {
MElement *el = *it;
_el[iEl] = el;
- const polynomialBasis *lagrange = el->getFunctionSpace();
+ const polynomialBasis *lagrange = (polynomialBasis*) el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_lag2Bez.find(lagrange->type) == _lag2Bez.end()) {
_gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);
@@ -236,8 +236,8 @@ void Mesh::distSqToStraight(std::vector<double> &dSq)
std::vector<SPoint3> sxyz(nVert());
for (int iEl = 0; iEl < nEl(); iEl++) {
MElement *el = _el[iEl];
- const polynomialBasis *lagrange = el->getFunctionSpace();
- const polynomialBasis *lagrange1 = el->getFunctionSpace(1);
+ const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
+ const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
int nV = lagrange->points.size1();
int nV1 = lagrange1->points.size1();
for (int i = 0; i < nV1; ++i) {
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 4eb27f3f4dcdbef2449f3f608f232d7f3c5fb23c..610a41b858b5f6ffb0095e267cd824bde472627a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -51,9 +51,8 @@ void OptHomMessage (const char *s, ...) {
double distMaxStraight (MElement *el)
{
-
- const polynomialBasis *lagrange = el->getFunctionSpace();
- const polynomialBasis *lagrange1 = el->getFunctionSpace(1);
+ const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
+ const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
int nV = lagrange->points.size1();
int nV1 = lagrange1->points.size1();
SPoint3 sxyz[256];