diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 000efdbd6ddbafbf8b7ce82f89a60bdffe31bc52..999913e9991d453898b915d36269f3aea945c736 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -96,6 +96,31 @@
#define MSH_QUA_12 39
#define MSH_QUA_16I 40
#define MSH_QUA_20 41
+#define MSH_TRI_28 42
+#define MSH_TRI_36 43
+#define MSH_TRI_45 44
+#define MSH_TRI_55 45
+#define MSH_TRI_66 46
+#define MSH_QUA_49 47
+#define MSH_QUA_64 48
+#define MSH_QUA_81 49
+#define MSH_QUA_100 50
+#define MSH_QUA_121 51
+#define MSH_TRI_18 52
+#define MSH_TRI_21I 53
+#define MSH_TRI_24 54
+#define MSH_TRI_27 55
+#define MSH_TRI_30 56
+#define MSH_QUA_24 57
+#define MSH_QUA_28 58
+#define MSH_QUA_32 59
+#define MSH_QUA_36I 60
+#define MSH_QUA_40 61
+#define MSH_LIN_7 62
+#define MSH_LIN_8 63
+#define MSH_LIN_9 64
+#define MSH_LIN_10 65
+#define MSH_LIN_11 66
#define MSH_NUM_TYPE 35
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 0913b8d2266ced55a2b5de5e0a6bc345a71d834f..d4cdaeb0be59be90b1968a6732fbfc8ce89dbc3b 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -385,6 +385,8 @@ bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
bool GFaceCompound::checkOrientation(int iter) const
{
+ return true;
+
//Only check orientation for stl files (1 patch)
// if(_compound.size() > 1.0) return true;
@@ -496,7 +498,7 @@ bool GFaceCompound::parametrize() const
//Multiscale Laplace parametrization
//-----------------
else if (_mapping == MULTISCALE){
- printf("multiscale exit \n");
+ // printf("multiscale exit \n");
std::vector<MElement*> _elements;
for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
@@ -523,6 +525,8 @@ bool GFaceCompound::parametrize() const
buildOct();
+ printStuff();
+
if (!checkOrientation(0)){
Msg::Info("--- Parametrization switched to convex combination map");
coordinates.clear();
@@ -534,7 +538,6 @@ bool GFaceCompound::parametrize() const
buildOct();
}
- printStuff();
double AR = checkAspectRatio();
if (floor(AR) > AR_MAX){
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 45002386829f1cbec8486eda2dcb524f7c0b2670..f0d2970503dd53630ad9e443d1527a014744abd7 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -348,15 +348,25 @@ int GModel::readMSH(const std::string &name)
if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary,
&numVertices) != 5)
return 0;
- if(numVertices != MElement::getInfoMSH(type)) return 0;
+ if(numVertices != MElement::getInfoMSH(type)) {
+ // printf("coucou1\n");
+ return 0;
+ }
}
else{
int numTags;
- if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
+ if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) {
+ // printf("coucou2\n");
+ return 0;
+ }
+ // printf("%d %d %d\n",num,type,numTags);
int numPartitions = 0;
for(int j = 0; j < numTags; j++){
int tag;
- if(fscanf(fp, "%d", &tag) != 1) return 0;
+ if(fscanf(fp, "%d", &tag) != 1) {
+ // printf("coucou3\n");
+ return 0;
+ }
if(j == 0) physical = tag;
else if(j == 1) elementary = tag;
else if(version < 2.2 && j == 2) partition = tag;
@@ -366,19 +376,35 @@ int GModel::readMSH(const std::string &name)
else if(j == numTags - 1) parent = tag;
}
if(!(numVertices = MElement::getInfoMSH(type))) {
- if(type != MSH_POLYG_ && type != MSH_POLYH_) return 0;
- if(fscanf(fp, "%d", &numVertices) != 1) return 0;
+ if(type != MSH_POLYG_ && type != MSH_POLYH_) {
+ // printf("coucou4\n");
+ return 0;
+ }
+ if(fscanf(fp, "%d", &numVertices) != 1){
+ // printf("coucou5\n");
+ return 0;
+ }
}
+
}
int *indices = new int[numVertices];
for(int j = 0; j < numVertices; j++)
- if(fscanf(fp, "%d", &indices[j]) != 1) return 0;
+ if(fscanf(fp, "%d", &indices[j]) != 1) {
+ // printf("coucou6\n");
+ return 0;
+ }
std::vector<MVertex*> vertices;
if(vertexVector.size()){
- if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0;
+ if(!getVertices(numVertices, indices, vertexVector, vertices)) {
+ // printf("coucou7\n");
+ return 0;
+ }
}
else{
- if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
+ if(!getVertices(numVertices, indices, vertexMap, vertices)) {
+ // printf("coucou8\n");
+ return 0;
+ }
}
MElement *e = createElementMSH(this, num, type, physical, elementary,
partition, vertices, elements, physicals);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 76f0a5b08feb9c420a71df5e5680113b5f21d36e..421d93724aa22dfef1115f37d0f9383c417d51a4 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -750,6 +750,11 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_LIN_4 : if(name) *name = "Line 4"; return 2 + 2;
case MSH_LIN_5 : if(name) *name = "Line 5"; return 2 + 3;
case MSH_LIN_6 : if(name) *name = "Line 6"; return 2 + 4;
+ case MSH_LIN_7 : if(name) *name = "Line 7"; return 2 + 5;
+ case MSH_LIN_8 : if(name) *name = "Line 8"; return 2 + 6;
+ case MSH_LIN_9 : if(name) *name = "Line 9"; return 2 + 7;
+ case MSH_LIN_10 : if(name) *name = "Line 10"; return 2 + 8;
+ case MSH_LIN_11 : if(name) *name = "Line 11"; return 2 + 9;
case MSH_TRI_3 : if(name) *name = "Triangle 3"; return 3;
case MSH_TRI_6 : if(name) *name = "Triangle 6"; return 3 + 3;
case MSH_TRI_9 : if(name) *name = "Triangle 9"; return 3 + 6;
@@ -758,12 +763,27 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_TRI_15 : if(name) *name = "Triangle 15"; return 3 + 9 + 3;
case MSH_TRI_15I: if(name) *name = "Triangle 15I"; return 3 + 12;
case MSH_TRI_21 : if(name) *name = "Triangle 21"; return 3 + 12 + 6;
+ case MSH_TRI_28 : if(name) *name = "Triangle 28"; return 3 + 15 + 10;
+ case MSH_TRI_36 : if(name) *name = "Triangle 36"; return 3 + 18 + 15;
+ case MSH_TRI_45 : if(name) *name = "Triangle 45"; return 3 + 21 + 21;
+ case MSH_TRI_55 : if(name) *name = "Triangle 55"; return 3 + 24 + 28;
+ case MSH_TRI_66 : if(name) *name = "Triangle 66"; return 3 + 27 + 36;
+ case MSH_TRI_18 : if(name) *name = "Triangle 18"; return 3 + 15;
+ case MSH_TRI_21I : if(name) *name = "Triangle 21I"; return 3 + 18;
+ case MSH_TRI_24 : if(name) *name = "Triangle 24"; return 3 + 21;
+ case MSH_TRI_27 : if(name) *name = "Triangle 27"; return 3 + 24;
+ case MSH_TRI_30 : if(name) *name = "Triangle 30"; return 3 + 27;
case MSH_QUA_4 : if(name) *name = "Quadrilateral 4"; return 4;
case MSH_QUA_8 : if(name) *name = "Quadrilateral 8"; return 4 + 4;
case MSH_QUA_9 : if(name) *name = "Quadrilateral 9"; return 9;
case MSH_QUA_16 : if(name) *name = "Quadrilateral 16"; return 16;
case MSH_QUA_25 : if(name) *name = "Quadrilateral 25"; return 25;
case MSH_QUA_36 : if(name) *name = "Quadrilateral 36"; return 36;
+ case MSH_QUA_49 : if(name) *name = "Quadrilateral 49"; return 49;
+ case MSH_QUA_64 : if(name) *name = "Quadrilateral 64"; return 64;
+ case MSH_QUA_81 : if(name) *name = "Quadrilateral 81"; return 81;
+ case MSH_QUA_100 : if(name) *name = "Quadrilateral 100"; return 100;
+ case MSH_QUA_121 : if(name) *name = "Quadrilateral 121"; return 121;
case MSH_QUA_12 : if(name) *name = "Quadrilateral 12"; return 12;
case MSH_QUA_16I : if(name) *name = "Quadrilateral 16I"; return 16;
case MSH_QUA_20 : if(name) *name = "Quadrilateral 20"; return 20;
@@ -811,6 +831,11 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_LIN_4: return new MLineN(v, num, part);
case MSH_LIN_5: return new MLineN(v, num, part);
case MSH_LIN_6: return new MLineN(v, num, part);
+ case MSH_LIN_7: return new MLineN(v, num, part);
+ case MSH_LIN_8: return new MLineN(v, num, part);
+ case MSH_LIN_9: return new MLineN(v, num, part);
+ case MSH_LIN_10: return new MLineN(v, num, part);
+ case MSH_LIN_11: return new MLineN(v, num, part);
case MSH_TRI_3: return new MTriangle(v, num, part);
case MSH_TRI_6: return new MTriangle6(v, num, part);
case MSH_TRI_9: return new MTriangleN(v, 3, num, part);
@@ -819,12 +844,22 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_TRI_15: return new MTriangleN(v, 4, num, part);
case MSH_TRI_15I:return new MTriangleN(v, 5, num, part);
case MSH_TRI_21: return new MTriangleN(v, 5, num, part);
+ case MSH_TRI_28: return new MTriangleN(v, 6, num, part);
+ case MSH_TRI_36: return new MTriangleN(v, 7, num, part);
+ case MSH_TRI_45: return new MTriangleN(v, 8, num, part);
+ case MSH_TRI_55: return new MTriangleN(v, 9, num, part);
+ case MSH_TRI_66: return new MTriangleN(v,10, num, part);
case MSH_QUA_4: return new MQuadrangle(v, num, part);
case MSH_QUA_8: return new MQuadrangle8(v, num, part);
case MSH_QUA_9: return new MQuadrangle9(v, num, part);
case MSH_QUA_16: return new MQuadrangleN(v, 3, num, part);
case MSH_QUA_25: return new MQuadrangleN(v, 4, num, part);
case MSH_QUA_36: return new MQuadrangleN(v, 5, num, part);
+ case MSH_QUA_49: return new MQuadrangleN(v, 6, num, part);
+ case MSH_QUA_64: return new MQuadrangleN(v, 7, num, part);
+ case MSH_QUA_81: return new MQuadrangleN(v, 8, num, part);
+ case MSH_QUA_100: return new MQuadrangleN(v, 9, num, part);
+ case MSH_QUA_121: return new MQuadrangleN(v, 10, num, part);
case MSH_POLYG_: return new MPolygon(v, num, part);
case MSH_TET_4: return new MTetrahedron(v, num, part);
case MSH_TET_10: return new MTetrahedron10(v, num, part);
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 2349f0a8324a3888c31b9a440514b9701a8b2486..b6cdc2fc7f795d6bdc6f854f789e6c500cab18f7 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -19,6 +19,11 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_LIN_4);
case 4: return &polynomialBases::find(MSH_LIN_5);
case 5: return &polynomialBases::find(MSH_LIN_6);
+ case 6: return &polynomialBases::find(MSH_LIN_7);
+ case 7: return &polynomialBases::find(MSH_LIN_8);
+ case 8: return &polynomialBases::find(MSH_LIN_9);
+ case 9: return &polynomialBases::find(MSH_LIN_10);
+ case 10: return &polynomialBases::find(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 9b8198958174fede4054711bd879406daba45f66..de78ae3daea133e0b000dd15610d260fd05dfedf 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -178,6 +178,11 @@ class MLineN : public MLine {
if(_vs.size() == 2) return MSH_LIN_4;
if(_vs.size() == 3) return MSH_LIN_5;
if(_vs.size() == 4) return MSH_LIN_6;
+ if(_vs.size() == 5) return MSH_LIN_7;
+ if(_vs.size() == 6) return MSH_LIN_8;
+ if(_vs.size() == 7) return MSH_LIN_9;
+ if(_vs.size() == 8) return MSH_LIN_10;
+ if(_vs.size() == 9) return MSH_LIN_11;
return 0;
}
};
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 80aaf08ee69b4cae5737a1d541f3bad6e8af342d..098e061b41ee72a5a046b441d785a5e5183d37b6 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -29,6 +29,11 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_QUA_12);
case 4: return &polynomialBases::find(MSH_QUA_16I);
case 5: return &polynomialBases::find(MSH_QUA_20);
+ case 6: return &polynomialBases::find(MSH_QUA_24);
+ case 7: return &polynomialBases::find(MSH_QUA_28);
+ case 8: return &polynomialBases::find(MSH_QUA_32);
+ case 9: return &polynomialBases::find(MSH_QUA_36I);
+ case 10: return &polynomialBases::find(MSH_QUA_40);
}
}
switch (order) {
@@ -37,6 +42,11 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_QUA_16);
case 4: return &polynomialBases::find(MSH_QUA_25);
case 5: return &polynomialBases::find(MSH_QUA_36);
+ case 6: return &polynomialBases::find(MSH_QUA_49);
+ case 7: return &polynomialBases::find(MSH_QUA_64);
+ case 8: return &polynomialBases::find(MSH_QUA_81);
+ case 9: return &polynomialBases::find(MSH_QUA_100);
+ case 10: return &polynomialBases::find(MSH_QUA_121);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
return 0;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 4cc751606b99a47f1a238924458175a50aa23cda..584ee6aa3d9385ce00786af2dddb001cb4d1e81d 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -375,6 +375,11 @@ class MQuadrangleN : public MQuadrangle {
if(_order==4 && _vs.size()+4==25) return MSH_QUA_25;
if(_order==5 && _vs.size()+4==20) return MSH_QUA_20;
if(_order==5 && _vs.size()+4==36) return MSH_QUA_36;
+ if(_order==6 && _vs.size()+4==49) return MSH_QUA_49;
+ if(_order==7 && _vs.size()+4==64) return MSH_QUA_64;
+ if(_order==8 && _vs.size()+4==81) return MSH_QUA_81;
+ if(_order==9 && _vs.size()+4==100) return MSH_QUA_100;
+ if(_order==10 && _vs.size()+4==121) return MSH_QUA_121;
return 0;
}
virtual void revert()
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 3aa43df13b83cc253dfbae7c88362d8405274076..3365ed4250fa82c9f62bdc6fff1e85238ae27be8 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -78,7 +78,12 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_TRI_9);
case 4: return &polynomialBases::find(MSH_TRI_12);
case 5: return &polynomialBases::find(MSH_TRI_15I);
- default: Msg::Error("Order %d triangle function space not implemented", order);
+ case 6: return &polynomialBases::find(MSH_TRI_18);
+ case 7: return &polynomialBases::find(MSH_TRI_21I);
+ case 8: return &polynomialBases::find(MSH_TRI_24);
+ case 9: return &polynomialBases::find(MSH_TRI_27);
+ case 10: return &polynomialBases::find(MSH_TRI_30);
+ default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
@@ -88,6 +93,11 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
case 3: return &polynomialBases::find(MSH_TRI_10);
case 4: return &polynomialBases::find(MSH_TRI_15);
case 5: return &polynomialBases::find(MSH_TRI_21);
+ case 6: return &polynomialBases::find(MSH_TRI_28);
+ case 7: return &polynomialBases::find(MSH_TRI_36);
+ case 8: return &polynomialBases::find(MSH_TRI_45);
+ case 9: return &polynomialBases::find(MSH_TRI_55);
+ case 10: return &polynomialBases::find(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 5568ac8c565dd6e025b2452edd84408f3d766c7e..08b124c9a9572e00c52cd79675ebf5bbf0325c92 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -266,7 +266,17 @@ class MTriangleN : public MTriangle {
if(_order == 4 && _vs.size() == 9) return 0;
if(_order == 4 && _vs.size() == 12) return 3;
if(_order == 5 && _vs.size() == 12) return 0;
- if(_order == 5 && _vs.size() == 18) return 6;
+ if(_order == 5 && _vs.size() == 18) return 6;
+ if(_order == 6 && _vs.size() == 25) return 10;
+ if(_order == 7 && _vs.size() == 33) return 12;
+ if(_order == 8 && _vs.size() == 42) return 15;
+ if(_order == 9 && _vs.size() == 52) return 21;
+ if(_order == 10 && _vs.size() == 63) return 28;
+ if(_order == 6 && _vs.size() == 15) return 0;
+ if(_order == 7 && _vs.size() == 18) return 0;
+ if(_order == 8 && _vs.size() == 21) return 0;
+ if(_order == 9 && _vs.size() == 24) return 0;
+ if(_order == 10 && _vs.size() == 27) return 0;
return 0;
}
virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
@@ -297,6 +307,11 @@ class MTriangleN : public MTriangle {
if(_order == 4 && _vs.size() == 12) return MSH_TRI_15;
if(_order == 5 && _vs.size() == 12) return MSH_TRI_15I;
if(_order == 5 && _vs.size() == 18) return MSH_TRI_21;
+ if(_order == 6 && _vs.size() == 25) return MSH_TRI_28;
+ if(_order == 7 && _vs.size() == 33) return MSH_TRI_36;
+ if(_order == 8 && _vs.size() == 42) return MSH_TRI_45;
+ if(_order == 9 && _vs.size() == 52) return MSH_TRI_55;
+ if(_order ==10 && _vs.size() == 63) return MSH_TRI_66;
return 0;
}
virtual void revert()
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 61a846627520c9d6c6fd3cba1a8d3e1983169a55..c1af7e7db62a766587074aaecd024d2fb34e14c3 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -416,7 +416,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
}
else{
std::vector<MVertex*> &vtcs = faceVertices[face];
- SPoint2 pts[20];
+ SPoint2 pts[100];
bool reparamOK = true;
if(!linear){
for(int k = 0; k < incomplete->getNumVertices(); k++)
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index c5270b06dfeef79ea92f152f36061814d0b000d5..fe0011e2b5e2dd622f8589cef12fe5034a49e1cc 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -270,12 +270,10 @@ void meshGEdge::operator() (GEdge *ge)
ge->setLength(length);
Points.clear();
- if(length == 0.){
- if(CTX::instance()->mesh.toleranceEdgeLength == 0.)
- Msg::Error("Curve %d has a zero length", ge->tag());
- else
- Msg::Debug("Curve %d has a zero length", ge->tag());
- }
+ if(length == 0. && CTX::instance()->mesh.toleranceEdgeLength == 0.)
+ Msg::Error("Curve %d has a zero length", ge->tag());
+ else if (length == 0.)
+ Msg::Debug("Curve %d has a zero length", ge->tag());
if(length < CTX::instance()->mesh.toleranceEdgeLength)
ge->setTooSmall(true);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 537117d1f59d2bf304a3b99a5715cc184730ca45..8217d95b0bca7dbd957e1a835a74d84e34c85bd4 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1235,7 +1235,7 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
-int debugSurface = -1;
+int debugSurface = -100;
void meshGFace::operator() (GFace *gf)
{
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 703320d751e97b0f2e9e4c9ebe0a2ba189638aac..54b150ecf14ecc813a703df959dc8cfdf6d5d7a7 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -375,9 +375,9 @@ double qmTriangleAngles (MTriangle *e) {
double c;
prosca(v1,v2,&c);
double x = acos(c)-M_PI/3;
- printf("Angle %g ", (x+M_PI/3)/M_PI*180);
+ // printf("Angle %g ", (x+M_PI/3)/M_PI*180);
double quality = (atan(a*(x+M_PI/9)) + atan(a*(M_PI/9-x)))/den;
- printf("Quality %g\n",quality);
+ // printf("Quality %g\n",quality);
worst_quality = std::min(worst_quality, quality);
}
//printf("\n");
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 3705ace1f3c5d9d64fe5e6faaaba7977417acea4..43fed498c7353ba87a73e17a213f6e726eaa8db2 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -716,7 +716,9 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
(const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{
if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
+ Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
+ monomial.size1(),point.size1(),
+ monomial.size2(),point.size2() );
return fullMatrix<double>(1, 1);
}
@@ -934,6 +936,36 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = generate1DPoints(5);
generate1dVertexClosure(F.closures);
break;
+ case MSH_LIN_7:
+ F.numFaces = 2;
+ F.monomials = generate1DMonomials(6);
+ F.points = generate1DPoints(6);
+ generate1dVertexClosure(F.closures);
+ break;
+ case MSH_LIN_8:
+ F.numFaces = 2;
+ F.monomials = generate1DMonomials(7);
+ F.points = generate1DPoints(7);
+ generate1dVertexClosure(F.closures);
+ break;
+ case MSH_LIN_9:
+ F.numFaces = 2;
+ F.monomials = generate1DMonomials(8);
+ F.points = generate1DPoints(8);
+ generate1dVertexClosure(F.closures);
+ break;
+ case MSH_LIN_10:
+ F.numFaces = 2;
+ F.monomials = generate1DMonomials(9);
+ F.points = generate1DPoints(9);
+ generate1dVertexClosure(F.closures);
+ break;
+ case MSH_LIN_11:
+ F.numFaces = 2;
+ F.monomials = generate1DMonomials(10);
+ F.points = generate1DPoints(10);
+ generate1dVertexClosure(F.closures);
+ break;
case MSH_TRI_3 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(1);
@@ -982,6 +1014,66 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(5, false);
generate2dEdgeClosure(F.closures, 5);
break;
+ case MSH_TRI_28 :
+ F.numFaces = 3;
+ F.monomials = generatePascalTriangle(6);
+ F.points = gmshGeneratePointsTriangle(6, false);
+ generate2dEdgeClosure(F.closures, 6);
+ break;
+ case MSH_TRI_36 :
+ F.numFaces = 3;
+ F.monomials = generatePascalTriangle(7);
+ F.points = gmshGeneratePointsTriangle(7, false);
+ generate2dEdgeClosure(F.closures, 7);
+ break;
+ case MSH_TRI_45 :
+ F.numFaces = 3;
+ F.monomials = generatePascalTriangle(8);
+ F.points = gmshGeneratePointsTriangle(8, false);
+ generate2dEdgeClosure(F.closures, 8);
+ break;
+ case MSH_TRI_55 :
+ F.numFaces = 3;
+ F.monomials = generatePascalTriangle(9);
+ F.points = gmshGeneratePointsTriangle(9, false);
+ generate2dEdgeClosure(F.closures, 9);
+ break;
+ case MSH_TRI_66 :
+ F.numFaces = 3;
+ F.monomials = generatePascalTriangle(10);
+ F.points = gmshGeneratePointsTriangle(10, false);
+ generate2dEdgeClosure(F.closures, 10);
+ break;
+ case MSH_TRI_18 :
+ F.numFaces = 3;
+ F.monomials = generatePascalSerendipityTriangle(6);
+ F.points = gmshGeneratePointsTriangle(6, true);
+ generate2dEdgeClosure(F.closures, 6);
+ break;
+ case MSH_TRI_21I :
+ F.numFaces = 3;
+ F.monomials = generatePascalSerendipityTriangle(7);
+ F.points = gmshGeneratePointsTriangle(7, true);
+ generate2dEdgeClosure(F.closures, 7);
+ break;
+ case MSH_TRI_24 :
+ F.numFaces = 3;
+ F.monomials = generatePascalSerendipityTriangle(8);
+ F.points = gmshGeneratePointsTriangle(8, true);
+ generate2dEdgeClosure(F.closures, 8);
+ break;
+ case MSH_TRI_27 :
+ F.numFaces = 3;
+ F.monomials = generatePascalSerendipityTriangle(9);
+ F.points = gmshGeneratePointsTriangle(9, true);
+ generate2dEdgeClosure(F.closures, 9);
+ break;
+ case MSH_TRI_30 :
+ F.numFaces = 3;
+ F.monomials = generatePascalSerendipityTriangle(10);
+ F.points = gmshGeneratePointsTriangle(10, true);
+ generate2dEdgeClosure(F.closures, 10);
+ break;
case MSH_TET_4 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(1);
@@ -1054,6 +1146,36 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(5,false);
generate2dEdgeClosure(F.closures, 5, 4);
break;
+ case MSH_QUA_49 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuad(6);
+ F.points = gmshGeneratePointsQuad(6,false);
+ generate2dEdgeClosure(F.closures, 6, 4);
+ break;
+ case MSH_QUA_64 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuad(7);
+ F.points = gmshGeneratePointsQuad(7,false);
+ generate2dEdgeClosure(F.closures, 7, 4);
+ break;
+ case MSH_QUA_81 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuad(8);
+ F.points = gmshGeneratePointsQuad(8,false);
+ generate2dEdgeClosure(F.closures, 8, 4);
+ break;
+ case MSH_QUA_100 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuad(9);
+ F.points = gmshGeneratePointsQuad(9,false);
+ generate2dEdgeClosure(F.closures, 9, 4);
+ break;
+ case MSH_QUA_121 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuad(10);
+ F.points = gmshGeneratePointsQuad(10,false);
+ generate2dEdgeClosure(F.closures, 10, 4);
+ break;
case MSH_QUA_8 :
F.numFaces = 4;
F.monomials = generatePascalQuadSerendip(2);
@@ -1078,6 +1200,36 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(5,true);
generate2dEdgeClosure(F.closures, 5, 4);
break;
+ case MSH_QUA_24 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuadSerendip(6);
+ F.points = gmshGeneratePointsQuad(6,true);
+ generate2dEdgeClosure(F.closures, 6, 4);
+ break;
+ case MSH_QUA_28 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuadSerendip(7);
+ F.points = gmshGeneratePointsQuad(7,true);
+ generate2dEdgeClosure(F.closures, 7, 4);
+ break;
+ case MSH_QUA_32 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuadSerendip(8);
+ F.points = gmshGeneratePointsQuad(8,true);
+ generate2dEdgeClosure(F.closures, 8, 4);
+ break;
+ case MSH_QUA_36I :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuadSerendip(9);
+ F.points = gmshGeneratePointsQuad(9,true);
+ generate2dEdgeClosure(F.closures, 9, 4);
+ break;
+ case MSH_QUA_40 :
+ F.numFaces = 4;
+ F.monomials = generatePascalQuadSerendip(10);
+ F.points = gmshGeneratePointsQuad(10,true);
+ generate2dEdgeClosure(F.closures, 10, 4);
+ break;
case MSH_PRI_6 : // first order
F.numFaces = 5;
F.monomials = generatePascalPrism(1);
@@ -1099,6 +1251,7 @@ const polynomialBasis &polynomialBases::find(int tag)
generate3dFaceClosure(F.closures, 1);
break;
}
+
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
// printf("Case: %d coeffs:\n",tag);
// for (int i = 0; i<F.coefficients.size1(); i++) {
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 3c4acab11ad53d4520fa7d60676fa8bfa0e8ef21..80cece67066ae08efb2630a70f1d26d208fbcc21 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -12,6 +12,57 @@
#include "fullMatrix.h"
#include <iostream>
+inline double pow_int (const double &a , const int &n) {
+ switch (n) {
+ case 0 : return 1.0;
+ case 1 : return a;
+ case 2 : return a*a;
+ case 3 : return a*a*a;
+ case 4 :
+ {
+ const double a2 = a*a;
+ return a2*a2;
+ }
+ case 5 :
+ {
+ const double a2 = a*a;
+ const double a3 = a*a*a;
+ return a2*a3;
+ }
+ case 6 :
+ {
+ const double a3 = a*a*a;
+ return a3*a3;
+ }
+ case 7 :
+ {
+ const double a3 = a*a*a;
+ return a3*a3*a;
+ }
+ case 8 :
+ {
+ const double a2 = a*a;
+ const double a4 = a2*a2;
+ return a4*a4;
+ }
+ case 9 :
+ {
+ const double a2 = a*a;
+ const double a4 = a2*a2;
+ return a4*a4*a;
+ }
+ case 10 :
+ {
+ const double a2 = a*a;
+ const double a4 = a2*a2;
+ return a4*a4*a2;
+ }
+ default :
+ return pow_int(a,n-1)*a;
+ }
+}
+
+
// presently those function spaces are only for simplices and quads;
// should be extended to other elements like hexes
class polynomialBasis
@@ -35,9 +86,9 @@ class polynomialBasis
inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
for (int j = 0; j < monomials.size1(); j++) {
- p[j] = pow(u, (int)monomials(j, 0));
- if (monomials.size2() > 1) p[j] *= pow(v, (int)monomials(j, 1));
- if (monomials.size2() > 2) p[j] *= pow(w, (int)monomials(j, 2));
+ p[j] = pow_int(u, (int)monomials(j, 0));
+ if (monomials.size2() > 1) p[j] *= pow_int(v, (int)monomials(j, 1));
+ if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
}
}
inline void f(double u, double v, double w, double *sf) const
@@ -62,7 +113,7 @@ class polynomialBasis
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0);
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0);
}
}
break;
@@ -74,12 +125,12 @@ class polynomialBasis
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0) *
- pow(v, monomials(j, 1));
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1));
if (monomials(j, 1) > 0)
grads[i][1] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1) - 1) * monomials(j, 1);
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1);
}
}
break;
@@ -91,19 +142,19 @@ class polynomialBasis
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 0)
grads[i][0] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0) *
- pow(v, monomials(j, 1)) *
- pow(w, monomials(j, 2));
+ pow_int(u, (int)(monomials(j, 0) - 1)) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2));
if (monomials(j, 1) > 0)
grads[i][1] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1) - 1) * monomials(j, 1) *
- pow(w, monomials(j, 2));
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)(monomials(j, 1) - 1)) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2));
if (monomials(j, 2) > 0)
grads[i][2] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1)) *
- pow(w, monomials(j, 2) - 1) * monomials(j, 2);
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)(monomials(j, 2) - 1)) * monomials(j, 2);
}
}
break;
@@ -120,7 +171,7 @@ class polynomialBasis
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1) // second derivative !=0
- hess[i][0][0] += coefficients(i, j) * pow(u, monomials(j, 0) - 2) *
+ hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
monomials(j, 0) * (monomials(j, 0) - 1);
}
}
@@ -132,16 +183,16 @@ class polynomialBasis
hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1) // second derivative !=0
- hess[i][0][0] += coefficients(i, j) * pow(u, monomials(j, 0) - 2) *
- monomials(j, 0) * (monomials(j, 0) - 1) * pow(v, monomials(j, 1));
+ hess[i][0][0] += coefficients(i, j) * pow_int(u, (int)(monomials(j, 0) - 2)) *
+ monomials(j, 0) * (monomials(j, 0) - 1) * pow_int(v, (int)monomials(j, 1));
if ((monomials(j, 1) > 0) && (monomials(j, 0) > 0))
hess[i][0][1] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0) *
- pow(v, monomials(j, 1) - 1) * monomials(j, 1);
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1);
if (monomials(j, 1) > 1)
hess[i][1][1] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1) - 1);
}
hess[i][1][0] = hess[i][0][1];
}
@@ -154,35 +205,35 @@ class polynomialBasis
for(int j = 0; j < coefficients.size2(); j++){
if (monomials(j, 0) > 1)
hess[i][0][0] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
- pow(v, monomials(j, 1)) *
- pow(w, monomials(j, 2));
+ pow_int(u, (int)monomials(j, 0) - 2) * monomials(j, 0) * (monomials(j, 0)-1) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2));
if ((monomials(j, 0) > 0) && (monomials(j, 1) > 0))
hess[i][0][1] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0) *
- pow(v, monomials(j, 1) - 1) * monomials(j, 1) *
- pow(w, monomials(j, 2));
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2));
if ((monomials(j, 0) > 0) && (monomials(j, 2) > 0))
hess[i][0][2] += coefficients(i, j) *
- pow(u, monomials(j, 0) - 1) * monomials(j, 0) *
- pow(v, monomials(j, 1)) *
- pow(w, monomials(j, 2) - 1) * monomials(j, 2);
+ pow_int(u, (int)monomials(j, 0) - 1) * monomials(j, 0) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
if (monomials(j, 1) > 1)
hess[i][1][1] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
- pow(w, monomials(j, 2));
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 2) * monomials(j, 1) * (monomials(j, 1)-1) *
+ pow_int(w, (int)monomials(j, 2));
if ((monomials(j, 1) > 0) && (monomials(j, 2) > 0))
hess[i][1][2] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1) - 1) * monomials(j, 1) *
- pow(w, monomials(j, 2) - 1) * monomials(j, 2);
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1) - 1) * monomials(j, 1) *
+ pow_int(w, (int)monomials(j, 2) - 1) * monomials(j, 2);
if (monomials(j, 2) > 1)
hess[i][2][2] += coefficients(i, j) *
- pow(u, monomials(j, 0)) *
- pow(v, monomials(j, 1)) *
- pow(w, monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
+ pow_int(u, (int)monomials(j, 0)) *
+ pow_int(v, (int)monomials(j, 1)) *
+ pow_int(w, (int)monomials(j, 2) - 2) * monomials(j, 2) * (monomials(j, 2) - 1);
}
hess[i][1][0] = hess[i][0][1];
hess[i][2][0] = hess[i][0][2];
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 6d878a2434caa7c1f9543a142270ae4b9b00a6c2..b966904ddd25df3a053001234111e7c0ba4c6c60 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -101,6 +101,7 @@ class dofManager{
linearSystem<dataMat> *_current;
public:
+ dofManager() : _current(0) { }
dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) {
_linearSystems.insert(std::make_pair("A", l1));