diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 3c96b885ede5df82be23f3e23c42cce89b926b2b..99f8a6d7ca515d617afa68c1b6a82f06deb2e9c7 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,4 +1,4 @@
-#include "polynomialBasis.h"
+#include "BasisFactory.h"
#include "Exception.h"
#include "TriLagrangeBasis.h"
@@ -57,7 +57,7 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
break;
}
- l = polynomialBases::find(tag);
+ l = (polynomialBasis*)BasisFactory::create(tag);
point = new fullMatrix<double>(triPoint(order));
// Set Basis Type //
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 1d76cb7a6bad51c1aff08fed0b5c02ceb8d4390f..bdc5c5132e16a4f74b2dde6ae891676cc2072f64 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,7 +112,7 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
{
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
- if (getType() == TYPE_PYR) return;
+// if (getType() == TYPE_PYR) return;
extern double mesh_functional_distorsion(MElement*,double,double);
if (getPolynomialOrder() == 1) return;
const bezierBasis *jac = getJacobianFuncSpace()->bezier;
@@ -1416,6 +1416,138 @@ int MElement::ParentTypeFromTag(int tag)
}
}
+// Gives the order corresponding to any element type.
+int MElement::OrderFromTag(int tag)
+{
+
+ switch (tag) {
+ case MSH_PNT : return 0;
+ case MSH_LIN_1 : return 0;
+ case MSH_LIN_2 : return 1;
+ case MSH_LIN_3 : return 2;
+ case MSH_LIN_4 : return 3;
+ case MSH_LIN_5 : return 4;
+ case MSH_LIN_6 : return 5;
+ case MSH_LIN_7 : return 6;
+ case MSH_LIN_8 : return 7;
+ case MSH_LIN_9 : return 8;
+ case MSH_LIN_10 : return 9;
+ case MSH_LIN_11 : return 10;
+ case MSH_TRI_1 : return 0;
+ case MSH_TRI_3 : return 1;
+ case MSH_TRI_6 : return 2;
+ case MSH_TRI_10 : return 3;
+ case MSH_TRI_15 : return 4;
+ case MSH_TRI_21 : return 5;
+ case MSH_TRI_28 : return 6;
+ case MSH_TRI_36 : return 7;
+ case MSH_TRI_45 : return 8;
+ case MSH_TRI_55 : return 9;
+ case MSH_TRI_66 : return 10;
+ case MSH_TRI_9 : return 3;
+ case MSH_TRI_12 : return 4;
+ case MSH_TRI_15I : return 5;
+ case MSH_TRI_18 : return 6;
+ case MSH_TRI_21I : return 7;
+ case MSH_TRI_24 : return 8;
+ case MSH_TRI_27 : return 9;
+ case MSH_TRI_30 : return 10;
+ case MSH_TET_1 : return 0;
+ case MSH_TET_4 : return 1;
+ case MSH_TET_10 : return 2;
+ case MSH_TET_20 : return 3;
+ case MSH_TET_35 : return 4;
+ case MSH_TET_56 : return 5;
+ case MSH_TET_84 : return 6;
+ case MSH_TET_120 : return 7;
+ case MSH_TET_165 : return 8;
+ case MSH_TET_220 : return 9;
+ case MSH_TET_286 : return 10;
+ case MSH_TET_34 : return 4;
+ case MSH_TET_52 : return 5;
+ case MSH_TET_74 : return 6;
+ case MSH_TET_100 : return 7;
+ case MSH_TET_130 : return 8;
+ case MSH_TET_164 : return 9;
+ case MSH_TET_202 : return 10;
+ case MSH_QUA_1 : return 0;
+ case MSH_QUA_4 : return 1;
+ case MSH_QUA_9 : return 2;
+ case MSH_QUA_16 : return 3;
+ case MSH_QUA_25 : return 4;
+ case MSH_QUA_36 : return 5;
+ case MSH_QUA_49 : return 6;
+ case MSH_QUA_64 : return 7;
+ case MSH_QUA_81 : return 8;
+ case MSH_QUA_100 : return 9;
+ case MSH_QUA_121 : return 10;
+ case MSH_QUA_8 : return 2;
+ case MSH_QUA_12 : return 3;
+ case MSH_QUA_16I : return 4;
+ case MSH_QUA_20 : return 5;
+ case MSH_QUA_24 : return 6;
+ case MSH_QUA_28 : return 7;
+ case MSH_QUA_32 : return 8;
+ case MSH_QUA_36I : return 9;
+ case MSH_QUA_40 : return 10;
+ case MSH_PRI_1 : return 0;
+ case MSH_PRI_6 : return 1;
+ case MSH_PRI_18 : return 2;
+ case MSH_PRI_40 : return 3;
+ case MSH_PRI_75 : return 4;
+ case MSH_PRI_126 : return 5;
+ case MSH_PRI_196 : return 6;
+ case MSH_PRI_288 : return 7;
+ case MSH_PRI_405 : return 8;
+ case MSH_PRI_550 : return 9;
+ case MSH_PRI_15 : return 2;
+ case MSH_PRI_38 : return 3;
+ case MSH_PRI_66 : return 4;
+ case MSH_PRI_102 : return 5;
+ case MSH_PRI_146 : return 6;
+ case MSH_PRI_198 : return 7;
+ case MSH_PRI_258 : return 8;
+ case MSH_PRI_326 : return 9;
+ case MSH_HEX_8 : return 1;
+ case MSH_HEX_27 : return 2;
+ case MSH_HEX_64 : return 3;
+ case MSH_HEX_125 : return 4;
+ case MSH_HEX_216 : return 5;
+ case MSH_HEX_343 : return 6;
+ case MSH_HEX_512 : return 7;
+ case MSH_HEX_729 : return 8;
+ case MSH_HEX_1000: return 9;
+ case MSH_HEX_20 : return 2;
+ case MSH_HEX_56 : return 3;
+ case MSH_HEX_98 : return 4;
+ case MSH_HEX_152 : return 5;
+ case MSH_HEX_222 : return 6;
+ case MSH_HEX_296 : return 7;
+ case MSH_HEX_386 : return 8;
+ case MSH_HEX_488 : return 9;
+ case MSH_PYR_5 : return 1;
+ case MSH_PYR_14 : return 2;
+ case MSH_PYR_30 : return 3;
+ case MSH_PYR_55 : return 4;
+ case MSH_PYR_91 : return 5;
+ case MSH_PYR_140 : return 6;
+ case MSH_PYR_204 : return 7;
+ case MSH_PYR_285 : return 8;
+ case MSH_PYR_385 : return 9;
+ case MSH_PYR_13 : return 2;
+ case MSH_PYR_29 : return 3;
+ case MSH_PYR_50 : return 4;
+ case MSH_PYR_77 : return 5;
+ case MSH_PYR_110 : return 6;
+ case MSH_PYR_149 : return 7;
+ case MSH_PYR_194 : return 8;
+ case MSH_PYR_245 : return 9;
+ default :
+ Msg::Error("Unknown element type %d: reverting to order 1",tag);
+ return 1;
+ }
+
+}
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7abe28e7bbb911cc9e7b3ccb3c503fa9081f2830..9029b64772a471e223c2cb82f842d81f1b9cf73c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -355,6 +355,7 @@ class MElement
std::map<MElement*, MElement*> &newDomains);
static int ParentTypeFromTag(int tag);
+ static int OrderFromTag(int tag);
};
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 38dbab818806c213a7b42583929e9425158e89c9..b6dc0cedc3d76bfc0b92ae4607bd6b0dc47b4dcb 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -184,7 +184,7 @@ const nodalBasis* MHexahedron::getFunctionSpace(int o) const
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);
+ default: Msg::Error("Order %d hex function space not implemented", order); break;
}
}
else {
@@ -199,7 +199,44 @@ const nodalBasis* MHexahedron::getFunctionSpace(int o) const
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);
+ default: Msg::Error("Order %d hex function space not implemented", order); break;
+ }
+ }
+ return 0;
+}
+
+const JacobianBasis* MHexahedron::getJacobianFuncSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ int nv = getNumVolumeVertices();
+
+ if ((nv == 0) && (o == -1)) {
+ switch (order) {
+ case 1: return JacobianBasis::find(MSH_HEX_8);
+ case 2: return JacobianBasis::find(MSH_HEX_20);
+ case 3: return JacobianBasis::find(MSH_HEX_56);
+ case 4: return JacobianBasis::find(MSH_HEX_98);
+ case 5: return JacobianBasis::find(MSH_HEX_152);
+ case 6: return JacobianBasis::find(MSH_HEX_222);
+ case 7: return JacobianBasis::find(MSH_HEX_296);
+ case 8: return JacobianBasis::find(MSH_HEX_386);
+ case 9: return JacobianBasis::find(MSH_HEX_488);
+ default: Msg::Error("Order %d hex incomplete Jacobian function space not implemented", order); break;
+ }
+ }
+ else {
+ switch (order) {
+ case 1: return JacobianBasis::find(MSH_HEX_8);
+ case 2: return JacobianBasis::find(MSH_HEX_27);
+ case 3: return JacobianBasis::find(MSH_HEX_64);
+ case 4: return JacobianBasis::find(MSH_HEX_125);
+ case 5: return JacobianBasis::find(MSH_HEX_216);
+ case 6: return JacobianBasis::find(MSH_HEX_343);
+ case 7: return JacobianBasis::find(MSH_HEX_512);
+ case 8: return JacobianBasis::find(MSH_HEX_729);
+ case 9: return JacobianBasis::find(MSH_HEX_1000);
+ default: Msg::Error("Order %d hex Jacobian function space not implemented", order); break;
}
}
return 0;
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index d089dbae6141ed459a5c29144bdfd73f26d98860..d5984c1ef576f87da06504d1a9baa291aa999ea9 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -58,6 +58,7 @@ class MHexahedron : public MElement {
virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual const nodalBasis* getFunctionSpace(int o=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(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/MPyramid.cpp b/Geo/MPyramid.cpp
index 5eb5db0e489b93df5439cc58ffc5f1f2fa4e1954..f4e42ef6ca7a4fb60502618e5476fce80723cfd7 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -66,6 +66,20 @@ const nodalBasis* MPyramid::getFunctionSpace(int o) const
return 0;
}
+const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
+{
+ int order = (o == -1) ? getPolynomialOrder() : o;
+
+ switch (order) {
+ case 1: return JacobianBasis::find(MSH_PYR_5);
+ case 2: return JacobianBasis::find(MSH_PYR_14);
+ case 3: return JacobianBasis::find(MSH_PYR_30);
+ default: Msg::Error("Order %d pyramid function space not implemented", order); break;
+ }
+
+ return 0;
+}
+
MPyramidN::~MPyramidN() {}
double MPyramidN::distoShapeMeasure()
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 1277b1d75533e7867ba7605d3c01cc6964c196df..36bbca031974586a96d7e87ba858bf76c62b6101 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -68,6 +68,7 @@ class MPyramid : public MElement {
virtual int getNumVertices() const { return 5; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual const nodalBasis* getFunctionSpace(int o=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
virtual int getNumEdges(){ return 8; }
virtual MEdge getEdge(int num)
{
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index a7333333b68aa4d8041482933675a36f327b8346..150cb77043ff073019103f8632aa3ec310269b8e 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -21,7 +21,6 @@
#include "Numeric.h"
#include "Context.h"
#include "fullMatrix.h"
-#include "polynomialBasis.h"
#include "BasisFactory.h"
#define SQU(a) ((a)*(a))
@@ -598,14 +597,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
switch (nPts){
case 0 : break;
case 1 : break;
- case 2 : points = polynomialBases::find(MSH_TRI_10)->points; start = 9; break;
- case 3 : points = polynomialBases::find(MSH_TRI_15)->points; start = 12; break;
- case 4 : points = polynomialBases::find(MSH_TRI_21)->points; start = 15; break;
- case 5 : points = polynomialBases::find(MSH_TRI_28)->points; start = 18; break;
- case 6 : points = polynomialBases::find(MSH_TRI_36)->points; start = 21; break;
- case 7 : points = polynomialBases::find(MSH_TRI_45)->points; start = 24; break;
- case 8 : points = polynomialBases::find(MSH_TRI_55)->points; start = 27; break;
- case 9 : points = polynomialBases::find(MSH_TRI_66)->points; start = 30; break;
+ case 2 : points = BasisFactory::create(MSH_TRI_10)->points; start = 9; break;
+ case 3 : points = BasisFactory::create(MSH_TRI_15)->points; start = 12; break;
+ case 4 : points = BasisFactory::create(MSH_TRI_21)->points; start = 15; break;
+ case 5 : points = BasisFactory::create(MSH_TRI_28)->points; start = 18; break;
+ case 6 : points = BasisFactory::create(MSH_TRI_36)->points; start = 21; break;
+ case 7 : points = BasisFactory::create(MSH_TRI_45)->points; start = 24; break;
+ case 8 : points = BasisFactory::create(MSH_TRI_55)->points; start = 27; break;
+ case 9 : points = BasisFactory::create(MSH_TRI_66)->points; start = 30; break;
default :
Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
break;
@@ -614,14 +613,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
case 4:
switch (nPts){
case 0 : break;
- case 1 : points = polynomialBases::find(MSH_QUA_9)->points; break;
- case 2 : points = polynomialBases::find(MSH_QUA_16)->points; break;
- case 3 : points = polynomialBases::find(MSH_QUA_25)->points; break;
- case 4 : points = polynomialBases::find(MSH_QUA_36)->points; break;
- case 5 : points = polynomialBases::find(MSH_QUA_49)->points; break;
- case 6 : points = polynomialBases::find(MSH_QUA_64)->points; break;
- case 7 : points = polynomialBases::find(MSH_QUA_81)->points; break;
- case 8 : points = polynomialBases::find(MSH_QUA_100)->points; break;
+ case 1 : points = BasisFactory::create(MSH_QUA_9)->points; break;
+ case 2 : points = BasisFactory::create(MSH_QUA_16)->points; break;
+ case 3 : points = BasisFactory::create(MSH_QUA_25)->points; break;
+ case 4 : points = BasisFactory::create(MSH_QUA_36)->points; break;
+ case 5 : points = BasisFactory::create(MSH_QUA_49)->points; break;
+ case 6 : points = BasisFactory::create(MSH_QUA_64)->points; break;
+ case 7 : points = BasisFactory::create(MSH_QUA_81)->points; break;
+ case 8 : points = BasisFactory::create(MSH_QUA_100)->points; break;
default :
Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
break;
@@ -743,14 +742,14 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
case TYPE_HEX :
switch (nPts){
case 0: return;
- case 1: points = polynomialBases::find(MSH_HEX_27)->points; break;
- case 2: points = polynomialBases::find(MSH_HEX_64)->points; break;
- case 3: points = polynomialBases::find(MSH_HEX_125)->points; break;
- case 4: points = polynomialBases::find(MSH_HEX_216)->points; break;
- case 5: points = polynomialBases::find(MSH_HEX_343)->points; break;
- case 6: points = polynomialBases::find(MSH_HEX_512)->points; break;
- case 7: points = polynomialBases::find(MSH_HEX_729)->points; break;
- case 8: points = polynomialBases::find(MSH_HEX_1000)->points; break;
+ case 1: points = BasisFactory::create(MSH_HEX_27)->points; break;
+ case 2: points = BasisFactory::create(MSH_HEX_64)->points; break;
+ case 3: points = BasisFactory::create(MSH_HEX_125)->points; break;
+ case 4: points = BasisFactory::create(MSH_HEX_216)->points; break;
+ case 5: points = BasisFactory::create(MSH_HEX_343)->points; break;
+ case 6: points = BasisFactory::create(MSH_HEX_512)->points; break;
+ case 7: points = BasisFactory::create(MSH_HEX_729)->points; break;
+ case 8: points = BasisFactory::create(MSH_HEX_1000)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 4479b6d4eb0be7c76b6270f37f88f38a77e64c9a..ce341f915630702a417446451451f3d436fa4dbd 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -32,181 +32,15 @@ const nodalBasis* BasisFactory::create(int elementType) {
case(TYPE_PRI):
case(TYPE_TET):
case(TYPE_HEX):
- B = new polynomialBasis();
+ B = new polynomialBasis(elementType);
break;
case(TYPE_PYR):
- B = new pyramidalBasis();
+ B = new pyramidalBasis(elementType);
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);
- 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));
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index e361960798f64e78940a03ef197e3f4b28ce98b0..2638b5bd1327713c3cdb6558c7a98acd1b6fba8b 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -10,6 +10,8 @@
}*/
+
+
BergotBasis::BergotBasis(int p): order(p)
{
// allocate function information and fill
@@ -41,6 +43,7 @@ BergotBasis::BergotBasis(int p): order(p)
}
+
BergotBasis::~BergotBasis()
{
delete [] iOrder;
@@ -49,10 +52,7 @@ BergotBasis::~BergotBasis()
}
-int BergotBasis::size() const
-{
- return (2*order*order*order + 9*order*order + 13*order + 6)/6;
-}
+
void BergotBasis::f(double u, double v, double w, double* val) const
{
@@ -99,7 +99,9 @@ void BergotBasis::f(double u, double v, double w, double* val) const
}
}
-inline void BergotBasis::df(double u, double v, double w, double grads[][3]) const
+
+
+void BergotBasis::df(double u, double v, double w, double grads[][3]) const
{
std::vector<double> uFcts(order+1);
legendre.f(u,&(uFcts[0]));
diff --git a/Numeric/BergotBasis.h b/Numeric/BergotBasis.h
index 0a68e68ce46ae34cb3db150fd8799ae812c54c5f..f03d5082bd271b8278ca6d80ba63411bdb8f1faf 100644
--- a/Numeric/BergotBasis.h
+++ b/Numeric/BergotBasis.h
@@ -11,15 +11,17 @@
#include "jacobiPolynomials.h"
#include "legendrePolynomials.h"
-class BergotBasis : public nodalBasis {
+class BergotBasis {
public:
- BergotBasis() {};
BergotBasis(int p);
- ~BergotBasis();
+ virtual ~BergotBasis();
+
+ int size() const { const int n = order+1; return n*(n+1)*(2*n+1)/6; }
+
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 df(double u, double v, double w, double grads[][3]) const;
void initialize() {};
@@ -36,8 +38,6 @@ class BergotBasis : public nodalBasis {
//! list of Jacobi polynomials up to order p in function of index i (\f$ \alpha = 2*i + 2\f$)
std::map<int,JacobiPolynomials> jacobi;
- int size() const;
-
};
#endif
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index aa2bb2bb932f5a5e7a23a0dd2df366f98ef1419a..383a213684f25dd41036d695349ef2ba589e3562 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
fullMatrix.cpp
BasisFactory.cpp
+ nodalBasis.cpp
polynomialBasis.cpp
JacobianBasis.cpp
BergotBasis.cpp
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index c05d04d19e2689b39b7b83b38fe6ced8a3cdec51..e21b2a1736e189c95ad3d29dffa43c3823ee5184 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -8,6 +8,9 @@
#include "JacobianBasis.h"
#include <vector>
#include "polynomialBasis.h"
+#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
+#include "BasisFactory.h"
// Bezier Exponents
static fullMatrix<double> generate1DExponents(int order)
@@ -449,8 +452,8 @@ static fullMatrix<double> generateExponentsHex(int order)
fullMatrix<double> lineExp = generate1DExponents(order);
for (int i = 0; i < 2; ++i) {
- for (int j = 0; i < 2; ++j) {
- for (int k = 0; i < 2; ++k) {
+ for (int j = 0; j < 2; ++j) {
+ for (int k = 0; k < 2; ++k) {
exponents(index, 0) = lineExp(i, 0);
exponents(index, 1) = lineExp(j, 0);
exponents(index, 2) = lineExp(k, 0);
@@ -460,8 +463,8 @@ static fullMatrix<double> generateExponentsHex(int order)
}
for (int i = 2; i < lineExp.size1(); ++i) {
- for (int j = 0; i < 2; ++j) {
- for (int k = 0; i < 2; ++k) {
+ for (int j = 0; j < 2; ++j) {
+ for (int k = 0; k < 2; ++k) {
exponents(index, 0) = lineExp(i, 0);
exponents(index, 1) = lineExp(j, 0);
exponents(index, 2) = lineExp(k, 0);
@@ -470,8 +473,8 @@ static fullMatrix<double> generateExponentsHex(int order)
}
}
for (int i = 0; i < lineExp.size1(); ++i) {
- for (int j = 2; i < lineExp.size1(); ++j) {
- for (int k = 0; i < 2; ++k) {
+ for (int j = 2; j < lineExp.size1(); ++j) {
+ for (int k = 0; k < 2; ++k) {
exponents(index, 0) = lineExp(i, 0);
exponents(index, 1) = lineExp(j, 0);
exponents(index, 2) = lineExp(k, 0);
@@ -480,8 +483,8 @@ static fullMatrix<double> generateExponentsHex(int order)
}
}
for (int i = 0; i < lineExp.size1(); ++i) {
- for (int j = 0; i < lineExp.size1(); ++j) {
- for (int k = 2; i < lineExp.size1(); ++k) {
+ for (int j = 0; j < lineExp.size1(); ++j) {
+ for (int k = 2; k < lineExp.size1(); ++k) {
exponents(index, 0) = lineExp(i, 0);
exponents(index, 1) = lineExp(j, 0);
exponents(index, 2) = lineExp(k, 0);
@@ -507,64 +510,6 @@ static fullMatrix<double> generate1DPoints(int order)
return line;
}
-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> generatePointsQuadRecur(int order, bool serendip)
{
int nbPoints = serendip ? order*4 : (order+1)*(order+1);
@@ -618,143 +563,6 @@ static fullMatrix<double> generatePointsQuad(int order, bool serendip)
return point;
}
-
-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> generatePointsPrism(int order, bool serendip)
{
const double prism18Pts[18][3] = {
@@ -1126,9 +934,10 @@ static fullMatrix<double> generateSubDivisor
-static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points
- , const fullMatrix<double> &monomials, const fullMatrix<double> &coefficients)
+static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points,
+ const fullMatrix<double> &monomials, const fullMatrix<double> &coefficients)
{
+
int nbPts = points.size1();
int nbDofs = monomials.size1();
int dim = points.size2();
@@ -1208,10 +1017,14 @@ static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &poi
}
}
}
+ break;
}
return;
+
}
+
+
std::map<int, bezierBasis> bezierBasis::_bbs;
const bezierBasis *bezierBasis::find(int tag)
{
@@ -1220,80 +1033,96 @@ const bezierBasis *bezierBasis::find(int tag)
return &it->second;
bezierBasis &B = _bbs[tag];
- const polynomialBasis *F = polynomialBases::find(tag);
- int dimSimplex;
- std::vector< fullMatrix<double> > subPoints;
- switch (F->parentType) {
- case TYPE_PNT :
- B.numLagPts = 1;
- B.exponents = generate1DExponents(0);
- B.points = generate1DPoints(0);
- dimSimplex = 0;
- break;
- case TYPE_LIN : {
- B.numLagPts = 2;
- B.exponents = generate1DExponents(F->order);
- B.points = generate1DPoints(F->order);
- dimSimplex = 0;
- subPoints = generateSubPointsLine(0);
- break;
- }
- case TYPE_TRI : {
- B.numLagPts = 3;
- B.exponents = generateExponentsTriangle(F->order);
- B.points = gmshGeneratePointsTriangle(F->order,false);
- dimSimplex = 2;
- subPoints = generateSubPointsTriangle(F->order);
- break;
- }
- case TYPE_QUA : {
- B.numLagPts = 4;
- B.exponents = generateExponentsQuad(F->order);
- B.points = generatePointsQuad(F->order,false);
- dimSimplex = 0;
- subPoints = generateSubPointsQuad(F->order);
- // B.points.print("points");
- break;
- }
- case TYPE_TET : {
- B.numLagPts = 4;
- B.exponents = generateExponentsTetrahedron(F->order);
- B.points = gmshGeneratePointsTetrahedron(F->order,false);
- dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(F->order);
- break;
- }
- case TYPE_PRI : {
- B.numLagPts = 6;
- B.exponents = generateExponentsPrism(F->order);
- B.points = generatePointsPrism(F->order, false);
- dimSimplex = 2;
- subPoints = generateSubPointsPrism(F->order);
- break;
+ B.order = MElement::OrderFromTag(tag);
+
+ if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+ B.numLagPts = 5;
+ B.points = gmshGeneratePointsPyramid(B.order,false);
+ B.matrixLag2Bez.resize(B.points.size1(),B.points.size1(),0.);
+ B.matrixBez2Lag.resize(B.points.size1(),B.points.size1(),0.);
+ for (int i=0; i<B.matrixBez2Lag.size1(); ++i) {
+ B.matrixLag2Bez(i,i) = 1.;
+ B.matrixBez2Lag(i,i) = 1.;
}
- case TYPE_HEX : {
- B.numLagPts = 8;
- B.exponents = generateExponentsHex(F->order);
- B.points = generatePointsHex(F->order, false);
- dimSimplex = 0;
- subPoints = generateSubPointsHex(F->order);
- break;
- }
- default : {
- Msg::Error("Unknown function space %d: reverting to TET_1", tag);
- F = polynomialBases::find(MSH_TET_1);
- B.numLagPts = 4;
- B.exponents = generateExponentsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0, false);
- dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(0);
+// B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, F->order, dimSimplex);
+// B.numDivisions = (int) pow(2., (int) B.points.size2());
+ }
+ else {
+ const polynomialBasis *F = (polynomialBasis*)BasisFactory::create(tag);
+ int dimSimplex;
+ std::vector< fullMatrix<double> > subPoints;
+ switch (F->parentType) {
+ case TYPE_PNT :
+ B.numLagPts = 1;
+ B.exponents = generate1DExponents(0);
+ B.points = generate1DPoints(0);
+ dimSimplex = 0;
+ break;
+ case TYPE_LIN : {
+ B.numLagPts = 2;
+ B.exponents = generate1DExponents(F->order);
+ B.points = generate1DPoints(F->order);
+ dimSimplex = 0;
+ subPoints = generateSubPointsLine(0);
+ break;
+ }
+ case TYPE_TRI : {
+ B.numLagPts = 3;
+ B.exponents = generateExponentsTriangle(F->order);
+ B.points = gmshGeneratePointsTriangle(F->order,false);
+ dimSimplex = 2;
+ subPoints = generateSubPointsTriangle(F->order);
+ break;
+ }
+ case TYPE_QUA : {
+ B.numLagPts = 4;
+ B.exponents = generateExponentsQuad(F->order);
+ B.points = generatePointsQuad(F->order,false);
+ dimSimplex = 0;
+ subPoints = generateSubPointsQuad(F->order);
+ // B.points.print("points");
+ break;
+ }
+ case TYPE_TET : {
+ B.numLagPts = 4;
+ B.exponents = generateExponentsTetrahedron(F->order);
+ B.points = gmshGeneratePointsTetrahedron(F->order,false);
+ dimSimplex = 3;
+ subPoints = generateSubPointsTetrahedron(F->order);
+ break;
+ }
+ case TYPE_PRI : {
+ B.numLagPts = 6;
+ B.exponents = generateExponentsPrism(F->order);
+ B.points = generatePointsPrism(F->order, false);
+ dimSimplex = 2;
+ subPoints = generateSubPointsPrism(F->order);
+ break;
+ }
+ case TYPE_HEX : {
+ B.numLagPts = 8;
+ B.exponents = generateExponentsHex(F->order);
+ B.points = generatePointsHex(F->order, false);
+ dimSimplex = 0;
+ subPoints = generateSubPointsHex(F->order);
+ break;
+ }
+ default : {
+ Msg::Error("Unknown function space %d: reverting to TET_1", tag);
+ B.numLagPts = 4;
+ B.exponents = generateExponentsTetrahedron(0);
+ B.points = gmshGeneratePointsTetrahedron(0, false);
+ dimSimplex = 3;
+ subPoints = generateSubPointsTetrahedron(0);
+ break;
+ }
}
+ B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.points, F->order, dimSimplex);
+ B.matrixBez2Lag.invert(B.matrixLag2Bez);
+ B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, F->order, dimSimplex);
+ B.numDivisions = (int) pow(2., (int) B.points.size2());
}
- B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.points, F->order, dimSimplex);
- B.matrixBez2Lag.invert(B.matrixLag2Bez);
- B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, F->order, dimSimplex);
- B.numDivisions = (int) pow(2., (int) B.points.size2());
return &B;
}
@@ -1301,24 +1130,38 @@ std::map<int, JacobianBasis> JacobianBasis::_fs;
const JacobianBasis *JacobianBasis::find(int tag)
{
std::map<int, JacobianBasis>::const_iterator it = _fs.find(tag);
- if (it != _fs.end())
- return &it->second;
+ if (it != _fs.end()) return &it->second;
JacobianBasis &J = _fs[tag];
- const polynomialBasis *F = polynomialBases::find(tag);
- int jacobianOrder;
- switch (F->parentType) {
- case TYPE_PNT : jacobianOrder = 0; break;
- case TYPE_LIN : jacobianOrder = F->order - 1; break;
- case TYPE_TRI : jacobianOrder = 2 * (F->order - 1); break;
- case TYPE_QUA : jacobianOrder = 2 * F->order - 1; break;
- case TYPE_TET : jacobianOrder = 3 * (F->order - 1); break;
- case TYPE_PRI : jacobianOrder = 3 * F->order - 1; break;
- case TYPE_HEX : jacobianOrder = 3 * F->order - 1; break;
+ if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+ switch (tag) {
+ case MSH_PYR_5 : J.bezier = bezierBasis::find(MSH_PYR_14); break; // TODO: Order 1, Jac. "order" 2, check this
+ case MSH_PYR_14 : J.bezier = bezierBasis::find(MSH_PYR_91); break; // TODO: Order 2, Jac. "order" 5, check this
+ case MSH_PYR_30 : J.bezier = bezierBasis::find(MSH_PYR_285); break; // TODO: Order 3, Jac. "order" 8, check this
default :
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- jacobianOrder = 0;
+ Msg::Error("Unknown Jacobian function space for element type %d: reverting to PYR_5", tag);
+ J.bezier = bezierBasis::find(MSH_PYR_14);
+ break;
+ }
+ }
+ else {
+ const int parentType = MElement::ParentTypeFromTag(tag), order = MElement::OrderFromTag(tag);
+ int jacobianOrder;
+ switch (parentType) {
+ case TYPE_PNT : jacobianOrder = 0; break;
+ case TYPE_LIN : jacobianOrder = order - 1; break;
+ case TYPE_TRI : jacobianOrder = 2 * (order - 1); break;
+ case TYPE_QUA : jacobianOrder = 2 * order - 1; break;
+ case TYPE_TET : jacobianOrder = 3 * (order - 1); break;
+ case TYPE_PRI : jacobianOrder = 3 * order - 1; break;
+ case TYPE_HEX : jacobianOrder = 3 * order - 1; break;
+ default :
+ Msg::Error("Unknown Jacobian function space for element type %d: reverting to TET_4", tag);
+ jacobianOrder = 0;
+ break;
+ }
+ J.bezier = bezierBasis::find(polynomialBasis::getTag(parentType, jacobianOrder, false));
+ polynomialBasis *F = (polynomialBasis*)BasisFactory::create(tag);
+ generateGradShapes(J, J.bezier->points, F->monomials, F->coefficients);
}
- J.bezier = bezierBasis::find(polynomialBasis::getTag(F->parentType, jacobianOrder, false));
- generateGradShapes(J, J.bezier->points, F->monomials, F->coefficients);
return &J;
}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index d876be6ddaa9d725ccd37e7323f1f6dd0e0f20a7..af8326a968153c54ea6a901bd7db55b7b4995da2 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -13,6 +13,7 @@ class bezierBasis {
private :
static std::map<int, bezierBasis> _bbs;
public :
+ int order;
int numLagPts;
int numDivisions;
// the 'numLagPts' first exponents are related to 'Lagrangian' values
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e11af47cadf4e0572cdb990b374375b9b86364d6
--- /dev/null
+++ b/Numeric/nodalBasis.cpp
@@ -0,0 +1,1589 @@
+// 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 <limits>
+#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 and 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;
+}
+
+
+
+static void generate1dVertexClosure(nodalBasis::clCont &closure, int order)
+{
+ closure.clear();
+ closure.resize(2);
+ closure[0].push_back(0);
+ closure[1].push_back(order == 0 ? 0 : 1);
+ closure[0].type = MSH_PNT;
+ closure[1].type = MSH_PNT;
+}
+
+
+
+static void generate1dVertexClosureFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef, int order)
+{
+ closure.clear();
+ closure.resize(2);
+ closure[0].push_back(0);
+ if (order != 0) {
+ closure[0].push_back(1);
+ closure[1].push_back(1);
+ }
+ closure[1].push_back(0);
+ for (int i = 0; i < order - 1; i++) {
+ closure[0].push_back(2 + i);
+ closure[1].push_back(2 + order - 2 - i);
+ }
+ closureRef.resize(2);
+ closureRef[0] = 0;
+ closureRef[1] = 0;
+}
+
+
+
+static void getFaceClosureTet(int iFace, int iSign, int iRotate,
+ nodalBasis::closure &closure, int order)
+{
+ closure.clear();
+ closure.resize((order + 1) * (order + 2) / 2);
+ closure.type = nodalBasis::getTag(TYPE_TRI, order, false);
+
+ switch (order){
+ case 0:
+ closure[0] = 0;
+ break;
+ default:
+ int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
+ int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
+ for (int i = 0; i < 3; ++i){
+ int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
+ }
+ for (int i = 0;i < 3; ++i){
+ int edgenumber = iSign *
+ face[iFace][(6 + i * iSign + (-1 + iSign) / 2 + iRotate) % 3]; //- iSign * iRotate
+ for (int k = 0; k < (order - 1); k++){
+ if (edgenumber > 0)
+ closure[3 + i * (order - 1) + k] =
+ 4 + (edgenumber - 1) * (order - 1) + k;
+ else
+ closure[3 + i * (order - 1) + k] =
+ 4 + (-edgenumber) * (order - 1) - 1 - k;
+ }
+ }
+ int fi = 3 + 3 * (order - 1);
+ int ti = 4 + 6 * (order - 1);
+ int ndofff = (order - 3 + 2) * (order - 3 + 1) / 2;
+ ti = ti + iFace * ndofff;
+ for (int k = 0; k < order / 3; k++){
+ int orderint = order - 3 - k * 3;
+ if (orderint > 0){
+ for (int ci = 0; ci < 3 ; ci++){
+ int shift = (3 + iSign * ci + iRotate) % 3; //- iSign * iRotate
+ closure[fi + ci] = ti + shift;
+ }
+ fi = fi + 3; ti = ti + 3;
+ for (int l = 0; l < orderint - 1; l++){
+ for (int ei = 0; ei < 3; ei++){
+ int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3;
+ //- iSign * iRotate
+ if (iSign > 0)
+ closure[fi + ei * (orderint - 1) + l] =
+ ti + edgenumber * (orderint - 1) + l;
+ else
+ closure[fi + ei * (orderint - 1) + l] =
+ ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
+ }
+ }
+ fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
+ }
+ else {
+ closure[fi] = ti;
+ ti++;
+ fi++;
+ }
+ }
+ break;
+ }
+}
+
+
+
+static void generate2dEdgeClosureFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef,
+ int order, int nNod, bool serendip)
+{
+ closure.clear();
+ closure.resize(2*nNod);
+ closureRef.resize(2*nNod);
+ int shift = 0;
+ for (int corder = order; corder>=0; corder -= (nNod == 3 ? 3 : 2)) {
+ if (corder == 0) {
+ for (int r = 0; r < nNod ; r++){
+ closure[r].push_back(shift);
+ closure[r+nNod].push_back(shift);
+ }
+ break;
+ }
+ for (int r = 0; r < nNod ; r++){
+ for (int j = 0; j < nNod; j++) {
+ closure[r].push_back(shift + (r + j) % nNod);
+ closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
+ }
+ }
+ shift += nNod;
+ int n = nNod*(corder-1);
+ for (int r = 0; r < nNod ; r++){
+ for (int j = 0; j < n; j++){
+ closure[r].push_back(shift + (j + (corder - 1) * r) % n);
+ closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
+ }
+ }
+ shift += n;
+ if (serendip) break;
+ }
+ for (int r = 0; r < nNod*2 ; r++) {
+ closure[r].type = nodalBasis::getTag(TYPE_LIN, order);
+ closureRef[r] = 0;
+ }
+}
+
+
+
+static void addEdgeNodes(nodalBasis::clCont &closureFull, const int *edges, int order)
+{
+ if (order < 2)
+ return;
+ int numNodes = 0;
+ for (int i = 0; edges[i] >= 0; ++i) {
+ numNodes = std::max(numNodes, edges[i] + 1);
+ }
+ std::vector<std::vector<int> > nodes2edges(numNodes, std::vector<int>(numNodes, -1));
+ for (int i = 0; edges[i] >= 0; i += 2) {
+ nodes2edges[edges[i]][edges[i + 1]] = i;
+ nodes2edges[edges[i + 1]][edges[i]] = i + 1;
+ }
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ std::vector<int> &cl = closureFull[iClosure];
+ for (int iEdge = 0; edges[iEdge] >= 0; iEdge += 2) {
+ if (cl.empty())
+ continue;
+ int n0 = cl[edges[iEdge]];
+ int n1 = cl[edges[iEdge + 1]];
+ int oEdge = nodes2edges[n0][n1];
+ if (oEdge == -1)
+ Msg::Error("invalid p1 closure or invalid edges list");
+ for (int i = 0 ; i < order - 1; i++)
+ cl.push_back(numNodes + oEdge/2 * (order - 1) + ((oEdge % 2) ? order - 2 - i : i));
+ }
+ }
+}
+
+
+
+static void generateFaceClosureTet(nodalBasis::clCont &closure, int order)
+{
+ closure.clear();
+ for (int iRotate = 0; iRotate < 3; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 4; iFace++){
+ nodalBasis::closure closure_face;
+ getFaceClosureTet(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
+ }
+ }
+ }
+}
+
+
+
+static void generateFaceClosureTetFull(nodalBasis::clCont &closureFull,
+ std::vector<int> &closureRef,
+ int order, bool serendip)
+{
+ closureFull.clear();
+ //input :
+ static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
+ static const int edges[] = {0, 1, 1, 2, 2, 0, 3, 0, 3, 2, 3, 1, -1};
+ static const int faceOrientation[6] = {0,1,2,5,3,4};
+ closureFull.resize(24);
+ closureRef.resize(24);
+ for (int i = 0; i < 24; i ++)
+ closureRef[i] = 0;
+ if (order == 0) {
+ for (int i = 0; i < 24; i ++) {
+ closureFull[i].push_back(0);
+ }
+ return;
+ }
+ //Mapping for the p1 nodes
+ nodalBasis::clCont closure;
+ generateFaceClosureTet(closure, 1);
+ for (unsigned int i = 0; i < closureFull.size(); i++) {
+ std::vector<int> &clFull = closureFull[i];
+ std::vector<int> &cl = closure[i];
+ clFull.resize(4, -1);
+ closureRef[i] = 0;
+ for (unsigned int j = 0; j < cl.size(); j ++)
+ clFull[closure[0][j]] = cl[j];
+ for (int j = 0; j < 4; j ++)
+ if (clFull[j] == -1)
+ clFull[j] = (6 - clFull[(j + 1) % 4] - clFull[(j + 2) % 4] - clFull[(j + 3) % 4]);
+ }
+ int nodes2Faces[4][4][4];
+ for (int i = 0; i < 4; i++) {
+ for (int iRotate = 0; iRotate < 3; iRotate ++) {
+ short int n0 = faces[i][(3 - iRotate) % 3];
+ short int n1 = faces[i][(4 - iRotate) % 3];
+ short int n2 = faces[i][(5 - iRotate) % 3];
+ nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
+ nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
+ }
+ }
+ nodalBasis::clCont closureTriangles;
+ std::vector<int> closureTrianglesRef;
+ if (order >= 3)
+ generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
+ addEdgeNodes(closureFull, edges, order);
+ for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
+ //faces
+ std::vector<int> &cl = closureFull[iClosure];
+ if (order >= 3) {
+ for (int iFace = 0; iFace < 4; iFace ++) {
+ int n0 = cl[faces[iFace][0]];
+ int n1 = cl[faces[iFace][1]];
+ int n2 = cl[faces[iFace][2]];
+ short int id = nodes2Faces[n0][n1][n2];
+ short int iTriClosure = faceOrientation [id % 6];
+ short int idFace = id/6;
+ int nOnFace = closureTriangles[iTriClosure].size();
+ for (int i = 0; i < nOnFace; i++) {
+ cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace +
+ closureTriangles[iTriClosure][i]);
+ }
+ }
+ }
+ }
+ if (order >= 4 && !serendip) {
+ nodalBasis::clCont insideClosure;
+ std::vector<int> fakeClosureRef;
+ generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
+ for (unsigned int i = 0; i < closureFull.size(); i ++) {
+ unsigned int shift = closureFull[i].size();
+ for (unsigned int j = 0; j < insideClosure[i].size(); j++)
+ closureFull[i].push_back(insideClosure[i][j] + shift);
+ }
+ }
+}
+
+
+
+void rotateHex(int iFace, int iRot, int iSign, double uI, double vI,
+ double &uO, double &vO, double &wO)
+{
+ if (iSign<0){
+ double tmp=uI;
+ uI=vI;
+ vI=tmp;
+ }
+ for (int i=0; i < iRot; i++){
+ double tmp=uI;
+ uI=-vI;
+ vI=tmp;
+ }
+ switch (iFace) {
+ case 0: uO = vI; vO = uI; wO=-1; break;
+ case 1: uO = uI; vO = -1; wO=vI; break;
+ case 2: uO =-1; vO = vI; wO=uI; break;
+ case 3: uO = 1; vO = uI; wO=vI; break;
+ case 4: uO =-uI; vO = 1; wO=vI; break;
+ case 5: uO =uI; vO = vI; wO=1; break;
+ }
+}
+
+
+
+void rotateHexFull(int iFace, int iRot, int iSign, double uI, double vI, double wI, double &uO, double &vO, double &wO)
+{
+ switch (iFace) {
+ case 0: uO = uI; vO = vI; wO = wI; break;
+ case 1: uO = wI; vO = uI; wO = vI; break;
+ case 2: uO = vI; vO = wI; wO = uI; break;
+ case 3: uO = wI; vO = vI; wO =-uI; break;
+ case 4: uO = wI; vO =-uI; wO =-vI; break;
+ case 5: uO = vI; vO = uI; wO =-wI; break;
+ }
+ for (int i = 0; i < iRot; i++){
+ double tmp = uO;
+ uO = -vO;
+ vO = tmp;
+ }
+ if (iSign<0){
+ double tmp = uO;
+ uO = vO;
+ vO = tmp;
+ }
+}
+
+
+
+inline static double pow2(double x)
+{
+ return x * x;
+}
+
+/*
+static void checkClosure(int tag) {
+ printf("TYPE = %i\n", tag);
+ const nodalBasis &fs = *nodalBases::find(tag);
+ for (int i = 0; i < fs.closures.size(); ++i) {
+ const std::vector<int> &clRef = fs.closures[fs.closureRef[i]];
+ const std::vector<int> &cl = fs.closures[i];
+ const std::vector<int> &clFull = fs.fullClosures[i];
+ printf("i = %i\n", i);
+ for (int j = 0; j < cl.size(); ++j) {
+ printf("%3i ", clFull[clRef[j]]);
+ }
+ printf("\n");
+ for (int j = 0; j < cl.size(); ++j) {
+ printf("%3i ", cl[j]);
+ }
+ printf("\n");
+ }
+}
+*/
+
+
+
+static void generateFaceClosureHex(nodalBasis::clCont &closure, int order,
+ bool serendip, const fullMatrix<double> &points)
+{
+ closure.clear();
+ const nodalBasis &fsFace = *BasisFactory::create(nodalBasis::getTag(TYPE_QUA, order, serendip));
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 6; iFace++) {
+ nodalBasis::closure cl;
+ cl.type = fsFace.type;
+ cl.resize(fsFace.points.size1());
+ for (unsigned int iNode = 0; iNode < cl.size(); ++iNode) {
+ double u,v,w;
+ rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0),
+ fsFace.points(iNode, 1), u, v, w);
+ cl[iNode] = 0;
+ double D = std::numeric_limits<double>::max();
+ for (int jNode = 0; jNode < points.size1(); ++jNode) {
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
+ if (d < D) {
+ cl[iNode] = jNode;
+ D = d;
+ }
+ }
+ }
+ closure.push_back(cl);
+ }
+ }
+ }
+}
+
+
+
+static void generateFaceClosureHexFull(nodalBasis::clCont &closure,
+ std::vector<int> &closureRef,
+ int order, bool serendip,
+ const fullMatrix<double> &points)
+{
+ closure.clear();
+ int clId = 0;
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 6; iFace++) {
+ nodalBasis::closure cl;
+ cl.resize(points.size1());
+ for (int iNode = 0; iNode < points.size1(); ++iNode) {
+ double u,v,w;
+ rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1),
+ points(iNode, 2), u, v, w);
+ int J = 0;
+ double D = std::numeric_limits<double>::max();
+ for (int jNode = 0; jNode < points.size1(); ++jNode) {
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
+ if (d < D) {
+ J = jNode;
+ D = d;
+ }
+ }
+ cl[J] = iNode;
+ }
+ closure.push_back(cl);
+ closureRef.push_back(0);
+ clId ++;
+ }
+ }
+ }
+}
+
+
+
+static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
+ nodalBasis::closure &closure, int order)
+{
+ if (order > 2)
+ Msg::Error("FaceClosure not implemented for prisms of order %d",order);
+ bool isTriangle = iFace<2;
+ int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
+ closure.clear();
+ if (isTriangle && iRotate > 2) return;
+ closure.resize(nNodes);
+ if (order==0) {
+ closure[0] = 0;
+ return;
+ }
+ int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2},
+ {1, 2, 5, 4}};
+ int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15},
+ {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
+ // int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
+ // {8, 13, 11, 7}, {9, 11, 14, 10}};
+ int nVertex = isTriangle ? 3 : 4;
+ closure.type = nodalBasis::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
+ for (int i = 0; i < nVertex; ++i){
+ int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
+ }
+ if (order==2) {
+ for (int i = 0; i < nVertex; ++i){
+ int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;
+ //- iSign * iRotate
+ closure[nVertex+i] = order2node[iFace][k];
+ }
+ if (!isTriangle)
+ closure[nNodes-1] = order2node[iFace][4]; // center
+ }
+}
+
+
+
+static void generateFaceClosurePrism(nodalBasis::clCont &closure, int order)
+{
+ closure.clear();
+ for (int iRotate = 0; iRotate < 4; iRotate++){
+ for (int iSign = 1; iSign >= -1; iSign -= 2){
+ for (int iFace = 0; iFace < 5; iFace++){
+ nodalBasis::closure closure_face;
+ getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
+ }
+ }
+ }
+}
+
+
+
+static void generateFaceClosurePrismFull(nodalBasis::clCont &closureFull,
+ std::vector<int> &closureRef, int order)
+{
+ nodalBasis::clCont closure;
+ closureFull.clear();
+ closureFull.resize(40);
+ closureRef.resize(40);
+ generateFaceClosurePrism(closure, 1);
+ int ref3 = -1, ref4a = -1, ref4b = -1;
+ for (unsigned int i = 0; i < closure.size(); i++) {
+ std::vector<int> &clFull = closureFull[i];
+ std::vector<int> &cl = closure[i];
+ if (cl.size() == 0)
+ continue;
+ clFull.resize(6, -1);
+ int &ref = cl.size() == 3 ? ref3 : (cl[0] / 3 + cl[1] / 3) % 2 ? ref4b : ref4a;
+ if (ref == -1)
+ ref = i;
+ closureRef[i] = ref;
+ for (unsigned int j = 0; j < cl.size(); j ++)
+ clFull[closure[ref][j]] = cl[j];
+ for (int j = 0; j < 6; j ++) {
+ if (clFull[j] == -1) {
+ int k = ((j / 3) + 1) % 2 * 3;
+ int sum = (clFull[k + (j + 1) % 3] + clFull[k + (j + 2) % 3]);
+ clFull[j] = ((sum / 6 + 1) % 2) * 3 + (12 - sum) % 3;
+ }
+ }
+ }
+ static const int edges[] = {0, 1, 0, 2, 0, 3, 1, 2, 1, 4, 2, 5,
+ 3, 4, 3, 5, 4, 5, -1};
+ addEdgeNodes(closureFull, edges, order);
+ if ( order < 2 )
+ return;
+ // face center nodes for p2 prism
+ static const int faces[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3},
+ {0, 3, 5, 2}, {1, 2, 5, 4}};
+
+ if ( order == 2 ) {
+ int nextFaceNode = 15;
+ int numFaces = 5;
+ int numFaceNodes = 4;
+ std::map<int,int> nodeSum2Face;
+ for (int iFace = 0; iFace < numFaces ; iFace ++) {
+ int nodeSum = 0;
+ for (int iNode = 0; iNode < numFaceNodes; iNode++ ) {
+ nodeSum += faces[iFace][iNode];
+ }
+ nodeSum2Face[nodeSum] = iFace;
+ }
+ for (unsigned int i = 0; i < closureFull.size(); i++ ) {
+ if (closureFull[i].empty())
+ continue;
+ for (int iFace = 0; iFace < numFaces; iFace++ ) {
+ int nodeSum = 0;
+ for (int iNode = 0; iNode < numFaceNodes; iNode++)
+ nodeSum += faces[iFace][iNode] < 0 ? faces[iFace][iNode] :
+ closureFull[i][ faces[iFace][iNode] ];
+ std::map<int,int>::iterator it = nodeSum2Face.find(nodeSum);
+ if (it == nodeSum2Face.end() )
+ Msg::Error("Prism face not found");
+ int mappedFaceId = it->second;
+ if ( mappedFaceId > 1) {
+ closureFull[i].push_back(nextFaceNode + mappedFaceId - 2);
+ }
+ }
+ }
+ } else {
+ Msg::Error("FaceClosureFull not implemented for prisms of order %d",order);
+ }
+
+}
+
+
+
+static void generate2dEdgeClosure(nodalBasis::clCont &closure, int order, int nNod = 3)
+{
+ closure.clear();
+ closure.resize(2*nNod);
+ for (int j = 0; j < nNod ; j++){
+ closure[j].push_back(j);
+ closure[j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back(j);
+ for (int i=0; i < order-1; i++){
+ closure[j].push_back( nNod + (order-1)*j + i );
+ closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
+ }
+ closure[j].type = closure[nNod+j].type = nodalBasis::getTag(TYPE_LIN, order);
+ }
+}
+
+
+
+static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
+{
+ closure.clear();
+ closure.resize(nb);
+ for (int i=0; i < nb; i++) {
+ closure[i].push_back(0);
+ closure[i].type = MSH_PNT;
+ }
+}
+
+
+
+nodalBasis::nodalBasis(int tag)
+{
+
+ type = tag;
+
+ 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_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 = false; 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_222 : 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_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;
+ }
+
+
+ switch (parentType) {
+ case TYPE_PNT :
+ numFaces = 1;
+ dimension = 0;
+ points = generate1DPoints(0);
+ break;
+ case TYPE_LIN :
+ numFaces = 2;
+ dimension = 1;
+ points = generate1DPoints(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);
+ generateClosureOrder0(fullClosures, 6);
+ closureRef.resize(6, 0);
+ }
+ else {
+ generate2dEdgeClosure(closures, order);
+ generate2dEdgeClosureFull(fullClosures, closureRef, order, 3, serendip);
+ }
+ break;
+ case TYPE_QUA :
+ numFaces = 4;
+ dimension = 2;
+ points = gmshGeneratePointsQuad(order, serendip);
+ 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 :
+ numFaces = 4;
+ dimension = 3;
+ points = gmshGeneratePointsTetrahedron(order, serendip);
+ 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 :
+ numFaces = 5;
+ dimension = 3;
+ points = gmshGeneratePointsPrism(order, serendip);
+ 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 :
+ numFaces = 6;
+ dimension = 3;
+ points = gmshGeneratePointsHex(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 6a283f4eb754829054f4b00202543ec22a422672..da1bc0d723db5e373bfaf0d334dec90da7b45d28 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -3,8 +3,8 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef BASIS_H
-#define BASIS_H
+#ifndef NODALBASIS_H
+#define NODALBASIS_H
#include "fullMatrix.h"
#include "GmshDefines.h"
@@ -16,20 +16,21 @@ class nodalBasis {
bool serendip;
fullMatrix<double> points;
- virtual void initialize() = 0;
+ nodalBasis(int tag);
+ virtual ~nodalBasis() {}
// Basis functions evaluation
- inline virtual void f(double u, double v, double w, double *sf) const {Msg::Fatal("Not implemented");};
- inline virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const {Msg::Fatal("Not implemented");};
+ virtual void f(double u, double v, double w, double *sf) const {Msg::Fatal("Not implemented");};
+ virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const {Msg::Fatal("Not implemented");};
// Basis functions gradients evaluation
- inline virtual void df(double u, double v, double w, double grads[][3]) const {Msg::Fatal("Not implemented");};
- inline virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const {Msg::Fatal("Not implemented");};
+ virtual void df(double u, double v, double w, double grads[][3]) const {Msg::Fatal("Not implemented");};
+ virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const {Msg::Fatal("Not implemented");};
- inline virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
- inline virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
+ virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
+ virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
- inline virtual int getNumShapeFunctions() const {Msg::Fatal("Not implemented"); return -1;}
+ virtual int getNumShapeFunctions() const {Msg::Fatal("Not implemented"); return -1;}
class closure : public std::vector<int> {
public:
@@ -46,100 +47,155 @@ class nodalBasis {
// quad face)
clCont closures, fullClosures;
std::vector<int> closureRef;
- inline virtual int getClosureType(int id) const {Msg::Fatal("Not implemented"); return -1;}
- inline virtual const std::vector<int> &getClosure(int id) const {Msg::Fatal("Not implemented"); std::vector<int> *ret=NULL; return *ret;}
- inline virtual const std::vector<int> &getFullClosure(int id) const {Msg::Fatal("Not implemented"); std::vector<int> *ret=NULL; return *ret;}
- inline virtual int getClosureId(int iFace, int iSign=1, int iRot=0) const {Msg::Fatal("Not implemented"); return -1;}
- inline virtual void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const {Msg::Fatal("Not implemented");iFace=-1; iSign=-1; iRot=-1;}
-
- static int getTag(int parentTag, int order, bool serendip = false)
- {
- switch (parentTag) {
- case TYPE_PNT :
- return MSH_PNT;
- case TYPE_LIN :
- switch(order) {
- case 0 : return MSH_LIN_1;
- case 1 : return MSH_LIN_2;
- case 2 : return MSH_LIN_3;
- case 3 : return MSH_LIN_4;
- case 4 : return MSH_LIN_5;
- case 5 : return MSH_LIN_6;
- case 6 : return MSH_LIN_7;
- case 7 : return MSH_LIN_8;
- case 8 : return MSH_LIN_9;
- case 9 : return MSH_LIN_10;
- case 10: return MSH_LIN_11;
- default : Msg::Error("line order %i unknown", order); return 0;
- }
- case TYPE_TRI :
- switch(order) {
- case 0 : return MSH_TRI_1;
- case 1 : return MSH_TRI_3;
- case 2 : return MSH_TRI_6;
- case 3 : return serendip ? MSH_TRI_9 : MSH_TRI_10;
- case 4 : return serendip ? MSH_TRI_12 : MSH_TRI_15;
- case 5 : return serendip ? MSH_TRI_15I: MSH_TRI_21;
- case 6 : return serendip ? MSH_TRI_18 : MSH_TRI_28;
- case 7 : return serendip ? MSH_TRI_21I: MSH_TRI_36;
- case 8 : return serendip ? MSH_TRI_24 : MSH_TRI_45;
- case 9 : return serendip ? MSH_TRI_27 : MSH_TRI_55;
- case 10: return serendip ? MSH_TRI_30 : MSH_TRI_66;
- default : Msg::Error("triangle order %i unknown", order); return 0;
- }
- case TYPE_QUA :
- switch(order) {
- case 0 : return MSH_QUA_1;
- case 1 : return MSH_QUA_4;
- case 2 : return serendip ? MSH_QUA_8 : MSH_QUA_9;
- case 3 : return serendip ? MSH_QUA_12 : MSH_QUA_16;
- case 4 : return serendip ? MSH_QUA_16I: MSH_QUA_25;
- case 5 : return serendip ? MSH_QUA_20 : MSH_QUA_36;
- case 6 : return serendip ? MSH_QUA_24 : MSH_QUA_49;
- case 7 : return serendip ? MSH_QUA_28 : MSH_QUA_64;
- case 8 : return serendip ? MSH_QUA_32 : MSH_QUA_81;
- case 9 : return serendip ? MSH_QUA_36I: MSH_QUA_100;
- case 10: return serendip ? MSH_QUA_40 : MSH_QUA_121;
- default : Msg::Error("quad order %i unknown", order); return 0;
- }
- case TYPE_TET :
- switch(order) {
- case 0 : return MSH_TET_1;
- case 1 : return MSH_TET_4;
- case 2 : return MSH_TET_10;
- case 3 : return MSH_TET_20;
- case 4 : return serendip ? MSH_TET_34 : MSH_TET_35;
- case 5 : return serendip ? MSH_TET_52 : MSH_TET_56;
- case 6 : return serendip ? MSH_TET_74 : MSH_TET_84;
- case 7 : return serendip ? MSH_TET_100: MSH_TET_120;
- case 8 : return serendip ? MSH_TET_130: MSH_TET_165;
- case 9 : return serendip ? MSH_TET_164: MSH_TET_220;
- case 10: return serendip ? MSH_TET_202: MSH_TET_286;
- default : Msg::Error("terahedron order %i unknown", order); return 0;
- }
- case TYPE_HEX :
- switch(order) {
- case 1 : return MSH_HEX_8;
- case 2 : return serendip ? MSH_HEX_20 : MSH_HEX_27;
- case 3 : return serendip ? MSH_HEX_56 : MSH_HEX_64;
- case 4 : return serendip ? MSH_HEX_98 : MSH_HEX_125;
- case 5 : return serendip ? MSH_HEX_152: MSH_HEX_216;
- case 6 : return serendip ? MSH_HEX_222: MSH_HEX_343;
- case 7 : return serendip ? MSH_HEX_296: MSH_HEX_512;
- case 8 : return serendip ? MSH_HEX_386: MSH_HEX_729;
- case 9 : return serendip ? MSH_HEX_488: MSH_HEX_1000;
- default : Msg::Error("hexahedron order %i unknown", order); return 0;
- }
- case TYPE_PRI :
- switch(order) {
- case 0 : return MSH_PRI_1;
- case 1 : return MSH_PRI_6;
- case 2 : return MSH_PRI_18;
- default : Msg::Error("prism order %i unknown", order); return 0;
- }
- default : Msg::Error("unknown element type %i", parentTag); return 0;
+
+ // for a given face/edge, with both a sign and a rotation, give an
+ // ordered list of nodes on this face/edge
+ virtual int getClosureType(int id) const { return closures[id].type; }
+ virtual const std::vector<int> &getClosure(int id) const { return closures[id]; }
+ virtual const std::vector<int> &getFullClosure(int id) const { return fullClosures[id]; }
+ inline int getClosureId(int iFace, int iSign=1, int iRot=0) const;
+ inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const;
+
+ static inline int getTag(int parentTag, int order, bool serendip = false);
+
+};
+
+
+
+inline int nodalBasis::getClosureId(int iFace, int iSign, int iRot) const
+{
+ return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
+}
+
+
+
+inline void nodalBasis::breakClosureId(int i, int &iFace, int &iSign, int &iRot) const
+{
+ iFace = i % numFaces;
+ i = (i - iFace)/numFaces;
+ iSign = i % 2;
+ iRot = (i - iSign) / 2;
+}
+
+
+
+inline int nodalBasis::getTag(int parentTag, int order, bool serendip)
+{
+ switch (parentTag) {
+ case TYPE_PNT :
+ return MSH_PNT;
+ case TYPE_LIN :
+ switch(order) {
+ case 0 : return MSH_LIN_1;
+ case 1 : return MSH_LIN_2;
+ case 2 : return MSH_LIN_3;
+ case 3 : return MSH_LIN_4;
+ case 4 : return MSH_LIN_5;
+ case 5 : return MSH_LIN_6;
+ case 6 : return MSH_LIN_7;
+ case 7 : return MSH_LIN_8;
+ case 8 : return MSH_LIN_9;
+ case 9 : return MSH_LIN_10;
+ case 10: return MSH_LIN_11;
+ default : Msg::Error("line order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_TRI :
+ switch(order) {
+ case 0 : return MSH_TRI_1;
+ case 1 : return MSH_TRI_3;
+ case 2 : return MSH_TRI_6;
+ case 3 : return serendip ? MSH_TRI_9 : MSH_TRI_10;
+ case 4 : return serendip ? MSH_TRI_12 : MSH_TRI_15;
+ case 5 : return serendip ? MSH_TRI_15I: MSH_TRI_21;
+ case 6 : return serendip ? MSH_TRI_18 : MSH_TRI_28;
+ case 7 : return serendip ? MSH_TRI_21I: MSH_TRI_36;
+ case 8 : return serendip ? MSH_TRI_24 : MSH_TRI_45;
+ case 9 : return serendip ? MSH_TRI_27 : MSH_TRI_55;
+ case 10: return serendip ? MSH_TRI_30 : MSH_TRI_66;
+ default : Msg::Error("triangle order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_QUA :
+ switch(order) {
+ case 0 : return MSH_QUA_1;
+ case 1 : return MSH_QUA_4;
+ case 2 : return serendip ? MSH_QUA_8 : MSH_QUA_9;
+ case 3 : return serendip ? MSH_QUA_12 : MSH_QUA_16;
+ case 4 : return serendip ? MSH_QUA_16I: MSH_QUA_25;
+ case 5 : return serendip ? MSH_QUA_20 : MSH_QUA_36;
+ case 6 : return serendip ? MSH_QUA_24 : MSH_QUA_49;
+ case 7 : return serendip ? MSH_QUA_28 : MSH_QUA_64;
+ case 8 : return serendip ? MSH_QUA_32 : MSH_QUA_81;
+ case 9 : return serendip ? MSH_QUA_36I: MSH_QUA_100;
+ case 10: return serendip ? MSH_QUA_40 : MSH_QUA_121;
+ default : Msg::Error("quad order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_TET :
+ switch(order) {
+ case 0 : return MSH_TET_1;
+ case 1 : return MSH_TET_4;
+ case 2 : return MSH_TET_10;
+ case 3 : return MSH_TET_20;
+ case 4 : return serendip ? MSH_TET_34 : MSH_TET_35;
+ case 5 : return serendip ? MSH_TET_52 : MSH_TET_56;
+ case 6 : return serendip ? MSH_TET_74 : MSH_TET_84;
+ case 7 : return serendip ? MSH_TET_100: MSH_TET_120;
+ case 8 : return serendip ? MSH_TET_130: MSH_TET_165;
+ case 9 : return serendip ? MSH_TET_164: MSH_TET_220;
+ case 10: return serendip ? MSH_TET_202: MSH_TET_286;
+ default : Msg::Error("terahedron order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_HEX :
+ switch(order) {
+ case 1 : return MSH_HEX_8;
+ case 2 : return serendip ? MSH_HEX_20 : MSH_HEX_27;
+ case 3 : return serendip ? MSH_HEX_56 : MSH_HEX_64;
+ case 4 : return serendip ? MSH_HEX_98 : MSH_HEX_125;
+ case 5 : return serendip ? MSH_HEX_152: MSH_HEX_216;
+ case 6 : return serendip ? MSH_HEX_222: MSH_HEX_343;
+ case 7 : return serendip ? MSH_HEX_296: MSH_HEX_512;
+ case 8 : return serendip ? MSH_HEX_386: MSH_HEX_729;
+ case 9 : return serendip ? MSH_HEX_488: MSH_HEX_1000;
+ default : Msg::Error("hexahedron order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_PRI :
+ switch(order) {
+ case 0 : return MSH_PRI_1;
+ case 1 : return MSH_PRI_6;
+ case 2 : return serendip ? MSH_PRI_15 : MSH_PRI_18;
+ case 3 : return serendip ? MSH_PRI_38 : MSH_PRI_40;
+ case 4 : return serendip ? MSH_PRI_66 : MSH_PRI_75;
+ case 5 : return serendip ? MSH_PRI_102 : MSH_PRI_126;
+ case 6 : return serendip ? MSH_PRI_146 : MSH_PRI_196;
+ case 7 : return serendip ? MSH_PRI_198 : MSH_PRI_288;
+ case 8 : return serendip ? MSH_PRI_258 : MSH_PRI_405;
+ case 9 : return serendip ? MSH_PRI_326 : MSH_PRI_550;
+ default : Msg::Error("prism order %i unknown", order); return 0;
}
+ break;
+ case TYPE_PYR :
+ switch(order) {
+ case 0 : return MSH_PYR_1;
+ case 1 : return MSH_PYR_5;
+ case 2 : return serendip ? MSH_PYR_13 : MSH_PYR_14;
+ case 3: return serendip ? MSH_PYR_29 : MSH_PYR_30;
+ case 4: return serendip ? MSH_PYR_50 : MSH_PYR_55;
+ case 5: return serendip ? MSH_PYR_77 : MSH_PYR_91;
+ case 6: return serendip ? MSH_PYR_110 : MSH_PYR_140;
+ case 7: return serendip ? MSH_PYR_149 : MSH_PYR_204;
+ case 8: return serendip ? MSH_PYR_194 : MSH_PYR_285;
+ case 9: return serendip ? MSH_PYR_245 : MSH_PYR_385;
+ default : Msg::Error("pyramid order %i unknown", order); return 0;
+ }
+ break;
+ default : Msg::Error("unknown element type %i", parentTag); return 0;
}
-
-};
+}
+
+
+
#endif
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index f7d53bb6923e803a5e8a8e034fd4050630c6a019..c32c1adec670c0079e87130ad0246f33b8eacf65 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -43,6 +43,7 @@ static void printClosure(polynomialBasis::clCont &fullClosure,
*/
+
fullMatrix<double> generate1DMonomials(int order)
{
fullMatrix<double> monomials(order + 1, 1);
@@ -50,18 +51,7 @@ fullMatrix<double> generate1DMonomials(int order)
return monomials;
}
-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> generatePascalTriangle(int order)
{
@@ -77,6 +67,8 @@ static fullMatrix<double> generatePascalTriangle(int order)
return monomials;
}
+
+
// generate the exterior hull of the Pascal triangle for Serendipity element
static fullMatrix<double> generatePascalSerendipityTriangle(int order)
@@ -107,6 +99,8 @@ static fullMatrix<double> generatePascalSerendipityTriangle(int order)
return monomials;
}
+
+
// generate all monomials xi^m * eta^n with n and m <= order
static fullMatrix<double> generatePascalQuad(int order)
@@ -136,6 +130,8 @@ static fullMatrix<double> generatePascalQuad(int order)
â‹® â‹®
*/
+
+
// generate all monomials xi^m * eta^n with n and m <= order
static fullMatrix<double> generatePascalHex(int order, bool serendip)
@@ -159,6 +155,8 @@ static fullMatrix<double> generatePascalHex(int order, bool serendip)
return monomials;
}
+
+
static fullMatrix<double> generatePascalQuadSerendip(int order)
{
fullMatrix<double> monomials( (order)*4, 2);
@@ -188,6 +186,8 @@ static fullMatrix<double> generatePascalQuadSerendip(int order)
return monomials;
}
+
+
// generate the monomials subspace of all monomials of order exactly == p
static fullMatrix<double> generateMonomialSubspace(int dim, int p)
@@ -222,6 +222,8 @@ static fullMatrix<double> generateMonomialSubspace(int dim, int p)
return monomials;
}
+
+
// generate external hull of the Pascal tetrahedron
static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
@@ -250,6 +252,8 @@ static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
return monomials;
}
+
+
// generate Pascal tetrahedron
static fullMatrix<double> generatePascalTetrahedron(int order)
@@ -269,6 +273,8 @@ static fullMatrix<double> generatePascalTetrahedron(int order)
return monomials;
}
+
+
// generate Pascal prism
static fullMatrix<double> generatePascalPrism(int order)
@@ -307,570 +313,7 @@ static fullMatrix<double> generatePascalPrism(int order)
return monomials;
}
-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> 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> generateLagrangeMonomialCoefficients
(const fullMatrix<double>& monomial, const fullMatrix<double>& point)
@@ -899,845 +342,386 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
return coefficient;
}
-static void generate1dVertexClosure(polynomialBasis::clCont &closure, int order)
-{
- closure.clear();
- closure.resize(2);
- closure[0].push_back(0);
- closure[1].push_back(order == 0 ? 0 : 1);
- closure[0].type = MSH_PNT;
- closure[1].type = MSH_PNT;
-}
-static void generate1dVertexClosureFull(polynomialBasis::clCont &closure,
- std::vector<int> &closureRef, int order)
-{
- closure.clear();
- closure.resize(2);
- closure[0].push_back(0);
- if (order != 0) {
- closure[0].push_back(1);
- closure[1].push_back(1);
- }
- closure[1].push_back(0);
- for (int i = 0; i < order - 1; i++) {
- closure[0].push_back(2 + i);
- closure[1].push_back(2 + order - 2 - i);
- }
- closureRef.resize(2);
- closureRef[0] = 0;
- closureRef[1] = 0;
-}
-static void getFaceClosureTet(int iFace, int iSign, int iRotate,
- polynomialBasis::closure &closure, int order)
+polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
{
- closure.clear();
- closure.resize((order + 1) * (order + 2) / 2);
- closure.type = nodalBasis::getTag(TYPE_TRI, order, false);
- switch (order){
- case 0:
- closure[0] = 0;
+ switch (parentType) {
+ case TYPE_PNT :
+ monomials = generate1DMonomials(0);
break;
- default:
- int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
- int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
- for (int i = 0; i < 3; ++i){
- int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
- closure[i] = order1node[iFace][k];
- }
- for (int i = 0;i < 3; ++i){
- int edgenumber = iSign *
- face[iFace][(6 + i * iSign + (-1 + iSign) / 2 + iRotate) % 3]; //- iSign * iRotate
- for (int k = 0; k < (order - 1); k++){
- if (edgenumber > 0)
- closure[3 + i * (order - 1) + k] =
- 4 + (edgenumber - 1) * (order - 1) + k;
- else
- closure[3 + i * (order - 1) + k] =
- 4 + (-edgenumber) * (order - 1) - 1 - k;
- }
- }
- int fi = 3 + 3 * (order - 1);
- int ti = 4 + 6 * (order - 1);
- int ndofff = (order - 3 + 2) * (order - 3 + 1) / 2;
- ti = ti + iFace * ndofff;
- for (int k = 0; k < order / 3; k++){
- int orderint = order - 3 - k * 3;
- if (orderint > 0){
- for (int ci = 0; ci < 3 ; ci++){
- int shift = (3 + iSign * ci + iRotate) % 3; //- iSign * iRotate
- closure[fi + ci] = ti + shift;
- }
- fi = fi + 3; ti = ti + 3;
- for (int l = 0; l < orderint - 1; l++){
- for (int ei = 0; ei < 3; ei++){
- int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3;
- //- iSign * iRotate
- if (iSign > 0)
- closure[fi + ei * (orderint - 1) + l] =
- ti + edgenumber * (orderint - 1) + l;
- else
- closure[fi + ei * (orderint - 1) + l] =
- ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
- }
- }
- fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
- }
- else {
- closure[fi] = ti;
- ti++;
- fi++;
- }
- }
+ case TYPE_LIN :
+ monomials = generate1DMonomials(order);
+ break;
+ case TYPE_TRI :
+ monomials = serendip ? generatePascalSerendipityTriangle(order) :
+ generatePascalTriangle(order);
+ break;
+ case TYPE_QUA :
+ monomials = serendip ? generatePascalQuadSerendip(order) :
+ generatePascalQuad(order);
+ break;
+ case TYPE_TET :
+ monomials = serendip ? generatePascalSerendipityTetrahedron(order) :
+ generatePascalTetrahedron(order);
+ break;
+ case TYPE_PRI :
+ monomials = generatePascalPrism(order);
+ break;
+ case TYPE_HEX :
+ monomials = generatePascalHex(order, serendip);
break;
}
-}
-static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure,
- std::vector<int> &closureRef,
- int order, int nNod, bool serendip)
-{
- closure.clear();
- closure.resize(2*nNod);
- closureRef.resize(2*nNod);
- int shift = 0;
- for (int corder = order; corder>=0; corder -= (nNod == 3 ? 3 : 2)) {
- if (corder == 0) {
- for (int r = 0; r < nNod ; r++){
- closure[r].push_back(shift);
- closure[r+nNod].push_back(shift);
- }
- break;
- }
- for (int r = 0; r < nNod ; r++){
- for (int j = 0; j < nNod; j++) {
- closure[r].push_back(shift + (r + j) % nNod);
- closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
- }
- }
- shift += nNod;
- int n = nNod*(corder-1);
- for (int r = 0; r < nNod ; r++){
- for (int j = 0; j < n; j++){
- closure[r].push_back(shift + (j + (corder - 1) * r) % n);
- closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
- }
- }
- shift += n;
- if (serendip) break;
- }
- for (int r = 0; r < nNod*2 ; r++) {
- closure[r].type = nodalBasis::getTag(TYPE_LIN, order);
- closureRef[r] = 0;
- }
-}
+ coefficients = generateLagrangeMonomialCoefficients(monomials, points);
-static void addEdgeNodes(polynomialBasis::clCont &closureFull, const int *edges, int order)
-{
- if (order < 2)
- return;
- int numNodes = 0;
- for (int i = 0; edges[i] >= 0; ++i) {
- numNodes = std::max(numNodes, edges[i] + 1);
- }
- std::vector<std::vector<int> > nodes2edges(numNodes, std::vector<int>(numNodes, -1));
- for (int i = 0; edges[i] >= 0; i += 2) {
- nodes2edges[edges[i]][edges[i + 1]] = i;
- nodes2edges[edges[i + 1]][edges[i]] = i + 1;
- }
- for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
- std::vector<int> &cl = closureFull[iClosure];
- for (int iEdge = 0; edges[iEdge] >= 0; iEdge += 2) {
- if (cl.empty())
- continue;
- int n0 = cl[edges[iEdge]];
- int n1 = cl[edges[iEdge + 1]];
- int oEdge = nodes2edges[n0][n1];
- if (oEdge == -1)
- Msg::Error("invalid p1 closure or invalid edges list");
- for (int i = 0 ; i < order - 1; i++)
- cl.push_back(numNodes + oEdge/2 * (order - 1) + ((oEdge % 2) ? order - 2 - i : i));
- }
- }
}
-static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
-{
- closure.clear();
- for (int iRotate = 0; iRotate < 3; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 4; iFace++){
- polynomialBasis::closure closure_face;
- getFaceClosureTet(iFace, iSign, iRotate, closure_face, order);
- closure.push_back(closure_face);
- }
- }
- }
-}
-static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull,
- std::vector<int> &closureRef,
- int order, bool serendip)
-{
- closureFull.clear();
- //input :
- static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
- static const int edges[] = {0, 1, 1, 2, 2, 0, 3, 0, 3, 2, 3, 1, -1};
- static const int faceOrientation[6] = {0,1,2,5,3,4};
- closureFull.resize(24);
- closureRef.resize(24);
- for (int i = 0; i < 24; i ++)
- closureRef[i] = 0;
- if (order == 0) {
- for (int i = 0; i < 24; i ++) {
- closureFull[i].push_back(0);
- }
- return;
- }
- //Mapping for the p1 nodes
- polynomialBasis::clCont closure;
- generateFaceClosureTet(closure, 1);
- for (unsigned int i = 0; i < closureFull.size(); i++) {
- std::vector<int> &clFull = closureFull[i];
- std::vector<int> &cl = closure[i];
- clFull.resize(4, -1);
- closureRef[i] = 0;
- for (unsigned int j = 0; j < cl.size(); j ++)
- clFull[closure[0][j]] = cl[j];
- for (int j = 0; j < 4; j ++)
- if (clFull[j] == -1)
- clFull[j] = (6 - clFull[(j + 1) % 4] - clFull[(j + 2) % 4] - clFull[(j + 3) % 4]);
- }
- int nodes2Faces[4][4][4];
- for (int i = 0; i < 4; i++) {
- for (int iRotate = 0; iRotate < 3; iRotate ++) {
- short int n0 = faces[i][(3 - iRotate) % 3];
- short int n1 = faces[i][(4 - iRotate) % 3];
- short int n2 = faces[i][(5 - iRotate) % 3];
- nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
- nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
- }
- }
- polynomialBasis::clCont closureTriangles;
- std::vector<int> closureTrianglesRef;
- if (order >= 3)
- generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
- addEdgeNodes(closureFull, edges, order);
- for (unsigned int iClosure = 0; iClosure < closureFull.size(); iClosure++) {
- //faces
- std::vector<int> &cl = closureFull[iClosure];
- if (order >= 3) {
- for (int iFace = 0; iFace < 4; iFace ++) {
- int n0 = cl[faces[iFace][0]];
- int n1 = cl[faces[iFace][1]];
- int n2 = cl[faces[iFace][2]];
- short int id = nodes2Faces[n0][n1][n2];
- short int iTriClosure = faceOrientation [id % 6];
- short int idFace = id/6;
- int nOnFace = closureTriangles[iTriClosure].size();
- for (int i = 0; i < nOnFace; i++) {
- cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace +
- closureTriangles[iTriClosure][i]);
- }
- }
- }
- }
- if (order >= 4 && !serendip) {
- polynomialBasis::clCont insideClosure;
- std::vector<int> fakeClosureRef;
- generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
- for (unsigned int i = 0; i < closureFull.size(); i ++) {
- unsigned int shift = closureFull[i].size();
- for (unsigned int j = 0; j < insideClosure[i].size(); j++)
- closureFull[i].push_back(insideClosure[i][j] + shift);
- }
- }
-}
-void rotateHex(int iFace, int iRot, int iSign, double uI, double vI,
- double &uO, double &vO, double &wO)
+polynomialBasis::~polynomialBasis()
{
- if (iSign<0){
- double tmp=uI;
- uI=vI;
- vI=tmp;
- }
- for (int i=0; i < iRot; i++){
- double tmp=uI;
- uI=-vI;
- vI=tmp;
- }
- switch (iFace) {
- case 0: uO = vI; vO = uI; wO=-1; break;
- case 1: uO = uI; vO = -1; wO=vI; break;
- case 2: uO =-1; vO = vI; wO=uI; break;
- case 3: uO = 1; vO = uI; wO=vI; break;
- case 4: uO =-uI; vO = 1; wO=vI; break;
- case 5: uO =uI; vO = vI; wO=1; break;
- }
}
-void rotateHexFull(int iFace, int iRot, int iSign, double uI, double vI, double wI, double &uO, double &vO, double &wO)
-{
- switch (iFace) {
- case 0: uO = uI; vO = vI; wO = wI; break;
- case 1: uO = wI; vO = uI; wO = vI; break;
- case 2: uO = vI; vO = wI; wO = uI; break;
- case 3: uO = wI; vO = vI; wO =-uI; break;
- case 4: uO = wI; vO =-uI; wO =-vI; break;
- case 5: uO = vI; vO = uI; wO =-wI; break;
- }
- for (int i = 0; i < iRot; i++){
- double tmp = uO;
- uO = -vO;
- vO = tmp;
- }
- if (iSign<0){
- double tmp = uO;
- uO = vO;
- vO = tmp;
- }
-}
-inline static double pow2(double x)
-{
- return x * x;
-}
-/*
-static void checkClosure(int tag) {
- printf("TYPE = %i\n", tag);
- const polynomialBasis &fs = *polynomialBases::find(tag);
- for (int i = 0; i < fs.closures.size(); ++i) {
- const std::vector<int> &clRef = fs.closures[fs.closureRef[i]];
- const std::vector<int> &cl = fs.closures[i];
- const std::vector<int> &clFull = fs.fullClosures[i];
- printf("i = %i\n", i);
- for (int j = 0; j < cl.size(); ++j) {
- printf("%3i ", clFull[clRef[j]]);
- }
- printf("\n");
- for (int j = 0; j < cl.size(); ++j) {
- printf("%3i ", cl[j]);
- }
- printf("\n");
- }
-}
-*/
+int polynomialBasis::getNumShapeFunctions() const { return coefficients.size1(); }
-static void generateFaceClosureHex(polynomialBasis::clCont &closure, int order,
- bool serendip, const fullMatrix<double> &points)
-{
- closure.clear();
- const polynomialBasis &fsFace = *polynomialBases::find
- (nodalBasis::getTag(TYPE_QUA, order, serendip));
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 6; iFace++) {
- polynomialBasis::closure cl;
- cl.type = fsFace.type;
- cl.resize(fsFace.points.size1());
- for (unsigned int iNode = 0; iNode < cl.size(); ++iNode) {
- double u,v,w;
- rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0),
- fsFace.points(iNode, 1), u, v, w);
- cl[iNode] = 0;
- double D = std::numeric_limits<double>::max();
- for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
- pow2(points(jNode, 2) - w);
- if (d < D) {
- cl[iNode] = jNode;
- D = d;
- }
- }
- }
- closure.push_back(cl);
- }
- }
- }
-}
-static void generateFaceClosureHexFull(polynomialBasis::clCont &closure,
- std::vector<int> &closureRef,
- int order, bool serendip,
- const fullMatrix<double> &points)
-{
- closure.clear();
- int clId = 0;
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 6; iFace++) {
- polynomialBasis::closure cl;
- cl.resize(points.size1());
- for (int iNode = 0; iNode < points.size1(); ++iNode) {
- double u,v,w;
- rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1),
- points(iNode, 2), u, v, w);
- int J = 0;
- double D = std::numeric_limits<double>::max();
- for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
- pow2(points(jNode, 2) - w);
- if (d < D) {
- J = jNode;
- D = d;
- }
- }
- cl[J] = iNode;
- }
- closure.push_back(cl);
- closureRef.push_back(0);
- clId ++;
- }
- }
- }
-}
-static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
- polynomialBasis::closure &closure, int order)
+void polynomialBasis::f(double u, double v, double w, double *sf) const
{
- if (order > 2)
- Msg::Error("FaceClosure not implemented for prisms of order %d",order);
- bool isTriangle = iFace<2;
- int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
- closure.clear();
- if (isTriangle && iRotate > 2) return;
- closure.resize(nNodes);
- if (order==0) {
- closure[0] = 0;
- return;
- }
- int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2},
- {1, 2, 5, 4}};
- int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15},
- {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
- // int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
- // {8, 13, 11, 7}, {9, 11, 14, 10}};
- int nVertex = isTriangle ? 3 : 4;
- closure.type = nodalBasis::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
- for (int i = 0; i < nVertex; ++i){
- int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
- closure[i] = order1node[iFace][k];
- }
- if (order==2) {
- for (int i = 0; i < nVertex; ++i){
- int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;
- //- iSign * iRotate
- closure[nVertex+i] = order2node[iFace][k];
+ double p[1256];
+ evaluateMonomials(u, v, w, p);
+ for (int i = 0; i < coefficients.size1(); i++) {
+ sf[i] = 0.0;
+ for (int j = 0; j < coefficients.size2(); j++) {
+ sf[i] += coefficients(i, j) * p[j];
}
- if (!isTriangle)
- closure[nNodes-1] = order2node[iFace][4]; // center
}
}
-static void generateFaceClosurePrism(polynomialBasis::clCont &closure, int order)
+
+
+void polynomialBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf) const
{
- closure.clear();
- for (int iRotate = 0; iRotate < 4; iRotate++){
- for (int iSign = 1; iSign >= -1; iSign -= 2){
- for (int iFace = 0; iFace < 5; iFace++){
- polynomialBasis::closure closure_face;
- getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
- closure.push_back(closure_face);
- }
- }
+ double p[1256];
+ sf.resize (coord.size1(), coefficients.size1());
+ for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
+ evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
+ for (int i = 0; i < coefficients.size1(); i++)
+ for (int j = 0; j < coefficients.size2(); j++)
+ sf(iPoint,i) += coefficients(i, j) * p[j];
}
}
-static void generateFaceClosurePrismFull(polynomialBasis::clCont &closureFull,
- std::vector<int> &closureRef, int order)
-{
- polynomialBasis::clCont closure;
- closureFull.clear();
- closureFull.resize(40);
- closureRef.resize(40);
- generateFaceClosurePrism(closure, 1);
- int ref3 = -1, ref4a = -1, ref4b = -1;
- for (unsigned int i = 0; i < closure.size(); i++) {
- std::vector<int> &clFull = closureFull[i];
- std::vector<int> &cl = closure[i];
- if (cl.size() == 0)
- continue;
- clFull.resize(6, -1);
- int &ref = cl.size() == 3 ? ref3 : (cl[0] / 3 + cl[1] / 3) % 2 ? ref4b : ref4a;
- if (ref == -1)
- ref = i;
- closureRef[i] = ref;
- for (unsigned int j = 0; j < cl.size(); j ++)
- clFull[closure[ref][j]] = cl[j];
- for (int j = 0; j < 6; j ++) {
- if (clFull[j] == -1) {
- int k = ((j / 3) + 1) % 2 * 3;
- int sum = (clFull[k + (j + 1) % 3] + clFull[k + (j + 2) % 3]);
- clFull[j] = ((sum / 6 + 1) % 2) * 3 + (12 - sum) % 3;
- }
- }
- }
- static const int edges[] = {0, 1, 0, 2, 0, 3, 1, 2, 1, 4, 2, 5,
- 3, 4, 3, 5, 4, 5, -1};
- addEdgeNodes(closureFull, edges, order);
- if ( order < 2 )
- return;
- // face center nodes for p2 prism
- static const int faces[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3},
- {0, 3, 5, 2}, {1, 2, 5, 4}};
-
- if ( order == 2 ) {
- int nextFaceNode = 15;
- int numFaces = 5;
- int numFaceNodes = 4;
- std::map<int,int> nodeSum2Face;
- for (int iFace = 0; iFace < numFaces ; iFace ++) {
- int nodeSum = 0;
- for (int iNode = 0; iNode < numFaceNodes; iNode++ ) {
- nodeSum += faces[iFace][iNode];
- }
- nodeSum2Face[nodeSum] = iFace;
- }
- for (unsigned int i = 0; i < closureFull.size(); i++ ) {
- if (closureFull[i].empty())
- continue;
- for (int iFace = 0; iFace < numFaces; iFace++ ) {
- int nodeSum = 0;
- for (int iNode = 0; iNode < numFaceNodes; iNode++)
- nodeSum += faces[iFace][iNode] < 0 ? faces[iFace][iNode] :
- closureFull[i][ faces[iFace][iNode] ];
- std::map<int,int>::iterator it = nodeSum2Face.find(nodeSum);
- if (it == nodeSum2Face.end() )
- Msg::Error("Prism face not found");
- int mappedFaceId = it->second;
- if ( mappedFaceId > 1) {
- closureFull[i].push_back(nextFaceNode + mappedFaceId - 2);
- }
- }
- }
- } else {
- Msg::Error("Prisms of order %d not implemented",order);
- }
-}
-static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
+void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
{
- closure.clear();
- closure.resize(2*nNod);
- for (int j = 0; j < nNod ; j++){
- closure[j].push_back(j);
- closure[j].push_back((j+1)%nNod);
- closure[nNod+j].push_back((j+1)%nNod);
- closure[nNod+j].push_back(j);
- for (int i=0; i < order-1; i++){
- closure[j].push_back( nNod + (order-1)*j + i );
- closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
+ double dfv[1256][3];
+ dfm.resize (coefficients.size1(), coord.size1() * 3, false);
+ int ii = 0;
+ int dimCoord = coord.size2();
+ for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
+ df(coord(iPoint,0), dimCoord > 1 ? coord(iPoint, 1) : 0., dimCoord > 2 ? coord(iPoint, 2) : 0., dfv);
+ for (int i = 0; i < coefficients.size1(); i++) {
+ dfm(i, iPoint * 3 + 0) = dfv[i][0];
+ dfm(i, iPoint * 3 + 1) = dfv[i][1];
+ dfm(i, iPoint * 3 + 2) = dfv[i][2];
+ ++ii;
}
- closure[j].type = closure[nNod+j].type = nodalBasis::getTag(TYPE_LIN, order);
}
}
-static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
-{
- closure.clear();
- closure.resize(nb);
- for (int i=0; i < nb; i++) {
- closure[i].push_back(0);
- closure[i].type = MSH_PNT;
- }
-}
-std::map<int, polynomialBasis> polynomialBases::fs;
-
-void polynomialBasis::initialize()
+void polynomialBasis::df(double u, double v, double w, double grads[][3]) const
{
- 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);
+ switch (monomials.size2()) {
+ case 1:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 0)
+ grads[i][0] += coefficients(i, j) *
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
}
- break;
- case TYPE_PRI:
- monomials = generatePascalPrism(order);
- if (order == 0) {
- generateClosureOrder0(closures,48);
- generateClosureOrder0(fullClosures,48);
- closureRef.resize(48, 0);
+ }
+ break;
+ case 2:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 0)
+ grads[i][0] += coefficients(i, j) *
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1));
+ if (monomials(j, 1) > 0)
+ grads[i][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
}
- else {
- generateFaceClosurePrism(closures, order);
- generateFaceClosurePrismFull(fullClosures, closureRef, order);
+ }
+ break;
+ case 3:
+ for (int i = 0; i < coefficients.size1(); i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ grads[i][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 0)
+ grads[i][0] += coefficients(i, j) *
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2));
+ if (monomials(j, 1) > 0)
+ grads[i][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2));
+ if (monomials(j, 2) > 0)
+ grads[i][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
}
- break;
- case TYPE_HEX:
- monomials = generatePascalHex(order, serendip);
- generateFaceClosureHex(closures, order, serendip, points);
- generateFaceClosureHexFull(fullClosures, closureRef, order, serendip, points);
- break;
+ }
+ break;
}
-
- coefficients = generateLagrangeMonomialCoefficients(monomials, points);
- }
+}
-const polynomialBasis *polynomialBases::find(int tag)
+
+void polynomialBasis::ddf(double u, double v, double w, double hess[][3][3]) const
{
- std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
- if (it != fs.end()) return &it->second;
-
- polynomialBasis F;
- F.type = tag;
- switch (tag) {
- case MSH_PNT : F.parentType = TYPE_PNT; F.order = 0; F.serendip = false; break;
- case MSH_LIN_1 : F.parentType = TYPE_LIN; F.order = 0; F.serendip = false; break;
- case MSH_LIN_2 : F.parentType = TYPE_LIN; F.order = 1; F.serendip = false; break;
- case MSH_LIN_3 : F.parentType = TYPE_LIN; F.order = 2; F.serendip = false; break;
- case MSH_LIN_4 : F.parentType = TYPE_LIN; F.order = 3; F.serendip = false; break;
- case MSH_LIN_5 : F.parentType = TYPE_LIN; F.order = 4; F.serendip = false; break;
- case MSH_LIN_6 : F.parentType = TYPE_LIN; F.order = 5; F.serendip = false; break;
- case MSH_LIN_7 : F.parentType = TYPE_LIN; F.order = 6; F.serendip = false; break;
- case MSH_LIN_8 : F.parentType = TYPE_LIN; F.order = 7; F.serendip = false; break;
- case MSH_LIN_9 : F.parentType = TYPE_LIN; F.order = 8; F.serendip = false; break;
- case MSH_LIN_10 : F.parentType = TYPE_LIN; F.order = 9; F.serendip = false; break;
- case MSH_LIN_11 : F.parentType = TYPE_LIN; F.order = 10;F.serendip = false; break;
- case MSH_TRI_1 : F.parentType = TYPE_TRI; F.order = 0; F.serendip = false; break;
- case MSH_TRI_3 : F.parentType = TYPE_TRI; F.order = 1; F.serendip = false; break;
- case MSH_TRI_6 : F.parentType = TYPE_TRI; F.order = 2; F.serendip = false; break;
- case MSH_TRI_10 : F.parentType = TYPE_TRI; F.order = 3; F.serendip = false; break;
- case MSH_TRI_15 : F.parentType = TYPE_TRI; F.order = 4; F.serendip = false; break;
- case MSH_TRI_21 : F.parentType = TYPE_TRI; F.order = 5; F.serendip = false; break;
- case MSH_TRI_28 : F.parentType = TYPE_TRI; F.order = 6; F.serendip = false; break;
- case MSH_TRI_36 : F.parentType = TYPE_TRI; F.order = 7; F.serendip = false; break;
- case MSH_TRI_45 : F.parentType = TYPE_TRI; F.order = 8; F.serendip = false; break;
- case MSH_TRI_55 : F.parentType = TYPE_TRI; F.order = 9; F.serendip = false; break;
- case MSH_TRI_66 : F.parentType = TYPE_TRI; F.order = 10;F.serendip = false; break;
- case MSH_TRI_9 : F.parentType = TYPE_TRI; F.order = 3; F.serendip = true; break;
- case MSH_TRI_12 : F.parentType = TYPE_TRI; F.order = 4; F.serendip = true; break;
- case MSH_TRI_15I : F.parentType = TYPE_TRI; F.order = 5; F.serendip = true; break;
- case MSH_TRI_18 : F.parentType = TYPE_TRI; F.order = 6; F.serendip = true; break;
- case MSH_TRI_21I : F.parentType = TYPE_TRI; F.order = 7; F.serendip = true; break;
- case MSH_TRI_24 : F.parentType = TYPE_TRI; F.order = 8; F.serendip = true; break;
- case MSH_TRI_27 : F.parentType = TYPE_TRI; F.order = 9; F.serendip = true; break;
- case MSH_TRI_30 : F.parentType = TYPE_TRI; F.order = 10;F.serendip = true; break;
- case MSH_TET_1 : F.parentType = TYPE_TET; F.order = 0; F.serendip = false; break;
- case MSH_TET_4 : F.parentType = TYPE_TET; F.order = 1; F.serendip = false; break;
- case MSH_TET_10 : F.parentType = TYPE_TET; F.order = 2; F.serendip = false; break;
- case MSH_TET_20 : F.parentType = TYPE_TET; F.order = 3; F.serendip = false; break;
- case MSH_TET_35 : F.parentType = TYPE_TET; F.order = 4; F.serendip = false; break;
- case MSH_TET_56 : F.parentType = TYPE_TET; F.order = 5; F.serendip = false; break;
- case MSH_TET_84 : F.parentType = TYPE_TET; F.order = 6; F.serendip = false; break;
- case MSH_TET_120 : F.parentType = TYPE_TET; F.order = 7; F.serendip = false; break;
- case MSH_TET_165 : F.parentType = TYPE_TET; F.order = 8; F.serendip = false; break;
- case MSH_TET_220 : F.parentType = TYPE_TET; F.order = 9; F.serendip = false; break;
- case MSH_TET_286 : F.parentType = TYPE_TET; F.order = 10;F.serendip = false; break;
- case MSH_TET_34 : F.parentType = TYPE_TET; F.order = 4; F.serendip = true; break;
- case MSH_TET_52 : F.parentType = TYPE_TET; F.order = 5; F.serendip = true; break;
- case MSH_TET_74 : F.parentType = TYPE_TET; F.order = 6; F.serendip = true; break;
- case MSH_TET_100 : F.parentType = TYPE_TET; F.order = 7; F.serendip = true; break;
- case MSH_TET_130 : F.parentType = TYPE_TET; F.order = 8; F.serendip = true; break;
- case MSH_TET_164 : F.parentType = TYPE_TET; F.order = 9; F.serendip = true; break;
- case MSH_TET_202 : F.parentType = TYPE_TET; F.order = 10;F.serendip = true; break;
- case MSH_QUA_1 : F.parentType = TYPE_QUA; F.order = 0; F.serendip = false; break;
- case MSH_QUA_4 : F.parentType = TYPE_QUA; F.order = 1; F.serendip = false; break;
- case MSH_QUA_9 : F.parentType = TYPE_QUA; F.order = 2; F.serendip = false; break;
- case MSH_QUA_16 : F.parentType = TYPE_QUA; F.order = 3; F.serendip = false; break;
- case MSH_QUA_25 : F.parentType = TYPE_QUA; F.order = 4; F.serendip = false; break;
- case MSH_QUA_36 : F.parentType = TYPE_QUA; F.order = 5; F.serendip = false; break;
- case MSH_QUA_49 : F.parentType = TYPE_QUA; F.order = 6; F.serendip = false; break;
- case MSH_QUA_64 : F.parentType = TYPE_QUA; F.order = 7; F.serendip = false; break;
- case MSH_QUA_81 : F.parentType = TYPE_QUA; F.order = 8; F.serendip = false; break;
- case MSH_QUA_100 : F.parentType = TYPE_QUA; F.order = 9; F.serendip = false; break;
- case MSH_QUA_121 : F.parentType = TYPE_QUA; F.order = 10;F.serendip = false; break;
- case MSH_QUA_8 : F.parentType = TYPE_QUA; F.order = 2; F.serendip = true; break;
- case MSH_QUA_12 : F.parentType = TYPE_QUA; F.order = 3; F.serendip = true; break;
- case MSH_QUA_16I : F.parentType = TYPE_QUA; F.order = 4; F.serendip = true; break;
- case MSH_QUA_20 : F.parentType = TYPE_QUA; F.order = 5; F.serendip = true; break;
- case MSH_QUA_24 : F.parentType = TYPE_QUA; F.order = 6; F.serendip = true; break;
- case MSH_QUA_28 : F.parentType = TYPE_QUA; F.order = 7; F.serendip = true; break;
- case MSH_QUA_32 : F.parentType = TYPE_QUA; F.order = 8; F.serendip = true; break;
- case MSH_QUA_36I : F.parentType = TYPE_QUA; F.order = 9; F.serendip = true; break;
- case MSH_QUA_40 : F.parentType = TYPE_QUA; F.order = 10;F.serendip = true; break;
- case MSH_PRI_1 : F.parentType = TYPE_PRI; F.order = 0; F.serendip = false; break;
- case MSH_PRI_6 : F.parentType = TYPE_PRI; F.order = 1; F.serendip = false; break;
- case MSH_PRI_18 : F.parentType = TYPE_PRI; F.order = 2; F.serendip = false; break;
- case MSH_HEX_8 : F.parentType = TYPE_HEX; F.order = 1; F.serendip = false; break;
- case MSH_HEX_27 : F.parentType = TYPE_HEX; F.order = 2; F.serendip = false; break;
- case MSH_HEX_64 : F.parentType = TYPE_HEX; F.order = 3; F.serendip = false; break;
- case MSH_HEX_125 : F.parentType = TYPE_HEX; F.order = 4; F.serendip = false; break;
- case MSH_HEX_216 : F.parentType = TYPE_HEX; F.order = 5; F.serendip = false; break;
- case MSH_HEX_343 : F.parentType = TYPE_HEX; F.order = 6; F.serendip = false; break;
- case MSH_HEX_512 : F.parentType = TYPE_HEX; F.order = 7; F.serendip = false; break;
- case MSH_HEX_729 : F.parentType = TYPE_HEX; F.order = 8; F.serendip = false; break;
- case MSH_HEX_1000: F.parentType = TYPE_HEX; F.order = 9; F.serendip = false; break;
- case MSH_HEX_20 : F.parentType = TYPE_HEX; F.order = 2; F.serendip = false; break;
- case MSH_HEX_56 : F.parentType = TYPE_HEX; F.order = 3; F.serendip = true; break;
- case MSH_HEX_98 : F.parentType = TYPE_HEX; F.order = 4; F.serendip = true; break;
- case MSH_HEX_152 : F.parentType = TYPE_HEX; F.order = 5; F.serendip = true; break;
- case MSH_HEX_222 : F.parentType = TYPE_HEX; F.order = 6; F.serendip = true; break;
- case MSH_HEX_296 : F.parentType = TYPE_HEX; F.order = 7; F.serendip = true; break;
- case MSH_HEX_386 : F.parentType = TYPE_HEX; F.order = 8; F.serendip = true; break;
- case MSH_HEX_488 : F.parentType = TYPE_HEX; F.order = 9; F.serendip = true; break;
- default :
- Msg::Error("Unknown function space %d: reverting to TET_4", tag);
- F.parentType = TYPE_TET; F.order = 1; F.serendip = false; break;
- }
+ switch (monomials.size2()) {
+ case 1:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
- switch (F.parentType) {
- case TYPE_PNT :
- F.numFaces = 1;
- F.dimension = 0;
- F.monomials = generate1DMonomials(0);
- F.points = generate1DPoints(0);
- break;
- case TYPE_LIN :
- F.numFaces = 2;
- F.dimension = 1;
- F.monomials = generate1DMonomials(F.order);
- F.points = generate1DPoints(F.order);
- generate1dVertexClosure(F.closures, F.order);
- generate1dVertexClosureFull(F.fullClosures, F.closureRef, F.order);
- break;
- case TYPE_TRI :
- F.numFaces = 3;
- F.dimension = 2;
- F.monomials = F.serendip ? generatePascalSerendipityTriangle(F.order) :
- generatePascalTriangle(F.order);
- F.points = gmshGeneratePointsTriangle(F.order, F.serendip);
- if (F.order == 0) {
- generateClosureOrder0(F.closures, 6);
- generateClosureOrder0(F.fullClosures, 6);
- F.closureRef.resize(6, 0);
- }
- else {
- generate2dEdgeClosure(F.closures, F.order);
- generate2dEdgeClosureFull(F.fullClosures, F.closureRef, F.order, 3, F.serendip);
- }
- break;
- case TYPE_QUA :
- F.numFaces = 4;
- F.dimension = 2;
- F.monomials = F.serendip ? generatePascalQuadSerendip(F.order) :
- generatePascalQuad(F.order);
- F.points = gmshGeneratePointsQuad(F.order, F.serendip);
- if (F.order == 0) {
- generateClosureOrder0(F.closures, 8);
- generateClosureOrder0(F.fullClosures, 8);
- F.closureRef.resize(8, 0);
- }
- else {
- generate2dEdgeClosure(F.closures, F.order, 4);
- generate2dEdgeClosureFull(F.fullClosures, F.closureRef, F.order, 4, F.serendip);
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 1) // second derivative !=0
+ hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
+ monomials(j, 0) * (monomials(j, 0) - 1);
+ }
}
break;
- case TYPE_TET :
- F.numFaces = 4;
- F.dimension = 3;
- F.monomials = F.serendip ? generatePascalSerendipityTetrahedron(F.order) :
- generatePascalTetrahedron(F.order);
- F.points = gmshGeneratePointsTetrahedron(F.order, F.serendip);
- if (F.order == 0) {
- generateClosureOrder0(F.closures,24);
- generateClosureOrder0(F.fullClosures, 24);
- F.closureRef.resize(24, 0);
- }
- else {
- generateFaceClosureTet(F.closures, F.order);
- generateFaceClosureTetFull(F.fullClosures, F.closureRef, F.order, F.serendip);
+ case 2:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 1) // second derivative !=0
+ hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
+ monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
+ if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
+ hess[i][0][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
+ if (monomials(j, 1) > 1)
+ hess[i][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
+ }
+ hess[i][1][0] = hess[i][0][1];
}
break;
- case TYPE_PRI :
- F.numFaces = 5;
- F.dimension = 3;
- F.monomials = generatePascalPrism(F.order);
- F.points = gmshGeneratePointsPrism(F.order, F.serendip);
- if (F.order == 0) {
- generateClosureOrder0(F.closures,48);
- generateClosureOrder0(F.fullClosures,48);
- F.closureRef.resize(48, 0);
- }
- else {
- generateFaceClosurePrism(F.closures, F.order);
- generateFaceClosurePrismFull(F.fullClosures, F.closureRef, F.order);
+ case 3:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 1)
+ hess[i][0][0] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2));
+
+ if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
+ hess[i][0][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2));
+ if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
+ hess[i][0][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
+ if (monomials(j, 1) > 1)
+ hess[i][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
+ pow_int(w, (int)monomials(j, 2));
+ if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
+ hess[i][1][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
+ if (monomials(j, 2) > 1)
+ hess[i][2][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
+ }
+ hess[i][1][0] = hess[i][0][1];
+ hess[i][2][0] = hess[i][0][2];
+ hess[i][2][1] = hess[i][1][2];
}
break;
- case TYPE_HEX :
- F.numFaces = 6;
- F.dimension = 3;
- F.monomials = generatePascalHex(F.order, F.serendip);
- F.points = gmshGeneratePointsHex(F.order, F.serendip);
- generateFaceClosureHex(F.closures, F.order, F.serendip, F.points);
- generateFaceClosureHexFull(F.fullClosures, F.closureRef, F.order, F.serendip, F.points);
- break;
}
- F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
- fs.insert(std::make_pair(tag, F));
- return &fs[tag];
}
-std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
-
-const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
-{
- std::pair<int,int> key(tag1,tag2);
- std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it =
- injector.find(key);
- if (it != injector.end()) return it->second;
-
- const polynomialBasis& fs1 = *find(tag1);
- const polynomialBasis& fs2 = *find(tag2);
- fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
- double sf[1256];
+void polynomialBasis::dddf(double u, double v, double w, double third[][3][3][3]) const
+{
+ switch (monomials.size2()) {
+ case 1:
+ for (int i = 0; i < coefficients.size1(); i++){
+ for (int p=0; p<3; p++)
+ for (int q=0; q<3; q++)
+ for (int r=0; r<3; r++)
+ third[i][p][q][r] = 0.;
- for (int i = 0; i < fs1.points.size1(); i++) {
- fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
- for (int j = 0; j < fs2.points.size1(); j++) inj(i, j) = sf[j];
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 2) // third derivative !=0
+ third[i][0][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 3)) *
+ monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2);
+ }
+ }
+ break;
+ case 2:
+ for (int i = 0; i < coefficients.size1(); i++){
+ for (int p=0; p<3; p++)
+ for (int q=0; q<3; q++)
+ for (int r=0; r<3; r++)
+ third[i][p][q][r] = 0.;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 2) // second derivative !=0
+ third[i][0][0][0] += coefficients(i, j) *
+ pow_int(u, (int)(monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) *
+ pow_int(v, (int)monomials(j, 1));
+ if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
+ third[i][0][0][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
+ if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
+ third[i][0][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
+ if (monomials(j, 1) > 2)
+ third[i][1][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2);
+ }
+ third[i][0][1][0] = third[i][0][0][1];
+ third[i][1][0][0] = third[i][0][0][1];
+ third[i][1][0][1] = third[i][0][1][1];
+ third[i][1][1][0] = third[i][0][1][1];
+ }
+ break;
+ case 3:
+ for (int i = 0; i < coefficients.size1(); i++){
+ for (int p=0; p<3; p++)
+ for (int q=0; q<3; q++)
+ for (int r=0; r<3; r++)
+ third[i][p][q][r] = 0.;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if (monomials(j, 0) > 2) // second derivative !=0
+ third[i][0][0][0] += coefficients(i, j) *
+ pow_int(u, (int)(monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) *
+ pow_int(v, (int)monomials(j, 1))*
+ pow_int(w, (int)monomials(j, 2));
+
+ if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
+ third[i][0][0][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1)*
+ pow_int(w, (int)monomials(j, 2));
+
+ if ((monomials(j, 0) > 1) && (monomials(j, 2) > 0))
+ third[i][0][0][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2) - 1)* monomials(j, 2);
+
+ if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
+ third[i][0][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)*
+ pow_int(w, (int)monomials(j, 2));
+
+ if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)&& (monomials(j, 2) > 0))
+ third[i][0][1][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
+ pow_int(v, (int)monomials(j, 1) - 1) *monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2);
+
+ if ((monomials(j, 0) > 0) && (monomials(j, 2) > 1))
+ third[i][0][2][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
+ pow_int(v, (int)monomials(j, 1))*
+ pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
+ if ((monomials(j, 1) > 2))
+ third[i][1][1][1] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)*
+ pow_int(w, (int)monomials(j, 2));
+
+ if ((monomials(j, 1) > 1) && (monomials(j, 2) > 0))
+ third[i][1][1][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)*
+ pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2) ;
+
+ if ((monomials(j, 1) > 0) && (monomials(j, 2) > 1))
+ third[i][1][2][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)-1) *monomials(j, 1)*
+ pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ;
+
+ if ((monomials(j, 2) > 2))
+ third[i][2][2][2] += coefficients(i, j) *
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1))*
+ pow_int(w, (int)monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2);
+
+
+ }
+ third[i][0][1][0] = third[i][0][0][1];
+ third[i][1][0][0] = third[i][0][0][1];
+
+ third[i][2][0][0] = third[i][0][0][2];
+ third[i][0][2][0] = third[i][0][0][2];
+
+ third[i][1][0][1] = third[i][0][1][1];
+ third[i][1][1][0] = third[i][0][1][1];
+
+ third[i][0][2][1] = third[i][0][1][2];
+ third[i][1][0][2] = third[i][0][1][2];
+ third[i][1][2][0] = third[i][0][1][2];
+ third[i][2][1][0] = third[i][0][1][2];
+ third[i][2][0][1] = third[i][0][1][2];
+
+ third[i][2][2][0] = third[i][0][2][2];
+ third[i][2][0][2] = third[i][0][2][2];
+
+ third[i][1][2][1] = third[i][1][1][2];
+ third[i][2][1][1] = third[i][1][1][2];
+
+ third[i][2][2][1] = third[i][1][2][2];
+ third[i][2][1][2] = third[i][1][2][2];
+ }
+ break;
}
-
- injector.insert(std::make_pair(key, inj));
- return injector[key];
}
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 17078d4e3e9f800731861e1cc7b193f2ca694c1e..4444efbb899053ac4368494879d8fc25e20f0f17 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -74,7 +74,7 @@ fullMatrix<double> generate1DMonomials(int order);
class polynomialBasis : public nodalBasis
{
// integrationOrder, closureId => df/dXi
- mutable std::map<int,std::vector<fullMatrix<double> > > _dfAtFace;
+// 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
@@ -85,372 +85,33 @@ class polynomialBasis : public nodalBasis
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
- void initialize();
+ polynomialBasis(int tag);
+ ~polynomialBasis();
- 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
- inline int getClosureType(int id) const
- {
- return closures[id].type;
- }
- inline const std::vector<int> &getClosure(int id) const
- {
- return closures[id];
- }
- inline const std::vector<int> &getFullClosure(int id) const
- {
- return fullClosures[id];
- }
- inline int getClosureId(int iFace, int iSign=1, int iRot=0) const
- {
- return iFace + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
- }
- inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const
- {
- iFace = i % numFaces;
- i = (i - iFace)/numFaces;
- iSign = i % 2;
- iRot = (i - iSign) / 2;
- }
- 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));
- }
- }
- inline void f(double u, double v, double w, double *sf) const
- {
- double p[1256];
- evaluateMonomials(u, v, w, p);
- for (int i = 0; i < coefficients.size1(); i++) {
- sf[i] = 0.0;
- for (int j = 0; j < coefficients.size2(); j++) {
- sf[i] += coefficients(i, j) * p[j];
- }
- }
- }
- inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const
- {
- double p[1256];
- sf.resize (coord.size1(), coefficients.size1());
- for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
- evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
- for (int i = 0; i < coefficients.size1(); i++)
- for (int j = 0; j < coefficients.size2(); j++)
- sf(iPoint,i) += coefficients(i, j) * p[j];
- }
- }
- inline void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
- {
- double dfv[1256][3];
- dfm.resize (coefficients.size1(), coord.size1() * 3, false);
- int ii = 0;
- int dimCoord = coord.size2();
- for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
- df(coord(iPoint,0), dimCoord > 1 ? coord(iPoint, 1) : 0., dimCoord > 2 ? coord(iPoint, 2) : 0., dfv);
- for (int i = 0; i < coefficients.size1(); i++) {
- dfm(i, iPoint * 3 + 0) = dfv[i][0];
- dfm(i, iPoint * 3 + 1) = dfv[i][1];
- dfm(i, iPoint * 3 + 2) = dfv[i][2];
- ++ii;
- }
- }
- }
- inline void df(double u, double v, double w, double grads[][3]) const
- {
- switch (monomials.size2()) {
- case 1:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 0)
- grads[i][0] += coefficients(i, j) *
- pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
- }
- }
- break;
- case 2:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 0)
- grads[i][0] += coefficients(i, j) *
- pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
- pow_int(v, (int)monomials(j, 1));
- if (monomials(j, 1) > 0)
- grads[i][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
- }
- }
- break;
- case 3:
- for (int i = 0; i < coefficients.size1(); i++){
- grads[i][0] = 0;
- grads[i][1] = 0;
- grads[i][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 0)
- grads[i][0] += coefficients(i, j) *
- pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)monomials(j, 2));
- if (monomials(j, 1) > 0)
- grads[i][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
- pow_int(w, (int)monomials(j, 2));
- if (monomials(j, 2) > 0)
- grads[i][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
- }
- }
- break;
- }
- }
- inline void ddf(double u, double v, double w, double hess[][3][3]) const
- {
- switch (monomials.size2()) {
- case 1:
- for (int i = 0; i < coefficients.size1(); i++){
- hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
- hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
- hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
-
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 1) // second derivative !=0
- hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
- monomials(j, 0) * (monomials(j, 0) - 1);
- }
- }
- break;
- case 2:
- for (int i = 0; i < coefficients.size1(); i++){
- hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
- hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
- hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 1) // second derivative !=0
- hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
- monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
- if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
- hess[i][0][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
- pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
- if (monomials(j, 1) > 1)
- hess[i][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
- }
- hess[i][1][0] = hess[i][0][1];
- }
- break;
- case 3:
- for (int i = 0; i < coefficients.size1(); i++){
- hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
- hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
- hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 1)
- hess[i][0][0] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)monomials(j, 2));
-
- if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
- hess[i][0][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
- pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
- pow_int(w, (int)monomials(j, 2));
- if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
- hess[i][0][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
- if (monomials(j, 1) > 1)
- hess[i][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
- pow_int(w, (int)monomials(j, 2));
- if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
- hess[i][1][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
- pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
- if (monomials(j, 2) > 1)
- hess[i][2][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
- }
- hess[i][1][0] = hess[i][0][1];
- hess[i][2][0] = hess[i][0][2];
- hess[i][2][1] = hess[i][1][2];
- }
- break;
- }
- }
- inline void dddf(double u, double v, double w, double third[][3][3][3]) const
- {
- switch (monomials.size2()) {
- case 1:
- for (int i = 0; i < coefficients.size1(); i++){
- for (int p=0; p<3; p++)
- for (int q=0; q<3; q++)
- for (int r=0; r<3; r++)
- third[i][p][q][r] = 0.;
+ virtual void f(double u, double v, double w, double *sf) const;
+ virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const;
+ virtual void df(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;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 2) // third derivative !=0
- third[i][0][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 3)) *
- monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2);
- }
- }
- break;
- case 2:
- for (int i = 0; i < coefficients.size1(); i++){
- for (int p=0; p<3; p++)
- for (int q=0; q<3; q++)
- for (int r=0; r<3; r++)
- third[i][p][q][r] = 0.;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 2) // second derivative !=0
- third[i][0][0][0] += coefficients(i, j) *
- pow_int(u, (int)(monomials(j, 0) - 3))*monomials(j, 0) * (monomials(j, 0) - 1) *(monomials(j, 0) - 2) *
- pow_int(v, (int)monomials(j, 1));
- if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
- third[i][0][0][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
- pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
- if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
- third[i][0][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
- pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
- if (monomials(j, 1) > 2)
- third[i][1][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1) - 3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2);
- }
- third[i][0][1][0] = third[i][0][0][1];
- third[i][1][0][0] = third[i][0][0][1];
- third[i][1][0][1] = third[i][0][1][1];
- third[i][1][1][0] = third[i][0][1][1];
- }
- break;
- case 3:
- for (int i = 0; i < coefficients.size1(); i++){
- for (int p=0; p<3; p++)
- for (int q=0; q<3; q++)
- for (int r=0; r<3; r++)
- third[i][p][q][r] = 0.;
- for(int j = 0; j < coefficients.size2(); j++){
- if (monomials(j, 0) > 2) // second derivative !=0
- third[i][0][0][0] += coefficients(i, j) *
- pow_int(u, (int)(monomials(j, 0) - 3)) *monomials(j, 0) * (monomials(j, 0) - 1)*(monomials(j, 0) - 2) *
- pow_int(v, (int)monomials(j, 1))*
- pow_int(w, (int)monomials(j, 2));
+ virtual int getNumShapeFunctions() const;
- if ((monomials(j, 0) > 1) && (monomials(j, 1) > 0))
- third[i][0][0][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
- pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1)*
- pow_int(w, (int)monomials(j, 2));
+ inline void evaluateMonomials(double u, double v, double w, double p[]) const;
- if ((monomials(j, 0) > 1) && (monomials(j, 2) > 0))
- third[i][0][0][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0)*(monomials(j, 0) - 1) *
- pow_int(v, (int)monomials(j, 1)) *
- pow_int(w, (int)monomials(j, 2) - 1)* monomials(j, 2);
-
- if ((monomials(j, 0) > 0) && (monomials(j, 1) > 1))
- third[i][0][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
- pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1)*
- pow_int(w, (int)monomials(j, 2));
-
- if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0)&& (monomials(j, 2) > 0))
- third[i][0][1][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
- pow_int(v, (int)monomials(j, 1) - 1) *monomials(j, 1) *
- pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2);
-
- if ((monomials(j, 0) > 0) && (monomials(j, 2) > 1))
- third[i][0][2][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0) - 1) *monomials(j, 0)*
- pow_int(v, (int)monomials(j, 1))*
- pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
- if ((monomials(j, 1) > 2))
- third[i][1][1][1] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1)-3) * monomials(j, 1) * (monomials(j, 1) - 1)*(monomials(j, 1) - 2)*
- pow_int(w, (int)monomials(j, 2));
-
- if ((monomials(j, 1) > 1) && (monomials(j, 2) > 0))
- third[i][1][1][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1)-2) * monomials(j, 1) * (monomials(j, 1) - 1)*
- pow_int(w, (int)monomials(j, 2) - 1) *monomials(j, 2) ;
-
- if ((monomials(j, 1) > 0) && (monomials(j, 2) > 1))
- third[i][1][2][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1)-1) *monomials(j, 1)*
- pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1) ;
-
- if ((monomials(j, 2) > 2))
- third[i][2][2][2] += coefficients(i, j) *
- pow_int(u, (int)monomials(j, 0)) *
- pow_int(v, (int)monomials(j, 1))*
- pow_int(w, (int)monomials(j, 2) - 3) * monomials(j, 2) * (monomials(j, 2) - 1)*(monomials(j, 2) - 2);
-
-
- }
- third[i][0][1][0] = third[i][0][0][1];
- third[i][1][0][0] = third[i][0][0][1];
-
- third[i][2][0][0] = third[i][0][0][2];
- third[i][0][2][0] = third[i][0][0][2];
-
- third[i][1][0][1] = third[i][0][1][1];
- third[i][1][1][0] = third[i][0][1][1];
-
- third[i][0][2][1] = third[i][0][1][2];
- third[i][1][0][2] = third[i][0][1][2];
- third[i][1][2][0] = third[i][0][1][2];
- third[i][2][1][0] = third[i][0][1][2];
- third[i][2][0][1] = third[i][0][1][2];
+};
- third[i][2][2][0] = third[i][0][2][2];
- third[i][2][0][2] = third[i][0][2][2];
- third[i][1][2][1] = third[i][1][1][2];
- third[i][2][1][1] = third[i][1][1][2];
- third[i][2][2][1] = third[i][1][2][2];
- third[i][2][1][2] = third[i][1][2][2];
- }
- break;
- }
+inline void polynomialBasis::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));
}
-};
+}
+
-class polynomialBases
-{
- private:
- static std::map<int, polynomialBasis> fs;
- static std::map<std::pair<int, int>, fullMatrix<double> > injector;
- public :
- static const polynomialBasis *find(int);
- static const fullMatrix<double> &findInjector(int, int);
-};
#endif
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 5c0bbe25c0dd3d39386747873fed9ace04263ee8..24a382d29ad734a8399de1c2bd01ba0ba224c368 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -4,17 +4,14 @@
// bugs and problems to <gmsh@geuz.org>.
#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
-pyramidalBasis::pyramidalBasis()
-{
-}
-void pyramidalBasis::initialize()
+pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
{
- BergotBasis bb(order);
- bergot = bb;
+ bergot = new BergotBasis(order);
int n = order+1;
int num_points = n*(n+1)*(2*n+1)/6;
@@ -28,31 +25,98 @@ void pyramidalBasis::initialize()
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);
+ 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);
+
}
+
+
+pyramidalBasis::~pyramidalBasis()
+{
+ delete bergot;
+}
+
+
+
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);
+ const int N = bergot->size();
+
+ double f[N];
+ bergot->f(u, v, w, f);
+
+ for (int i = 0; i < N; i++) {
+ val[i] = 0.;
+ for (int j = 0; j < N; j++) val[i] += (*VDMinv)(i,j)*f[j];
+ }
+
+}
+
+
+
+inline void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
+{
+
+ const int N = bergot->size(), NPts = coord.size1();
+
+ sf.resize(NPts, N);
+ double f[N];
+ fullVector<double> fVect(f,N);
- 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];
+ for (int iPt=0; iPt<NPts; iPt++) {
+ bergot->f(coord(iPt,0), coord(iPt,0), coord(iPt,0), f);
+ for (int i = 0; i < N; i++) {
+ sf(iPt,i) = 0.;
+ for (int j = 0; j < N; j++) sf(iPt,i) += (*VDMinv)(i,j)*f[j];
}
}
+
}
+
+
inline void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
{
+ const int N = bergot->size();
+
+ double df[N][3];
+ bergot->df(u, v, w, df);
+
+ 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)*df[j][0];
+ grads[i][1] += (*VDMinv)(i,j)*df[j][1];
+ grads[i][2] += (*VDMinv)(i,j)*df[j][2];
+ }
+ }
+
+}
+
+
+
+inline void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
+{
+
+ const int N = bergot->size(), NPts = coord.size1();
+
+ double dfv[N][3];
+ dfm.resize (NPts, 3*N, false);
+
+ for (int iPt=0; iPt<NPts; iPt++) {
+ df(coord(iPt,0), coord(iPt,1), coord(iPt,2), dfv);
+ for (int i = 0; i < N; i++) {
+ dfm(i, 3*iPt+0) = dfv[i][0];
+ dfm(i, 3*iPt+1) = dfv[i][1];
+ dfm(i, 3*iPt+2) = dfv[i][2];
+ }
+ }
+
}
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index 2b5b52aa6a6ba49b0f594a08e274a66815bec1eb..929ff85f85739e1bee1d49f4538299cc81c20822 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -11,6 +11,8 @@
#include "nodalBasis.h"
#include "BergotBasis.h"
+
+
class pyramidalBasis: public nodalBasis
{
private:
@@ -18,16 +20,19 @@ class pyramidalBasis: public nodalBasis
fullMatrix<double>* VDMinv;
// Orthogonal basis for the pyramid
- BergotBasis bergot;
+ 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;
+ pyramidalBasis(int tag);
+ ~pyramidalBasis();
- void initialize();
+ virtual void f(double u, double v, double w, double *val) const;
+ virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf);
+ virtual void df(double u, double v, double w, double grads[][3]) const;
+ virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
};
+
+
#endif
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 35575f7369510447334ce69169796256d4aaacfe..335edf3da5c5f3b12f1e492742b1b4164f151e16 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -6,7 +6,7 @@
#include "PViewDataList.h"
#include "GmshMessage.h"
#include "GmshDefines.h"
-#include "polynomialBasis.h"
+#include "BasisFactory.h"
#include "Numeric.h"
#include "SmoothData.h"
#include "Context.h"
@@ -883,11 +883,11 @@ void PViewDataList::setOrder2(int type)
case TYPE_TET: typeMSH = MSH_TET_10; break;
case TYPE_HEX: typeMSH = MSH_HEX_27; break;
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
- case TYPE_PYR: typeMSH = MSH_PYR_14; break;
+// case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
- const polynomialBasis *fs = polynomialBases::find(typeMSH);
+ const polynomialBasis *fs = (polynomialBasis*)BasisFactory::create(typeMSH);
if(!fs){
- Msg::Error("Could not find function space for element type %d", typeMSH);
+ Msg::Error("Could not find polynomial function space for element type %d", typeMSH);
return;
}
setInterpolationMatrices(type, fs->coefficients, fs->monomials,
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 5a488939754d6b65e37fae30449eb95f621ca9bf..478fd0cace1de5f701bee190b07ffcaa56c202a0 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -620,7 +620,7 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
void DI_Element::getShapeFunctions(double u, double v, double w, double s[], int o) const
{
//printf("type elem =%d order =%d\n", type(), 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");
}
@@ -628,7 +628,7 @@ void DI_Element::getShapeFunctions(double u, double v, double w, double s[], int
void DI_Element::getGradShapeFunctions(double u, double v, double w, double s[][3],
int o) const
{
- 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");
}
@@ -642,7 +642,7 @@ void DI_Element::setPolynomialOrder (int o) {
if(o == 1) return;
- const polynomialBasis *fs = getFunctionSpace(o);
+ const nodalBasis *fs = getFunctionSpace(o);
if(!fs) Msg::Error("Function space not implemented for this type of element");
mid_ = new DI_Point[nbMid()];
@@ -668,7 +668,7 @@ void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vect
if(o == 1) return;
- const polynomialBasis *fs = getFunctionSpace(o);
+ const nodalBasis *fs = getFunctionSpace(o);
if(!fs) Msg::Error("Function space not implemented for this type of element");
mid_ = new DI_Point[nbMid()];
@@ -1042,10 +1042,10 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
}
// DI_Line methods --------------------------------------------------------------------------------
-const polynomialBasis* DI_Line::getFunctionSpace(int o) const{
+const nodalBasis* DI_Line::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_LIN, order);
- return polynomialBases::find(tag);
+ return BasisFactory::create(tag);
}
void DI_Line::computeIntegral() {
@@ -1062,11 +1062,11 @@ void DI_Line::computeIntegral() {
// }
}
// DI_Triangle methods ----------------------------------------------------------------------------
-const polynomialBasis* DI_Triangle::getFunctionSpace(int o) const
+const nodalBasis* DI_Triangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_TRI, order);
- return polynomialBases::find(tag);
+ return BasisFactory::create(tag);
}
void DI_Triangle::computeIntegral() {
integral_ = TriSurf(pt(0), pt(1), pt(2));
@@ -1087,10 +1087,10 @@ double DI_Triangle::quality() const {
}
// DI_Quad methods --------------------------------------------------------------------------------
-const polynomialBasis* DI_Quad::getFunctionSpace(int o) const{
+const nodalBasis* DI_Quad::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_QUA, order);
- return polynomialBases::find(tag);
+ return BasisFactory::create(tag);
}
void DI_Quad::computeIntegral() {
@@ -1111,10 +1111,10 @@ void DI_Quad::computeIntegral() {
}
// DI_Tetra methods -------------------------------------------------------------------------------
-const polynomialBasis* DI_Tetra::getFunctionSpace(int o) const{
+const nodalBasis* DI_Tetra::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_TET, order);
- return polynomialBases::find(tag);
+ return BasisFactory::create(tag);
}
void DI_Tetra::computeIntegral() {
@@ -1126,10 +1126,10 @@ double DI_Tetra::quality() const {
// Hexahedron methods -----------------------------------------------------------------------------
-const polynomialBasis* DI_Hexa::getFunctionSpace(int o) const{
+const nodalBasis* DI_Hexa::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_HEX, order);
- return polynomialBases::find(tag);
+ return BasisFactory::create(tag);
}
void DI_Hexa::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index fdc1b93f8a397733adb7a8ee4985633acdf25e4f..c791cfc1fba67c07781382aed74dd52d2e0e4ed4 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -9,7 +9,7 @@
#include <map>
#include <cmath>
#include "gmshLevelset.h"
-#include "polynomialBasis.h"
+#include "BasisFactory.h"
#include "GmshDefines.h"
// Element type
@@ -48,8 +48,8 @@ class DI_Point
DI_Point (double x, double y, double z, const DI_Element *e,
const std::vector<gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
virtual ~DI_Point(){}
- virtual const polynomialBasis* getFunctionSpace(int o) const
- { return polynomialBases::find(MSH_PNT); }
+ virtual const nodalBasis* getFunctionSpace(int o) const
+ { return BasisFactory::create(MSH_PNT); }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
s[0] = 1.;
@@ -245,7 +245,7 @@ class DI_Element
if(pts_) delete [] pts_;
if(mid_) delete [] mid_;
}
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
+ virtual const nodalBasis* getFunctionSpace(int order=-1) const { return 0; }
// return type
virtual int type() const = 0;
// return the dimension of the element
@@ -410,7 +410,7 @@ class DI_Line : public DI_Element
inline int nbVert() const {return 2;}
inline int nbMid() const {return polOrder_ - 1;}
inline int nbEdg() const {return 1;}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 2.;}
void getRefIntegrationPoints (const int polynomialOrder,
@@ -477,7 +477,7 @@ class DI_Triangle : public DI_Element
return 0.5 * (polOrder_ - 1) * (polOrder_ + 4);
}
inline int nbEdg() const {return 3;}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 0.5;}
void getRefIntegrationPoints (const int polynomialOrder,
@@ -561,7 +561,7 @@ class DI_Quad : public DI_Element
return (polOrder_ - 1) * (polOrder_ + 3);
}
inline int nbEdg() const {return 4;}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 4.;}
void getRefIntegrationPoints (const int polynomialOrder,
@@ -643,7 +643,7 @@ class DI_Tetra : public DI_Element
return (polOrder_ + 1) * (polOrder_ + 2) * (polOrder_ + 3) / 6 - 4;
}
inline int nbEdg() const {return 6;}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 1. / 6.;}
inline void getRefIntegrationPoints (const int polynomialOrder,
@@ -745,7 +745,7 @@ class DI_Hexa : public DI_Element
return (polOrder_ + 1) * (polOrder_ + 1) * (polOrder_ + 1) - 8;
}
inline int nbEdg() const {return 12;}
- virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
+ virtual const nodalBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 8.;}
void getRefIntegrationPoints (const int polynomialOrder,
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index e940e8b3951b390ea1331087a49b14e757a3cdad..a37d265549d5d42dd52a683a2c1775de6fae214c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -150,10 +150,6 @@ SVector3 Mesh::getNormalEl(int iEl)
return n;
break;
}
- case TYPE_TET: {
- return SVector3(0.);
- break;
- }
default:
std::cout << "ERROR: getNormalEl: Unknown element type" << std::endl;
break;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 610a41b858b5f6ffb0095e267cd824bde472627a..b713254b1e1bda8282e80ea24da9ee8e1f0d97c1 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -329,6 +329,9 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
fflush(stdout);
OptHOM temp(&entity, toOptimize[i].first, toOptimize[i].second, method);
+std::ostringstream ossI1;
+ossI1 << "initial_" << entity.tag() << "ITER_" << i << ".msh";
+temp.mesh.writeMSH(ossI1.str().c_str());
int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
OptHomMessage("Jacobian optimization succeed, starting svd optimization");