diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 18b578588860cc0c56ae236fadee15236003cdf1..8921c9a31560da0f5c3b52b695ee60d505f33add 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -709,7 +709,7 @@ bool GFace::buildSTLTriangulation()
// by default we assume that straight lines are geodesics
SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
{
- if(CTX.mesh.second_order_experimental && geomType() != GEntity::Plane){
+ if(CTX.mesh.second_order_experimental && geomType() != GEntity::Plane ){
// FIXME: this is buggy -- remove the CTX option once we do it in
// a robust manner
GPoint gp1 = point(pt1.x(), pt1.y());
@@ -719,7 +719,10 @@ SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
gp1.y() + t * (gp2.y() - gp1.y()),
gp1.z() + t * (gp2.z() - gp1.z())),
(double*)guess);
- return SPoint2(gp.u(), gp.v());
+ if (gp.g())
+ return SPoint2(gp.u(), gp.v());
+ else
+ return pt1 + (pt2 - pt1) * t;
}
else{
return pt1 + (pt2 - pt1) * t;
diff --git a/Geo/GPoint.h b/Geo/GPoint.h
index 7eed7c4d9f4db57f0a85cf1b21b09d2044a6c2bc..4024727b680787c6aa40fc456ba5e9034a22ac11 100644
--- a/Geo/GPoint.h
+++ b/Geo/GPoint.h
@@ -25,6 +25,7 @@ class GPoint
inline double &z() { return Z; }
inline double u() const { return par[0]; }
inline double v() const { return par[1]; }
+ inline const GEntity* g() const { return e; }
GPoint (double _x=0, double _y=0, double _z=0, const GEntity *onwhat=0)
: X(_x), Y(_y), Z(_z), e(onwhat)
{
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e500de751cb2fbff1b0439d5dbb00f9622b2240d..ee316e53952414639e0688cc9681355da97ce376 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1057,6 +1057,145 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
}
}
+
+int MTetrahedronN::getNumEdgesRep(){ return 6 * numSubEdges; }
+
+void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
+ static int ed [6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
+ int iEdge = num / numSubEdges;
+ int iSubEdge = num % numSubEdges;
+
+
+ int iVertex1 = ed [iEdge][0];
+ int iVertex2 = ed [iEdge][1];
+ double t1 = (double) iSubEdge / (double) numSubEdges;
+ double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+
+ double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+ double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+
+ SPoint3 pnt1, pnt2;
+ pnt(u1,v1,w1,pnt1);
+ pnt(u2,v2,w2,pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
+}
+
+
+int MTetrahedronN::getNumFacesRep(){ return 4 * numSubEdges * numSubEdges ; }
+
+void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
+ static int fak [4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}};
+ int iFace = num / (numSubEdges*numSubEdges);
+ int iSubFace = num % (numSubEdges*numSubEdges);
+
+ int iVertex1 = fak [iFace][0];
+ int iVertex2 = fak [iFace][1];
+ int iVertex3 = fak [iFace][2];
+
+ /*
+ 0
+ 0 1
+ 0 1 2
+ 0 1 2 3
+ 0 1 2 3 4
+ 0 1 2 3 4 5
+ */
+
+ // on the first layer, we have (numSubEdges-1) * 2 + 1 triangles
+ // on the second layer, we have (numSubEdges-2) * 2 + 1 triangles
+ // on the ith layer, we have (numSubEdges-1-i) * 2 + 1 triangles
+ int ix = 0, iy = 0;
+ int nbt = 0;
+ for (int i = 0; i < numSubEdges; i++){
+ int nbl = (numSubEdges - i - 1) * 2 + 1;
+ nbt += nbl;
+ if (nbt > iSubFace){
+ iy = i;
+ ix = nbl - (nbt - iSubFace);
+ break;
+ }
+ }
+
+ const double d = 1. / numSubEdges;
+
+ SPoint3 pnt1, pnt2, pnt3;
+ double J1[2][3], J2[2][3], J3[2][3];
+ double u1,v1,u2,v2,u3,v3;
+ if (ix % 2 == 0){
+ u1 = ix / 2 * d; v1= iy*d;
+ u2 = (ix / 2 + 1) * d ; v2 = iy * d;
+ u3 = ix / 2 * d ; v3 = (iy+1) * d;
+ }
+ else{
+ u1 = (ix / 2 + 1) * d; v1= iy * d;
+ u2 = (ix / 2 + 1) * d; v2= (iy + 1) * d;
+ u3 = ix / 2 * d ; v3 = (iy + 1) * d;
+ }
+
+ double U1 = pp[iVertex1][0] * (1.-u1-v1) + pp[iVertex2][0] * u1 + pp[iVertex3][0] * v1;
+ double U2 = pp[iVertex1][0] * (1.-u2-v2) + pp[iVertex2][0] * u2 + pp[iVertex3][0] * v2;
+ double U3 = pp[iVertex1][0] * (1.-u3-v3) + pp[iVertex2][0] * u3 + pp[iVertex3][0] * v3;
+
+ double V1 = pp[iVertex1][1] * (1.-u1-v1) + pp[iVertex2][1] * u1 + pp[iVertex3][1] * v1;
+ double V2 = pp[iVertex1][1] * (1.-u2-v2) + pp[iVertex2][1] * u2 + pp[iVertex3][1] * v2;
+ double V3 = pp[iVertex1][1] * (1.-u3-v3) + pp[iVertex2][1] * u3 + pp[iVertex3][1] * v3;
+
+ double W1 = pp[iVertex1][2] * (1.-u1-v1) + pp[iVertex2][2] * u1 + pp[iVertex3][2] * v1;
+ double W2 = pp[iVertex1][2] * (1.-u2-v2) + pp[iVertex2][2] * u2 + pp[iVertex3][2] * v2;
+ double W3 = pp[iVertex1][2] * (1.-u3-v3) + pp[iVertex2][2] * u3 + pp[iVertex3][2] * v3;
+
+ pnt(U1,V1,W1,pnt1);
+ pnt(U2,V2,W2,pnt2);
+ pnt(U3,V3,W3,pnt3);
+
+ x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
+
+ // facetted first
+
+ SVector3 d1(x[1]-x[0],y[1]-y[0],z[1]-z[0]);
+ SVector3 d2(x[2]-x[0],y[2]-y[0],z[2]-z[0]);
+ n[0] = crossprod(d1, d2);
+ n[0].normalize();
+ n[1] = n[0];
+ n[2] = n[0];
+
+ return;
+
+ {
+ SVector3 d1(J1[0][0], J1[0][1], J1[0][2]);
+ SVector3 d2(J1[1][0], J1[1][1], J1[1][2]);
+ n[0] = crossprod(d1, d2);
+ n[0].normalize();
+ }
+ {
+ SVector3 d1(J2[0][0], J2[0][1], J2[0][2]);
+ SVector3 d2(J2[1][0], J2[1][1], J2[1][2]);
+ n[1] = crossprod(d1, d2);
+ n[1].normalize();
+ }
+ {
+ SVector3 d1(J3[0][0], J3[0][1], J3[0][2]);
+ SVector3 d2(J3[1][0], J3[1][1], J3[1][2]);
+ n[2] = crossprod(d1, d2);
+ n[2].normalize();
+ }
+}
+
+
+
+
MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 806cba53e7316e0a81f50fa03655a9c802db4a43..8d66a06ad44765558c354f14a641c79292c76913 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -1619,6 +1619,12 @@ class MTetrahedronN : public MTetrahedron {
}
}
}
+
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumEdgesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumFacesRep();
+
virtual void jac(double u, double v, double w, double j[3][3])
{
MTetrahedron::jac(_order, _vs , u, v, w, j);
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 95802691d3549d1da14cc9f0fd4378240cb06169..6935469d6ed5c162615c837fe185e7986bc40554 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -240,7 +240,7 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
v.Pos.Z = qp.z();
bool result = ProjectPointOnSurface(s, v, u);
if (!result)
- printf("Project Point on surface %d\n",result);
+ return GPoint(-1.e22, -1.e22, -1.e22, 0, u);
return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 041d34142922a3f96b3fb35b6dd4ea5b0528b75e..a49e912c7dffa3919d0a7923260fd3ea45f6fc74 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -54,8 +54,9 @@ bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
double UU, VV;
if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
param = SPoint2(UU, VV);
- else
- param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ else
+ return false;
+ // param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
}
return true;
}
@@ -175,6 +176,8 @@ bool computeEquidistantParameters(GFace *gf, double u0, double uN, double v0, do
v[N] = vN;
t[N] = 1.0;
+ return true;
+
// create the tangent matrix
const int M = N - 2;
Double_Matrix J(M, M);
@@ -466,7 +469,10 @@ void getFaceVertices(GFace *gf,
}
if (reparamOK){
GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
- v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
+ if (gp.g())
+ v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
+ else
+ v = new MVertex(X,Y,Z, gf);
}
else{
v = new MVertex(X,Y,Z, gf);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 68f649ed82124068b0d3fc5b78c3a038f9550a1b..78b33d9cabe625d85627ea9d789eeb717ebc04e0 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -179,8 +179,11 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
Count++;
}
}
- guess[0]/=Count;
- guess[1]/=Count;
+ if (Count != 0){
+ guess[0]/=Count;
+ guess[1]/=Count;
+ }
+
for(int j = 0; j < 3; j++){
if(v[j]->onWhat()->dim() == 3){
v[j]->onWhat()->mesh_vertices.erase
@@ -192,16 +195,20 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// PARAMETRIC COORDINATES SHOULD BE SET for the vertex !!!!!!!!!!!!!
// This is not 100 % safe yet, so we reserve that operation for high order
// meshes.
- Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher, guess %g %g",gf->tag(),guess[0],guess[1]);
+ // Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher, guess %g %g",gf->tag(),guess[0],guess[1]);
GPoint gp = gf->closestPoint (SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
- Msg::Debug("The point has been projected back to the surface (%g %g %g) -> (%g %g %g), (%g,%g)",
- v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z(),gp.u(),gp.v());
// To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
//v1b = new MVertex(gp.x(),gp.y(),gp.z(),gf);
- v1b = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- //v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
-
+ if (gp.g()){
+ v1b = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
+ Msg::Info("The point has been projected back to the surface (%g %g %g) -> (%g %g %g)",
+ v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z());
+ }
+ else{
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
+ Msg::Warning("The point was not projected back to the surface (%g %g %g)", v[j]->x(), v[j]->y(), v[j]->z());
+ }
}
else{
v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index f3b45dfc8e8845d24a271b747ee24e4985ebdc20..75ea2cbb2797e4ed71a3491e87518c0b3ddbb456 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -220,7 +220,38 @@ double qmDistorsionOfMapping (MTriangle *e)
return dmin;
}
+static double mesh_functional_distorsion(MTetrahedron *t, double u, double v, double w)
+{
+ // compute uncurved element jacobian d_u x and d_v x
+ double mat[3][3];
+ t->jac(1, 0, 0, 0, 0, mat);
+ const double det1 = det3x3(mat);
+ t->jac(u, v, w, mat);
+ const double detN = det3x3(mat);
+
+ // printf("%g %g %g = %g %g\n",u,v,w,det1,detN);
+
+ if (det1 == 0 || detN == 0) return 0;
+ double dist = std::min(detN/det1, det1/detN);
+ return dist;
+}
+
+
double qmDistorsionOfMapping (MTetrahedron *e)
{
- return 1.0;
+ if (e->getPolynomialOrder() == 1)return 1.0;
+ IntPt *pts;
+ int npts;
+ e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+ double dmin;
+ for (int i=0;i<npts;i++){
+ const double u = pts[i].pt[0];
+ const double v = pts[i].pt[1];
+ const double w = pts[i].pt[2];
+ const double di = mesh_functional_distorsion (e,u,v,w);
+ dmin = (i==0)? di : std::min(dmin,di);
+ }
+ // printf("DMIN = %g\n\n",dmin);
+
+ return dmin< 0 ? 0 :dmin;
}