From 36ccac01c945b192d2b2ebf0d03d43c6eab50a06 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Thu, 6 Jun 2013 20:35:07 +0000
Subject: [PATCH] Structure Bases. All find/create functions are now in
BasisFactory
---
Geo/MHexahedron.cpp | 80 +-
Geo/MLine.cpp | 42 +-
Geo/MPoint.h | 4 +-
Geo/MPrism.cpp | 20 +-
Geo/MPyramid.cpp | 46 +-
Geo/MQuadrangle.cpp | 84 +-
Geo/MTetrahedron.cpp | 76 +-
Geo/MTriangle.cpp | 86 +-
Geo/discreteEdge.cpp | 2 +-
Mesh/HighOrder.cpp | 80 +-
Mesh/meshGFaceRecombine.h | 2 +-
Numeric/BasisFactory.cpp | 37 +-
Numeric/BasisFactory.h | 11 +-
Numeric/CMakeLists.txt | 11 +-
Numeric/JacobianBasis.cpp | 1087 +----------------
Numeric/JacobianBasis.h | 19 +-
Numeric/bezierBasis.cpp | 1069 ++++++++++++++++
Numeric/bezierBasis.h | 30 +
Numeric/jacobiPolynomials.h | 2 +-
Numeric/legendrePolynomials.h | 2 +-
Numeric/nodalBasis.cpp | 2 +-
Numeric/polynomialBasis.h | 10 +-
Plugin/AnalyseCurvedMesh.cpp | 2 +-
Post/PViewDataList.cpp | 2 +-
contrib/DiscreteIntegration/Integration3D.cpp | 10 +-
contrib/DiscreteIntegration/Integration3D.h | 2 +-
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 2 +-
27 files changed, 1426 insertions(+), 1394 deletions(-)
create mode 100644 Numeric/bezierBasis.cpp
create mode 100644 Numeric/bezierBasis.h
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 52d0a067f5..29cfc61496 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -198,31 +198,31 @@ const nodalBasis* MHexahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return BasisFactory::create(MSH_HEX_1);
- case 1: return BasisFactory::create(MSH_HEX_8);
- case 2: return BasisFactory::create(MSH_HEX_20);
- case 3: return BasisFactory::create(MSH_HEX_56);
- case 4: return BasisFactory::create(MSH_HEX_98);
- case 5: return BasisFactory::create(MSH_HEX_152);
- case 6: return BasisFactory::create(MSH_HEX_218);
- case 7: return BasisFactory::create(MSH_HEX_296);
- case 8: return BasisFactory::create(MSH_HEX_386);
- case 9: return BasisFactory::create(MSH_HEX_488);
+ case 0: return BasisFactory::getNodalBasis(MSH_HEX_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_HEX_8);
+ case 2: return BasisFactory::getNodalBasis(MSH_HEX_20);
+ case 3: return BasisFactory::getNodalBasis(MSH_HEX_56);
+ case 4: return BasisFactory::getNodalBasis(MSH_HEX_98);
+ case 5: return BasisFactory::getNodalBasis(MSH_HEX_152);
+ case 6: return BasisFactory::getNodalBasis(MSH_HEX_218);
+ case 7: return BasisFactory::getNodalBasis(MSH_HEX_296);
+ case 8: return BasisFactory::getNodalBasis(MSH_HEX_386);
+ case 9: return BasisFactory::getNodalBasis(MSH_HEX_488);
default: Msg::Error("Order %d hex function space not implemented", order); break;
}
}
else {
switch (order) {
- case 0: return BasisFactory::create(MSH_HEX_1);
- case 1: return BasisFactory::create(MSH_HEX_8);
- case 2: return BasisFactory::create(MSH_HEX_27);
- case 3: return BasisFactory::create(MSH_HEX_64);
- case 4: return BasisFactory::create(MSH_HEX_125);
- case 5: return BasisFactory::create(MSH_HEX_216);
- case 6: return BasisFactory::create(MSH_HEX_343);
- case 7: return BasisFactory::create(MSH_HEX_512);
- case 8: return BasisFactory::create(MSH_HEX_729);
- case 9: return BasisFactory::create(MSH_HEX_1000);
+ case 0: return BasisFactory::getNodalBasis(MSH_HEX_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_HEX_8);
+ case 2: return BasisFactory::getNodalBasis(MSH_HEX_27);
+ case 3: return BasisFactory::getNodalBasis(MSH_HEX_64);
+ case 4: return BasisFactory::getNodalBasis(MSH_HEX_125);
+ case 5: return BasisFactory::getNodalBasis(MSH_HEX_216);
+ case 6: return BasisFactory::getNodalBasis(MSH_HEX_343);
+ case 7: return BasisFactory::getNodalBasis(MSH_HEX_512);
+ case 8: return BasisFactory::getNodalBasis(MSH_HEX_729);
+ case 9: return BasisFactory::getNodalBasis(MSH_HEX_1000);
default: Msg::Error("Order %d hex function space not implemented", order); break;
}
}
@@ -237,31 +237,31 @@ const JacobianBasis* MHexahedron::getJacobianFuncSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return JacobianBasis::find(MSH_HEX_1);
- 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_218);
- case 7: return JacobianBasis::find(MSH_HEX_296);
- case 8: return JacobianBasis::find(MSH_HEX_386);
- case 9: return JacobianBasis::find(MSH_HEX_488);
+ case 0: return BasisFactory::getJacobianBasis(MSH_HEX_1);
+ case 1: return BasisFactory::getJacobianBasis(MSH_HEX_8);
+ case 2: return BasisFactory::getJacobianBasis(MSH_HEX_20);
+ case 3: return BasisFactory::getJacobianBasis(MSH_HEX_56);
+ case 4: return BasisFactory::getJacobianBasis(MSH_HEX_98);
+ case 5: return BasisFactory::getJacobianBasis(MSH_HEX_152);
+ case 6: return BasisFactory::getJacobianBasis(MSH_HEX_218);
+ case 7: return BasisFactory::getJacobianBasis(MSH_HEX_296);
+ case 8: return BasisFactory::getJacobianBasis(MSH_HEX_386);
+ case 9: return BasisFactory::getJacobianBasis(MSH_HEX_488);
default: Msg::Error("Order %d hex incomplete Jacobian function space not implemented", order); break;
}
}
else {
switch (order) {
- case 0: return JacobianBasis::find(MSH_HEX_1);
- 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);
+ case 0: return BasisFactory::getJacobianBasis(MSH_HEX_1);
+ case 1: return BasisFactory::getJacobianBasis(MSH_HEX_8);
+ case 2: return BasisFactory::getJacobianBasis(MSH_HEX_27);
+ case 3: return BasisFactory::getJacobianBasis(MSH_HEX_64);
+ case 4: return BasisFactory::getJacobianBasis(MSH_HEX_125);
+ case 5: return BasisFactory::getJacobianBasis(MSH_HEX_216);
+ case 6: return BasisFactory::getJacobianBasis(MSH_HEX_343);
+ case 7: return BasisFactory::getJacobianBasis(MSH_HEX_512);
+ case 8: return BasisFactory::getJacobianBasis(MSH_HEX_729);
+ case 9: return BasisFactory::getJacobianBasis(MSH_HEX_1000);
default: Msg::Error("Order %d hex Jacobian function space not implemented", order); break;
}
}
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 91d9404595..7b33b69a68 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -16,17 +16,17 @@ const nodalBasis* MLine::getFunctionSpace(int o) const
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 0: return BasisFactory::create(MSH_LIN_1);
- case 1: return BasisFactory::create(MSH_LIN_2);
- case 2: return BasisFactory::create(MSH_LIN_3);
- case 3: return BasisFactory::create(MSH_LIN_4);
- case 4: return BasisFactory::create(MSH_LIN_5);
- case 5: return BasisFactory::create(MSH_LIN_6);
- case 6: return BasisFactory::create(MSH_LIN_7);
- case 7: return BasisFactory::create(MSH_LIN_8);
- case 8: return BasisFactory::create(MSH_LIN_9);
- case 9: return BasisFactory::create(MSH_LIN_10);
- case 10: return BasisFactory::create(MSH_LIN_11);
+ case 0: return BasisFactory::getNodalBasis(MSH_LIN_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_LIN_2);
+ case 2: return BasisFactory::getNodalBasis(MSH_LIN_3);
+ case 3: return BasisFactory::getNodalBasis(MSH_LIN_4);
+ case 4: return BasisFactory::getNodalBasis(MSH_LIN_5);
+ case 5: return BasisFactory::getNodalBasis(MSH_LIN_6);
+ case 6: return BasisFactory::getNodalBasis(MSH_LIN_7);
+ case 7: return BasisFactory::getNodalBasis(MSH_LIN_8);
+ case 8: return BasisFactory::getNodalBasis(MSH_LIN_9);
+ case 9: return BasisFactory::getNodalBasis(MSH_LIN_10);
+ case 10: return BasisFactory::getNodalBasis(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
@@ -37,16 +37,16 @@ const JacobianBasis* MLine::getJacobianFuncSpace(int o) const
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1: return JacobianBasis::find(MSH_LIN_2);
- case 2: return JacobianBasis::find(MSH_LIN_3);
- case 3: return JacobianBasis::find(MSH_LIN_4);
- case 4: return JacobianBasis::find(MSH_LIN_5);
- case 5: return JacobianBasis::find(MSH_LIN_6);
- case 6: return JacobianBasis::find(MSH_LIN_7);
- case 7: return JacobianBasis::find(MSH_LIN_8);
- case 8: return JacobianBasis::find(MSH_LIN_9);
- case 9: return JacobianBasis::find(MSH_LIN_10);
- case 10: return JacobianBasis::find(MSH_LIN_11);
+ case 1: return BasisFactory::getJacobianBasis(MSH_LIN_2);
+ case 2: return BasisFactory::getJacobianBasis(MSH_LIN_3);
+ case 3: return BasisFactory::getJacobianBasis(MSH_LIN_4);
+ case 4: return BasisFactory::getJacobianBasis(MSH_LIN_5);
+ case 5: return BasisFactory::getJacobianBasis(MSH_LIN_6);
+ case 6: return BasisFactory::getJacobianBasis(MSH_LIN_7);
+ case 7: return BasisFactory::getJacobianBasis(MSH_LIN_8);
+ case 8: return BasisFactory::getJacobianBasis(MSH_LIN_9);
+ case 9: return BasisFactory::getJacobianBasis(MSH_LIN_10);
+ case 10: return BasisFactory::getJacobianBasis(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 66ae8582c1..08ccc7fbc1 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -63,11 +63,11 @@ class MPoint : public MElement {
}
virtual const nodalBasis* getFunctionSpace(int o) const
{
- return BasisFactory::create(MSH_PNT);
+ return BasisFactory::getNodalBasis(MSH_PNT);
}
virtual const JacobianBasis* getJacobianFuncSpace(int o) const
{
- return JacobianBasis::find(MSH_PNT);
+ return BasisFactory::getJacobianBasis(MSH_PNT);
}
virtual bool isInside(double u, double v, double w) const
{
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index f364c55674..508bf117d7 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -40,17 +40,17 @@ const nodalBasis* MPrism::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return BasisFactory::create(MSH_PRI_1);
- case 1: return BasisFactory::create(MSH_PRI_6);
- case 2: return BasisFactory::create(MSH_PRI_18);
+ case 0: return BasisFactory::getNodalBasis(MSH_PRI_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_PRI_6);
+ case 2: return BasisFactory::getNodalBasis(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return BasisFactory::create(MSH_PRI_1);
- case 1: return BasisFactory::create(MSH_PRI_6);
- case 2: return BasisFactory::create(MSH_PRI_18);
+ case 0: return BasisFactory::getNodalBasis(MSH_PRI_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_PRI_6);
+ case 2: return BasisFactory::getNodalBasis(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
@@ -65,15 +65,15 @@ const JacobianBasis* MPrism::getJacobianFuncSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBasis::find(MSH_PRI_6);
- case 2: return JacobianBasis::find(MSH_PRI_18);
+ case 1: return BasisFactory::getJacobianBasis(MSH_PRI_6);
+ case 2: return BasisFactory::getJacobianBasis(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBasis::find(MSH_PRI_6);
- case 2: return JacobianBasis::find(MSH_PRI_18);
+ case 1: return BasisFactory::getJacobianBasis(MSH_PRI_6);
+ case 2: return BasisFactory::getJacobianBasis(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index f6e3c41307..1a44f7c4b6 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -35,31 +35,31 @@ const nodalBasis* MPyramid::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return BasisFactory::create(MSH_PYR_5);
- case 2: return BasisFactory::create(MSH_PYR_14);
- case 0: return BasisFactory::create(MSH_PYR_1);
- case 3: return BasisFactory::create(MSH_PYR_29);
- case 4: return BasisFactory::create(MSH_PYR_50);
- case 5: return BasisFactory::create(MSH_PYR_77);
- case 6: return BasisFactory::create(MSH_PYR_110);
- case 7: return BasisFactory::create(MSH_PYR_149);
- case 8: return BasisFactory::create(MSH_PYR_194);
- case 9: return BasisFactory::create(MSH_PYR_245);
+ case 1: return BasisFactory::getNodalBasis(MSH_PYR_5);
+ case 2: return BasisFactory::getNodalBasis(MSH_PYR_14);
+ case 0: return BasisFactory::getNodalBasis(MSH_PYR_1);
+ case 3: return BasisFactory::getNodalBasis(MSH_PYR_29);
+ case 4: return BasisFactory::getNodalBasis(MSH_PYR_50);
+ case 5: return BasisFactory::getNodalBasis(MSH_PYR_77);
+ case 6: return BasisFactory::getNodalBasis(MSH_PYR_110);
+ case 7: return BasisFactory::getNodalBasis(MSH_PYR_149);
+ case 8: return BasisFactory::getNodalBasis(MSH_PYR_194);
+ case 9: return BasisFactory::getNodalBasis(MSH_PYR_245);
default: Msg::Error("Order %d pyramid function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return BasisFactory::create(MSH_PYR_1);
- case 1: return BasisFactory::create(MSH_PYR_5);
- case 2: return BasisFactory::create(MSH_PYR_14);
- case 3: return BasisFactory::create(MSH_PYR_30);
- case 4: return BasisFactory::create(MSH_PYR_55);
- case 5: return BasisFactory::create(MSH_PYR_91);
- case 6: return BasisFactory::create(MSH_PYR_140);
- case 7: return BasisFactory::create(MSH_PYR_204);
- case 8: return BasisFactory::create(MSH_PYR_285);
- case 9: return BasisFactory::create(MSH_PYR_385);
+ case 0: return BasisFactory::getNodalBasis(MSH_PYR_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_PYR_5);
+ case 2: return BasisFactory::getNodalBasis(MSH_PYR_14);
+ case 3: return BasisFactory::getNodalBasis(MSH_PYR_30);
+ case 4: return BasisFactory::getNodalBasis(MSH_PYR_55);
+ case 5: return BasisFactory::getNodalBasis(MSH_PYR_91);
+ case 6: return BasisFactory::getNodalBasis(MSH_PYR_140);
+ case 7: return BasisFactory::getNodalBasis(MSH_PYR_204);
+ case 8: return BasisFactory::getNodalBasis(MSH_PYR_285);
+ case 9: return BasisFactory::getNodalBasis(MSH_PYR_385);
default: Msg::Error("Order %d pyramid function space not implemented", order);
}
}
@@ -71,9 +71,9 @@ 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);
+ case 1: return BasisFactory::getJacobianBasis(MSH_PYR_5);
+ case 2: return BasisFactory::getJacobianBasis(MSH_PYR_14);
+ case 3: return BasisFactory::getJacobianBasis(MSH_PYR_30);
default: Msg::Error("Order %d pyramid function space not implemented", order); break;
}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 6808e69e23..0f6b54877c 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -25,31 +25,31 @@ const nodalBasis* MQuadrangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 0: return BasisFactory::create(MSH_QUA_1);
- case 1: return BasisFactory::create(MSH_QUA_4);
- case 2: return BasisFactory::create(MSH_QUA_8);
- case 3: return BasisFactory::create(MSH_QUA_12);
- case 4: return BasisFactory::create(MSH_QUA_16I);
- case 5: return BasisFactory::create(MSH_QUA_20);
- case 6: return BasisFactory::create(MSH_QUA_24);
- case 7: return BasisFactory::create(MSH_QUA_28);
- case 8: return BasisFactory::create(MSH_QUA_32);
- case 9: return BasisFactory::create(MSH_QUA_36I);
- case 10: return BasisFactory::create(MSH_QUA_40);
+ case 0: return BasisFactory::getNodalBasis(MSH_QUA_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_QUA_4);
+ case 2: return BasisFactory::getNodalBasis(MSH_QUA_8);
+ case 3: return BasisFactory::getNodalBasis(MSH_QUA_12);
+ case 4: return BasisFactory::getNodalBasis(MSH_QUA_16I);
+ case 5: return BasisFactory::getNodalBasis(MSH_QUA_20);
+ case 6: return BasisFactory::getNodalBasis(MSH_QUA_24);
+ case 7: return BasisFactory::getNodalBasis(MSH_QUA_28);
+ case 8: return BasisFactory::getNodalBasis(MSH_QUA_32);
+ case 9: return BasisFactory::getNodalBasis(MSH_QUA_36I);
+ case 10: return BasisFactory::getNodalBasis(MSH_QUA_40);
}
}
switch (order) {
- case 0: return BasisFactory::create(MSH_QUA_1);
- case 1: return BasisFactory::create(MSH_QUA_4);
- case 2: return BasisFactory::create(MSH_QUA_9);
- case 3: return BasisFactory::create(MSH_QUA_16);
- case 4: return BasisFactory::create(MSH_QUA_25);
- case 5: return BasisFactory::create(MSH_QUA_36);
- case 6: return BasisFactory::create(MSH_QUA_49);
- case 7: return BasisFactory::create(MSH_QUA_64);
- case 8: return BasisFactory::create(MSH_QUA_81);
- case 9: return BasisFactory::create(MSH_QUA_100);
- case 10: return BasisFactory::create(MSH_QUA_121);
+ case 0: return BasisFactory::getNodalBasis(MSH_QUA_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_QUA_4);
+ case 2: return BasisFactory::getNodalBasis(MSH_QUA_9);
+ case 3: return BasisFactory::getNodalBasis(MSH_QUA_16);
+ case 4: return BasisFactory::getNodalBasis(MSH_QUA_25);
+ case 5: return BasisFactory::getNodalBasis(MSH_QUA_36);
+ case 6: return BasisFactory::getNodalBasis(MSH_QUA_49);
+ case 7: return BasisFactory::getNodalBasis(MSH_QUA_64);
+ case 8: return BasisFactory::getNodalBasis(MSH_QUA_81);
+ case 9: return BasisFactory::getNodalBasis(MSH_QUA_100);
+ case 10: return BasisFactory::getNodalBasis(MSH_QUA_121);
default: Msg::Error("Order %d quadrangle function space not implemented", order);
}
return 0;
@@ -63,29 +63,29 @@ const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBasis::find(MSH_QUA_4);
- case 2: return JacobianBasis::find(MSH_QUA_8);
- case 3: return JacobianBasis::find(MSH_QUA_12);
- case 4: return JacobianBasis::find(MSH_QUA_16I);
- case 5: return JacobianBasis::find(MSH_QUA_20);
- case 6: return JacobianBasis::find(MSH_QUA_24);
- case 7: return JacobianBasis::find(MSH_QUA_28);
- case 8: return JacobianBasis::find(MSH_QUA_32);
- case 9: return JacobianBasis::find(MSH_QUA_36I);
- case 10: return JacobianBasis::find(MSH_QUA_40);
+ case 1: return BasisFactory::getJacobianBasis(MSH_QUA_4);
+ case 2: return BasisFactory::getJacobianBasis(MSH_QUA_8);
+ case 3: return BasisFactory::getJacobianBasis(MSH_QUA_12);
+ case 4: return BasisFactory::getJacobianBasis(MSH_QUA_16I);
+ case 5: return BasisFactory::getJacobianBasis(MSH_QUA_20);
+ case 6: return BasisFactory::getJacobianBasis(MSH_QUA_24);
+ case 7: return BasisFactory::getJacobianBasis(MSH_QUA_28);
+ case 8: return BasisFactory::getJacobianBasis(MSH_QUA_32);
+ case 9: return BasisFactory::getJacobianBasis(MSH_QUA_36I);
+ case 10: return BasisFactory::getJacobianBasis(MSH_QUA_40);
}
}
switch (order) {
- case 1: return JacobianBasis::find(MSH_QUA_4);
- case 2: return JacobianBasis::find(MSH_QUA_9);
- case 3: return JacobianBasis::find(MSH_QUA_16);
- case 4: return JacobianBasis::find(MSH_QUA_25);
- case 5: return JacobianBasis::find(MSH_QUA_36);
- case 6: return JacobianBasis::find(MSH_QUA_49);
- case 7: return JacobianBasis::find(MSH_QUA_64);
- case 8: return JacobianBasis::find(MSH_QUA_81);
- case 9: return JacobianBasis::find(MSH_QUA_100);
- case 10: return JacobianBasis::find(MSH_QUA_121);
+ case 1: return BasisFactory::getJacobianBasis(MSH_QUA_4);
+ case 2: return BasisFactory::getJacobianBasis(MSH_QUA_9);
+ case 3: return BasisFactory::getJacobianBasis(MSH_QUA_16);
+ case 4: return BasisFactory::getJacobianBasis(MSH_QUA_25);
+ case 5: return BasisFactory::getJacobianBasis(MSH_QUA_36);
+ case 6: return BasisFactory::getJacobianBasis(MSH_QUA_49);
+ case 7: return BasisFactory::getJacobianBasis(MSH_QUA_64);
+ case 8: return BasisFactory::getJacobianBasis(MSH_QUA_81);
+ case 9: return BasisFactory::getJacobianBasis(MSH_QUA_100);
+ case 10: return BasisFactory::getJacobianBasis(MSH_QUA_121);
default: Msg::Error("Order %d quadrangle function space not implemented", order);
}
return 0;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 5d1d29a0b1..af35188863 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -8,7 +8,7 @@
#include "Numeric.h"
#include "Context.h"
#include "BasisFactory.h"
-#include "JacobianBasis.h"
+
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
@@ -108,33 +108,33 @@ const nodalBasis* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 0: return BasisFactory::create(MSH_TET_1);
- case 1: return BasisFactory::create(MSH_TET_4);
- case 2: return BasisFactory::create(MSH_TET_10);
- case 3: return BasisFactory::create(MSH_TET_20);
- case 4: return BasisFactory::create(MSH_TET_34);
- case 5: return BasisFactory::create(MSH_TET_52);
- case 6: return BasisFactory::create(MSH_TET_74);
- case 7: return BasisFactory::create(MSH_TET_100);
- case 8: return BasisFactory::create(MSH_TET_130);
- case 9: return BasisFactory::create(MSH_TET_164);
- case 10: return BasisFactory::create(MSH_TET_202);
+ case 0: return BasisFactory::getNodalBasis(MSH_TET_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_TET_4);
+ case 2: return BasisFactory::getNodalBasis(MSH_TET_10);
+ case 3: return BasisFactory::getNodalBasis(MSH_TET_20);
+ case 4: return BasisFactory::getNodalBasis(MSH_TET_34);
+ case 5: return BasisFactory::getNodalBasis(MSH_TET_52);
+ case 6: return BasisFactory::getNodalBasis(MSH_TET_74);
+ case 7: return BasisFactory::getNodalBasis(MSH_TET_100);
+ case 8: return BasisFactory::getNodalBasis(MSH_TET_130);
+ case 9: return BasisFactory::getNodalBasis(MSH_TET_164);
+ case 10: return BasisFactory::getNodalBasis(MSH_TET_202);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return BasisFactory::create(MSH_TET_1);
- case 1: return BasisFactory::create(MSH_TET_4);
- case 2: return BasisFactory::create(MSH_TET_10);
- case 3: return BasisFactory::create(MSH_TET_20);
- case 4: return BasisFactory::create(MSH_TET_35);
- case 5: return BasisFactory::create(MSH_TET_56);
- case 6: return BasisFactory::create(MSH_TET_84);
- case 7: return BasisFactory::create(MSH_TET_120);
- case 8: return BasisFactory::create(MSH_TET_165);
- case 9: return BasisFactory::create(MSH_TET_220);
- case 10: return BasisFactory::create(MSH_TET_286);
+ case 0: return BasisFactory::getNodalBasis(MSH_TET_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_TET_4);
+ case 2: return BasisFactory::getNodalBasis(MSH_TET_10);
+ case 3: return BasisFactory::getNodalBasis(MSH_TET_20);
+ case 4: return BasisFactory::getNodalBasis(MSH_TET_35);
+ case 5: return BasisFactory::getNodalBasis(MSH_TET_56);
+ case 6: return BasisFactory::getNodalBasis(MSH_TET_84);
+ case 7: return BasisFactory::getNodalBasis(MSH_TET_120);
+ case 8: return BasisFactory::getNodalBasis(MSH_TET_165);
+ case 9: return BasisFactory::getNodalBasis(MSH_TET_220);
+ case 10: return BasisFactory::getNodalBasis(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
@@ -149,26 +149,26 @@ const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBasis::find(MSH_TET_4);
- case 2: return JacobianBasis::find(MSH_TET_10);
- case 3: return JacobianBasis::find(MSH_TET_20);
- case 4: return JacobianBasis::find(MSH_TET_34);
- case 5: return JacobianBasis::find(MSH_TET_52);
+ case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
+ case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
+ case 3: return BasisFactory::getJacobianBasis(MSH_TET_20);
+ case 4: return BasisFactory::getJacobianBasis(MSH_TET_34);
+ case 5: return BasisFactory::getJacobianBasis(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBasis::find(MSH_TET_4);
- case 2: return JacobianBasis::find(MSH_TET_10);
- case 3: return JacobianBasis::find(MSH_TET_20);
- case 4: return JacobianBasis::find(MSH_TET_35);
- case 5: return JacobianBasis::find(MSH_TET_56);
- case 6: return JacobianBasis::find(MSH_TET_84);
- case 7: return JacobianBasis::find(MSH_TET_120);
- case 8: return JacobianBasis::find(MSH_TET_165);
- case 9: return JacobianBasis::find(MSH_TET_220);
- case 10: return JacobianBasis::find(MSH_TET_286);
+ case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
+ case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
+ case 3: return BasisFactory::getJacobianBasis(MSH_TET_20);
+ case 4: return BasisFactory::getJacobianBasis(MSH_TET_35);
+ case 5: return BasisFactory::getJacobianBasis(MSH_TET_56);
+ case 6: return BasisFactory::getJacobianBasis(MSH_TET_84);
+ case 7: return BasisFactory::getJacobianBasis(MSH_TET_120);
+ case 8: return BasisFactory::getJacobianBasis(MSH_TET_165);
+ case 9: return BasisFactory::getJacobianBasis(MSH_TET_220);
+ case 10: return BasisFactory::getJacobianBasis(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 49d68dee7b..3c73eb8de4 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -127,33 +127,33 @@ const nodalBasis* MTriangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 0: return BasisFactory::create(MSH_TRI_1);
- case 1: return BasisFactory::create(MSH_TRI_3);
- case 2: return BasisFactory::create(MSH_TRI_6);
- case 3: return BasisFactory::create(MSH_TRI_9);
- case 4: return BasisFactory::create(MSH_TRI_12);
- case 5: return BasisFactory::create(MSH_TRI_15I);
- case 6: return BasisFactory::create(MSH_TRI_18);
- case 7: return BasisFactory::create(MSH_TRI_21I);
- case 8: return BasisFactory::create(MSH_TRI_24);
- case 9: return BasisFactory::create(MSH_TRI_27);
- case 10: return BasisFactory::create(MSH_TRI_30);
+ case 0: return BasisFactory::getNodalBasis(MSH_TRI_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_TRI_3);
+ case 2: return BasisFactory::getNodalBasis(MSH_TRI_6);
+ case 3: return BasisFactory::getNodalBasis(MSH_TRI_9);
+ case 4: return BasisFactory::getNodalBasis(MSH_TRI_12);
+ case 5: return BasisFactory::getNodalBasis(MSH_TRI_15I);
+ case 6: return BasisFactory::getNodalBasis(MSH_TRI_18);
+ case 7: return BasisFactory::getNodalBasis(MSH_TRI_21I);
+ case 8: return BasisFactory::getNodalBasis(MSH_TRI_24);
+ case 9: return BasisFactory::getNodalBasis(MSH_TRI_27);
+ case 10: return BasisFactory::getNodalBasis(MSH_TRI_30);
default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
switch (order) {
- case 0: return BasisFactory::create(MSH_TRI_1);
- case 1: return BasisFactory::create(MSH_TRI_3);
- case 2: return BasisFactory::create(MSH_TRI_6);
- case 3: return BasisFactory::create(MSH_TRI_10);
- case 4: return BasisFactory::create(MSH_TRI_15);
- case 5: return BasisFactory::create(MSH_TRI_21);
- case 6: return BasisFactory::create(MSH_TRI_28);
- case 7: return BasisFactory::create(MSH_TRI_36);
- case 8: return BasisFactory::create(MSH_TRI_45);
- case 9: return BasisFactory::create(MSH_TRI_55);
- case 10: return BasisFactory::create(MSH_TRI_66);
+ case 0: return BasisFactory::getNodalBasis(MSH_TRI_1);
+ case 1: return BasisFactory::getNodalBasis(MSH_TRI_3);
+ case 2: return BasisFactory::getNodalBasis(MSH_TRI_6);
+ case 3: return BasisFactory::getNodalBasis(MSH_TRI_10);
+ case 4: return BasisFactory::getNodalBasis(MSH_TRI_15);
+ case 5: return BasisFactory::getNodalBasis(MSH_TRI_21);
+ case 6: return BasisFactory::getNodalBasis(MSH_TRI_28);
+ case 7: return BasisFactory::getNodalBasis(MSH_TRI_36);
+ case 8: return BasisFactory::getNodalBasis(MSH_TRI_45);
+ case 9: return BasisFactory::getNodalBasis(MSH_TRI_55);
+ case 10: return BasisFactory::getNodalBasis(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
@@ -169,32 +169,32 @@ const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 0: return JacobianBasis::find(MSH_TRI_1);
- case 1: return JacobianBasis::find(MSH_TRI_3);
- case 2: return JacobianBasis::find(MSH_TRI_6);
- case 3: return JacobianBasis::find(MSH_TRI_9);
- case 4: return JacobianBasis::find(MSH_TRI_12);
- case 5: return JacobianBasis::find(MSH_TRI_15I);
- case 6: return JacobianBasis::find(MSH_TRI_18);
- case 7: return JacobianBasis::find(MSH_TRI_21I);
- case 8: return JacobianBasis::find(MSH_TRI_24);
- case 9: return JacobianBasis::find(MSH_TRI_27);
- case 10: return JacobianBasis::find(MSH_TRI_30);
+ case 0: return BasisFactory::getJacobianBasis(MSH_TRI_1);
+ case 1: return BasisFactory::getJacobianBasis(MSH_TRI_3);
+ case 2: return BasisFactory::getJacobianBasis(MSH_TRI_6);
+ case 3: return BasisFactory::getJacobianBasis(MSH_TRI_9);
+ case 4: return BasisFactory::getJacobianBasis(MSH_TRI_12);
+ case 5: return BasisFactory::getJacobianBasis(MSH_TRI_15I);
+ case 6: return BasisFactory::getJacobianBasis(MSH_TRI_18);
+ case 7: return BasisFactory::getJacobianBasis(MSH_TRI_21I);
+ case 8: return BasisFactory::getJacobianBasis(MSH_TRI_24);
+ case 9: return BasisFactory::getJacobianBasis(MSH_TRI_27);
+ case 10: return BasisFactory::getJacobianBasis(MSH_TRI_30);
default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBasis::find(MSH_TRI_3);
- case 2: return JacobianBasis::find(MSH_TRI_6);
- case 3: return JacobianBasis::find(MSH_TRI_10);
- case 4: return JacobianBasis::find(MSH_TRI_15);
- case 5: return JacobianBasis::find(MSH_TRI_21);
- case 6: return JacobianBasis::find(MSH_TRI_28);
- case 7: return JacobianBasis::find(MSH_TRI_36);
- case 8: return JacobianBasis::find(MSH_TRI_45);
- case 9: return JacobianBasis::find(MSH_TRI_55);
- case 10: return JacobianBasis::find(MSH_TRI_66);
+ case 1: return BasisFactory::getJacobianBasis(MSH_TRI_3);
+ case 2: return BasisFactory::getJacobianBasis(MSH_TRI_6);
+ case 3: return BasisFactory::getJacobianBasis(MSH_TRI_10);
+ case 4: return BasisFactory::getJacobianBasis(MSH_TRI_15);
+ case 5: return BasisFactory::getJacobianBasis(MSH_TRI_21);
+ case 6: return BasisFactory::getJacobianBasis(MSH_TRI_28);
+ case 7: return BasisFactory::getJacobianBasis(MSH_TRI_36);
+ case 8: return BasisFactory::getJacobianBasis(MSH_TRI_45);
+ case 9: return BasisFactory::getJacobianBasis(MSH_TRI_55);
+ case 10: return BasisFactory::getJacobianBasis(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c79bd36304..7ce3f119ff 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -501,7 +501,7 @@ double discreteEdge::curvature(double par) const
double c0, c1;
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
- std::cout << "Need to compute discrete curvature (in discreteEdge)" << std::endl;
+ Msg::Warning("Need to compute discrete curvature (in discreteEdge)");
Curvature::typeOfCurvature type = Curvature::RUSIN; //RUSIN; //RBF
curvature.computeCurvature(model(), type);
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index dbc93f24e1..59f54a2bb0 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -472,14 +472,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
switch (nPts){
case 0 : break;
case 1 : 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;
+ case 2 : points = BasisFactory::getNodalBasis(MSH_TRI_10)->points; start = 9; break;
+ case 3 : points = BasisFactory::getNodalBasis(MSH_TRI_15)->points; start = 12; break;
+ case 4 : points = BasisFactory::getNodalBasis(MSH_TRI_21)->points; start = 15; break;
+ case 5 : points = BasisFactory::getNodalBasis(MSH_TRI_28)->points; start = 18; break;
+ case 6 : points = BasisFactory::getNodalBasis(MSH_TRI_36)->points; start = 21; break;
+ case 7 : points = BasisFactory::getNodalBasis(MSH_TRI_45)->points; start = 24; break;
+ case 8 : points = BasisFactory::getNodalBasis(MSH_TRI_55)->points; start = 27; break;
+ case 9 : points = BasisFactory::getNodalBasis(MSH_TRI_66)->points; start = 30; break;
default :
Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
break;
@@ -488,14 +488,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
case 4:
switch (nPts){
case 0 : 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;
+ case 1 : points = BasisFactory::getNodalBasis(MSH_QUA_9)->points; break;
+ case 2 : points = BasisFactory::getNodalBasis(MSH_QUA_16)->points; break;
+ case 3 : points = BasisFactory::getNodalBasis(MSH_QUA_25)->points; break;
+ case 4 : points = BasisFactory::getNodalBasis(MSH_QUA_36)->points; break;
+ case 5 : points = BasisFactory::getNodalBasis(MSH_QUA_49)->points; break;
+ case 6 : points = BasisFactory::getNodalBasis(MSH_QUA_64)->points; break;
+ case 7 : points = BasisFactory::getNodalBasis(MSH_QUA_81)->points; break;
+ case 8 : points = BasisFactory::getNodalBasis(MSH_QUA_100)->points; break;
default :
Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
break;
@@ -599,15 +599,15 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
case 0: return;
case 1: return;
case 2:
- BasisFactory::create(MSH_TET_20)->points.print();
- points = BasisFactory::create(MSH_TET_20)->points; break;
- case 3: points = BasisFactory::create(MSH_TET_35)->points; break;
- case 4: points = BasisFactory::create(MSH_TET_56)->points; break;
- case 5: points = BasisFactory::create(MSH_TET_84)->points; break;
- case 6: points = BasisFactory::create(MSH_TET_120)->points; break;
- case 7: points = BasisFactory::create(MSH_TET_165)->points; break;
- case 8: points = BasisFactory::create(MSH_TET_220)->points; break;
- case 9: points = BasisFactory::create(MSH_TET_286)->points; break;
+ BasisFactory::getNodalBasis(MSH_TET_20)->points.print();
+ points = BasisFactory::getNodalBasis(MSH_TET_20)->points; break;
+ case 3: points = BasisFactory::getNodalBasis(MSH_TET_35)->points; break;
+ case 4: points = BasisFactory::getNodalBasis(MSH_TET_56)->points; break;
+ case 5: points = BasisFactory::getNodalBasis(MSH_TET_84)->points; break;
+ case 6: points = BasisFactory::getNodalBasis(MSH_TET_120)->points; break;
+ case 7: points = BasisFactory::getNodalBasis(MSH_TET_165)->points; break;
+ case 8: points = BasisFactory::getNodalBasis(MSH_TET_220)->points; break;
+ case 9: points = BasisFactory::getNodalBasis(MSH_TET_286)->points; break;
default:
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
@@ -617,14 +617,14 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
case TYPE_HEX :
switch (nPts){
case 0: return;
- 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;
+ case 1: points = BasisFactory::getNodalBasis(MSH_HEX_27)->points; break;
+ case 2: points = BasisFactory::getNodalBasis(MSH_HEX_64)->points; break;
+ case 3: points = BasisFactory::getNodalBasis(MSH_HEX_125)->points; break;
+ case 4: points = BasisFactory::getNodalBasis(MSH_HEX_216)->points; break;
+ case 5: points = BasisFactory::getNodalBasis(MSH_HEX_343)->points; break;
+ case 6: points = BasisFactory::getNodalBasis(MSH_HEX_512)->points; break;
+ case 7: points = BasisFactory::getNodalBasis(MSH_HEX_729)->points; break;
+ case 8: points = BasisFactory::getNodalBasis(MSH_HEX_1000)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
@@ -635,13 +635,13 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
switch (nPts){
case 0:
case 1: return;
- case 2: points = BasisFactory::create(MSH_PYR_30)->points; break;
- case 3: points = BasisFactory::create(MSH_PYR_55)->points; break;
- case 4: points = BasisFactory::create(MSH_PYR_91)->points; break;
- case 5: points = BasisFactory::create(MSH_PYR_140)->points; break;
- case 6: points = BasisFactory::create(MSH_PYR_204)->points; break;
- case 7: points = BasisFactory::create(MSH_PYR_285)->points; break;
- case 8: points = BasisFactory::create(MSH_PYR_385)->points; break;
+ case 2: points = BasisFactory::getNodalBasis(MSH_PYR_30)->points; break;
+ case 3: points = BasisFactory::getNodalBasis(MSH_PYR_55)->points; break;
+ case 4: points = BasisFactory::getNodalBasis(MSH_PYR_91)->points; break;
+ case 5: points = BasisFactory::getNodalBasis(MSH_PYR_140)->points; break;
+ case 6: points = BasisFactory::getNodalBasis(MSH_PYR_204)->points; break;
+ case 7: points = BasisFactory::getNodalBasis(MSH_PYR_285)->points; break;
+ case 8: points = BasisFactory::getNodalBasis(MSH_PYR_385)->points; break;
default :
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
diff --git a/Mesh/meshGFaceRecombine.h b/Mesh/meshGFaceRecombine.h
index 486f70f366..94c7cc8f61 100644
--- a/Mesh/meshGFaceRecombine.h
+++ b/Mesh/meshGFaceRecombine.h
@@ -416,7 +416,7 @@ class Rec2DChange {
void *_info;
public :
- Rec2DChange() {Msg::Error("[Rec2DChange] I should not be created in this manner");}
+ Rec2DChange() {Msg::Fatal("[Rec2DChange] I should not be created in this manner");}
Rec2DChange(int);
Rec2DChange(Rec2DEdge*, bool toHide = false);
Rec2DChange(Rec2DVertex*, bool toHide = false);
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 46adb2e19f..36730f5950 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -11,9 +11,11 @@
#include "BasisFactory.h"
std::map<int, nodalBasis*> BasisFactory::fs;
+std::map<int, JacobianBasis*> BasisFactory::js;
+std::map<int, bezierBasis*> BasisFactory::bs;
-const nodalBasis* BasisFactory::create(int elementType) {
-
+const nodalBasis* BasisFactory::getNodalBasis(int elementType)
+{
// If the Basis has already been built, return it.
std::map<int, nodalBasis*>::const_iterator it = fs.find(elementType);
if (it != fs.end()) {
@@ -22,7 +24,7 @@ const nodalBasis* BasisFactory::create(int elementType) {
// Get the parent type to see which kind of basis
// we want to create
int parentType = MElement::ParentTypeFromTag(elementType);
- nodalBasis* B = 0;
+ nodalBasis* F = NULL;
switch(parentType) {
case(TYPE_PNT):
@@ -32,19 +34,40 @@ const nodalBasis* BasisFactory::create(int elementType) {
case(TYPE_PRI):
case(TYPE_TET):
case(TYPE_HEX):
- B = new polynomialBasis(elementType);
+ F = new polynomialBasis(elementType);
break;
case(TYPE_PYR):
- B = new pyramidalBasis(elementType);
+ F = new pyramidalBasis(elementType);
break;
default:
Msg::Error("Unknown type of element.");
- return 0;
+ return NULL;
}
// FIXME: check if already exists to deallocate if necessary
- fs.insert(std::make_pair(elementType, B));
+ fs.insert(std::make_pair(elementType, F));
return fs[elementType];
+}
+
+const bezierBasis* BasisFactory::getBezierBasis(int elementType)
+{
+ std::map<int, bezierBasis*>::const_iterator it = bs.find(elementType);
+ if (it != bs.end())
+ return it->second;
+
+ bezierBasis* B = new bezierBasis(elementType);
+ if (B) bs.insert(std::make_pair(elementType, B));
+ return B;
+}
+
+const JacobianBasis* BasisFactory::getJacobianBasis(int elementType)
+{
+ std::map<int, JacobianBasis*>::const_iterator it = js.find(elementType);
+ if (it != js.end())
+ return it->second;
+ JacobianBasis* J = new JacobianBasis(elementType);
+ if (J) js.insert(std::make_pair(elementType, J));
+ return J;
}
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 7932c3eae3..af4f6d7e98 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -9,16 +9,21 @@
#include "MElement.h"
#include "MPyramid.h"
#include "nodalBasis.h"
+#include "JacobianBasis.h"
+
class BasisFactory
{
private:
static std::map<int, nodalBasis*> fs;
+ static std::map<int, bezierBasis*> bs;
+ static std::map<int, JacobianBasis*> js;
public:
- // Caution: the returned pointer should be
- // checked against 0.
- static const nodalBasis* create(int elementType);
+ // Caution: the returned pointer can be NULL
+ static const nodalBasis* getNodalBasis(int elementType);
+ static const bezierBasis* getBezierBasis(int elementType);
+ static const JacobianBasis* getJacobianBasis(int elementType);
};
#endif
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 0931e65645..231bfe602c 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -8,12 +8,13 @@ set(SRC
fullMatrix.cpp
BasisFactory.cpp
nodalBasis.cpp
- polynomialBasis.cpp
+ polynomialBasis.cpp
+ pyramidalBasis.cpp
+ BergotBasis.cpp
+ jacobiPolynomials.cpp
+ legendrePolynomials.cpp
+ bezierBasis.cpp
JacobianBasis.cpp
- BergotBasis.cpp
- pyramidalBasis.cpp
- jacobiPolynomials.cpp
- legendrePolynomials.cpp
pointsGenerators.cpp
GaussIntegration.cpp
GaussQuadratureLin.cpp
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 8ac87aa574..4b4e661032 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -5,7 +5,7 @@
#include "GmshDefines.h"
#include "GmshMessage.h"
-#include "JacobianBasis.h"
+
#include <vector>
#include "polynomialBasis.h"
#include "pyramidalBasis.h"
@@ -13,1068 +13,6 @@
#include "BasisFactory.h"
#include "Numeric.h"
-// Bezier Exponents
-static fullMatrix<double> generate1DExponents(int order)
-{
- fullMatrix<double> exponents(order + 1, 1);
- exponents(0, 0) = 0;
- if (order > 0) {
- exponents(1, 0) = order;
- for (int i = 2; i < order + 1; i++) exponents(i, 0) = i - 1;
- }
-
- return exponents;
-}
-
-static fullMatrix<double> generateExponentsTriangle(int order)
-{
- int nbPoints = (order + 1) * (order + 2) / 2;
- fullMatrix<double> exponents(nbPoints, 2);
-
- exponents(0, 0) = 0;
- exponents(0, 1) = 0;
-
- if (order > 0) {
- exponents(1, 0) = order;
- exponents(1, 1) = 0;
- exponents(2, 0) = 0;
- exponents(2, 1) = order;
-
- if (order > 1) {
- int index = 3, ksi = 0, eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi++;
- exponents(index, 0) = ksi;
- exponents(index, 1) = eta;
- }
-
- ksi = order;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi--;
- eta++;
- exponents(index, 0) = ksi;
- exponents(index, 1) = eta;
- }
-
- eta = order;
- ksi = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta--;
- exponents(index, 0) = ksi;
- exponents(index, 1) = eta;
- }
-
- if (order > 2) {
- fullMatrix<double> inner = generateExponentsTriangle(order - 3);
- inner.add(1);
- exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
-
- return exponents;
-}
-static fullMatrix<double> generateExponentsQuad(int order)
-{
- int nbPoints = (order+1)*(order+1);
- fullMatrix<double> exponents(nbPoints, 2);
-
- exponents(0, 0) = 0;
- exponents(0, 1) = 0;
-
- if (order > 0) {
- exponents(1, 0) = order;
- exponents(1, 1) = 0;
- exponents(2, 0) = order;
- exponents(2, 1) = order;
- exponents(3, 0) = 0;
- exponents(3, 1) = order;
-
- 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++) {
- exponents(index, 0) = exponents(p0, 0) + i*(exponents(p1,0)-exponents(p0,0))/order;
- exponents(index, 1) = exponents(p0, 1) + i*(exponents(p1,1)-exponents(p0,1))/order;
- }
- }
- if (order > 2) {
- fullMatrix<double> inner = generateExponentsQuad(order - 2);
- inner.add(1);
- exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
-
- // exponents.print("expq");
-
- return exponents;
-}
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{ // uv surface - orientation v0-v2-v1
- 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;
- }
-}
-
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{ // uw surface - orientation v0-v1-v3
- 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;
- }
-}
-
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{ // vw surface - orientation v0-v3-v2
- 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;
- }
-}
-
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{ // uvw surface - orientation v3-v1-v2
- 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> generateExponentsTetrahedron(int order)
-{
- int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
-
- fullMatrix<double> exponents(nbPoints, 3);
-
- exponents(0, 0) = 0;
- exponents(0, 1) = 0;
- exponents(0, 2) = 0;
-
- if (order > 0) {
- exponents(1, 0) = order;
- exponents(1, 1) = 0;
- exponents(1, 2) = 0;
-
- exponents(2, 0) = 0;
- exponents(2, 1) = order;
- exponents(2, 2) = 0;
-
- exponents(3, 0) = 0;
- exponents(3, 1) = 0;
- exponents(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++) {
- exponents(4 + k, 0) = k + 1;
- exponents(4 + order - 1 + k, 0) = order - 1 - k;
- exponents(4 + 2 * (order - 1) + k, 0) = 0;
- exponents(4 + 3 * (order - 1) + k, 0) = 0;
- // exponents(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
- // exponents(4 + 5 * (order - 1) + k, 0) = 0.;
- exponents(4 + 4 * (order - 1) + k, 0) = 0;
- exponents(4 + 5 * (order - 1) + k, 0) = k+1;
-
- exponents(4 + k, 1) = 0;
- exponents(4 + order - 1 + k, 1) = k + 1;
- exponents(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- exponents(4 + 3 * (order - 1) + k, 1) = 0;
- // exponents(4 + 4 * (order - 1) + k, 1) = 0.;
- // exponents(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- exponents(4 + 4 * (order - 1) + k, 1) = k+1;
- exponents(4 + 5 * (order - 1) + k, 1) = 0;
-
- exponents(4 + k, 2) = 0;
- exponents(4 + order - 1 + k, 2) = 0;
- exponents(4 + 2 * (order - 1) + k, 2) = 0;
- exponents(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- exponents(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- exponents(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++){
- exponents(ns + i, 0) = u[i] * (order - 3) + 1;
- exponents(ns + i, 1) = v[i] * (order - 3) + 1;
- exponents(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++){
- exponents(ns + i, 0) = u[i] * (order - 3) + 1;
- exponents(ns + i, 1) = v[i] * (order - 3) ;
- exponents(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++){
- exponents(ns + i, 0) = u[i] * (order - 3);
- exponents(ns + i, 1) = v[i] * (order - 3) + 1;
- exponents(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++){
- exponents(ns + i, 0) = u[i] * (order - 3) + 1;
- exponents(ns + i, 1) = v[i] * (order - 3) + 1;
- exponents(ns + i, 2) = w[i] * (order - 3) + 1;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (order > 3) {
-
- fullMatrix<double> interior = generateExponentsTetrahedron(order - 4);
- for (int k = 0; k < interior.size1(); k++) {
- exponents(ns + k, 0) = 1 + interior(k, 0);
- exponents(ns + k, 1) = 1 + interior(k, 1);
- exponents(ns + k, 2) = 1 + interior(k, 2);
- }
- }
- }
- }
- }
-
- return exponents;
-}
-
-static fullMatrix<double> generateExponentsPrism(int order)
-{
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> exponents(nbPoints, 3);
-
- int index = 0;
- fullMatrix<double> triExp = generateExponentsTriangle(order);
- fullMatrix<double> lineExp = generate1DExponents(order);
- // First exponents (i < 3) must relate to corner
- for (int i = 0; i < triExp.size1(); i++) {
- exponents(index,0) = triExp(i,0);
- exponents(index,1) = triExp(i,1);
- exponents(index,2) = lineExp(0,0);
- index ++;
- exponents(index,0) = triExp(i,0);
- exponents(index,1) = triExp(i,1);
- exponents(index,2) = lineExp(1,0);
- index ++;
- }
- for (int j = 2; j <lineExp.size1() ; j++) {
- for (int i = 0; i < triExp.size1(); i++) {
- exponents(index,0) = triExp(i,0);
- exponents(index,1) = triExp(i,1);
- exponents(index,2) = lineExp(j,0);
- index ++;
- }
- }
-
- return exponents;
-}
-
-static fullMatrix<double> generateExponentsHex(int order)
-{
- int nbPoints = (order+1) * (order+1) * (order+1);
- fullMatrix<double> exponents(nbPoints, 3);
-
- if (order == 0) {
- exponents(0, 0) = .0;
- return exponents;
- }
-
- int index = 0;
- fullMatrix<double> lineExp = generate1DExponents(order);
-
- for (int i = 0; i < 2; ++i) {
- 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);
- ++index;
- }
- }
- }
-
- for (int i = 2; i < lineExp.size1(); ++i) {
- 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);
- ++index;
- }
- }
- }
- for (int i = 0; i < lineExp.size1(); ++i) {
- 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);
- ++index;
- }
- }
- }
- for (int i = 0; i < lineExp.size1(); ++i) {
- 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);
- ++index;
- }
- }
- }
-
- return exponents;
-}
-
-// Sampling Points
-static fullMatrix<double> generate1DPoints(int order)
-{
- fullMatrix<double> line(order + 1, 1);
- line(0,0) = .0;
- if (order > 0) {
- line(1, 0) = 1.;
- double dd = 1. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
- }
-
- return line;
-}
-
-static fullMatrix<double> generatePointsQuadRecur(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;
- }
- }
- if (order >= 2 && !serendip) {
- fullMatrix<double> inner = generatePointsQuadRecur(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> generatePointsQuad(int order, bool serendip)
-{
- fullMatrix<double> point = generatePointsQuadRecur(order, serendip);
- // rescale to [0,1] x [0,1]
- for (int i=0;i<point.size1();i++){
- point(i,0) = (1.+point(i,0))*.5;
- point(i,1) = (1.+point(i,1))*.5;
- }
- return point;
-}
-
-static fullMatrix<double> generatePointsPrism(int order, bool serendip)
-{
- const double prism18Pts[18][3] = {
- {0, 0, 0}, // 0
- {1, 0, 0}, // 1
- {0, 1, 0}, // 2
- {0, 0, 1}, // 3
- {1, 0, 1}, // 4
- {0, 1, 1}, // 5
- {.5, 0, 0}, // 6
- {0, .5, 0}, // 7
- {0, 0, .5}, // 8
- {.5, .5, 0}, // 9
- {1, 0, .5}, // 10
- {0, 1, .5}, // 11
- {.5, 0, 1}, // 12
- {0, .5, 1}, // 13
- {.5, .5, 1}, // 14
- {.5, 0, .5}, // 15
- {0, .5, .5}, // 16
- {.5, .5, .5}, // 17
- };
-
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> point(nbPoints, 3);
-
- int index = 0;
-
- if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- point(i,j) = prism18Pts[i][j];
- else {
- fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
- fullMatrix<double> linePoints = generate1DPoints(order);
- 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> generatePointsHex(int order, bool serendip)
-{
- int nbPoints = (order+1) * (order+1) * (order+1);
- fullMatrix<double> point(nbPoints, 3);
-
- fullMatrix<double> linePoints = generate1DPoints(order);
- int index = 0;
-
- for (int i = 0; i < linePoints.size1(); ++i) {
- for (int j = 0; j < linePoints.size1(); ++j) {
- for (int k = 0; k < linePoints.size1(); ++k) {
- point(index, 0) = linePoints(i, 0);
- point(index, 1) = linePoints(j, 0);
- point(index, 2) = linePoints(k, 0);
- ++index;
- }
- }
- }
-
- return point;
-}
-
-// Sub Control Points
-static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
-{
- std::vector< fullMatrix<double> > subPoints(2);
-
- subPoints[0] = generate1DPoints(order);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- subPoints[1].add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
-{
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
-
- subPoints[0] = gmshGeneratePointsTriangle(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- subPoints[3].scale(-1.);
- subPoints[3].add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
-{
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
-
- subPoints[0] = generatePointsQuad(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox1;
- fullMatrix<double> prox2;
-
- subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[1], 0, 1);
- prox1.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[2], 1, 1);
- prox1.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[3], 2, 1);
- prox1.add(.5);
-
- // u := .5-u-w
- // v := .5-v-w
- // w := w
- subPoints[4].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[4], 0, 2);
- prox1.scale(-1.);
- prox1.add(.5);
- prox1.setAsProxy(subPoints[4], 0, 1);
- prox2.setAsProxy(subPoints[4], 2, 1);
- prox1.add(prox2, -1.);
- prox1.setAsProxy(subPoints[4], 1, 1);
- prox1.add(prox2, -1.);
-
- // u := u
- // v := .5-v
- // w := w+v
- subPoints[5].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[5], 2, 1);
- prox2.setAsProxy(subPoints[5], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := .5-u
- // v := v
- // w := w+u
- subPoints[6].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[6], 2, 1);
- prox2.setAsProxy(subPoints[6], 0, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := u+w
- // v := v+w
- // w := .5-w
- subPoints[7].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[7], 0, 1);
- prox2.setAsProxy(subPoints[7], 2, 1);
- prox1.add(prox2);
- prox1.setAsProxy(subPoints[7], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
-
- subPoints[0] = generatePointsPrism(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[0]);
- prox.setAsProxy(subPoints[3], 0, 2);
- prox.scale(-1.);
- prox.add(.5);
-
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
-
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
-
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
-
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
-
- return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
-{
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
-
- subPoints[0] = generatePointsHex(order, false);
- subPoints[0].scale(.5);
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
-
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
-
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
-
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
-
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
-
- return subPoints;
-}
-
-// Matrices generation
-static int nChoosek(int n, int k)
-{
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination.");
- return 1;
- }
-
- if (k > n/2) k = n-k;
- if (k == 1)
- return n;
- if (k == 0)
- return 1;
-
- int c = 1;
- for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
- return c;
-}
-
-
-static fullMatrix<double> generateBez2LagMatrix
- (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order, int dimSimplex)
-{
-
- if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
-
- int ndofs = exponent.size1();
- int dim = exponent.size2();
-
- fullMatrix<double> bez2Lag(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
-
- {
- double pointCompl = 1.;
- int exponentCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exponentCompl, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k));
- pointCompl -= point(j, k);
- exponentCompl -= (int) exponent(i, k);
- }
- dd *= pow(pointCompl, exponentCompl);
- }
-
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k))
- * pow(1. - point(j, k), order - exponent(i, k));
-
- bez2Lag(j, i) = dd;
- }
- }
- return bez2Lag;
-}
-
-
-
-static fullMatrix<double> generateSubDivisor
- (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
-{
- if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exponents.size1(), lag2Bez.size1(),
- exponents.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
-
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
-
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> subDivisor(nbSubPts, nbPts);
-
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
- lag2Bez.mult(intermediate1, intermediate2);
- subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
- }
- return subDivisor;
-}
-
-
-
-// Convert Bezier points to Lagrange points
-static fullMatrix<double> bez2LagPoints(bool dim1, bool dim2, bool dim3, const fullMatrix<double> &bezierPoints)
-{
-
- // FIXME BUG for 2nd order quads we seem to try to access bezierPoints(i, 2);
- // but bezierPoints only has 2 columns
-
- const int numPt = bezierPoints.size1();
- fullMatrix<double> lagPoints(numPt,3);
-
- for (int i=0; i<numPt; i++) {
- lagPoints(i,0) = dim1 ? -1. + 2*bezierPoints(i,0) : bezierPoints(i,0);
- lagPoints(i,1) = dim2 ? -1. + 2*bezierPoints(i,1) : bezierPoints(i,1);
- lagPoints(i,2) = dim3 ? -1. + 2*bezierPoints(i,2) : bezierPoints(i,2);
- }
-
- return lagPoints;
-
-}
-
-
-
-std::map<int, bezierBasis> bezierBasis::_bbs;
-const bezierBasis *bezierBasis::find(int tag)
-{
- std::map<int, bezierBasis>::const_iterator it = _bbs.find(tag);
- if (it != _bbs.end())
- return &it->second;
-
- bezierBasis &B = _bbs[tag];
- B.order = MElement::OrderFromTag(tag);
-
- if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
- B.dim = 3;
- B.numLagPts = 5;
- B.bezierPoints = gmshGeneratePointsPyramid(B.order,false);
- B.lagPoints = B.bezierPoints;
- B.matrixLag2Bez.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
- B.matrixBez2Lag.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
- for (int i=0; i<B.matrixBez2Lag.size1(); ++i) {
- B.matrixLag2Bez(i,i) = 1.;
- B.matrixBez2Lag(i,i) = 1.;
- }
- // TODO: Implement subdidivsor
- }
- else {
- int dimSimplex;
- std::vector< fullMatrix<double> > subPoints;
- switch (MElement::ParentTypeFromTag(tag)) {
- case TYPE_PNT :
- B.dim = 0;
- B.numLagPts = 1;
- B.exponents = generate1DExponents(0);
- B.bezierPoints = generate1DPoints(0);
- B.lagPoints = B.bezierPoints;
- dimSimplex = 0;
- break;
- case TYPE_LIN : {
- B.dim = 1;
- B.numLagPts = 2;
- B.exponents = generate1DExponents(B.order);
- B.bezierPoints = generate1DPoints(B.order);
- B.lagPoints = bez2LagPoints(true,false,false,B.bezierPoints);
- dimSimplex = 0;
- subPoints = generateSubPointsLine(0);
- break;
- }
- case TYPE_TRI : {
- B.dim = 2;
- B.numLagPts = 3;
- B.exponents = generateExponentsTriangle(B.order);
- B.bezierPoints = gmshGeneratePointsTriangle(B.order,false);
- B.lagPoints = B.bezierPoints;
- dimSimplex = 2;
- subPoints = generateSubPointsTriangle(B.order);
- break;
- }
- case TYPE_QUA : {
- B.dim = 2;
- B.numLagPts = 4;
- B.exponents = generateExponentsQuad(B.order);
- B.bezierPoints = generatePointsQuad(B.order,false);
- B.lagPoints = bez2LagPoints(true,true,false,B.bezierPoints);
- dimSimplex = 0;
- subPoints = generateSubPointsQuad(B.order);
- // B.points.print("points");
- break;
- }
- case TYPE_TET : {
- B.dim = 3;
- B.numLagPts = 4;
- B.exponents = generateExponentsTetrahedron(B.order);
- B.bezierPoints = gmshGeneratePointsTetrahedron(B.order,false);
- B.lagPoints = B.bezierPoints;
- dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(B.order);
- break;
- }
- case TYPE_PRI : {
- B.dim = 3;
- B.numLagPts = 6;
- B.exponents = generateExponentsPrism(B.order);
- B.bezierPoints = generatePointsPrism(B.order, false);
- B.lagPoints = bez2LagPoints(false,false,true,B.bezierPoints);
- dimSimplex = 2;
- subPoints = generateSubPointsPrism(B.order);
- break;
- }
- case TYPE_HEX : {
- B.dim = 3;
- B.numLagPts = 8;
- B.exponents = generateExponentsHex(B.order);
- B.bezierPoints = generatePointsHex(B.order, false);
- B.lagPoints = bez2LagPoints(true,true,true,B.bezierPoints);
- dimSimplex = 0;
- subPoints = generateSubPointsHex(B.order);
- break;
- }
- default : {
- Msg::Error("Unknown function space %d: reverting to TET_1", tag);
- B.dim = 3;
- B.numLagPts = 4;
- B.exponents = generateExponentsTetrahedron(0);
- B.bezierPoints = gmshGeneratePointsTetrahedron(0, false);
- B.lagPoints = B.bezierPoints;
- dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(0);
- break;
- }
- }
- B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.bezierPoints, B.order, dimSimplex);
- B.matrixBez2Lag.invert(B.matrixLag2Bez);
- B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, B.order, dimSimplex);
- B.numDivisions = (int) pow(2., (int) B.bezierPoints.size2());
- }
-
- return &B;
-}
-
-
JacobianBasis::JacobianBasis(int tag)
{
@@ -1113,10 +51,10 @@ JacobianBasis::JacobianBasis(int tag)
}
// Store Bezier basis
- bezier = bezierBasis::find(jacType);
+ bezier = BasisFactory::getBezierBasis(jacType);
// Store shape function gradients of mapping at Jacobian nodes
- const nodalBasis *mapBasis = BasisFactory::create(tag);
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
fullMatrix<double> allDPsi;
mapBasis->df(getPoints(), allDPsi);
numJacNodes = getPoints().size1();
@@ -1132,7 +70,7 @@ JacobianBasis::JacobianBasis(int tag)
}
// Compute matrix for lifting from primary Jacobian basis to Jacobian basis
- const nodalBasis *primJacBasis = BasisFactory::create(primJacType);
+ const nodalBasis *primJacBasis = BasisFactory::getNodalBasis(primJacType);
numPrimJacNodes = primJacBasis->getNumShapeFunctions();
matrixPrimJac2Jac.resize(numJacNodes,numPrimJacNodes);
primJacBasis->f(getPoints(),matrixPrimJac2Jac);
@@ -1140,7 +78,7 @@ JacobianBasis::JacobianBasis(int tag)
// Compute shape function gradients of primary mapping at barycenter,
// in order to compute normal to straight element
const int primMapType = polynomialBasis::getTag(parentType, 1, false);
- const nodalBasis *primMapBasis = BasisFactory::create(primMapType);
+ const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
numPrimMapNodes = primMapBasis->getNumShapeFunctions();
double xBar = 0., yBar = 0., zBar = 0.;
for (int i=0; i<numPrimMapNodes; i++) {
@@ -1168,21 +106,6 @@ JacobianBasis::JacobianBasis(int tag)
-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;
-
- JacobianBasis *B = new JacobianBasis(tag);
- _fs.insert(std::make_pair(tag, B));
- return B;
-
-}
-
-
-
// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
{
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 8901c1d9a5..2e8f52b23c 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -8,28 +8,12 @@
#include <map>
#include <vector>
+#include "bezierBasis.h"
#include "fullMatrix.h"
-class MElement;
-
-class bezierBasis {
- private :
- static std::map<int, bezierBasis> _bbs;
- public :
- int dim, order;
- int numLagPts;
- int numDivisions;
- // the 'numLagPts' first exponents are related to 'Lagrangian' values
- fullMatrix<double> exponents; // exponents of Bezier FF
- fullMatrix<double> bezierPoints, lagPoints; // sampling point
- fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
- fullMatrix<double> subDivisor;
- static const bezierBasis *find(int);
-};
class JacobianBasis {
private:
- static std::map<int, JacobianBasis*> _fs;
const bezierBasis *bezier;
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
fullVector<double> primGradShapeBarX, primGradShapeBarY, primGradShapeBarZ;
@@ -37,7 +21,6 @@ class JacobianBasis {
int numJacNodes, numMapNodes, numPrimJacNodes, numPrimMapNodes;
inline const fullMatrix<double> &getPoints() const { return bezier->lagPoints; }
public :
- static const JacobianBasis *find(int);
JacobianBasis(int tag);
inline int getNumJacNodes() const { return numJacNodes; }
inline int getNumMapNodes() const { return numMapNodes; }
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
new file mode 100644
index 0000000000..9f0592dc93
--- /dev/null
+++ b/Numeric/bezierBasis.cpp
@@ -0,0 +1,1069 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include <vector>
+#include "polynomialBasis.h"
+#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
+#include "BasisFactory.h"
+#include "Numeric.h"
+
+// Bezier Exponents
+static fullMatrix<double> generate1DExponents(int order)
+{
+ fullMatrix<double> exponents(order + 1, 1);
+ exponents(0, 0) = 0;
+ if (order > 0) {
+ exponents(1, 0) = order;
+ for (int i = 2; i < order + 1; i++) exponents(i, 0) = i - 1;
+ }
+
+ return exponents;
+}
+
+static fullMatrix<double> generateExponentsTriangle(int order)
+{
+ int nbPoints = (order + 1) * (order + 2) / 2;
+ fullMatrix<double> exponents(nbPoints, 2);
+
+ exponents(0, 0) = 0;
+ exponents(0, 1) = 0;
+
+ if (order > 0) {
+ exponents(1, 0) = order;
+ exponents(1, 1) = 0;
+ exponents(2, 0) = 0;
+ exponents(2, 1) = order;
+
+ if (order > 1) {
+ int index = 3, ksi = 0, eta = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi++;
+ exponents(index, 0) = ksi;
+ exponents(index, 1) = eta;
+ }
+
+ ksi = order;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ ksi--;
+ eta++;
+ exponents(index, 0) = ksi;
+ exponents(index, 1) = eta;
+ }
+
+ eta = order;
+ ksi = 0;
+
+ for (int i = 0; i < order - 1; i++, index++) {
+ eta--;
+ exponents(index, 0) = ksi;
+ exponents(index, 1) = eta;
+ }
+
+ if (order > 2) {
+ fullMatrix<double> inner = generateExponentsTriangle(order - 3);
+ inner.add(1);
+ exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+
+ return exponents;
+}
+static fullMatrix<double> generateExponentsQuad(int order)
+{
+ int nbPoints = (order+1)*(order+1);
+ fullMatrix<double> exponents(nbPoints, 2);
+
+ exponents(0, 0) = 0;
+ exponents(0, 1) = 0;
+
+ if (order > 0) {
+ exponents(1, 0) = order;
+ exponents(1, 1) = 0;
+ exponents(2, 0) = order;
+ exponents(2, 1) = order;
+ exponents(3, 0) = 0;
+ exponents(3, 1) = order;
+
+ 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++) {
+ exponents(index, 0) = exponents(p0, 0) + i*(exponents(p1,0)-exponents(p0,0))/order;
+ exponents(index, 1) = exponents(p0, 1) + i*(exponents(p1,1)-exponents(p0,1))/order;
+ }
+ }
+ if (order > 2) {
+ fullMatrix<double> inner = generateExponentsQuad(order - 2);
+ inner.add(1);
+ exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+ }
+ }
+ }
+
+ // exponents.print("expq");
+
+ return exponents;
+}
+static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
+
+static void nodepositionface0(int order, double *u, double *v, double *w)
+{ // uv surface - orientation v0-v2-v1
+ 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;
+ }
+}
+
+static void nodepositionface1(int order, double *u, double *v, double *w)
+{ // uw surface - orientation v0-v1-v3
+ 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;
+ }
+}
+
+static void nodepositionface2(int order, double *u, double *v, double *w)
+{ // vw surface - orientation v0-v3-v2
+ 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;
+ }
+}
+
+static void nodepositionface3(int order, double *u, double *v, double *w)
+{ // uvw surface - orientation v3-v1-v2
+ 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> generateExponentsTetrahedron(int order)
+{
+ int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
+
+ fullMatrix<double> exponents(nbPoints, 3);
+
+ exponents(0, 0) = 0;
+ exponents(0, 1) = 0;
+ exponents(0, 2) = 0;
+
+ if (order > 0) {
+ exponents(1, 0) = order;
+ exponents(1, 1) = 0;
+ exponents(1, 2) = 0;
+
+ exponents(2, 0) = 0;
+ exponents(2, 1) = order;
+ exponents(2, 2) = 0;
+
+ exponents(3, 0) = 0;
+ exponents(3, 1) = 0;
+ exponents(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++) {
+ exponents(4 + k, 0) = k + 1;
+ exponents(4 + order - 1 + k, 0) = order - 1 - k;
+ exponents(4 + 2 * (order - 1) + k, 0) = 0;
+ exponents(4 + 3 * (order - 1) + k, 0) = 0;
+ // exponents(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
+ // exponents(4 + 5 * (order - 1) + k, 0) = 0.;
+ exponents(4 + 4 * (order - 1) + k, 0) = 0;
+ exponents(4 + 5 * (order - 1) + k, 0) = k+1;
+
+ exponents(4 + k, 1) = 0;
+ exponents(4 + order - 1 + k, 1) = k + 1;
+ exponents(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ exponents(4 + 3 * (order - 1) + k, 1) = 0;
+ // exponents(4 + 4 * (order - 1) + k, 1) = 0.;
+ // exponents(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
+ exponents(4 + 4 * (order - 1) + k, 1) = k+1;
+ exponents(4 + 5 * (order - 1) + k, 1) = 0;
+
+ exponents(4 + k, 2) = 0;
+ exponents(4 + order - 1 + k, 2) = 0;
+ exponents(4 + 2 * (order - 1) + k, 2) = 0;
+ exponents(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+ exponents(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ exponents(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++){
+ exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+ exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+ exponents(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++){
+ exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+ exponents(ns + i, 1) = v[i] * (order - 3) ;
+ exponents(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++){
+ exponents(ns + i, 0) = u[i] * (order - 3);
+ exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+ exponents(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++){
+ exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+ exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+ exponents(ns + i, 2) = w[i] * (order - 3) + 1;
+ }
+
+ ns = ns + nbdofface;
+
+ delete [] u;
+ delete [] v;
+ delete [] w;
+
+ if (order > 3) {
+
+ fullMatrix<double> interior = generateExponentsTetrahedron(order - 4);
+ for (int k = 0; k < interior.size1(); k++) {
+ exponents(ns + k, 0) = 1 + interior(k, 0);
+ exponents(ns + k, 1) = 1 + interior(k, 1);
+ exponents(ns + k, 2) = 1 + interior(k, 2);
+ }
+ }
+ }
+ }
+ }
+
+ return exponents;
+}
+
+static fullMatrix<double> generateExponentsPrism(int order)
+{
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<double> exponents(nbPoints, 3);
+
+ int index = 0;
+ fullMatrix<double> triExp = generateExponentsTriangle(order);
+ fullMatrix<double> lineExp = generate1DExponents(order);
+ // First exponents (i < 3) must relate to corner
+ for (int i = 0; i < triExp.size1(); i++) {
+ exponents(index,0) = triExp(i,0);
+ exponents(index,1) = triExp(i,1);
+ exponents(index,2) = lineExp(0,0);
+ index ++;
+ exponents(index,0) = triExp(i,0);
+ exponents(index,1) = triExp(i,1);
+ exponents(index,2) = lineExp(1,0);
+ index ++;
+ }
+ for (int j = 2; j <lineExp.size1() ; j++) {
+ for (int i = 0; i < triExp.size1(); i++) {
+ exponents(index,0) = triExp(i,0);
+ exponents(index,1) = triExp(i,1);
+ exponents(index,2) = lineExp(j,0);
+ index ++;
+ }
+ }
+
+ return exponents;
+}
+
+static fullMatrix<double> generateExponentsHex(int order)
+{
+ int nbPoints = (order+1) * (order+1) * (order+1);
+ fullMatrix<double> exponents(nbPoints, 3);
+
+ if (order == 0) {
+ exponents(0, 0) = .0;
+ return exponents;
+ }
+
+ int index = 0;
+ fullMatrix<double> lineExp = generate1DExponents(order);
+
+ for (int i = 0; i < 2; ++i) {
+ 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);
+ ++index;
+ }
+ }
+ }
+
+ for (int i = 2; i < lineExp.size1(); ++i) {
+ 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);
+ ++index;
+ }
+ }
+ }
+ for (int i = 0; i < lineExp.size1(); ++i) {
+ 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);
+ ++index;
+ }
+ }
+ }
+ for (int i = 0; i < lineExp.size1(); ++i) {
+ 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);
+ ++index;
+ }
+ }
+ }
+
+ return exponents;
+}
+
+// Sampling Points
+static fullMatrix<double> generate1DPoints(int order)
+{
+ fullMatrix<double> line(order + 1, 1);
+ line(0,0) = .0;
+ if (order > 0) {
+ line(1, 0) = 1.;
+ double dd = 1. / order;
+ for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
+ }
+
+ return line;
+}
+
+static fullMatrix<double> generatePointsQuadRecur(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;
+ }
+ }
+ if (order >= 2 && !serendip) {
+ fullMatrix<double> inner = generatePointsQuadRecur(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> generatePointsQuad(int order, bool serendip)
+{
+ fullMatrix<double> point = generatePointsQuadRecur(order, serendip);
+ // rescale to [0,1] x [0,1]
+ for (int i=0;i<point.size1();i++){
+ point(i,0) = (1.+point(i,0))*.5;
+ point(i,1) = (1.+point(i,1))*.5;
+ }
+ return point;
+}
+
+static fullMatrix<double> generatePointsPrism(int order, bool serendip)
+{
+ const double prism18Pts[18][3] = {
+ {0, 0, 0}, // 0
+ {1, 0, 0}, // 1
+ {0, 1, 0}, // 2
+ {0, 0, 1}, // 3
+ {1, 0, 1}, // 4
+ {0, 1, 1}, // 5
+ {.5, 0, 0}, // 6
+ {0, .5, 0}, // 7
+ {0, 0, .5}, // 8
+ {.5, .5, 0}, // 9
+ {1, 0, .5}, // 10
+ {0, 1, .5}, // 11
+ {.5, 0, 1}, // 12
+ {0, .5, 1}, // 13
+ {.5, .5, 1}, // 14
+ {.5, 0, .5}, // 15
+ {0, .5, .5}, // 16
+ {.5, .5, .5}, // 17
+ };
+
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ fullMatrix<double> point(nbPoints, 3);
+
+ int index = 0;
+
+ if (order == 2)
+ for (int i =0; i<18; i++)
+ for (int j=0; j<3;j++)
+ point(i,j) = prism18Pts[i][j];
+ else {
+ fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
+ fullMatrix<double> linePoints = generate1DPoints(order);
+ 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> generatePointsHex(int order, bool serendip)
+{
+ int nbPoints = (order+1) * (order+1) * (order+1);
+ fullMatrix<double> point(nbPoints, 3);
+
+ fullMatrix<double> linePoints = generate1DPoints(order);
+ int index = 0;
+
+ for (int i = 0; i < linePoints.size1(); ++i) {
+ for (int j = 0; j < linePoints.size1(); ++j) {
+ for (int k = 0; k < linePoints.size1(); ++k) {
+ point(index, 0) = linePoints(i, 0);
+ point(index, 1) = linePoints(j, 0);
+ point(index, 2) = linePoints(k, 0);
+ ++index;
+ }
+ }
+ }
+
+ return point;
+}
+
+// Sub Control Points
+static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(2);
+
+ subPoints[0] = generate1DPoints(order);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ subPoints[1].add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
+
+ subPoints[0] = gmshGeneratePointsTriangle(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ subPoints[3].scale(-1.);
+ subPoints[3].add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generatePointsQuad(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox1;
+ fullMatrix<double> prox2;
+
+ subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[1], 0, 1);
+ prox1.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[2], 1, 1);
+ prox1.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[3], 2, 1);
+ prox1.add(.5);
+
+ // u := .5-u-w
+ // v := .5-v-w
+ // w := w
+ subPoints[4].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[4], 0, 2);
+ prox1.scale(-1.);
+ prox1.add(.5);
+ prox1.setAsProxy(subPoints[4], 0, 1);
+ prox2.setAsProxy(subPoints[4], 2, 1);
+ prox1.add(prox2, -1.);
+ prox1.setAsProxy(subPoints[4], 1, 1);
+ prox1.add(prox2, -1.);
+
+ // u := u
+ // v := .5-v
+ // w := w+v
+ subPoints[5].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[5], 2, 1);
+ prox2.setAsProxy(subPoints[5], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := .5-u
+ // v := v
+ // w := w+u
+ subPoints[6].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[6], 2, 1);
+ prox2.setAsProxy(subPoints[6], 0, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := u+w
+ // v := v+w
+ // w := .5-w
+ subPoints[7].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[7], 0, 1);
+ prox2.setAsProxy(subPoints[7], 2, 1);
+ prox1.add(prox2);
+ prox1.setAsProxy(subPoints[7], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generatePointsPrism(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[3], 0, 2);
+ prox.scale(-1.);
+ prox.add(.5);
+
+ subPoints[4].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[4], 2, 1);
+ prox.add(.5);
+
+ subPoints[5].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[5], 2, 1);
+ prox.add(.5);
+
+ subPoints[6].copy(subPoints[2]);
+ prox.setAsProxy(subPoints[6], 2, 1);
+ prox.add(.5);
+
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
+
+ subPoints[0] = generatePointsHex(order, false);
+ subPoints[0].scale(.5);
+
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
+
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
+
+ subPoints[4].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[4], 2, 1);
+ prox.add(.5);
+
+ subPoints[5].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[5], 2, 1);
+ prox.add(.5);
+
+ subPoints[6].copy(subPoints[2]);
+ prox.setAsProxy(subPoints[6], 2, 1);
+ prox.add(.5);
+
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+// Matrices generation
+static int nChoosek(int n, int k)
+{
+ if (n < k || k < 0) {
+ Msg::Error("Wrong argument for combination.");
+ return 1;
+ }
+
+ if (k > n/2) k = n-k;
+ if (k == 1)
+ return n;
+ if (k == 0)
+ return 1;
+
+ int c = 1;
+ for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+ return c;
+}
+
+
+static fullMatrix<double> generateBez2LagMatrix
+ (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+ int order, int dimSimplex)
+{
+
+ if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
+ exponent.size1(),point.size1(),
+ exponent.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
+
+ int ndofs = exponent.size1();
+ int dim = exponent.size2();
+
+ fullMatrix<double> bez2Lag(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+
+ {
+ double pointCompl = 1.;
+ int exponentCompl = order;
+ for (int k = 0; k < dimSimplex; k++) {
+ dd *= nChoosek(exponentCompl, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k));
+ pointCompl -= point(j, k);
+ exponentCompl -= (int) exponent(i, k);
+ }
+ dd *= pow(pointCompl, exponentCompl);
+ }
+
+ for (int k = dimSimplex; k < dim; k++)
+ dd *= nChoosek(order, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k))
+ * pow(1. - point(j, k), order - exponent(i, k));
+
+ bez2Lag(j, i) = dd;
+ }
+ }
+ return bez2Lag;
+}
+
+
+
+static fullMatrix<double> generateSubDivisor
+ (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
+{
+ if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exponents.size1(), lag2Bez.size1(),
+ exponents.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
+
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
+
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> subDivisor(nbSubPts, nbPts);
+
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
+ lag2Bez.mult(intermediate1, intermediate2);
+ subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+ }
+ return subDivisor;
+}
+
+
+
+// Convert Bezier points to Lagrange points
+static fullMatrix<double> bez2LagPoints(bool dim1, bool dim2, bool dim3, const fullMatrix<double> &bezierPoints)
+{
+
+ // FIXME BUG for 2nd order quads we seem to try to access bezierPoints(i, 2);
+ // but bezierPoints only has 2 columns
+
+ const int numPt = bezierPoints.size1();
+ fullMatrix<double> lagPoints(numPt,3);
+
+ for (int i=0; i<numPt; i++) {
+ lagPoints(i,0) = dim1 ? -1. + 2*bezierPoints(i,0) : bezierPoints(i,0);
+ lagPoints(i,1) = dim2 ? -1. + 2*bezierPoints(i,1) : bezierPoints(i,1);
+ lagPoints(i,2) = dim3 ? -1. + 2*bezierPoints(i,2) : bezierPoints(i,2);
+ }
+
+ return lagPoints;
+
+}
+
+
+bezierBasis::bezierBasis(int tag)
+{
+ order = MElement::OrderFromTag(tag);
+
+ if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+ dim = 3;
+ numLagPts = 5;
+ bezierPoints = gmshGeneratePointsPyramid(order,false);
+ lagPoints = bezierPoints;
+ matrixLag2Bez.resize(bezierPoints.size1(), bezierPoints.size1(), false);
+ matrixBez2Lag.resize(bezierPoints.size1(), bezierPoints.size1(), false);
+ for (int i=0; i<matrixBez2Lag.size1(); ++i) {
+ matrixLag2Bez(i,i) = 1.;
+ matrixBez2Lag(i,i) = 1.;
+ }
+ // TODO: Implement subdidivsor
+ }
+ else {
+ int dimSimplex;
+ std::vector< fullMatrix<double> > subPoints;
+ switch (MElement::ParentTypeFromTag(tag)) {
+ case TYPE_PNT :
+ dim = 0;
+ numLagPts = 1;
+ exponents = generate1DExponents(0);
+ bezierPoints = generate1DPoints(0);
+ lagPoints = bezierPoints;
+ dimSimplex = 0;
+ break;
+ case TYPE_LIN : {
+ dim = 1;
+ numLagPts = 2;
+ exponents = generate1DExponents(order);
+ bezierPoints = generate1DPoints(order);
+ lagPoints = bez2LagPoints(true,false,false,bezierPoints);
+ dimSimplex = 0;
+ subPoints = generateSubPointsLine(0);
+ break;
+ }
+ case TYPE_TRI : {
+ dim = 2;
+ numLagPts = 3;
+ exponents = generateExponentsTriangle(order);
+ bezierPoints = gmshGeneratePointsTriangle(order,false);
+ lagPoints = bezierPoints;
+ dimSimplex = 2;
+ subPoints = generateSubPointsTriangle(order);
+ break;
+ }
+ case TYPE_QUA : {
+ dim = 2;
+ numLagPts = 4;
+ exponents = generateExponentsQuad(order);
+ bezierPoints = generatePointsQuad(order,false);
+ lagPoints = bez2LagPoints(true,true,false,bezierPoints);
+ dimSimplex = 0;
+ subPoints = generateSubPointsQuad(order);
+ // points.print("points");
+ break;
+ }
+ case TYPE_TET : {
+ dim = 3;
+ numLagPts = 4;
+ exponents = generateExponentsTetrahedron(order);
+ bezierPoints = gmshGeneratePointsTetrahedron(order,false);
+ lagPoints = bezierPoints;
+ dimSimplex = 3;
+ subPoints = generateSubPointsTetrahedron(order);
+ break;
+ }
+ case TYPE_PRI : {
+ dim = 3;
+ numLagPts = 6;
+ exponents = generateExponentsPrism(order);
+ bezierPoints = generatePointsPrism(order, false);
+ lagPoints = bez2LagPoints(false,false,true,bezierPoints);
+ dimSimplex = 2;
+ subPoints = generateSubPointsPrism(order);
+ break;
+ }
+ case TYPE_HEX : {
+ dim = 3;
+ numLagPts = 8;
+ exponents = generateExponentsHex(order);
+ bezierPoints = generatePointsHex(order, false);
+ lagPoints = bez2LagPoints(true,true,true,bezierPoints);
+ dimSimplex = 0;
+ subPoints = generateSubPointsHex(order);
+ break;
+ }
+ default : {
+ Msg::Error("Unknown function space %d: reverting to TET_1", tag);
+ dim = 3;
+ numLagPts = 4;
+ exponents = generateExponentsTetrahedron(0);
+ bezierPoints = gmshGeneratePointsTetrahedron(0, false);
+ lagPoints = bezierPoints;
+ dimSimplex = 3;
+ subPoints = generateSubPointsTetrahedron(0);
+ break;
+ }
+ }
+ Msg::Info("%d %d %d %d %d %d", exponents.size1(), exponents.size2(), bezierPoints.size1(), bezierPoints.size2(), order, dimSimplex);
+ matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+ Msg::Info("%d %d %d %d", matrixBez2Lag.size1(), matrixBez2Lag.size2(), matrixLag2Bez.size1(), matrixLag2Bez.size2());
+ Msg::Info("%f %f %f %f", bezierPoints(0, 0), bezierPoints(0, 1), exponents(0, 0), exponents(0, 1));
+ matrixBez2Lag.invert(matrixLag2Bez);
+ subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
+ numDivisions = (int) pow(2., (int) bezierPoints.size2());
+ }
+}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
new file mode 100644
index 0000000000..e8241ff6b7
--- /dev/null
+++ b/Numeric/bezierBasis.h
@@ -0,0 +1,30 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _BEZIER_BASIS_H_
+#define _BEZIER_BASIS_H_
+
+#include <map>
+#include <vector>
+#include "fullMatrix.h"
+
+class bezierBasis {
+ private :
+ public :
+ int dim, order;
+ int numLagPts;
+ int numDivisions;
+ // the 'numLagPts' first exponents are related to 'Lagrangian' values
+ bezierBasis() {Msg::Fatal("Should not be created this way");};
+ bezierBasis(int tag);
+ fullMatrix<double> exponents; // exponents of Bezier FF
+ fullMatrix<double> bezierPoints, lagPoints; // sampling point
+ fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
+ fullMatrix<double> subDivisor;
+};
+
+#endif
+
+
diff --git a/Numeric/jacobiPolynomials.h b/Numeric/jacobiPolynomials.h
index 9d83f5077f..b3d510c341 100644
--- a/Numeric/jacobiPolynomials.h
+++ b/Numeric/jacobiPolynomials.h
@@ -11,7 +11,7 @@
class JacobiPolynomials {
public:
- JacobiPolynomials() {};
+ JacobiPolynomials() : n(0) {};
JacobiPolynomials(double a, double b, int o);
~JacobiPolynomials();
diff --git a/Numeric/legendrePolynomials.h b/Numeric/legendrePolynomials.h
index f7fd118ddf..381e62e249 100644
--- a/Numeric/legendrePolynomials.h
+++ b/Numeric/legendrePolynomials.h
@@ -11,7 +11,7 @@
class LegendrePolynomials {
public:
- LegendrePolynomials() {};
+ LegendrePolynomials() : n(0) {};
LegendrePolynomials(int o);
~LegendrePolynomials();
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 81675224b7..38126f89c9 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -1141,7 +1141,7 @@ 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));
+ const nodalBasis &fsFace = *BasisFactory::getNodalBasis(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++) {
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 4516507d35..0c09ef633c 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -30,8 +30,7 @@ inline double pow_int(const double &a, const int &n)
case 5 :
{
const double a2 = a*a;
- const double a3 = a*a*a;
- return a2*a3;
+ return a2*a2*a;
}
case 6 :
{
@@ -51,9 +50,8 @@ inline double pow_int(const double &a, const int &n)
}
case 9 :
{
- const double a2 = a*a;
- const double a4 = a2*a2;
- return a4*a4*a;
+ const double a3 = a*a*a;
+ return a3*a3*a3;
}
case 10 :
{
@@ -93,7 +91,7 @@ class polynomialBasis : public nodalBasis
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void ddf(double u, double v, double w, double hess[][3][3]) const;
- virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
+ virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
virtual int getNumShapeFunctions() const;
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 05e5653bde..1e4f5438ce 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -6,7 +6,7 @@
#include "Gmsh.h"
#include "GModel.h"
#include "GmshMessage.h"
-#include "JacobianBasis.h"
+
#include "polynomialBasis.h"
#include "AnalyseCurvedMesh.h"
#include "Context.h"
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 01feaae979..211bb22467 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -903,7 +903,7 @@ void PViewDataList::setOrder2(int type)
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
// case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
- const polynomialBasis *fs = (polynomialBasis*)BasisFactory::create(typeMSH);
+ const polynomialBasis *fs = (polynomialBasis*)BasisFactory::getNodalBasis(typeMSH);
if(!fs){
Msg::Error("Could not find polynomial function space for element type %d",
typeMSH);
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 478fd0cace..3d2dbaba46 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -1045,7 +1045,7 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
const nodalBasis* DI_Line::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_LIN, order);
- return BasisFactory::create(tag);
+ return BasisFactory::getNodalBasis(tag);
}
void DI_Line::computeIntegral() {
@@ -1066,7 +1066,7 @@ const nodalBasis* DI_Triangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_TRI, order);
- return BasisFactory::create(tag);
+ return BasisFactory::getNodalBasis(tag);
}
void DI_Triangle::computeIntegral() {
integral_ = TriSurf(pt(0), pt(1), pt(2));
@@ -1090,7 +1090,7 @@ double DI_Triangle::quality() const {
const nodalBasis* DI_Quad::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_QUA, order);
- return BasisFactory::create(tag);
+ return BasisFactory::getNodalBasis(tag);
}
void DI_Quad::computeIntegral() {
@@ -1114,7 +1114,7 @@ void DI_Quad::computeIntegral() {
const nodalBasis* DI_Tetra::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_TET, order);
- return BasisFactory::create(tag);
+ return BasisFactory::getNodalBasis(tag);
}
void DI_Tetra::computeIntegral() {
@@ -1129,7 +1129,7 @@ double DI_Tetra::quality() const {
const nodalBasis* DI_Hexa::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
int tag = polynomialBasis::getTag(TYPE_HEX, order);
- return BasisFactory::create(tag);
+ return BasisFactory::getNodalBasis(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 c791cfc1fb..132405c539 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -49,7 +49,7 @@ class DI_Point
const std::vector<gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
virtual ~DI_Point(){}
virtual const nodalBasis* getFunctionSpace(int o) const
- { return BasisFactory::create(MSH_PNT); }
+ { return BasisFactory::getNodalBasis(MSH_PNT); }
virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
s[0] = 1.;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 481d6a2b8e..e8f526c165 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -8,7 +8,7 @@
#include "MVertex.h"
#include "ParamCoord.h"
#include "polynomialBasis.h"
-#include "JacobianBasis.h"
+
--
GitLab