diff --git a/FunctionSpace/LineLagrangeBasis.cpp b/FunctionSpace/LineLagrangeBasis.cpp
index c9ab665b93f28cd9ab80b6ab638d188075dd57b8..e6023784c3c946bbb85bbeb8a83fae22efd1ae4c 100644
--- a/FunctionSpace/LineLagrangeBasis.cpp
+++ b/FunctionSpace/LineLagrangeBasis.cpp
@@ -33,7 +33,7 @@ LineLagrangeBasis::~LineLagrangeBasis(void){
}
unsigned int LineLagrangeBasis::getTag(unsigned int order){
- unsigned int tag = MElement::getTag(TYPE_LIN, order, false);
+ unsigned int tag = ElementType::getTag(TYPE_LIN, order, false);
if(tag)
return tag;
diff --git a/FunctionSpace/QuadLagrangeBasis.cpp b/FunctionSpace/QuadLagrangeBasis.cpp
index 78c5cfa220825a5f2864b169e659a87b9b2d3149..9e8c4ab3fc13c3eaf7c91a0cde26e82ea2997ffc 100644
--- a/FunctionSpace/QuadLagrangeBasis.cpp
+++ b/FunctionSpace/QuadLagrangeBasis.cpp
@@ -33,7 +33,7 @@ QuadLagrangeBasis::~QuadLagrangeBasis(void){
}
unsigned int QuadLagrangeBasis::getTag(unsigned int order){
- unsigned int tag = MElement::getTag(TYPE_QUA, order, false);
+ unsigned int tag = ElementType::getTag(TYPE_QUA, order, false);
if(tag)
return tag;
diff --git a/FunctionSpace/TetLagrangeBasis.cpp b/FunctionSpace/TetLagrangeBasis.cpp
index b4cae893d67be596f2aa31ea14426387b6e39ed7..0c6f3eeb4825b6185f7d2b893ebead696277f733 100644
--- a/FunctionSpace/TetLagrangeBasis.cpp
+++ b/FunctionSpace/TetLagrangeBasis.cpp
@@ -33,7 +33,7 @@ TetLagrangeBasis::~TetLagrangeBasis(void){
}
unsigned int TetLagrangeBasis::getTag(unsigned int order){
- unsigned int tag = MElement::getTag(TYPE_TET, order, false);
+ unsigned int tag = ElementType::getTag(TYPE_TET, order, false);
if(tag)
return tag;
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index bf41e9765a2504bc1c7607a7796f2606008bdbc9..ccf6db229795f3755df4f4d767536ce9f7998d0a 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -41,7 +41,7 @@ TriLagrangeBasis::~TriLagrangeBasis(void){
}
unsigned int TriLagrangeBasis::getTag(unsigned int order){
- unsigned int tag = MElement::getTag(TYPE_TRI, order, false);
+ unsigned int tag = ElementType::getTag(TYPE_TRI, order, false);
if(tag)
return tag;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index abee8693e7906c39c60d77db40c0d1045ca38c19..759f3bb38219c6c10594e8a98a0fa61a23e5cef9 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1384,544 +1384,6 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
return newEl;
}
-// Gives the parent type corresponding to
-// any element type.
-// Ex. : MSH_TRI_3 -> TYPE_TRI
-int MElement::ParentTypeFromTag(int tag)
-{
- switch(tag) {
- case(MSH_PNT):
- return TYPE_PNT;
- case(MSH_LIN_2): case(MSH_LIN_3):
- case(MSH_LIN_4): case(MSH_LIN_5):
- case(MSH_LIN_6): case(MSH_LIN_7):
- case(MSH_LIN_8): case(MSH_LIN_9):
- case(MSH_LIN_10): case(MSH_LIN_11):
- case(MSH_LIN_B): case(MSH_LIN_C):
- case(MSH_LIN_1):
- return TYPE_LIN;
- case(MSH_TRI_3): case(MSH_TRI_6):
- case(MSH_TRI_9): case(MSH_TRI_10):
- case(MSH_TRI_12): case(MSH_TRI_15):
- case(MSH_TRI_15I): case(MSH_TRI_21):
- case(MSH_TRI_28): case(MSH_TRI_36):
- case(MSH_TRI_45): case(MSH_TRI_55):
- case(MSH_TRI_66): case(MSH_TRI_18):
- case(MSH_TRI_21I): case(MSH_TRI_24):
- case(MSH_TRI_27): case(MSH_TRI_30):
- case(MSH_TRI_B): case(MSH_TRI_1):
- return TYPE_TRI;
- case(MSH_QUA_4): case(MSH_QUA_9):
- case(MSH_QUA_8): case(MSH_QUA_16):
- case(MSH_QUA_25): case(MSH_QUA_36):
- case(MSH_QUA_12): case(MSH_QUA_16I):
- case(MSH_QUA_20): case(MSH_QUA_49):
- case(MSH_QUA_64): case(MSH_QUA_81):
- case(MSH_QUA_100): case(MSH_QUA_121):
- case(MSH_QUA_24): case(MSH_QUA_28):
- case(MSH_QUA_32): case(MSH_QUA_36I):
- case(MSH_QUA_40): case(MSH_QUA_1):
- return TYPE_QUA;
- case(MSH_TET_4): case(MSH_TET_10):
- case(MSH_TET_20): case(MSH_TET_35):
- case(MSH_TET_56): case(MSH_TET_34):
- case(MSH_TET_52): case(MSH_TET_84):
- case(MSH_TET_120): case(MSH_TET_165):
- case(MSH_TET_220): case(MSH_TET_286):
- case(MSH_TET_74): case(MSH_TET_100):
- case(MSH_TET_130): case(MSH_TET_164):
- case(MSH_TET_202): case(MSH_TET_1):
- return TYPE_TET;
- case(MSH_PYR_5): case(MSH_PYR_14):
- case(MSH_PYR_13): case(MSH_PYR_30):
- case(MSH_PYR_55): case(MSH_PYR_91):
- case(MSH_PYR_140): case(MSH_PYR_204):
- case(MSH_PYR_285): case(MSH_PYR_385):
- case(MSH_PYR_21): case(MSH_PYR_29):
- case(MSH_PYR_37): case(MSH_PYR_45):
- case(MSH_PYR_53): case(MSH_PYR_61):
- case(MSH_PYR_69): case(MSH_PYR_1):
- return TYPE_PYR;
- case(MSH_PRI_6): case(MSH_PRI_18):
- case(MSH_PRI_15): case(MSH_PRI_1):
- case(MSH_PRI_40): case(MSH_PRI_75):
- case(MSH_PRI_126): case(MSH_PRI_196):
- case(MSH_PRI_288): case(MSH_PRI_405):
- case(MSH_PRI_550): case(MSH_PRI_24):
- case(MSH_PRI_33): case(MSH_PRI_42):
- case(MSH_PRI_51): case(MSH_PRI_60):
- case(MSH_PRI_69): case(MSH_PRI_78):
- return TYPE_PRI;
- case(MSH_HEX_8): case(MSH_HEX_27):
- case(MSH_HEX_20): case(MSH_HEX_1):
- case(MSH_HEX_64): case(MSH_HEX_125):
- case(MSH_HEX_216): case(MSH_HEX_343):
- case(MSH_HEX_512): case(MSH_HEX_729):
- case(MSH_HEX_1000): case(MSH_HEX_32):
- case(MSH_HEX_44): case(MSH_HEX_56):
- case(MSH_HEX_68): case(MSH_HEX_80):
- case(MSH_HEX_92): case(MSH_HEX_104):
- return TYPE_HEX;
- case(MSH_POLYG_): case(MSH_POLYG_B):
- return TYPE_POLYG;
- case(MSH_POLYH_):
- return TYPE_POLYH;
- case(MSH_PNT_SUB): case(MSH_LIN_SUB):
- case(MSH_TRI_SUB): case(MSH_TET_SUB):
- return TYPE_XFEM;
- default:
- Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
- return TYPE_TET;
- }
-}
-
-// Gives the order corresponding to any element type.
-int MElement::OrderFromTag(int tag)
-{
-
- switch (tag) {
- case MSH_PNT : return 0;
- case MSH_LIN_1 : return 0;
- case MSH_LIN_2 : return 1;
- case MSH_LIN_3 : return 2;
- case MSH_LIN_4 : return 3;
- case MSH_LIN_5 : return 4;
- case MSH_LIN_6 : return 5;
- case MSH_LIN_7 : return 6;
- case MSH_LIN_8 : return 7;
- case MSH_LIN_9 : return 8;
- case MSH_LIN_10 : return 9;
- case MSH_LIN_11 : return 10;
- case MSH_TRI_1 : return 0;
- case MSH_TRI_3 : return 1;
- case MSH_TRI_6 : return 2;
- case MSH_TRI_10 : return 3;
- case MSH_TRI_15 : return 4;
- case MSH_TRI_21 : return 5;
- case MSH_TRI_28 : return 6;
- case MSH_TRI_36 : return 7;
- case MSH_TRI_45 : return 8;
- case MSH_TRI_55 : return 9;
- case MSH_TRI_66 : return 10;
- case MSH_TRI_9 : return 3;
- case MSH_TRI_12 : return 4;
- case MSH_TRI_15I : return 5;
- case MSH_TRI_18 : return 6;
- case MSH_TRI_21I : return 7;
- case MSH_TRI_24 : return 8;
- case MSH_TRI_27 : return 9;
- case MSH_TRI_30 : return 10;
- case MSH_TET_1 : return 0;
- case MSH_TET_4 : return 1;
- case MSH_TET_10 : return 2;
- case MSH_TET_20 : return 3;
- case MSH_TET_35 : return 4;
- case MSH_TET_56 : return 5;
- case MSH_TET_84 : return 6;
- case MSH_TET_120 : return 7;
- case MSH_TET_165 : return 8;
- case MSH_TET_220 : return 9;
- case MSH_TET_286 : return 10;
- case MSH_TET_34 : return 4;
- case MSH_TET_52 : return 5;
- case MSH_TET_74 : return 6;
- case MSH_TET_100 : return 7;
- case MSH_TET_130 : return 8;
- case MSH_TET_164 : return 9;
- case MSH_TET_202 : return 10;
- case MSH_QUA_1 : return 0;
- case MSH_QUA_4 : return 1;
- case MSH_QUA_9 : return 2;
- case MSH_QUA_16 : return 3;
- case MSH_QUA_25 : return 4;
- case MSH_QUA_36 : return 5;
- case MSH_QUA_49 : return 6;
- case MSH_QUA_64 : return 7;
- case MSH_QUA_81 : return 8;
- case MSH_QUA_100 : return 9;
- case MSH_QUA_121 : return 10;
- case MSH_QUA_8 : return 2;
- case MSH_QUA_12 : return 3;
- case MSH_QUA_16I : return 4;
- case MSH_QUA_20 : return 5;
- case MSH_QUA_24 : return 6;
- case MSH_QUA_28 : return 7;
- case MSH_QUA_32 : return 8;
- case MSH_QUA_36I : return 9;
- case MSH_QUA_40 : return 10;
- case MSH_PRI_1 : return 0;
- case MSH_PRI_6 : return 1;
- case MSH_PRI_18 : return 2;
- case MSH_PRI_40 : return 3;
- case MSH_PRI_75 : return 4;
- case MSH_PRI_126 : return 5;
- case MSH_PRI_196 : return 6;
- case MSH_PRI_288 : return 7;
- case MSH_PRI_405 : return 8;
- case MSH_PRI_550 : return 9;
- case MSH_PRI_15 : return 2;
- case MSH_PRI_24 : return 3;
- case MSH_PRI_33 : return 4;
- case MSH_PRI_42 : return 5;
- case MSH_PRI_51 : return 6;
- case MSH_PRI_60 : return 7;
- case MSH_PRI_69 : return 8;
- case MSH_PRI_78 : return 9;
- case MSH_HEX_1 : return 0;
- case MSH_HEX_8 : return 1;
- case MSH_HEX_27 : return 2;
- case MSH_HEX_64 : return 3;
- case MSH_HEX_125 : return 4;
- case MSH_HEX_216 : return 5;
- case MSH_HEX_343 : return 6;
- case MSH_HEX_512 : return 7;
- case MSH_HEX_729 : return 8;
- case MSH_HEX_1000: return 9;
- case MSH_HEX_20 : return 2;
- case MSH_HEX_32 : return 3;
- case MSH_HEX_44 : return 4;
- case MSH_HEX_56 : return 5;
- case MSH_HEX_68 : return 6;
- case MSH_HEX_80 : return 7;
- case MSH_HEX_92 : return 8;
- case MSH_HEX_104 : return 9;
- case MSH_PYR_1 : return 0;
- case MSH_PYR_5 : return 1;
- case MSH_PYR_14 : return 2;
- case MSH_PYR_30 : return 3;
- case MSH_PYR_55 : return 4;
- case MSH_PYR_91 : return 5;
- case MSH_PYR_140 : return 6;
- case MSH_PYR_204 : return 7;
- case MSH_PYR_285 : return 8;
- case MSH_PYR_385 : return 9;
- case MSH_PYR_13 : return 2;
- case MSH_PYR_21 : return 3;
- case MSH_PYR_29 : return 4;
- case MSH_PYR_37 : return 5;
- case MSH_PYR_45 : return 6;
- case MSH_PYR_53 : return 7;
- case MSH_PYR_61 : return 8;
- case MSH_PYR_69 : return 9;
- default :
- Msg::Error("Unknown element type %d: reverting to order 1",tag);
- return 1;
- }
-
-}
-
-// Gives > 0 if element type is in Serendipity Family.
-// Gives < 2 if element type is in Not Serendipity Family.
-int MElement::SerendipityFromTag(int tag)
-{
- switch (tag) {
- case MSH_PNT : case MSH_LIN_1 :
- case MSH_LIN_2 : case MSH_LIN_3 :
- case MSH_LIN_4 : case MSH_LIN_5 :
- case MSH_LIN_6 : case MSH_LIN_7 :
- case MSH_LIN_8 : case MSH_LIN_9 :
- case MSH_LIN_10 : case MSH_LIN_11 :
-
- case MSH_TRI_1 : case MSH_TRI_3 :
- case MSH_TRI_6 :
-
- case MSH_QUA_1 : case MSH_QUA_4 :
-
- case MSH_TET_1 : case MSH_TET_4 :
- case MSH_TET_10 : case MSH_TET_20 :
-
- case MSH_PRI_1 : case MSH_PRI_6 :
-
- case MSH_HEX_1 : case MSH_HEX_8 :
-
- case MSH_PYR_1 : case MSH_PYR_5 :
-
- return 1; // Serendipity or not
-
-
- case MSH_TRI_10 : case MSH_TRI_15 :
- case MSH_TRI_21 : case MSH_TRI_28 :
- case MSH_TRI_36 : case MSH_TRI_45 :
- case MSH_TRI_55 : case MSH_TRI_66 :
-
- case MSH_QUA_9 : case MSH_QUA_16 :
- case MSH_QUA_25 : case MSH_QUA_36 :
- case MSH_QUA_49 : case MSH_QUA_64 :
- case MSH_QUA_81 : case MSH_QUA_100 :
- case MSH_QUA_121 :
-
- case MSH_TET_35 : case MSH_TET_56 :
- case MSH_TET_84 : case MSH_TET_120 :
- case MSH_TET_165 : case MSH_TET_220 :
- case MSH_TET_286 :
-
- case MSH_PRI_18 : case MSH_PRI_40 :
- case MSH_PRI_75 : case MSH_PRI_126 :
- case MSH_PRI_196 : case MSH_PRI_288 :
- case MSH_PRI_405 : case MSH_PRI_550 :
-
- case MSH_HEX_27 : case MSH_HEX_64 :
- case MSH_HEX_125 : case MSH_HEX_216 :
- case MSH_HEX_343 : case MSH_HEX_512 :
- case MSH_HEX_729 : case MSH_HEX_1000:
-
- case MSH_PYR_14 : case MSH_PYR_30 :
- case MSH_PYR_55 : case MSH_PYR_91 :
- case MSH_PYR_140 : case MSH_PYR_204 :
- case MSH_PYR_285 : case MSH_PYR_385 :
-
- return 0; // Not Serendipity
-
-
- case MSH_TRI_9 : case MSH_TRI_12 :
- case MSH_TRI_15I : case MSH_TRI_18 :
- case MSH_TRI_21I : case MSH_TRI_24 :
- case MSH_TRI_27 : case MSH_TRI_30 :
-
- case MSH_QUA_8 : case MSH_QUA_12 :
- case MSH_QUA_16I : case MSH_QUA_20 :
- case MSH_QUA_24 : case MSH_QUA_28 :
- case MSH_QUA_32 : case MSH_QUA_36I :
- case MSH_QUA_40 :
-
- case MSH_TET_34 : case MSH_TET_52 :
- case MSH_TET_74 : case MSH_TET_100 :
- case MSH_TET_130 : case MSH_TET_164 :
- case MSH_TET_202 :
-
- case MSH_PRI_15 : case MSH_PRI_24 :
- case MSH_PRI_33 : case MSH_PRI_42 :
- case MSH_PRI_51 : case MSH_PRI_60 :
- case MSH_PRI_69 : case MSH_PRI_78 :
-
- case MSH_HEX_20 : case MSH_HEX_32 :
- case MSH_HEX_44 : case MSH_HEX_56 :
- case MSH_HEX_68 : case MSH_HEX_80 :
- case MSH_HEX_92 : case MSH_HEX_104 :
-
- case MSH_PYR_13 : case MSH_PYR_21 :
- case MSH_PYR_29 : case MSH_PYR_37 :
- case MSH_PYR_45 : case MSH_PYR_53 :
- case MSH_PYR_61 : case MSH_PYR_69 :
-
- return 2; // Only Serendipity
-
- default :
- Msg::Error("Unknown element type %d: assuming not serendipity",tag);
- return 0;
- }
-}
-
-// Gives the dimension corresponding to any element type.
-int MElement::DimensionFromTag(int tag)
-{
- switch(tag) {
- case(MSH_PNT): case(MSH_PNT_SUB):
- return 0;
-
- case(MSH_LIN_2): case(MSH_LIN_3):
- case(MSH_LIN_4): case(MSH_LIN_5):
- case(MSH_LIN_6): case(MSH_LIN_7):
- case(MSH_LIN_8): case(MSH_LIN_9):
- case(MSH_LIN_10): case(MSH_LIN_11):
- case(MSH_LIN_B): case(MSH_LIN_C):
- case(MSH_LIN_1): case(MSH_LIN_SUB):
- return 1;
-
- case(MSH_TRI_3): case(MSH_TRI_6):
- case(MSH_TRI_9): case(MSH_TRI_10):
- case(MSH_TRI_12): case(MSH_TRI_15):
- case(MSH_TRI_15I): case(MSH_TRI_21):
- case(MSH_TRI_28): case(MSH_TRI_36):
- case(MSH_TRI_45): case(MSH_TRI_55):
- case(MSH_TRI_66): case(MSH_TRI_18):
- case(MSH_TRI_21I): case(MSH_TRI_24):
- case(MSH_TRI_27): case(MSH_TRI_30):
- case(MSH_TRI_B): case(MSH_TRI_1):
- case(MSH_TRI_SUB):
-
- case(MSH_QUA_4): case(MSH_QUA_9):
- case(MSH_QUA_8): case(MSH_QUA_16):
- case(MSH_QUA_25): case(MSH_QUA_36):
- case(MSH_QUA_12): case(MSH_QUA_16I):
- case(MSH_QUA_20): case(MSH_QUA_49):
- case(MSH_QUA_64): case(MSH_QUA_81):
- case(MSH_QUA_100): case(MSH_QUA_121):
- case(MSH_QUA_24): case(MSH_QUA_28):
- case(MSH_QUA_32): case(MSH_QUA_36I):
- case(MSH_QUA_40): case(MSH_QUA_1):
-
- case(MSH_POLYG_): case(MSH_POLYG_B):
- return 2;
-
- case(MSH_TET_4): case(MSH_TET_10):
- case(MSH_TET_20): case(MSH_TET_35):
- case(MSH_TET_56): case(MSH_TET_34):
- case(MSH_TET_52): case(MSH_TET_84):
- case(MSH_TET_120): case(MSH_TET_165):
- case(MSH_TET_220): case(MSH_TET_286):
- case(MSH_TET_74): case(MSH_TET_100):
- case(MSH_TET_130): case(MSH_TET_164):
- case(MSH_TET_202): case(MSH_TET_1):
- case(MSH_TET_SUB):
-
- case(MSH_PYR_5): case(MSH_PYR_14):
- case(MSH_PYR_13): case(MSH_PYR_30):
- case(MSH_PYR_55): case(MSH_PYR_91):
- case(MSH_PYR_140): case(MSH_PYR_204):
- case(MSH_PYR_285): case(MSH_PYR_385):
- case(MSH_PYR_21): case(MSH_PYR_29):
- case(MSH_PYR_37): case(MSH_PYR_45):
- case(MSH_PYR_53): case(MSH_PYR_61):
- case(MSH_PYR_69): case(MSH_PYR_1):
-
- case(MSH_PRI_6): case(MSH_PRI_18):
- case(MSH_PRI_15): case(MSH_PRI_1):
- case(MSH_PRI_40): case(MSH_PRI_75):
- case(MSH_PRI_126): case(MSH_PRI_196):
- case(MSH_PRI_288): case(MSH_PRI_405):
- case(MSH_PRI_550): case(MSH_PRI_24):
- case(MSH_PRI_33): case(MSH_PRI_42):
- case(MSH_PRI_51): case(MSH_PRI_60):
- case(MSH_PRI_69): case(MSH_PRI_78):
-
- case(MSH_HEX_8): case(MSH_HEX_27):
- case(MSH_HEX_20): case(MSH_HEX_1):
- case(MSH_HEX_64): case(MSH_HEX_125):
- case(MSH_HEX_216): case(MSH_HEX_343):
- case(MSH_HEX_512): case(MSH_HEX_729):
- case(MSH_HEX_1000): case(MSH_HEX_32):
- case(MSH_HEX_44): case(MSH_HEX_56):
- case(MSH_HEX_68): case(MSH_HEX_80):
- case(MSH_HEX_92): case(MSH_HEX_104):
-
- case(MSH_POLYH_):
- return 3;
-
- default:
- Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
- return 3;
- }
-}
-
-int MElement::getTag(int parentTag, int order, bool serendip)
-{
- switch (parentTag) {
- case TYPE_PNT :
- return MSH_PNT;
- case TYPE_LIN :
- switch(order) {
- case 0 : return MSH_LIN_1;
- case 1 : return MSH_LIN_2;
- case 2 : return MSH_LIN_3;
- case 3 : return MSH_LIN_4;
- case 4 : return MSH_LIN_5;
- case 5 : return MSH_LIN_6;
- case 6 : return MSH_LIN_7;
- case 7 : return MSH_LIN_8;
- case 8 : return MSH_LIN_9;
- case 9 : return MSH_LIN_10;
- case 10: return MSH_LIN_11;
- default : Msg::Error("line order %i unknown", order); return 0;
- }
- break;
- case TYPE_TRI :
- switch(order) {
- case 0 : return MSH_TRI_1;
- case 1 : return MSH_TRI_3;
- case 2 : return MSH_TRI_6;
- case 3 : return serendip ? MSH_TRI_9 : MSH_TRI_10;
- case 4 : return serendip ? MSH_TRI_12 : MSH_TRI_15;
- case 5 : return serendip ? MSH_TRI_15I: MSH_TRI_21;
- case 6 : return serendip ? MSH_TRI_18 : MSH_TRI_28;
- case 7 : return serendip ? MSH_TRI_21I: MSH_TRI_36;
- case 8 : return serendip ? MSH_TRI_24 : MSH_TRI_45;
- case 9 : return serendip ? MSH_TRI_27 : MSH_TRI_55;
- case 10: return serendip ? MSH_TRI_30 : MSH_TRI_66;
- default : Msg::Error("triangle order %i unknown", order); return 0;
- }
- break;
- case TYPE_QUA :
- switch(order) {
- case 0 : return MSH_QUA_1;
- case 1 : return MSH_QUA_4;
- case 2 : return serendip ? MSH_QUA_8 : MSH_QUA_9;
- case 3 : return serendip ? MSH_QUA_12 : MSH_QUA_16;
- case 4 : return serendip ? MSH_QUA_16I: MSH_QUA_25;
- case 5 : return serendip ? MSH_QUA_20 : MSH_QUA_36;
- case 6 : return serendip ? MSH_QUA_24 : MSH_QUA_49;
- case 7 : return serendip ? MSH_QUA_28 : MSH_QUA_64;
- case 8 : return serendip ? MSH_QUA_32 : MSH_QUA_81;
- case 9 : return serendip ? MSH_QUA_36I: MSH_QUA_100;
- case 10: return serendip ? MSH_QUA_40 : MSH_QUA_121;
- default : Msg::Error("quad order %i unknown", order); return 0;
- }
- break;
- case TYPE_TET :
- switch(order) {
- case 0 : return MSH_TET_1;
- case 1 : return MSH_TET_4;
- case 2 : return MSH_TET_10;
- case 3 : return MSH_TET_20;
- case 4 : return serendip ? MSH_TET_34 : MSH_TET_35;
- case 5 : return serendip ? MSH_TET_52 : MSH_TET_56;
- case 6 : return serendip ? MSH_TET_74 : MSH_TET_84;
- case 7 : return serendip ? MSH_TET_100: MSH_TET_120;
- case 8 : return serendip ? MSH_TET_130: MSH_TET_165;
- case 9 : return serendip ? MSH_TET_164: MSH_TET_220;
- case 10: return serendip ? MSH_TET_202: MSH_TET_286;
- default : Msg::Error("terahedron order %i unknown", order); return 0;
- }
- break;
- case TYPE_HEX :
- switch(order) {
- case 0 : return MSH_HEX_1;
- case 1 : return MSH_HEX_8;
- case 2 : return serendip ? MSH_HEX_20 : MSH_HEX_27;
- case 3 : return serendip ? MSH_HEX_32 : MSH_HEX_64;
- case 4 : return serendip ? MSH_HEX_44 : MSH_HEX_125;
- case 5 : return serendip ? MSH_HEX_56: MSH_HEX_216;
- case 6 : return serendip ? MSH_HEX_68: MSH_HEX_343;
- case 7 : return serendip ? MSH_HEX_80: MSH_HEX_512;
- case 8 : return serendip ? MSH_HEX_92: MSH_HEX_729;
- case 9 : return serendip ? MSH_HEX_104: MSH_HEX_1000;
- default : Msg::Error("hexahedron order %i unknown", order); return 0;
- }
- break;
- case TYPE_PRI :
- switch(order) {
- case 0 : return MSH_PRI_1;
- case 1 : return MSH_PRI_6;
- case 2 : return serendip ? MSH_PRI_15 : MSH_PRI_18;
- case 3 : return serendip ? MSH_PRI_24 : MSH_PRI_40;
- case 4 : return serendip ? MSH_PRI_33 : MSH_PRI_75;
- case 5 : return serendip ? MSH_PRI_42 : MSH_PRI_126;
- case 6 : return serendip ? MSH_PRI_51 : MSH_PRI_196;
- case 7 : return serendip ? MSH_PRI_60 : MSH_PRI_288;
- case 8 : return serendip ? MSH_PRI_69 : MSH_PRI_405;
- case 9 : return serendip ? MSH_PRI_78 : MSH_PRI_550;
- default : Msg::Error("prism order %i unknown", order); return 0;
- }
- break;
- case TYPE_PYR :
- switch(order) {
- case 0 : return MSH_PYR_1;
- case 1 : return MSH_PYR_5;
- case 2 : return serendip ? MSH_PYR_13 : MSH_PYR_14;
- case 3: return serendip ? MSH_PYR_21 : MSH_PYR_30;
- case 4: return serendip ? MSH_PYR_29 : MSH_PYR_55;
- case 5: return serendip ? MSH_PYR_37 : MSH_PYR_91;
- case 6: return serendip ? MSH_PYR_45 : MSH_PYR_140;
- case 7: return serendip ? MSH_PYR_53 : MSH_PYR_204;
- case 8: return serendip ? MSH_PYR_61 : MSH_PYR_285;
- case 9: return serendip ? MSH_PYR_69 : MSH_PYR_385;
- default : Msg::Error("pyramid order %i unknown", order); return 0;
- }
- break;
- default : Msg::Error("unknown element type %i", parentTag); return 0;
- }
-}
-
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part, bool owner, MElement *parent,
MElement *d1, MElement *d2)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 910bab5210499a91d197576741f9b2cca4263f0b..328c3adf612e6d6ad558cfcaff327f680c72f209 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -364,12 +364,6 @@ class MElement
virtual MElement *copy(std::map<int, MVertex*> &vertexMap,
std::map<MElement*, MElement*> &newParents,
std::map<MElement*, MElement*> &newDomains);
-
- static int ParentTypeFromTag(int tag);
- static int OrderFromTag(int tag);
- static int SerendipityFromTag(int tag);
- static int DimensionFromTag(int tag);
- static int getTag(int parentTag, int order, bool serendip = false);
};
class MElementFactory{
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index b276baaa4ee0584b0cfdee2be5986ea15dfbdc3a..60f8ea8aa10ea9847393232f2d583a63621fdff2 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -9,21 +9,22 @@
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
#include "BasisFactory.h"
+#include "MElement.h"
std::map<int, nodalBasis*> BasisFactory::fs;
std::map<int, JacobianBasis*> BasisFactory::js;
-std::map<int, bezierBasis*> BasisFactory::bs;
+BasisFactory::Cont_bezierBasis BasisFactory::bs;
-const nodalBasis* BasisFactory::getNodalBasis(int elementType)
+const nodalBasis* BasisFactory::getNodalBasis(int tag)
{
// If the Basis has already been built, return it.
- std::map<int, nodalBasis*>::const_iterator it = fs.find(elementType);
+ std::map<int, nodalBasis*>::const_iterator it = fs.find(tag);
if (it != fs.end()) {
return it->second;
}
// Get the parent type to see which kind of basis
// we want to create
- int parentType = MElement::ParentTypeFromTag(elementType);
+ int parentType = ElementType::ParentTypeFromTag(tag);
nodalBasis* F = NULL;
switch(parentType) {
@@ -34,42 +35,42 @@ const nodalBasis* BasisFactory::getNodalBasis(int elementType)
case(TYPE_PRI):
case(TYPE_TET):
case(TYPE_HEX):
- F = new polynomialBasis(elementType);
+ F = new polynomialBasis(tag);
break;
case(TYPE_PYR):
- F = new pyramidalBasis(elementType);
+ F = new pyramidalBasis(tag);
break;
default:
- Msg::Error("Unknown type of element %d (in BasisFactory)", elementType);
+ Msg::Error("Unknown type of element %d (in BasisFactory)", tag);
return NULL;
}
// FIXME: check if already exists to deallocate if necessary
- fs.insert(std::make_pair(elementType, F));
+ fs.insert(std::make_pair(tag, F));
- return fs[elementType];
+ return F;
}
-const bezierBasis* BasisFactory::getBezierBasis(int elementType)
+const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
{
- std::map<int, bezierBasis*>::const_iterator it = bs.find(elementType);
- if (it != bs.end())
+ std::map<int, JacobianBasis*>::const_iterator it = js.find(tag);
+ if (it != js.end())
return it->second;
- bezierBasis* B = new bezierBasis(elementType);
- if (B) bs.insert(std::make_pair(elementType, B));
- return B;
+ JacobianBasis* J = new JacobianBasis(tag);
+ js.insert(std::make_pair(tag, J));
+ return J;
}
-const JacobianBasis* BasisFactory::getJacobianBasis(int elementType)
+const bezierBasis* BasisFactory::getBezierBasis(int parentType, int order)
{
- std::map<int, JacobianBasis*>::const_iterator it = js.find(elementType);
- if (it != js.end())
+ Cont_bezierBasis::iterator it = bs.find(std::make_pair(parentType, order));
+ if (it != bs.end())
return it->second;
- JacobianBasis* J = new JacobianBasis(elementType);
- if (J) js.insert(std::make_pair(elementType, J));
- return J;
+ bezierBasis* B = new bezierBasis(parentType, order);
+ bs.insert(std::make_pair(std::make_pair(parentType, order), B));
+ return B;
}
void BasisFactory::clearAll()
@@ -88,7 +89,7 @@ void BasisFactory::clearAll()
}
js.clear();
- std::map<int, bezierBasis*>::iterator itB = bs.begin();
+ Cont_bezierBasis::iterator itB = bs.begin();
while (itB != bs.end()) {
delete itB->second;
itB++;
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 61d34fd548870d82994883763100cb6f7fe46c61..0e441643f586711ae16bc652e4c58569b78699c1 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -11,21 +11,27 @@
#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;
+ typedef std::map<std::pair<int, int>, bezierBasis*> Cont_bezierBasis;
+
+ private:
+ static std::map<int, nodalBasis*> fs;
+ static std::map<int, JacobianBasis*> js;
+ static Cont_bezierBasis bs;
+ // store bezier bases by parentType and order (no serendipity..)
- public:
- // 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);
+ public:
+ // Caution: the returned pointer can be NULL
+ static const nodalBasis* getNodalBasis(int tag);
+ static const JacobianBasis* getJacobianBasis(int tag);
+ static const bezierBasis* getBezierBasis(int parentTag, int order);
+ static inline const bezierBasis* getBezierBasis(int tag) {
+ return getBezierBasis(ElementType::ParentTypeFromTag(tag),
+ ElementType::OrderFromTag(tag) );
+ }
- static void clearAll();
+ static void clearAll();
};
#endif
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 231bfe602cdffc0a6525889903140f41d5362a9e..e4f78f93a1e821ce0e1759873d97a2fd7e39abf0 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -15,7 +15,8 @@ set(SRC
legendrePolynomials.cpp
bezierBasis.cpp
JacobianBasis.cpp
- pointsGenerators.cpp
+ pointsGenerators.cpp
+ ElementType.cpp
GaussIntegration.cpp
GaussQuadratureLin.cpp
GaussQuadratureTri.cpp
diff --git a/Numeric/ElementType.cpp b/Numeric/ElementType.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..588995930f0b191a0390411801268bdd780365d1
--- /dev/null
+++ b/Numeric/ElementType.cpp
@@ -0,0 +1,539 @@
+// 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 "ElementType.h"
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+
+int ElementType::ParentTypeFromTag(int tag)
+{
+ switch(tag) {
+ case(MSH_PNT):
+ return TYPE_PNT;
+ case(MSH_LIN_2): case(MSH_LIN_3):
+ case(MSH_LIN_4): case(MSH_LIN_5):
+ case(MSH_LIN_6): case(MSH_LIN_7):
+ case(MSH_LIN_8): case(MSH_LIN_9):
+ case(MSH_LIN_10): case(MSH_LIN_11):
+ case(MSH_LIN_B): case(MSH_LIN_C):
+ case(MSH_LIN_1):
+ return TYPE_LIN;
+ case(MSH_TRI_3): case(MSH_TRI_6):
+ case(MSH_TRI_9): case(MSH_TRI_10):
+ case(MSH_TRI_12): case(MSH_TRI_15):
+ case(MSH_TRI_15I): case(MSH_TRI_21):
+ case(MSH_TRI_28): case(MSH_TRI_36):
+ case(MSH_TRI_45): case(MSH_TRI_55):
+ case(MSH_TRI_66): case(MSH_TRI_18):
+ case(MSH_TRI_21I): case(MSH_TRI_24):
+ case(MSH_TRI_27): case(MSH_TRI_30):
+ case(MSH_TRI_B): case(MSH_TRI_1):
+ return TYPE_TRI;
+ case(MSH_QUA_4): case(MSH_QUA_9):
+ case(MSH_QUA_8): case(MSH_QUA_16):
+ case(MSH_QUA_25): case(MSH_QUA_36):
+ case(MSH_QUA_12): case(MSH_QUA_16I):
+ case(MSH_QUA_20): case(MSH_QUA_49):
+ case(MSH_QUA_64): case(MSH_QUA_81):
+ case(MSH_QUA_100): case(MSH_QUA_121):
+ case(MSH_QUA_24): case(MSH_QUA_28):
+ case(MSH_QUA_32): case(MSH_QUA_36I):
+ case(MSH_QUA_40): case(MSH_QUA_1):
+ return TYPE_QUA;
+ case(MSH_TET_4): case(MSH_TET_10):
+ case(MSH_TET_20): case(MSH_TET_35):
+ case(MSH_TET_56): case(MSH_TET_34):
+ case(MSH_TET_52): case(MSH_TET_84):
+ case(MSH_TET_120): case(MSH_TET_165):
+ case(MSH_TET_220): case(MSH_TET_286):
+ case(MSH_TET_74): case(MSH_TET_100):
+ case(MSH_TET_130): case(MSH_TET_164):
+ case(MSH_TET_202): case(MSH_TET_1):
+ return TYPE_TET;
+ case(MSH_PYR_5): case(MSH_PYR_14):
+ case(MSH_PYR_13): case(MSH_PYR_30):
+ case(MSH_PYR_55): case(MSH_PYR_91):
+ case(MSH_PYR_140): case(MSH_PYR_204):
+ case(MSH_PYR_285): case(MSH_PYR_385):
+ case(MSH_PYR_21): case(MSH_PYR_29):
+ case(MSH_PYR_37): case(MSH_PYR_45):
+ case(MSH_PYR_53): case(MSH_PYR_61):
+ case(MSH_PYR_69): case(MSH_PYR_1):
+ return TYPE_PYR;
+ case(MSH_PRI_6): case(MSH_PRI_18):
+ case(MSH_PRI_15): case(MSH_PRI_1):
+ case(MSH_PRI_40): case(MSH_PRI_75):
+ case(MSH_PRI_126): case(MSH_PRI_196):
+ case(MSH_PRI_288): case(MSH_PRI_405):
+ case(MSH_PRI_550): case(MSH_PRI_24):
+ case(MSH_PRI_33): case(MSH_PRI_42):
+ case(MSH_PRI_51): case(MSH_PRI_60):
+ case(MSH_PRI_69): case(MSH_PRI_78):
+ return TYPE_PRI;
+ case(MSH_HEX_8): case(MSH_HEX_27):
+ case(MSH_HEX_20): case(MSH_HEX_1):
+ case(MSH_HEX_64): case(MSH_HEX_125):
+ case(MSH_HEX_216): case(MSH_HEX_343):
+ case(MSH_HEX_512): case(MSH_HEX_729):
+ case(MSH_HEX_1000): case(MSH_HEX_32):
+ case(MSH_HEX_44): case(MSH_HEX_56):
+ case(MSH_HEX_68): case(MSH_HEX_80):
+ case(MSH_HEX_92): case(MSH_HEX_104):
+ return TYPE_HEX;
+ case(MSH_POLYG_): case(MSH_POLYG_B):
+ return TYPE_POLYG;
+ case(MSH_POLYH_):
+ return TYPE_POLYH;
+ case(MSH_PNT_SUB): case(MSH_LIN_SUB):
+ case(MSH_TRI_SUB): case(MSH_TET_SUB):
+ return TYPE_XFEM;
+ default:
+ Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
+ return TYPE_TET;
+ }
+}
+
+int ElementType::OrderFromTag(int tag)
+{
+
+ switch (tag) {
+ case MSH_PNT : return 0;
+ case MSH_LIN_1 : return 0;
+ case MSH_LIN_2 : return 1;
+ case MSH_LIN_3 : return 2;
+ case MSH_LIN_4 : return 3;
+ case MSH_LIN_5 : return 4;
+ case MSH_LIN_6 : return 5;
+ case MSH_LIN_7 : return 6;
+ case MSH_LIN_8 : return 7;
+ case MSH_LIN_9 : return 8;
+ case MSH_LIN_10 : return 9;
+ case MSH_LIN_11 : return 10;
+ case MSH_TRI_1 : return 0;
+ case MSH_TRI_3 : return 1;
+ case MSH_TRI_6 : return 2;
+ case MSH_TRI_10 : return 3;
+ case MSH_TRI_15 : return 4;
+ case MSH_TRI_21 : return 5;
+ case MSH_TRI_28 : return 6;
+ case MSH_TRI_36 : return 7;
+ case MSH_TRI_45 : return 8;
+ case MSH_TRI_55 : return 9;
+ case MSH_TRI_66 : return 10;
+ case MSH_TRI_9 : return 3;
+ case MSH_TRI_12 : return 4;
+ case MSH_TRI_15I : return 5;
+ case MSH_TRI_18 : return 6;
+ case MSH_TRI_21I : return 7;
+ case MSH_TRI_24 : return 8;
+ case MSH_TRI_27 : return 9;
+ case MSH_TRI_30 : return 10;
+ case MSH_TET_1 : return 0;
+ case MSH_TET_4 : return 1;
+ case MSH_TET_10 : return 2;
+ case MSH_TET_20 : return 3;
+ case MSH_TET_35 : return 4;
+ case MSH_TET_56 : return 5;
+ case MSH_TET_84 : return 6;
+ case MSH_TET_120 : return 7;
+ case MSH_TET_165 : return 8;
+ case MSH_TET_220 : return 9;
+ case MSH_TET_286 : return 10;
+ case MSH_TET_34 : return 4;
+ case MSH_TET_52 : return 5;
+ case MSH_TET_74 : return 6;
+ case MSH_TET_100 : return 7;
+ case MSH_TET_130 : return 8;
+ case MSH_TET_164 : return 9;
+ case MSH_TET_202 : return 10;
+ case MSH_QUA_1 : return 0;
+ case MSH_QUA_4 : return 1;
+ case MSH_QUA_9 : return 2;
+ case MSH_QUA_16 : return 3;
+ case MSH_QUA_25 : return 4;
+ case MSH_QUA_36 : return 5;
+ case MSH_QUA_49 : return 6;
+ case MSH_QUA_64 : return 7;
+ case MSH_QUA_81 : return 8;
+ case MSH_QUA_100 : return 9;
+ case MSH_QUA_121 : return 10;
+ case MSH_QUA_8 : return 2;
+ case MSH_QUA_12 : return 3;
+ case MSH_QUA_16I : return 4;
+ case MSH_QUA_20 : return 5;
+ case MSH_QUA_24 : return 6;
+ case MSH_QUA_28 : return 7;
+ case MSH_QUA_32 : return 8;
+ case MSH_QUA_36I : return 9;
+ case MSH_QUA_40 : return 10;
+ case MSH_PRI_1 : return 0;
+ case MSH_PRI_6 : return 1;
+ case MSH_PRI_18 : return 2;
+ case MSH_PRI_40 : return 3;
+ case MSH_PRI_75 : return 4;
+ case MSH_PRI_126 : return 5;
+ case MSH_PRI_196 : return 6;
+ case MSH_PRI_288 : return 7;
+ case MSH_PRI_405 : return 8;
+ case MSH_PRI_550 : return 9;
+ case MSH_PRI_15 : return 2;
+ case MSH_PRI_24 : return 3;
+ case MSH_PRI_33 : return 4;
+ case MSH_PRI_42 : return 5;
+ case MSH_PRI_51 : return 6;
+ case MSH_PRI_60 : return 7;
+ case MSH_PRI_69 : return 8;
+ case MSH_PRI_78 : return 9;
+ case MSH_HEX_1 : return 0;
+ case MSH_HEX_8 : return 1;
+ case MSH_HEX_27 : return 2;
+ case MSH_HEX_64 : return 3;
+ case MSH_HEX_125 : return 4;
+ case MSH_HEX_216 : return 5;
+ case MSH_HEX_343 : return 6;
+ case MSH_HEX_512 : return 7;
+ case MSH_HEX_729 : return 8;
+ case MSH_HEX_1000: return 9;
+ case MSH_HEX_20 : return 2;
+ case MSH_HEX_32 : return 3;
+ case MSH_HEX_44 : return 4;
+ case MSH_HEX_56 : return 5;
+ case MSH_HEX_68 : return 6;
+ case MSH_HEX_80 : return 7;
+ case MSH_HEX_92 : return 8;
+ case MSH_HEX_104 : return 9;
+ case MSH_PYR_1 : return 0;
+ case MSH_PYR_5 : return 1;
+ case MSH_PYR_14 : return 2;
+ case MSH_PYR_30 : return 3;
+ case MSH_PYR_55 : return 4;
+ case MSH_PYR_91 : return 5;
+ case MSH_PYR_140 : return 6;
+ case MSH_PYR_204 : return 7;
+ case MSH_PYR_285 : return 8;
+ case MSH_PYR_385 : return 9;
+ case MSH_PYR_13 : return 2;
+ case MSH_PYR_21 : return 3;
+ case MSH_PYR_29 : return 4;
+ case MSH_PYR_37 : return 5;
+ case MSH_PYR_45 : return 6;
+ case MSH_PYR_53 : return 7;
+ case MSH_PYR_61 : return 8;
+ case MSH_PYR_69 : return 9;
+ default :
+ Msg::Error("Unknown element type %d: reverting to order 1",tag);
+ return 1;
+ }
+
+}
+
+int ElementType::DimensionFromTag(int tag)
+{
+ switch(tag) {
+ case(MSH_PNT): case(MSH_PNT_SUB):
+ return 0;
+
+ case(MSH_LIN_2): case(MSH_LIN_3):
+ case(MSH_LIN_4): case(MSH_LIN_5):
+ case(MSH_LIN_6): case(MSH_LIN_7):
+ case(MSH_LIN_8): case(MSH_LIN_9):
+ case(MSH_LIN_10): case(MSH_LIN_11):
+ case(MSH_LIN_B): case(MSH_LIN_C):
+ case(MSH_LIN_1): case(MSH_LIN_SUB):
+ return 1;
+
+ case(MSH_TRI_3): case(MSH_TRI_6):
+ case(MSH_TRI_9): case(MSH_TRI_10):
+ case(MSH_TRI_12): case(MSH_TRI_15):
+ case(MSH_TRI_15I): case(MSH_TRI_21):
+ case(MSH_TRI_28): case(MSH_TRI_36):
+ case(MSH_TRI_45): case(MSH_TRI_55):
+ case(MSH_TRI_66): case(MSH_TRI_18):
+ case(MSH_TRI_21I): case(MSH_TRI_24):
+ case(MSH_TRI_27): case(MSH_TRI_30):
+ case(MSH_TRI_B): case(MSH_TRI_1):
+ case(MSH_TRI_SUB):
+
+ case(MSH_QUA_4): case(MSH_QUA_9):
+ case(MSH_QUA_8): case(MSH_QUA_16):
+ case(MSH_QUA_25): case(MSH_QUA_36):
+ case(MSH_QUA_12): case(MSH_QUA_16I):
+ case(MSH_QUA_20): case(MSH_QUA_49):
+ case(MSH_QUA_64): case(MSH_QUA_81):
+ case(MSH_QUA_100): case(MSH_QUA_121):
+ case(MSH_QUA_24): case(MSH_QUA_28):
+ case(MSH_QUA_32): case(MSH_QUA_36I):
+ case(MSH_QUA_40): case(MSH_QUA_1):
+
+ case(MSH_POLYG_): case(MSH_POLYG_B):
+ return 2;
+
+ case(MSH_TET_4): case(MSH_TET_10):
+ case(MSH_TET_20): case(MSH_TET_35):
+ case(MSH_TET_56): case(MSH_TET_34):
+ case(MSH_TET_52): case(MSH_TET_84):
+ case(MSH_TET_120): case(MSH_TET_165):
+ case(MSH_TET_220): case(MSH_TET_286):
+ case(MSH_TET_74): case(MSH_TET_100):
+ case(MSH_TET_130): case(MSH_TET_164):
+ case(MSH_TET_202): case(MSH_TET_1):
+ case(MSH_TET_SUB):
+
+ case(MSH_PYR_5): case(MSH_PYR_14):
+ case(MSH_PYR_13): case(MSH_PYR_30):
+ case(MSH_PYR_55): case(MSH_PYR_91):
+ case(MSH_PYR_140): case(MSH_PYR_204):
+ case(MSH_PYR_285): case(MSH_PYR_385):
+ case(MSH_PYR_21): case(MSH_PYR_29):
+ case(MSH_PYR_37): case(MSH_PYR_45):
+ case(MSH_PYR_53): case(MSH_PYR_61):
+ case(MSH_PYR_69): case(MSH_PYR_1):
+
+ case(MSH_PRI_6): case(MSH_PRI_18):
+ case(MSH_PRI_15): case(MSH_PRI_1):
+ case(MSH_PRI_40): case(MSH_PRI_75):
+ case(MSH_PRI_126): case(MSH_PRI_196):
+ case(MSH_PRI_288): case(MSH_PRI_405):
+ case(MSH_PRI_550): case(MSH_PRI_24):
+ case(MSH_PRI_33): case(MSH_PRI_42):
+ case(MSH_PRI_51): case(MSH_PRI_60):
+ case(MSH_PRI_69): case(MSH_PRI_78):
+
+ case(MSH_HEX_8): case(MSH_HEX_27):
+ case(MSH_HEX_20): case(MSH_HEX_1):
+ case(MSH_HEX_64): case(MSH_HEX_125):
+ case(MSH_HEX_216): case(MSH_HEX_343):
+ case(MSH_HEX_512): case(MSH_HEX_729):
+ case(MSH_HEX_1000): case(MSH_HEX_32):
+ case(MSH_HEX_44): case(MSH_HEX_56):
+ case(MSH_HEX_68): case(MSH_HEX_80):
+ case(MSH_HEX_92): case(MSH_HEX_104):
+
+ case(MSH_POLYH_):
+ return 3;
+
+ default:
+ Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
+ return 3;
+ }
+}
+
+int ElementType::SerendipityFromTag(int tag)
+{
+ switch (tag) {
+ case MSH_PNT : case MSH_LIN_1 :
+ case MSH_LIN_2 : case MSH_LIN_3 :
+ case MSH_LIN_4 : case MSH_LIN_5 :
+ case MSH_LIN_6 : case MSH_LIN_7 :
+ case MSH_LIN_8 : case MSH_LIN_9 :
+ case MSH_LIN_10 : case MSH_LIN_11 :
+
+ case MSH_TRI_1 : case MSH_TRI_3 :
+ case MSH_TRI_6 :
+
+ case MSH_QUA_1 : case MSH_QUA_4 :
+
+ case MSH_TET_1 : case MSH_TET_4 :
+ case MSH_TET_10 : case MSH_TET_20 :
+
+ case MSH_PRI_1 : case MSH_PRI_6 :
+
+ case MSH_HEX_1 : case MSH_HEX_8 :
+
+ case MSH_PYR_1 : case MSH_PYR_5 :
+
+ return 1; // Serendipity or not
+
+
+ case MSH_TRI_10 : case MSH_TRI_15 :
+ case MSH_TRI_21 : case MSH_TRI_28 :
+ case MSH_TRI_36 : case MSH_TRI_45 :
+ case MSH_TRI_55 : case MSH_TRI_66 :
+
+ case MSH_QUA_9 : case MSH_QUA_16 :
+ case MSH_QUA_25 : case MSH_QUA_36 :
+ case MSH_QUA_49 : case MSH_QUA_64 :
+ case MSH_QUA_81 : case MSH_QUA_100 :
+ case MSH_QUA_121 :
+
+ case MSH_TET_35 : case MSH_TET_56 :
+ case MSH_TET_84 : case MSH_TET_120 :
+ case MSH_TET_165 : case MSH_TET_220 :
+ case MSH_TET_286 :
+
+ case MSH_PRI_18 : case MSH_PRI_40 :
+ case MSH_PRI_75 : case MSH_PRI_126 :
+ case MSH_PRI_196 : case MSH_PRI_288 :
+ case MSH_PRI_405 : case MSH_PRI_550 :
+
+ case MSH_HEX_27 : case MSH_HEX_64 :
+ case MSH_HEX_125 : case MSH_HEX_216 :
+ case MSH_HEX_343 : case MSH_HEX_512 :
+ case MSH_HEX_729 : case MSH_HEX_1000:
+
+ case MSH_PYR_14 : case MSH_PYR_30 :
+ case MSH_PYR_55 : case MSH_PYR_91 :
+ case MSH_PYR_140 : case MSH_PYR_204 :
+ case MSH_PYR_285 : case MSH_PYR_385 :
+
+ return 0; // Not Serendipity
+
+
+ case MSH_TRI_9 : case MSH_TRI_12 :
+ case MSH_TRI_15I : case MSH_TRI_18 :
+ case MSH_TRI_21I : case MSH_TRI_24 :
+ case MSH_TRI_27 : case MSH_TRI_30 :
+
+ case MSH_QUA_8 : case MSH_QUA_12 :
+ case MSH_QUA_16I : case MSH_QUA_20 :
+ case MSH_QUA_24 : case MSH_QUA_28 :
+ case MSH_QUA_32 : case MSH_QUA_36I :
+ case MSH_QUA_40 :
+
+ case MSH_TET_34 : case MSH_TET_52 :
+ case MSH_TET_74 : case MSH_TET_100 :
+ case MSH_TET_130 : case MSH_TET_164 :
+ case MSH_TET_202 :
+
+ case MSH_PRI_15 : case MSH_PRI_24 :
+ case MSH_PRI_33 : case MSH_PRI_42 :
+ case MSH_PRI_51 : case MSH_PRI_60 :
+ case MSH_PRI_69 : case MSH_PRI_78 :
+
+ case MSH_HEX_20 : case MSH_HEX_32 :
+ case MSH_HEX_44 : case MSH_HEX_56 :
+ case MSH_HEX_68 : case MSH_HEX_80 :
+ case MSH_HEX_92 : case MSH_HEX_104 :
+
+ case MSH_PYR_13 : case MSH_PYR_21 :
+ case MSH_PYR_29 : case MSH_PYR_37 :
+ case MSH_PYR_45 : case MSH_PYR_53 :
+ case MSH_PYR_61 : case MSH_PYR_69 :
+
+ return 2; // Only Serendipity
+
+ default :
+ Msg::Error("Unknown element type %d: assuming not serendipity",tag);
+ return 0;
+ }
+}
+
+int ElementType::getTag(int parentTag, int order, bool serendip)
+{
+ switch (parentTag) {
+ case TYPE_PNT :
+ return MSH_PNT;
+ case TYPE_LIN :
+ switch(order) {
+ case 0 : return MSH_LIN_1;
+ case 1 : return MSH_LIN_2;
+ case 2 : return MSH_LIN_3;
+ case 3 : return MSH_LIN_4;
+ case 4 : return MSH_LIN_5;
+ case 5 : return MSH_LIN_6;
+ case 6 : return MSH_LIN_7;
+ case 7 : return MSH_LIN_8;
+ case 8 : return MSH_LIN_9;
+ case 9 : return MSH_LIN_10;
+ case 10: return MSH_LIN_11;
+ default : Msg::Error("line order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_TRI :
+ switch(order) {
+ case 0 : return MSH_TRI_1;
+ case 1 : return MSH_TRI_3;
+ case 2 : return MSH_TRI_6;
+ case 3 : return serendip ? MSH_TRI_9 : MSH_TRI_10;
+ case 4 : return serendip ? MSH_TRI_12 : MSH_TRI_15;
+ case 5 : return serendip ? MSH_TRI_15I: MSH_TRI_21;
+ case 6 : return serendip ? MSH_TRI_18 : MSH_TRI_28;
+ case 7 : return serendip ? MSH_TRI_21I: MSH_TRI_36;
+ case 8 : return serendip ? MSH_TRI_24 : MSH_TRI_45;
+ case 9 : return serendip ? MSH_TRI_27 : MSH_TRI_55;
+ case 10: return serendip ? MSH_TRI_30 : MSH_TRI_66;
+ default : Msg::Error("triangle order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_QUA :
+ switch(order) {
+ case 0 : return MSH_QUA_1;
+ case 1 : return MSH_QUA_4;
+ case 2 : return serendip ? MSH_QUA_8 : MSH_QUA_9;
+ case 3 : return serendip ? MSH_QUA_12 : MSH_QUA_16;
+ case 4 : return serendip ? MSH_QUA_16I: MSH_QUA_25;
+ case 5 : return serendip ? MSH_QUA_20 : MSH_QUA_36;
+ case 6 : return serendip ? MSH_QUA_24 : MSH_QUA_49;
+ case 7 : return serendip ? MSH_QUA_28 : MSH_QUA_64;
+ case 8 : return serendip ? MSH_QUA_32 : MSH_QUA_81;
+ case 9 : return serendip ? MSH_QUA_36I: MSH_QUA_100;
+ case 10: return serendip ? MSH_QUA_40 : MSH_QUA_121;
+ default : Msg::Error("quad order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_TET :
+ switch(order) {
+ case 0 : return MSH_TET_1;
+ case 1 : return MSH_TET_4;
+ case 2 : return MSH_TET_10;
+ case 3 : return MSH_TET_20;
+ case 4 : return serendip ? MSH_TET_34 : MSH_TET_35;
+ case 5 : return serendip ? MSH_TET_52 : MSH_TET_56;
+ case 6 : return serendip ? MSH_TET_74 : MSH_TET_84;
+ case 7 : return serendip ? MSH_TET_100: MSH_TET_120;
+ case 8 : return serendip ? MSH_TET_130: MSH_TET_165;
+ case 9 : return serendip ? MSH_TET_164: MSH_TET_220;
+ case 10: return serendip ? MSH_TET_202: MSH_TET_286;
+ default : Msg::Error("terahedron order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_HEX :
+ switch(order) {
+ case 0 : return MSH_HEX_1;
+ case 1 : return MSH_HEX_8;
+ case 2 : return serendip ? MSH_HEX_20 : MSH_HEX_27;
+ case 3 : return serendip ? MSH_HEX_32 : MSH_HEX_64;
+ case 4 : return serendip ? MSH_HEX_44 : MSH_HEX_125;
+ case 5 : return serendip ? MSH_HEX_56: MSH_HEX_216;
+ case 6 : return serendip ? MSH_HEX_68: MSH_HEX_343;
+ case 7 : return serendip ? MSH_HEX_80: MSH_HEX_512;
+ case 8 : return serendip ? MSH_HEX_92: MSH_HEX_729;
+ case 9 : return serendip ? MSH_HEX_104: MSH_HEX_1000;
+ default : Msg::Error("hexahedron order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_PRI :
+ switch(order) {
+ case 0 : return MSH_PRI_1;
+ case 1 : return MSH_PRI_6;
+ case 2 : return serendip ? MSH_PRI_15 : MSH_PRI_18;
+ case 3 : return serendip ? MSH_PRI_24 : MSH_PRI_40;
+ case 4 : return serendip ? MSH_PRI_33 : MSH_PRI_75;
+ case 5 : return serendip ? MSH_PRI_42 : MSH_PRI_126;
+ case 6 : return serendip ? MSH_PRI_51 : MSH_PRI_196;
+ case 7 : return serendip ? MSH_PRI_60 : MSH_PRI_288;
+ case 8 : return serendip ? MSH_PRI_69 : MSH_PRI_405;
+ case 9 : return serendip ? MSH_PRI_78 : MSH_PRI_550;
+ default : Msg::Error("prism order %i unknown", order); return 0;
+ }
+ break;
+ case TYPE_PYR :
+ switch(order) {
+ case 0 : return MSH_PYR_1;
+ case 1 : return MSH_PYR_5;
+ case 2 : return serendip ? MSH_PYR_13 : MSH_PYR_14;
+ case 3: return serendip ? MSH_PYR_21 : MSH_PYR_30;
+ case 4: return serendip ? MSH_PYR_29 : MSH_PYR_55;
+ case 5: return serendip ? MSH_PYR_37 : MSH_PYR_91;
+ case 6: return serendip ? MSH_PYR_45 : MSH_PYR_140;
+ case 7: return serendip ? MSH_PYR_53 : MSH_PYR_204;
+ case 8: return serendip ? MSH_PYR_61 : MSH_PYR_285;
+ case 9: return serendip ? MSH_PYR_69 : MSH_PYR_385;
+ default : Msg::Error("pyramid order %i unknown", order); return 0;
+ }
+ break;
+ default : Msg::Error("unknown element type %i", parentTag); return 0;
+ }
+}
diff --git a/Numeric/ElementType.h b/Numeric/ElementType.h
new file mode 100644
index 0000000000000000000000000000000000000000..3338abc6d8e6a0e38e326f05a9fb28107a7b9c83
--- /dev/null
+++ b/Numeric/ElementType.h
@@ -0,0 +1,26 @@
+// 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 _ELEMENTTYPE_H_
+#define _ELEMENTTYPE_H_
+
+class ElementType
+{
+ public :
+ // Give parent type, order & dimension
+ // corresponding to any element type.
+ static int ParentTypeFromTag(int tag);
+ static int OrderFromTag(int tag);
+ static int DimensionFromTag(int tag);
+
+ // Gives > 0 if element type is in Serendipity Family.
+ // Gives < 2 if element type is in 'Normal' Family.
+ // 1 is for element that is either Serendipity or not !
+ static int SerendipityFromTag(int tag);
+
+ static int getTag(int parentTag, int order, bool serendip = false);
+};
+
+#endif
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 1acc387780c66bffb9079c26342313bd850daaf2..95656b01a429afd4813eb07de7aedd295cc74673 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -16,11 +16,12 @@
JacobianBasis::JacobianBasis(int tag)
{
+ const int parentType = ElementType::ParentTypeFromTag(tag);
+ const int order = ElementType::OrderFromTag(tag);
- const int parentType = MElement::ParentTypeFromTag(tag);
int jacType = 0, primJacType = 0;
- if (parentType == TYPE_PYR) {
+ /*if (parentType == TYPE_PYR) {
switch (tag) {
case MSH_PYR_5 : jacType = MSH_PYR_14; primJacType = MSH_PYR_14; break; // TODO: Order 1, Jac. "order" 2, check this
case MSH_PYR_14 : jacType = MSH_PYR_91; primJacType = MSH_PYR_14; break; // TODO: Order 2, Jac. "order" 5, check this
@@ -30,8 +31,7 @@ JacobianBasis::JacobianBasis(int tag)
break;
}
}
- else {
- const int order = MElement::OrderFromTag(tag);
+ else {*/
int jacobianOrder, primJacobianOrder;
switch (parentType) {
case TYPE_PNT : jacobianOrder = 0; primJacobianOrder = 0; break;
@@ -46,9 +46,9 @@ JacobianBasis::JacobianBasis(int tag)
jacobianOrder = 0;
break;
}
- jacType = MElement::getTag(parentType, jacobianOrder, false);
- primJacType = MElement::getTag(parentType, primJacobianOrder, false);
- }
+ jacType = ElementType::getTag(parentType, jacobianOrder, false);
+ primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
+ //}
// Store Bezier basis
bezier = BasisFactory::getBezierBasis(jacType);
@@ -77,7 +77,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 = MElement::getTag(parentType, 1, false);
+ const int primMapType = ElementType::getTag(parentType, 1, false);
const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
numPrimMapNodes = primMapBasis->getNumShapeFunctions();
double xBar = 0., yBar = 0., zBar = 0.;
@@ -104,8 +104,6 @@ JacobianBasis::JacobianBasis(int tag)
}
-
-
// 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
{
@@ -137,8 +135,6 @@ double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullM
}
-
-
// Computes (unit) normal to straight surface element at barycenter (with norm as return value)
double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
{
@@ -163,8 +159,6 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
}
-
-
static inline double calcDet3D(double dxdX, double dydX, double dzdX,
double dxdY, double dydY, double dzdY,
double dxdZ, double dydZ, double dzdZ)
@@ -173,8 +167,6 @@ static inline double calcDet3D(double dxdX, double dydX, double dzdX,
- dxdZ*dydY*dzdX - dxdY*dydX*dzdZ - dydZ*dzdY*dxdX;
}
-
-
// Returns absolute value of Jacobian of straight volume element at barycenter
double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
{
@@ -196,8 +188,6 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
}
-
-
// Calculate (signed) Jacobian at mapping's nodes for one element, with normal vectors to
// straight element for regularization
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
@@ -230,8 +220,6 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVe
}
-
-
// Calculate scaled (signed) Jacobian at mapping's nodes for one element, with normal vectors to
// straight element for regularization and scaling
void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
@@ -268,8 +256,6 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
}
-
-
// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization
// of line and surface Jacobians in 3D
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
@@ -326,8 +312,6 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
}
-
-
// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization
// of line and surface Jacobians in 3D
void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
@@ -430,8 +414,6 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
}
-
-
void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const
@@ -485,8 +467,6 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
}
-
-
// Calculate (signed) Jacobian at mapping's nodes for several elements
// TODO: Optimize and test 1D & 2D
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 2e8f52b23c84600c3576e7f5bdfe34e2495bde03..d5f34792bfe5d0fe36e62ba3b22cfead7b63729c 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -20,8 +20,11 @@ class JacobianBasis {
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numMapNodes, numPrimJacNodes, numPrimMapNodes;
inline const fullMatrix<double> &getPoints() const { return bezier->lagPoints; }
+
public :
JacobianBasis(int tag);
+
+ // get methods
inline int getNumJacNodes() const { return numJacNodes; }
inline int getNumMapNodes() const { return numMapNodes; }
inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
@@ -29,6 +32,8 @@ class JacobianBasis {
inline int getNumDivisions() const { return bezier->numDivisions; }
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
inline int getNumLagPts() const { return bezier->numLagPts; }
+
+ //
double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
double getPrimJac3D(const fullMatrix<double> &nodesXYZ) const;
@@ -43,6 +48,8 @@ class JacobianBasis {
void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+
+ //
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
bezier->matrixLag2Bez.mult(jac,bez);
}
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 9f0592dc93802fe471eb8190332e45f609e7f7bc..f8884abebae5ff040d7676ba00167dea36d8d742 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -13,633 +13,14 @@
#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[0] = gmshGenerateMonomialsLine(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
subPoints[1].add(.5);
@@ -652,8 +33,8 @@ 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[0] = gmshGenerateMonomialsTriangle(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -675,8 +56,8 @@ 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[0] = gmshGenerateMonomialsQuadrangle(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -699,8 +80,8 @@ static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
fullMatrix<double> prox1;
fullMatrix<double> prox2;
- subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
- subPoints[0].scale(.5);
+ subPoints[0] = gmshGenerateMonomialsTetrahedron(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox1.setAsProxy(subPoints[1], 0, 1);
@@ -768,8 +149,8 @@ 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[0] = gmshGenerateMonomialsPrism(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -808,8 +189,8 @@ 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[0] = gmshGenerateMonomialsHexahedron(order);
+ subPoints[0].scale(.5/order);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -861,7 +242,6 @@ static int nChoosek(int n, int k)
return c;
}
-
static fullMatrix<double> generateBez2LagMatrix
(const fullMatrix<double> &exponent, const fullMatrix<double> &point,
int order, int dimSimplex)
@@ -905,8 +285,6 @@ static fullMatrix<double> generateBez2LagMatrix
return bez2Lag;
}
-
-
static fullMatrix<double> generateSubDivisor
(const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
@@ -933,34 +311,11 @@ static fullMatrix<double> generateSubDivisor
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)
+void bezierBasis::_construct(int parentType, int p)
{
- order = MElement::OrderFromTag(tag);
+ order = p;
- if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+ /*if (parentType == TYPE_PYR) {
dim = 3;
numLagPts = 5;
bezierPoints = gmshGeneratePointsPyramid(order,false);
@@ -973,97 +328,93 @@ bezierBasis::bezierBasis(int tag)
}
// TODO: Implement subdidivsor
}
- else {
+ else {*/
int dimSimplex;
+ fullMatrix<double> exponents;
std::vector< fullMatrix<double> > subPoints;
- switch (MElement::ParentTypeFromTag(tag)) {
+
+ switch (parentType) {
case TYPE_PNT :
dim = 0;
numLagPts = 1;
- exponents = generate1DExponents(0);
- bezierPoints = generate1DPoints(0);
- lagPoints = bezierPoints;
dimSimplex = 0;
+ exponents = gmshGenerateMonomialsLine(0);
+ lagPoints = gmshGeneratePointsLine(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);
+ exponents = gmshGenerateMonomialsLine(order);
+ lagPoints = gmshGeneratePointsLine(order);
+ subPoints = generateSubPointsLine(order);
break;
}
case TYPE_TRI : {
dim = 2;
numLagPts = 3;
- exponents = generateExponentsTriangle(order);
- bezierPoints = gmshGeneratePointsTriangle(order,false);
- lagPoints = bezierPoints;
dimSimplex = 2;
+ exponents = gmshGenerateMonomialsTriangle(order);
+ lagPoints = gmshGeneratePointsTriangle(order);
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;
+ exponents = gmshGenerateMonomialsQuadrangle(order);
+ lagPoints = gmshGeneratePointsQuadrangle(order);
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;
+ exponents = gmshGenerateMonomialsTetrahedron(order);
+ lagPoints = gmshGeneratePointsTetrahedron(order);
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;
+ exponents = gmshGenerateMonomialsPrism(order);
+ lagPoints = gmshGeneratePointsPrism(order);
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;
+ exponents = gmshGenerateMonomialsHexahedron(order);
+ lagPoints = gmshGeneratePointsHexahedron(order);
subPoints = generateSubPointsHex(order);
break;
}
default : {
- Msg::Error("Unknown function space %d: reverting to TET_1", tag);
+ Msg::Error("Unknown function space of parentType %d : "
+ "reverting to TET_1", parentType);
dim = 3;
+ order = 0;
numLagPts = 4;
- exponents = generateExponentsTetrahedron(0);
- bezierPoints = gmshGeneratePointsTetrahedron(0, false);
- lagPoints = bezierPoints;
dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(0);
+ exponents = gmshGenerateMonomialsTetrahedron(order);
+ lagPoints = gmshGeneratePointsTetrahedron(order);
+ subPoints = generateSubPointsTetrahedron(order);
break;
}
}
- Msg::Info("%d %d %d %d %d %d", exponents.size1(), exponents.size2(), bezierPoints.size1(), bezierPoints.size2(), order, dimSimplex);
+
+ fullMatrix<double> bezierPoints = exponents;
+ bezierPoints.scale(1./order);
+
+ numDivisions = static_cast<int>(subPoints.size());
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
index e8241ff6b7881c360c2d4f6a6bb519b731873876..62aca79d0810d1c243f81c3ec0f832a119d1bab8 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -9,20 +9,38 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
+#include "ElementType.h"
+
+class MElement;
class bezierBasis {
- private :
- public :
+ public:
+ //private :
+ // the 'numLagPts' first exponents are related to 'Lagrangian' values
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;
+
+ public :
+ fullMatrix<double> lagPoints; // sampling point
+ fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
+
+ inline bezierBasis(int tag) {
+ _construct(ElementType::ParentTypeFromTag(tag), ElementType::OrderFromTag(tag));
+ }
+ inline bezierBasis(int parendtType, int order) {
+ _construct(parendtType, order);
+ }
+
+ // get methods
+ inline int getDim() {return dim;}
+ inline int getOrder() {return order;}
+ inline int getNumLagPts() {return numLagPts;}
+ inline int getNumDivision() {return numDivisions;}
+
+ private :
+ void _construct(int parendtType, int order);
};
#endif
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 874d8d1e13ec18b7a109aa4b51ab9dfa3f5bb54d..c6c97ea1e936c0609a11636483572e21d819f398 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -49,7 +49,7 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate,
{
closure.clear();
closure.resize((order + 1) * (order + 2) / 2);
- closure.type = MElement::getTag(TYPE_TRI, order, false);
+ closure.type = ElementType::getTag(TYPE_TRI, order, false);
switch (order){
case 0:
@@ -146,7 +146,7 @@ static void generate2dEdgeClosureFull(nodalBasis::clCont &closure,
if (serendip) break;
}
for (int r = 0; r < nNod*2 ; r++) {
- closure[r].type = MElement::getTag(TYPE_LIN, order);
+ closure[r].type = ElementType::getTag(TYPE_LIN, order);
closureRef[r] = 0;
}
}
@@ -362,7 +362,7 @@ static void generateFaceClosureHex(nodalBasis::clCont &closure, int order,
bool serendip, const fullMatrix<double> &points)
{
closure.clear();
- const nodalBasis &fsFace = *BasisFactory::getNodalBasis(MElement::getTag(TYPE_QUA, order, serendip));
+ const nodalBasis &fsFace = *BasisFactory::getNodalBasis(ElementType::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++) {
@@ -451,7 +451,7 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
// int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
// {8, 13, 11, 7}, {9, 11, 14, 10}};
int nVertex = isTriangle ? 3 : 4;
- closure.type = MElement::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
+ closure.type = ElementType::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
for (int i = 0; i < nVertex; ++i){
int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
closure[i] = order1node[iFace][k];
@@ -575,7 +575,7 @@ static void generate2dEdgeClosure(nodalBasis::clCont &closure, int order, int nN
closure[j].push_back( nNod + (order-1)*j + i );
closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
}
- closure[j].type = closure[nNod+j].type = MElement::getTag(TYPE_LIN, order);
+ closure[j].type = closure[nNod+j].type = ElementType::getTag(TYPE_LIN, order);
}
}
@@ -596,10 +596,10 @@ static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
nodalBasis::nodalBasis(int tag)
{
type = tag;
- parentType = MElement::ParentTypeFromTag(tag);
- order = MElement::OrderFromTag(tag);
- serendip = MElement::SerendipityFromTag(tag) > 1;
- dimension = MElement::DimensionFromTag(tag);
+ parentType = ElementType::ParentTypeFromTag(tag);
+ order = ElementType::OrderFromTag(tag);
+ serendip = ElementType::SerendipityFromTag(tag) > 1;
+ dimension = ElementType::DimensionFromTag(tag);
switch (parentType) {
case TYPE_PNT :
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index d2c2b50571caa54061ed6c77f7a85b0bf2cf2bcc..f5eca95871aec45e5cba614b37b31af77a03a07f 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -11,29 +11,12 @@
#include "MPrism.h"
#include "MPyramid.h"
-void getDoubleMonomials(fullMatrix<double>& doubleMon,
- fullMatrix<int> (*genIntMonomials)(int, bool),
- int order, bool serendip )
-{
- fullMatrix<int> mon;
- mon = (*genIntMonomials)(order, serendip);
-
- doubleMon.resize(mon.size1(), mon.size2());
-
- for (int i = 0; i < mon.size1(); ++i) {
- for (int j = 0; j < mon.size2(); ++j) {
- doubleMon(i, j) = static_cast<double>(mon(i, j));
- }
- }
-}
// Points Generators
fullMatrix<double> gmshGeneratePointsLine(int order)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsLine, order, false);
-
+ fullMatrix<double> points = gmshGenerateMonomialsLine(order);
if (order == 0) return points;
points.scale(2./order);
@@ -43,9 +26,7 @@ fullMatrix<double> gmshGeneratePointsLine(int order)
fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsTriangle, order, serendip);
-
+ fullMatrix<double> points = gmshGenerateMonomialsTriangle(order, serendip);
if (order == 0) return points;
points.scale(1./order);
@@ -54,8 +35,7 @@ fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsQuadrangle, order, serendip);
+ fullMatrix<double> points = gmshGenerateMonomialsQuadrangle(order, serendip);
if (order == 0) return points;
@@ -66,8 +46,7 @@ fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip)
fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsTetrahedron, order, serendip);
+ fullMatrix<double> points = gmshGenerateMonomialsTetrahedron(order, serendip);
if (order == 0) return points;
@@ -77,8 +56,7 @@ fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsHexahedron, order, serendip);
+ fullMatrix<double> points = gmshGenerateMonomialsHexahedron(order, serendip);
if (order == 0) return points;
@@ -89,8 +67,7 @@ fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip)
fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsPrism, order, serendip);
+ fullMatrix<double> points = gmshGenerateMonomialsPrism(order, serendip);
if (order == 0) return points;
@@ -107,8 +84,7 @@ fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
{
- fullMatrix<double> points;
- getDoubleMonomials(points, gmshGenerateMonomialsPyramid, order, serendip);
+ fullMatrix<double> points = gmshGenerateMonomialsPyramid(order, serendip);
if (order == 0) return points;
@@ -123,9 +99,9 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
// Monomials Generators
-fullMatrix<int> gmshGenerateMonomialsLine(int order, bool serendip)
+fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip)
{
- fullMatrix<int> monomials(order + 1, 1);
+ fullMatrix<double> monomials(order + 1, 1);
monomials(0,0) = 0;
if (order > 0) {
monomials(1, 0) = order;
@@ -134,11 +110,11 @@ fullMatrix<int> gmshGenerateMonomialsLine(int order, bool serendip)
return monomials;
}
-fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
+fullMatrix<double> gmshGenerateMonomialsTriangle(int order, bool serendip)
{
int nbMonomials = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
if (serendip && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 2);
+ fullMatrix<double> monomials(nbMonomials, 2);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -165,7 +141,7 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
}
}
if (!serendip) {
- fullMatrix<int> inner = gmshGenerateMonomialsTriangle(order-3);
+ fullMatrix<double> inner = gmshGenerateMonomialsTriangle(order-3);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
}
@@ -174,11 +150,11 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
return monomials;
}
-fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints)
+fullMatrix<double> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? order*4 : (order+1)*(order+1);
if (forSerendipPoints && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 2);
+ fullMatrix<double> monomials(nbMonomials, 2);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -209,7 +185,7 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoint
}
if (!forSerendipPoints) {
- fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2);
+ fullMatrix<double> inner = gmshGenerateMonomialsQuadrangle(order-2);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
}
@@ -227,10 +203,10 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoint
�. �.
*/
-fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order)
+fullMatrix<double> gmshGenerateMonomialsQuadSerendipity(int order)
{
int nbMonomials = order ? order*4 : 1;
- fullMatrix<int> monomials(nbMonomials, 2);
+ fullMatrix<double> monomials(nbMonomials, 2);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -267,7 +243,7 @@ fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order)
//KH : caveat : node coordinates are not yet coherent with node numbering associated
// to numbering of principal vertices of face !!!!
-fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
+fullMatrix<double> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
{
int nbMonomials =
(serendip ?
@@ -275,7 +251,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
(order + 1) * (order + 2) * (order + 3) / 6);
if (serendip && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -315,7 +291,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
}
if (order > 2) {
- fullMatrix<int> dudv = gmshGenerateMonomialsTriangle(order - 3);
+ fullMatrix<double> dudv = gmshGenerateMonomialsTriangle(order - 3);
dudv.add(1);
for (int iface = 0; iface < 4; ++iface) {
@@ -340,7 +316,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
}
if (!serendip && order > 3) {
- fullMatrix<int> inner = gmshGenerateMonomialsTetrahedron(order - 4);
+ fullMatrix<double> inner = gmshGenerateMonomialsTetrahedron(order - 4);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
}
@@ -350,12 +326,12 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
return monomials;
}
-fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
+fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 :
(order + 1) * (order + 1)*(order + 2)/2;
if (forSerendipPoints && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -400,10 +376,10 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
}
if (!forSerendipPoints) {
- fullMatrix<int> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
+ fullMatrix<double> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
dudvQ.add(1);
- fullMatrix<int> dudvT;
+ fullMatrix<double> dudvT;
if (order > 2) dudvT = gmshGenerateMonomialsTriangle(order - 3);
dudvT.add(1);
@@ -411,7 +387,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
int i0, i1, i2;
i0 = MPrism::faces_prism(iface, 0);
i1 = MPrism::faces_prism(iface, 1);
- fullMatrix<int> dudv;
+ fullMatrix<double> dudv;
if (MPrism::faces_prism(iface, 3) != -1) {
i2 = MPrism::faces_prism(iface, 3);
dudv.setAsProxy(dudvQ);
@@ -439,8 +415,8 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
}
if (order > 2) {
- fullMatrix<int> triMonomials = gmshGenerateMonomialsTriangle(order - 3);
- fullMatrix<int> lineMonomials = gmshGenerateMonomialsLine(order - 2);
+ fullMatrix<double> triMonomials = gmshGenerateMonomialsTriangle(order - 3);
+ fullMatrix<double> lineMonomials = gmshGenerateMonomialsLine(order - 2);
for (int i = 0; i < triMonomials.size1(); ++i) {
for (int j = 0; j < lineMonomials.size1(); ++j, ++index) {
@@ -456,10 +432,10 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
return monomials;
}
-fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order)
+fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order)
{
int nbMonomials = order ? 6 + (order-1) * 9 : 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -528,13 +504,12 @@ fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order)
return monomials;
}
-
-fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
+fullMatrix<double> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 :
(order+1)*(order+1)*(order+1);
if (forSerendipPoints && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -587,7 +562,7 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
}
if (!forSerendipPoints) {
- fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
+ fullMatrix<double> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
dudv.add(1);
for (int iface = 0; iface < 6; ++iface) {
@@ -611,7 +586,7 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
}
}
- fullMatrix<int> inner = gmshGenerateMonomialsHexahedron(order - 2);
+ fullMatrix<double> inner = gmshGenerateMonomialsHexahedron(order - 2);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
}
@@ -620,11 +595,10 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
return monomials;
}
-
-fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order)
+fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order)
{
int nbMonomials = order ? 8 + (order-1) * 12 : 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -691,13 +665,12 @@ fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order)
return monomials;
}
-
-fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
+fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 5 + (order-1)*8 :
(order+1)*((order+1)+1)*(2*(order+1)+1)/6;
if (forSerendipPoints && !order) nbMonomials = 1;
- fullMatrix<int> monomials(nbMonomials, 3);
+ fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
@@ -738,10 +711,10 @@ fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
}
if (!forSerendipPoints) {
- fullMatrix<int> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
+ fullMatrix<double> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
dudvQ.add(1);
- fullMatrix<int> dudvT;
+ fullMatrix<double> dudvT;
if (order > 2) dudvT = gmshGenerateMonomialsTriangle(order - 3);
dudvT.add(1);
@@ -749,7 +722,7 @@ fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
int i0, i1, i2;
i0 = MPyramid::faces_pyramid(iface, 0);
i1 = MPyramid::faces_pyramid(iface, 1);
- fullMatrix<int> dudv;
+ fullMatrix<double> dudv;
if (MPyramid::faces_pyramid(iface, 3) != -1) {
i2 = MPyramid::faces_pyramid(iface, 3);
dudv.setAsProxy(dudvQ);
@@ -777,7 +750,7 @@ fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
}
if (order > 2) {
- fullMatrix<int> inner = gmshGenerateMonomialsPyramid(order - 3);
+ fullMatrix<double> inner = gmshGenerateMonomialsPyramid(order - 3);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
}
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 15a4501d97fcca3ed589e4b3e777e1b127100fa6..59814e24099259095ed8c1df74033f0e7f01f7ae 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -20,32 +20,32 @@
fullMatrix<double> gmshGeneratePointsLine(int order);
-fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip);
-fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip);
+fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip = false);
+fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip = false);
-fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip);
-fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip);
-fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip);
+fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip = false);
+fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip = false);
+fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip = false);
-fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip);
+fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip = false);
// Monomial exponents
-fullMatrix<int> gmshGenerateMonomialsLine(int order, bool serendip = false);
+fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip = false);
-fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
-fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints = false);
-fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order);
+fullMatrix<double> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
+fullMatrix<double> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints = false);
+fullMatrix<double> gmshGenerateMonomialsQuadSerendipity(int order);
-fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false);
-fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints = false);
-fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order);
-fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
-fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order);
+fullMatrix<double> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false);
+fullMatrix<double> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints = false);
+fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order);
+fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
+fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order);
-fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
-//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //TODO
+fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
+//fullMatrix<double> gmshGenerateMonomialsPyramidSerendipity(int order); //TODO
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index dadc221bcab6023822461bd6ba0b64bf7164881b..3f985904f18b423f3c382839b1fd4c2d9fd7019f 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -75,41 +75,33 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
{
- fullMatrix<int> monom;
switch (parentType) {
case TYPE_PNT :
- monom = gmshGenerateMonomialsLine(0);
+ monomials = gmshGenerateMonomialsLine(0);
break;
case TYPE_LIN :
- monom = gmshGenerateMonomialsLine(order);
+ monomials = gmshGenerateMonomialsLine(order);
break;
case TYPE_TRI :
- monom = gmshGenerateMonomialsTriangle(order, serendip);
+ monomials = gmshGenerateMonomialsTriangle(order, serendip);
break;
case TYPE_QUA :
- monom = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
+ monomials = serendip ? gmshGenerateMonomialsQuadSerendipity(order) :
gmshGenerateMonomialsQuadrangle(order);
break;
case TYPE_TET :
- monom = gmshGenerateMonomialsTetrahedron(order, serendip);
+ monomials = gmshGenerateMonomialsTetrahedron(order, serendip);
break;
case TYPE_PRI :
- monom = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
+ monomials = serendip ? gmshGenerateMonomialsPrismSerendipity(order) :
gmshGenerateMonomialsPrism(order);
break;
case TYPE_HEX :
- monom = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
+ monomials = serendip ? gmshGenerateMonomialsHexaSerendipity(order) :
gmshGenerateMonomialsHexahedron(order);
break;
}
- monomials.resize(monom.size1(), monom.size2());
- for (int i = 0; i < monom.size1(); ++i) {
- for (int j = 0; j < monom.size2(); ++j) {
- monomials(i, j) = static_cast<double>(monom(i, j));
- }
- }
-
coefficients = generateLagrangeMonomialCoefficients(monomials, points);
}
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 1716ca129f20afaa4765759d540e6882e99433fa..c91fdfc94919ab157d456d7f891b57c01445f8f7 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -1044,7 +1044,7 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
// DI_Line methods --------------------------------------------------------------------------------
const nodalBasis* DI_Line::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
- int tag = MElement::getTag(TYPE_LIN, order);
+ int tag = ElementType::getTag(TYPE_LIN, order);
return BasisFactory::getNodalBasis(tag);
}
@@ -1065,7 +1065,7 @@ void DI_Line::computeIntegral() {
const nodalBasis* DI_Triangle::getFunctionSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
- int tag = MElement::getTag(TYPE_TRI, order);
+ int tag = ElementType::getTag(TYPE_TRI, order);
return BasisFactory::getNodalBasis(tag);
}
void DI_Triangle::computeIntegral() {
@@ -1089,7 +1089,7 @@ double DI_Triangle::quality() const {
// DI_Quad methods --------------------------------------------------------------------------------
const nodalBasis* DI_Quad::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
- int tag = MElement::getTag(TYPE_QUA, order);
+ int tag = ElementType::getTag(TYPE_QUA, order);
return BasisFactory::getNodalBasis(tag);
}
@@ -1113,7 +1113,7 @@ void DI_Quad::computeIntegral() {
// DI_Tetra methods -------------------------------------------------------------------------------
const nodalBasis* DI_Tetra::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
- int tag = MElement::getTag(TYPE_TET, order);
+ int tag = ElementType::getTag(TYPE_TET, order);
return BasisFactory::getNodalBasis(tag);
}
@@ -1128,7 +1128,7 @@ double DI_Tetra::quality() const {
// Hexahedron methods -----------------------------------------------------------------------------
const nodalBasis* DI_Hexa::getFunctionSpace(int o) const{
int order = (o == -1) ? getPolynomialOrder() : o;
- int tag = MElement::getTag(TYPE_HEX, order);
+ int tag = ElementType::getTag(TYPE_HEX, order);
return BasisFactory::getNodalBasis(tag);
}
void DI_Hexa::computeIntegral() {