diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d53b76d3d9bcab162908de654fcc37b6502d2995..650a2fbd0f15984b167cc51d3d06ea3e926c0f9c 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.31 2007-02-27 17:15:46 remacle Exp $
+// $Id: MElement.cpp,v 1.32 2007-03-14 15:47:49 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -375,6 +375,7 @@ void MTriangle::circumcenterXY(double *res) const
}
int MTriangleN::getNumFacesRep(){ return 1; }
+
MFace MTriangleN::getFaceRep(int num)
{
return MFace(_v[0],_v[1],_v[2]);
@@ -389,3 +390,222 @@ int MTriangleN::getNumFaceVertices(){
if (_order == 5 && _vs.size() == 18) return 6;
throw;
}
+
+int P1[3][2] = {
+ {0,0},
+ {1,0},
+ {0,1}
+};
+
+int P2[6][2] = {
+ {0,0},
+ {1,0},
+ {0,1},
+ {2,0},
+ {0,2},
+ {1,1}
+};
+
+int P3[9][2] = {
+ {0,0},
+ {1,0},
+ {0,1},
+ {2,0},
+ {0,2},
+ {3,0},
+ {0,3},
+ {2,1},
+ {1,2}
+};
+int P4[12][2] = {
+ {0,0},
+ {1,0},
+ {0,1},
+ {2,0},
+ {0,2},
+ {3,0},
+ {0,3},
+ {4,0},
+ {0,4},
+ {3,1},
+ {1,3},
+ {2,2}
+};
+
+int P5[15][2] = {
+ {0,0},
+ {1,0},
+ {0,1},
+ {2,0},
+ {0,2},
+ {3,0},
+ {0,3},
+ {4,0},
+ {0,4},
+ {5,0},
+ {0,5},
+ {4,1},
+ {3,2},
+ {2,3},
+ {1,4}
+};
+
+double coef1[3][3]={
+{ 1.00000000, -1.00000000, -1.00000000},
+{ 0.00000000, 1.00000000, 0.00000000},
+{ 0.00000000, 0.00000000, 1.00000000}
+};
+
+double coef2[6][6]={
+{ 1.00000000, -3.00000000, -3.00000000, 2.00000000, 2.00000000, 4.00000000},
+{ 0.00000000, -1.00000000, 0.00000000, 2.00000000, -0.00000000, -0.00000000},
+{ 0.00000000, 0.00000000, -1.00000000, -0.00000000, 2.00000000, -0.00000000},
+{ 0.00000000, 4.00000000, 0.00000000, -4.00000000, -0.00000000, -4.00000000},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, 4.00000000},
+{ 0.00000000, 0.00000000, 4.00000000, -0.00000000, -4.00000000, -4.00000000}
+};
+
+double coef3[9][9]={
+{ 1.00000000, -5.50000000, -5.50000000, 9.00000000, 9.00000000, -4.50000000, -4.50000000, 4.50000000, 4.50000000},
+{ 0.00000000, 1.00000000, 0.00000000, -4.50000000, -0.00000000, 4.50000000, -0.00000000, -0.00000000, -0.00000000},
+{ 0.00000000, 0.00000000, 1.00000000, -0.00000000, -4.50000000, -0.00000000, 4.50000000, -0.00000000, -0.00000000},
+{ 0.00000000, 9.00000000, 0.00000000, -22.50000000, -0.00000000, 13.50000000, 0.00000000, 4.50000000, -9.00000000},
+{ 0.00000000, -4.50000000, -0.00000000, 18.00000000, 0.00000000, -13.50000000, -0.00000000, -9.00000000, 4.50000000},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, 9.00000000, -4.50000000},
+{ 0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, -4.50000000, 9.00000000},
+{ 0.00000000, 0.00000000, -4.50000000, -0.00000000, 18.00000000, -0.00000000, -13.50000000, 4.50000000, -9.00000000},
+{ 0.00000000, 0.00000000, 9.00000000, -0.00000000, -22.50000000, -0.00000000, 13.50000000, -9.00000000, 4.50000000}
+};
+double coef4[12][12]={
+{ 1.00000000, -8.33333333, -8.33333333, 23.33333333, 23.33333333, -26.66666667, -26.66666667, 10.66666667, 10.66666667, 9.33333333, 9.33333333, -2.66666667},
+{ 0.00000000, -1.00000000, 0.00000000, 7.33333333, -0.00000000, -16.00000000, 0.00000000, 10.66666667, -0.00000000, -0.00000000, -0.00000000, -0.00000000},
+{ 0.00000000, 0.00000000, -1.00000000, -0.00000000, 7.33333333, -0.00000000, -16.00000000, -0.00000000, 10.66666667, -0.00000000, -0.00000000, -0.00000000},
+{ 0.00000000, 16.00000000, -0.00000000, -69.33333333, 0.00000000, 96.00000000, -0.00000000, -42.66666667, 0.00000000, -5.33333333, -16.00000000, 21.33333333},
+{ 0.00000000, -12.00000000, 0.00000000, 76.00000000, -0.00000000, -128.00000000, 0.00000000, 64.00000000, -0.00000000, 12.00000000, 12.00000000, -40.00000000},
+{ 0.00000000, 5.33333333, -0.00000000, -37.33333333, 0.00000000, 74.66666667, -0.00000000, -42.66666667, 0.00000000, -16.00000000, -5.33333333, 21.33333333},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, 16.00000000, 5.33333333, -21.33333333},
+{ 0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, -12.00000000, -12.00000000, 40.00000000},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, 5.33333333, 16.00000000, -21.33333333},
+{ 0.00000000, 0.00000000, 5.33333333, -0.00000000, -37.33333333, -0.00000000, 74.66666667, -0.00000000, -42.66666667, -5.33333333, -16.00000000, 21.33333333},
+{ 0.00000000, 0.00000000, -12.00000000, -0.00000000, 76.00000000, -0.00000000, -128.00000000, -0.00000000, 64.00000000, 12.00000000, 12.00000000, -40.00000000},
+{ 0.00000000, 0.00000000, 16.00000000, -0.00000000, -69.33333333, -0.00000000, 96.00000000, -0.00000000, -42.66666667, -16.00000000, -5.33333333, 21.33333333}
+};
+double coef5[15][15]={
+{ 1.00000000, -11.41666667, -11.41666667, 46.87500000, 46.87500000, -88.54166667, -88.54166667, 78.12500000, 78.12500000, -26.04166667, -26.04166667, 10.41666667, 5.20833333, 5.20833333, 10.41666667},
+{ 0.00000000, 1.00000000, -0.00000000, -10.41666667, 0.00000000, 36.45833333, -0.00000000, -52.08333333, 0.00000000, 26.04166667, 0.00000000, -0.00000000, 0.00000000, -0.00000000, 0.00000000},
+{ 0.00000000, 0.00000000, 1.00000000, -0.00000000, -10.41666667, -0.00000000, 36.45833333, -0.00000000, -52.08333333, 0.00000000, 26.04166667, -0.00000000, 0.00000000, -0.00000000, 0.00000000},
+{ 0.00000000, 25.00000000, -0.00000000, -160.41666667, 0.00000000, 369.79166667, -0.00000000, -364.58333333, 0.00000000, 130.20833333, -0.00000000, 6.25000000, -38.54166667, 60.41666667, -25.00000000},
+{ 0.00000000, -25.00000000, 0.00000000, 222.91666667, -0.00000000, -614.58333333, 0.00000000, 677.08333333, -0.00000000, -260.41666667, 0.00000000, -16.66666667, 95.83333333, -122.91666667, 25.00000000},
+{ 0.00000000, 16.66666667, -0.00000000, -162.50000000, 0.00000000, 510.41666667, -0.00000000, -625.00000000, 0.00000000, 260.41666667, -0.00000000, 25.00000000, -122.91666667, 95.83333333, -16.66666667},
+{ 0.00000000, -6.25000000, 0.00000000, 63.54166667, -0.00000000, -213.54166667, 0.00000000, 286.45833333, -0.00000000, -130.20833333, 0.00000000, -25.00000000, 60.41666667, -38.54166667, 6.25000000},
+{ 0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, 0.00000000, -0.00000000, 25.00000000, -60.41666667, 38.54166667, -6.25000000},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, 0.00000000, -25.00000000, 122.91666667, -95.83333333, 16.66666667},
+{ 0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, 0.00000000, -0.00000000, 16.66666667, -95.83333333, 122.91666667, -25.00000000},
+{ 0.00000000, 0.00000000, 0.00000000, -0.00000000, -0.00000000, -0.00000000, 0.00000000, -0.00000000, -0.00000000, 0.00000000, 0.00000000, -6.25000000, 38.54166667, -60.41666667, 25.00000000},
+{ 0.00000000, 0.00000000, -6.25000000, -0.00000000, 63.54166667, -0.00000000, -213.54166667, -0.00000000, 286.45833333, 0.00000000, -130.20833333, 6.25000000, -38.54166667, 60.41666667, -25.00000000},
+{ 0.00000000, 0.00000000, 16.66666667, -0.00000000, -162.50000000, -0.00000000, 510.41666667, -0.00000000, -625.00000000, 0.00000000, 260.41666667, -16.66666667, 95.83333333, -122.91666667, 25.00000000},
+{ 0.00000000, 0.00000000, -25.00000000, -0.00000000, 222.91666667, -0.00000000, -614.58333333, -0.00000000, 677.08333333, 0.00000000, -260.41666667, 25.00000000, -122.91666667, 95.83333333, -16.66666667},
+{ 0.00000000, 0.00000000, 25.00000000, -0.00000000, -160.41666667, -0.00000000, 369.79166667, -0.00000000, -364.58333333, 0.00000000, 130.20833333, -25.00000000, 60.41666667, -38.54166667, 6.25000000}
+};
+
+
+void GradGeomShapeFunctionP1 (double u, double v, double grads[6][2])
+{
+ for (int i=0;i<3;i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ for (int j=0;j<3;j++){
+ if (P1[j][0] > 0)grads[i][0] += coef1[i][j] * pow(u,P1[j][0] - 1 ) * pow(v,P1[j][1] ) ;
+ if (P1[j][1] > 0)grads[i][1] += coef1[i][j] * pow(u,P1[j][0] ) * pow(v,P1[j][1] -1 ) ;
+ }
+ }
+}
+void GradGeomShapeFunctionP2 (double u, double v, double grads[6][2])
+{
+ for (int i=0;i<6;i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ for (int j=0;j<6;j++){
+ if (P2[j][0] > 0)grads[i][0] += coef2[i][j] * pow(u,P2[j][0] - 1 ) * pow(v,P2[j][1] ) ;
+ if (P2[j][1] > 0)grads[i][1] += coef2[i][j] * pow(u,P2[j][0] ) * pow(v,P2[j][1] -1 ) ;
+ }
+ }
+}
+void GradGeomShapeFunctionP3 (double u, double v, double grads[9][2])
+{
+ for (int i=0;i<9;i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ for (int j=0;j<9;j++){
+ if (P3[j][0] > 0)grads[i][0] += coef3[i][j] * pow(u,P3[j][0] - 1 ) * pow(v,P3[j][1] ) ;
+ if (P3[j][1] > 0)grads[i][1] += coef3[i][j] * pow(u,P3[j][0] ) * pow(v,P3[j][1] -1 ) ;
+ }
+ }
+}
+void GradGeomShapeFunctionP4 (double u, double v, double grads[12][2])
+{
+ for (int i=0;i<12;i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ for (int j=0;j<12;j++){
+ if (P4[j][0] > 0)grads[i][0] += coef3[i][j] * pow(u,P4[j][0] - 1 ) * pow(v,P4[j][1] ) ;
+ if (P4[j][1] > 0)grads[i][1] += coef3[i][j] * pow(u,P4[j][0] ) * pow(v,P4[j][1] -1 ) ;
+ }
+ }
+}
+
+void GradGeomShapeFunctionP5 (double u, double v, double grads[15][2])
+{
+ for (int i=0;i<15;i++){
+ grads[i][0] = 0;
+ grads[i][1] = 0;
+ for (int j=0;j<15;j++){
+ if (P5[j][0] > 0)grads[i][0] += coef3[i][j] * pow(u,P5[j][0] - 1 ) * pow(v,P5[j][1] ) ;
+ if (P5[j][1] > 0)grads[i][1] += coef3[i][j] * pow(u,P5[j][0] ) * pow(v,P5[j][1] -1 ) ;
+ }
+ }
+}
+
+
+void MTriangle::jac ( int ord, MVertex *vs[] , double uu, double vv , double j[2][2])
+{
+ double grads[256][2];
+ switch (ord)
+ {
+ case 1:
+ GradGeomShapeFunctionP1 ( uu , vv , grads );break;
+ case 2:
+ GradGeomShapeFunctionP2 ( uu , vv , grads );break;
+ case 3:
+ GradGeomShapeFunctionP3 ( uu , vv , grads );break;
+ case 4:
+ GradGeomShapeFunctionP4 ( uu , vv , grads );break;
+ case 5:
+ GradGeomShapeFunctionP5 ( uu , vv , grads );break;
+ default:
+ throw;
+ }
+ j[0][0] = 0 ; for (int i=0;i<3;i++)j[0][0] += grads [i][0] * _v[i] -> x() ;
+ j[1][0] = 0 ; for (int i=0;i<3;i++)j[1][0] += grads [i][1] * _v[i] -> x() ;
+ j[0][1] = 0 ; for (int i=0;i<3;i++)j[0][1] += grads [i][0] * _v[i] -> y() ;
+ j[1][1] = 0 ; for (int i=0;i<3;i++)j[1][1] += grads [i][1] * _v[i] -> y() ;
+ for (int i=3;i<3*ord;i++)j[0][0] += grads [i][0] * vs[i-3] -> x() ;
+ for (int i=3;i<3*ord;i++)j[1][0] += grads [i][1] * vs[i-3] -> x() ;
+ for (int i=3;i<3*ord;i++)j[0][1] += grads [i][0] * vs[i-3] -> y() ;
+ for (int i=3;i<3*ord;i++)j[1][1] += grads [i][1] * vs[i-3] -> y() ;
+}
+
+void MTriangleN::jac ( double uu, double vv , double j[2][2])
+{
+ MTriangle::jac (_order,&(*(_vs.begin())),uu,vv,j);
+}
+
+void MTriangle6::jac ( double uu, double vv , double j[2][2])
+{
+ MTriangle::jac (2,_vs,uu,vv,j);
+}
+void MTriangle::jac ( double uu, double vv , double j[2][2])
+{
+ jac (1,0,uu,vv,j);
+}
+
diff --git a/Geo/MElement.h b/Geo/MElement.h
index a640b6b6eb3bc5ff3f119d237fca5e1d6eda0b34..73971a5302a7bcf7e49554b51600b7edb26ba311 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -288,6 +288,7 @@ class MLineN : public MLine {
class MTriangle : public MElement {
protected:
MVertex *_v[3];
+ virtual void jac ( int order, MVertex *v[], double u, double v , double jac[2][2]) ;
public :
MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
: MElement(num, part)
@@ -342,6 +343,7 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
+ virtual void jac ( double u, double v , double j[2][2]) ;
};
class MTriangle6 : public MTriangle {
@@ -401,6 +403,7 @@ class MTriangle6 : public MTriangle {
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
+ virtual void jac ( double u, double v , double j[2][2]) ;
};
class MTriangleN : public MTriangle {
@@ -464,6 +467,7 @@ class MTriangleN : public MTriangle {
inv.insert (inv.begin(),_vs.rbegin(),_vs.rend());
_vs = inv;
}
+ virtual void jac ( double u, double v , double jac[2][2]) ;
};
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 67aea10db72f3ab6f4c0ae2cc4eb9cba9594cd32..4d76f08a20a3172b555054b2bceae2c6dbbd43a3 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.8 2007-03-13 16:11:43 remacle Exp $
+// $Id: HighOrder.cpp,v 1.9 2007-03-14 15:47:49 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -626,6 +626,20 @@ bool straightLine ( std::vector<MVertex *> &l , MVertex *n1, MVertex *n2)
}
+void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
+{
+ double u[7] = {0,1,0,.5,0,.5,.3333};
+ double v[7] = {0,0,1,.5,.5,0,.3333};
+ double j[2][2];
+ for (int i=0;i<7;i++)
+ {
+ t->jac ( u[i], v[i] , j );
+ double det = det2x2 ( j );
+ minJ = std::min ( det , minJ );
+ maxJ = std::max ( det , maxJ );
+ }
+}
+
void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
{
typedef std::map<std::pair<MVertex*, MVertex*> , std::vector<MElement*> > edge2tris;
@@ -638,6 +652,7 @@ void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
e2t[p].push_back(t);
}
}
+
for ( edge2tris::iterator it = e2t.begin() ; it != e2t.end() ; ++it)
{
std::pair<MVertex*, MVertex*> edge = it->first;
@@ -650,6 +665,7 @@ void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
MVertex *n4 = edge.second;
MTriangle *t1 = (MTriangle*)triangles[0];
MTriangle *t2 = (MTriangle*)triangles[1];
+
MVertex *n1 = t1->getOtherVertex (n2,n4);
MVertex *n3 = t2->getOtherVertex (n2,n4);
if (n1<n2)
@@ -768,12 +784,48 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
- if (CTX.mesh.smooth_internal_edges)
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- {
- for (int i=0;i<10;i++)smoothInternalEdges ( *it , edgeVertices);
- }
-
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ {
+ if (CTX.mesh.smooth_internal_edges)
+ for (int i=0;i<10;i++)
+ smoothInternalEdges ( *it , edgeVertices);
+
+ bool twod = true;
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ MTriangle *t = (*it)->triangles[i];
+ if(t->getVertex(0)->z() != 0.0 || t->getVertex(1)->z() != 0.0 ||t->getVertex(2)->z() != 0.0)twod = false;
+ }
+
+
+ if (twod && order > 3)
+ {
+ Msg(INFO,"Validity check of higher order meshes is still bugged for orders > 3");
+ }
+ else if (twod)
+ {
+ double minJ = 1.e22;
+ double maxJ = -1.e22;
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ double minJloc = 1.e22;
+ double maxJloc = -1.e22;
+
+ MTriangle *t = (*it)->triangles[i];
+ getMinMaxJac (t,minJloc,maxJloc);
+ minJ = std::min(minJ,minJloc);
+ maxJ = std::max(maxJ,maxJloc);
+ if (minJloc*maxJloc < 0)
+ Msg(GERROR,"Triangle %d %d %d has negative jacobians in GFace %d",t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),t->getVertex(2)->getNum(),(*it)->tag());
+ }
+
+ if (minJ*maxJ < 0)
+ Msg(GERROR,"There exists triangles with negative jacobians in GFace %d",(*it)->tag());
+ else
+ Msg(INFO,"Jacobian range of face %d is %g %g",(*it)->tag(),minJ,maxJ);
+
+ }
+ }
+
double t2 = Cpu();
Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
Msg(STATUS1, "Mesh");