From ac640a8567e348561617924fbbfd7e1a34acc135 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 10 Jun 2008 08:37:34 +0000
Subject: [PATCH] Fixed bug in 2D mesh, added high order in-outs for tets,
added high order mesh visu
---
Common/GmshDefines.h | 3 +
Geo/GFace.cpp | 4 +-
Geo/GFace.h | 2 +-
Geo/GModelIO_Mesh.cpp | 5 +-
Geo/Geo.cpp | 4 +-
Geo/MElement.cpp | 150 +++++++-
Geo/MElement.h | 69 +++-
Geo/OCCFace.cpp | 18 +-
Geo/OCCFace.h | 2 +-
Geo/gmshFace.cpp | 6 +-
Geo/gmshFace.h | 2 +-
Mesh/BDS.cpp | 203 +++++++---
Mesh/BDS.h | 34 +-
Mesh/HighOrder.cpp | 19 +-
Mesh/Makefile | 3 +-
Mesh/meshGFace.cpp | 59 +--
Mesh/meshGFaceBDS.h | 1 +
Mesh/meshGFaceQuadrilateralize.cpp | 586 +++++++++++++++++++++++++++++
Mesh/meshGFaceQuadrilateralize.h | 25 ++
Mesh/meshGRegion.cpp | 43 ++-
Numeric/FunctionSpace.cpp | 22 +-
21 files changed, 1140 insertions(+), 120 deletions(-)
create mode 100644 Mesh/meshGFaceQuadrilateralize.cpp
create mode 100644 Mesh/meshGFaceQuadrilateralize.h
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 3992b51936..c89d9e04ab 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -79,6 +79,9 @@
#define MSH_LIN_4 26
#define MSH_LIN_5 27
#define MSH_LIN_6 28
+#define MSH_TET_20 29
+#define MSH_TET_35 30
+#define MSH_TET_56 31
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 5026a837d6..0116c247ef 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.64 2008-06-03 21:39:01 geuzaine Exp $
+// $Id: GFace.cpp,v 1.65 2008-06-10 08:37:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -575,7 +575,7 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p) const
return SPoint2(U, V);
}
-GPoint GFace::closestPoint(const SPoint3 & queryPoint) const
+GPoint GFace::closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const
{
Msg::Error("Closet point not implemented for this type of surface");
return GPoint(0, 0, 0);
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 8d1810e89b..18c8825b87 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -124,7 +124,7 @@ class GFace : public GEntity
virtual int containsParam(const SPoint2 &pt) const;
// Return the point on the face closest to the given point.
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
// Return the normal to the face at the given parameter location.
virtual SVector3 normal(const SPoint2 ¶m) const = 0;
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 94d6715e4c..7be95c325d 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.55 2008-06-07 17:20:46 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.56 2008-06-10 08:37:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -139,6 +139,9 @@ static int getNumVerticesForElementTypeMSH(int type)
case MSH_QUA_9 : return 4 + 4 + 1;
case MSH_TET_4 : return 4;
case MSH_TET_10 : return 4 + 6;
+ case MSH_TET_20 : return 4 + 12 + 4;
+ case MSH_TET_35 : return 4 + 16 + 12 + 3;
+ case MSH_TET_56 : return 4 + 20 + 24 + 8;
case MSH_HEX_8 : return 8;
case MSH_HEX_20 : return 8 + 12;
case MSH_HEX_27 : return 8 + 12 + 6 + 1;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index a068b31f59..a800505f5e 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.112 2008-06-07 17:20:46 geuzaine Exp $
+// $Id: Geo.cpp,v 1.113 2008-06-10 08:37:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -2954,7 +2954,7 @@ bool ProjectPointOnCurve(Curve *c, Vertex *v, Vertex *RES, Vertex *DER)
bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
{
- double x[3] = { 0.5, 0.5, 0.5 };
+ double x[3] = { 0.5, u[0], u[1] };
Vertex vv;
int check;
SURFACE = s;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e4ef0e2de4..e185b08086 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.70 2008-06-05 11:52:49 samtech Exp $
+// $Id: MElement.cpp,v 1.71 2008-06-10 08:37:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -718,6 +718,78 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
#endif
}
+void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
+{
+#if !defined(HAVE_GMSH_EMBEDDED)
+ double sf[256];
+ switch(ord){
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, sf); break;
+ default : throw;
+ }
+
+ double x = 0 ; for(int i = 0; i < 4; i++) x += sf[i] * _v[i]->x();
+ double y = 0 ; for(int i = 0; i < 4; i++) y += sf[i] * _v[i]->y();
+ double z = 0 ; for(int i = 0; i < 4; i++) z += sf[i] * _v[i]->z();
+
+ const int N = (ord+1)*(ord+2)*(ord+3)/6;
+
+ for(int i = 4; i < N+4; i++) x += sf[i] * vs[i - 4]->x();
+ for(int i = 4; i < N+4; i++) y += sf[i] * vs[i - 4]->y();
+ for(int i = 4; i < N+4; i++) z += sf[i] * vs[i - 4]->z();
+
+ p = SPoint3(x,y,z);
+#endif
+}
+
+void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double j[3][3])
+{
+#if defined(HAVE_GMSH_EMBEDDED)
+ return;
+#else
+ double grads[256][2];
+ switch(ord){
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).df(uu, vv, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).df(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[2][0] = 0 ; for(int i = 0; i < 3; i++) j[2][0] += grads [i][2] * _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();
+ j[2][1] = 0 ; for(int i = 0; i < 3; i++) j[2][1] += grads [i][2] * _v[i]->y();
+ j[0][2] = 0 ; for(int i = 0; i < 3; i++) j[0][2] += grads [i][0] * _v[i]->z();
+ j[1][2] = 0 ; for(int i = 0; i < 3; i++) j[1][2] += grads [i][1] * _v[i]->z();
+ j[2][2] = 0 ; for(int i = 0; i < 3; i++) j[2][2] += grads [i][2] * _v[i]->z();
+
+ if (ord == 1) return;
+
+ const int N = (ord+1)*(ord+2)*(ord+3)/6;
+
+ for(int i = 3; i < N+4; i++) j[0][0] += grads[i][0] * vs[i - 4]->x();
+ for(int i = 3; i < N+4; i++) j[1][0] += grads[i][1] * vs[i - 4]->x();
+ for(int i = 3; i < N+4; i++) j[2][0] += grads[i][2] * vs[i - 4]->x();
+
+ for(int i = 3; i < N+4; i++) j[0][1] += grads[i][0] * vs[i - 4]->y();
+ for(int i = 3; i < N+4; i++) j[1][1] += grads[i][1] * vs[i - 4]->y();
+ for(int i = 3; i < N+4; i++) j[2][1] += grads[i][2] * vs[i - 4]->y();
+
+ for(int i = 3; i < N+4; i++) j[0][2] += grads[i][0] * vs[i - 4]->z();
+ for(int i = 3; i < N+4; i++) j[1][2] += grads[i][1] * vs[i - 4]->z();
+ for(int i = 3; i < N+4; i++) j[2][2] += grads[i][2] * vs[i - 4]->z();
+#endif
+}
+
+
+
int MTriangle6::getNumEdgesRep(){ return 3 * 6; }
void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -754,7 +826,78 @@ void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
}
}
-int MTriangleN::getNumEdgesRep(){ return 3 * 12; }
+const int numSubEdges = 6;
+
+int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; }
+
+void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){
+
+ // 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, iy;
+ int nbt = 0;
+ for (int i=0;i<numSubEdges;i++){
+ int nbl = (numSubEdges-i-1)*2 + 1;
+ nbt += nbl;
+ if (nbt > num){
+ iy = i;
+ ix = nbl-(nbt-num);
+ break;
+ }
+ }
+
+ const double d = 1./numSubEdges;
+
+ SPoint3 pnt1, pnt2, pnt3;
+ double J1[2][3],J2[2][3],J3[2][3];
+ if (ix %2 == 0){
+ pnt(ix/2*d, iy*d, pnt1);
+ pnt((ix/2+1)*d, iy*d, pnt2);
+ pnt(ix/2*d, (iy+1)*d, pnt3);
+ jac(ix/2*d, iy*d, J1);
+ jac((ix/2+1)*d, iy*d, J2);
+ jac(ix/2*d, (iy+1)*d, J3);
+ }
+ else{
+ pnt((ix/2+1)*d, iy*d, pnt1);
+ pnt((ix/2+1)*d, (iy+1)*d, pnt2);
+ pnt(ix/2*d, (iy+1)*d, pnt3);
+ jac((ix/2+1)*d, iy*d, J1);
+ jac((ix/2+1)*d, (iy+1)*d, J2);
+ jac(ix/2*d, (iy+1)*d, J3);
+ }
+ {
+ 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();
+ }
+
+ 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();
+
+
+
+
+}
+
+
+int MTriangleN::getNumEdgesRep(){ return 3 * numSubEdges; }
void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -822,6 +965,9 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_PYR_5: return new MPyramid(v, num, part);
case MSH_PYR_13: return new MPyramid13(v, num, part);
case MSH_PYR_14: return new MPyramid14(v, num, part);
+ case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
+ case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
+ case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
default: return 0;
}
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 950ca87286..e25d730ff5 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -595,12 +595,9 @@ class MTriangleN : public MTriangle {
}
virtual int getNumEdgeVertices(){ return _order - 1; }
virtual int getNumEdgesRep();
+ virtual int getNumFacesRep();
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
- virtual int getNumFacesRep(){ return 0; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- throw;
- }
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getTypeForMSH()
{
if(_order == 3 && _vs.size() == 6) return MSH_TRI_9;
@@ -936,6 +933,8 @@ class MTetrahedron : public MElement {
default : s = 0.; break;
}
}
+ virtual void jac(int order, MVertex *verts[], double u, double v, double jac[3][3]);
+ virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
@@ -1033,6 +1032,66 @@ class MTetrahedron10 : public MTetrahedron {
}
};
+
+class MTetrahedronN : public MTetrahedron {
+ protected:
+ std::vector<MVertex *> _vs;
+ const short _order;
+ public:
+ MTetrahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
+ std::vector<MVertex*> &v, int order, int num=0, int part=0)
+ : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order)
+ {
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ MTetrahedronN(std::vector<MVertex*> &v, int order, int num=0, int part=0)
+ : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order)
+ {
+ for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ ~MTetrahedronN(){}
+ virtual int getPolynomialOrder(){ return _order; }
+ virtual int getNumVertices(){ return 4 + _vs.size(); }
+ virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual int getNumFaceVertices()
+ {
+ switch(getTypeForMSH()){
+ case MSH_TET_20 : return 1;
+ case MSH_TET_35 : return 3;
+ case MSH_TET_56 : return 6;
+ default :
+ throw;
+ }
+ }
+ virtual int getNumEdgeVertices(){ return _order - 1; }
+ virtual int getTypeForMSH()
+ {
+ // (p+1)*(p+2)*(p+3)/6
+ if(_order == 3 && _vs.size() + 4 == 20) return MSH_TET_20;
+ if(_order == 4 && _vs.size() + 4 == 35) return MSH_TET_35;
+ if(_order == 5 && _vs.size() + 4 == 56) return MSH_TET_56;
+ throw;
+ }
+ virtual void revert()
+ {
+ MVertex *tmp;
+ tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
+ std::vector<MVertex*> inv;
+ inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
+ _vs = inv;
+ }
+ virtual void jac(double u, double v, double j[3][3])
+ {
+ MTetrahedron::jac(_order, &(*(_vs.begin())), u, v, j);
+ }
+ virtual void pnt(double u, double v, SPoint3 &p)
+ {
+ MTetrahedron::pnt(_order, &(*(_vs.begin())), u, v, p);
+ }
+};
+
+
class MHexahedron : public MElement {
protected:
MVertex *_v[8];
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 4e1f3226f7..e95367c9dc 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.41 2008-05-25 07:10:57 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.42 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -143,17 +143,17 @@ GPoint OCCFace::point(double par1, double par2) const
return GPoint(val.X(), val.Y(), val.Z(), this, pp);
}
-GPoint OCCFace::closestPoint(const SPoint3 & qp) const
+GPoint OCCFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
{
gp_Pnt pnt(qp.x(), qp.y(), qp.z());
GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
-
- if(!proj.NbPoints()){
- Msg::Error("OCC Project Point on Surface FAIL");
- return GPoint(0, 0);
- }
-
- double pp[2];
+
+ if(!proj.NbPoints()){
+ Msg::Error("OCC Project Point on Surface FAIL");
+ return GPoint(0, 0);
+ }
+
+ double pp[2] = {initialGuess[0],initialGuess[1]};
proj.LowerDistanceParameters(pp[0], pp[1]);
Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 6d634af9d3..b3d2f3f5a5 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -40,7 +40,7 @@ class OCCFace : public GFace {
virtual ~OCCFace(){}
Range<double> parBounds(int i) const;
virtual GPoint point(double par1, double par2) const;
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
virtual int containsPoint(const SPoint3 &pt) const;
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index c15d0979da..6c86c92fb9 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.58 2008-06-05 11:52:49 samtech Exp $
+// $Id: gmshFace.cpp,v 1.59 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -196,10 +196,10 @@ GPoint gmshFace::point(double par1, double par2) const
}
}
-GPoint gmshFace::closestPoint(const SPoint3 & qp) const
+GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
{
Vertex v;
- double u[2];
+ double u[2] = {initialGuess[0],initialGuess[1]};
v.Pos.X = qp.x();
v.Pos.Y = qp.y();
v.Pos.Z = qp.z();
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index e7736899b8..80f8ba5efe 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -33,7 +33,7 @@ class gmshFace : public GFace {
Range<double> parBounds(int i) const;
void setModelEdges(std::list<GEdge*>&);
virtual GPoint point(double par1, double par2) const;
- virtual GPoint closestPoint(const SPoint3 & queryPoint) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
virtual int containsPoint(const SPoint3 &pt) const;
virtual double getMetricEigenvalue(const SPoint2 &);
virtual SVector3 normal(const SPoint2 ¶m) const;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index e78ef4b8ba..c57d6d8e83 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.106 2008-05-04 08:31:15 geuzaine Exp $
+// $Id: BDS.cpp,v 1.107 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -29,6 +29,7 @@
#include "qualityMeasures.h"
bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t);
+bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t);
void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf)
{
@@ -39,28 +40,55 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
std::list<BDS_Face*>::iterator tite = t.end();
while(tit != tite) {
BDS_Point *pts[4];
- (*tit)->getNodes(pts);
- if(param)
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- pts[0]->u, pts[0]->v, 0.0,
- pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,
- pts[0]->iD, pts[1]->iD, pts[2]->iD);
- else{
- if(!gf)
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- pts[0]->X, pts[0]->Y, pts[0]->Z,
- pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
- pts[0]->iD, pts[1]->iD, pts[2]->iD);
- else
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- pts[0]->X, pts[0]->Y, pts[0]->Z,
- pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
- gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
- gf->curvature(SPoint2(pts[2]->u, pts[2]->v)));
+ if (!(*tit)->deleted){
+ (*tit)->getNodes(pts);
+ if(param)
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ pts[0]->u, pts[0]->v, 0.0,
+ pts[1]->u, pts[1]->v, 0.0,
+ pts[2]->u, pts[2]->v, 0.0,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD);
+ else{
+ if(!gf){
+ if (pts[3]){
+ fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z,
+ pts[1]->X, pts[1]->Y, pts[1]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[3]->X, pts[3]->Y, pts[3]->Z,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD, pts[3]->iD);
+ }
+ else{
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z,
+ pts[1]->X, pts[1]->Y, pts[1]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD);
+ }
+ }
+ else{
+ if (pts[3]){
+ fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z,
+ pts[1]->X, pts[1]->Y, pts[1]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[3]->X, pts[3]->Y, pts[3]->Z,
+ gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
+ gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvature(SPoint2(pts[2]->u, pts[2]->v)),
+ gf->curvature(SPoint2(pts[3]->u, pts[3]->v)));
+ }
+ else{
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z,
+ pts[1]->X, pts[1]->Y, pts[1]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ gf->curvature(SPoint2(pts[0]->u, pts[0]->v)),
+ gf->curvature(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvature(SPoint2(pts[2]->u, pts[2]->v)));
+ }
+ }
+ }
}
++tit;
}
@@ -192,7 +220,7 @@ BDS_Edge *BDS_Mesh::find_edge(int num1, int num2)
}
int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
- double x3, double y3, double x4, double y4)
+ double x3, double y3, double x4, double y4,double x[2])
{
// double p1[2] = {x1,y1};
@@ -214,7 +242,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
// return 0;
double mat[2][2];
- double rhs[2], x[2];
+ double rhs[2];
mat[0][0] = (x2 - x1);
mat[0][1] = -(x4 - x3);
mat[1][0] = (y2 - y1);
@@ -228,6 +256,32 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
return 0;
}
+BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
+
+ std::list<BDS_Face*> ts;
+ p1->getTriangles(ts);
+ std::list<BDS_Face*>::iterator it = ts.begin();
+ std::list<BDS_Face*>::iterator ite = ts.end();
+ while(it != ite) {
+ BDS_Face *t = *it;
+ if (!t->e4){
+ BDS_Edge *e= t->oppositeEdge (p1);
+ BDS_Face *f= e->otherFace (t);
+ if (!f->e4){
+ BDS_Point *p2b = f->oppositeVertex(e);
+ if (p2 == p2b){
+ if (swap_edge(e, BDS_SwapEdgeTestQuality(false,false))){
+ printf("coucou\n");
+ return find_edge (p1,p2->iD);
+ }
+ }
+ }
+ }
+ ++it;
+ }
+ return 0;
+}
+
BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered)
{
@@ -244,6 +298,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
int ix = 0;
int ixMax = 300;
+ double x[2];
while(1){
std::vector<BDS_Edge*> intersected;
std::list<BDS_Edge*>::iterator it = edges.begin();
@@ -256,7 +311,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
if(Intersect_Edges_2d(e->p1->u, e->p1->v,
e->p2->u, e->p2->v,
p1->u, p1->v,
- p2->u, p2->v)){
+ p2->u, p2->v,x)){
// intersect
if(e2r && e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0)) != e2r->end()){
std::set<EdgeToRecover>::iterator itr1 =
@@ -271,6 +326,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
not_recovered->insert(EdgeToRecover(e->p1->iD, e->p2->iD, itr1->ge));
selfIntersection = true;
}
+ if (e->numfaces() != e->numTriangles())return 0;
intersected.push_back(e);
}
++it;
@@ -1102,8 +1158,11 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
// use robust predicates for not allowing to revert a triangle by
// moving one of its vertices
+
bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t)
{
+ if (t->e4)
+ return test_move_point_parametric_quad(p,u,v,t);
BDS_Point *pts[4];
t->getNodes(pts);
@@ -1142,6 +1201,46 @@ bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_F
return ori_init*ori_final > 0;
}
+bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t)
+{
+ BDS_Point *pts[4];
+ t->getNodes(pts);
+
+ double pa[2] = {pts[0]->u, pts[0]->v};
+ double pb[2] = {pts[1]->u, pts[1]->v};
+ double pc[2] = {pts[2]->u, pts[2]->v};
+ double pd[2] = {pts[3]->u, pts[3]->v};
+
+ double ori_init1 = gmsh::orient2d(pa, pb, pc);
+ double ori_init2 = gmsh::orient2d(pc, pd, pa);
+
+ if(p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
+ }
+ else if(p == pts[1]){
+ pb[0] = u;
+ pb[1] = v;
+ }
+ else if(p == pts[2]){
+ pc[0] = u;
+ pc[1] = v;
+ }
+ else if(p == pts[3]){
+ pd[0] = u;
+ pd[1] = v;
+ }
+ else
+ throw;
+
+ double ori_final1 = gmsh::orient2d(pa, pb, pc);
+ double ori_final2 = gmsh::orient2d(pc, pd, pa);
+ // allow to move a point when a triangle was flat
+ return ori_init1*ori_final1 > 0 && ori_init2*ori_final2 > 0 ;
+}
+
+
+
// d^2_i = (x^2_i - x)^T M (x_i - x)
// = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2
@@ -1265,18 +1364,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
while(ited != itede) {
BDS_Edge *e = *ited;
BDS_Point *n = e->othervertex(p);
-// double uv[2] = {(n->u + p->u)/2.0,(n->v + p->v)/2.0};
-// double metric[3];
-// buildMetric ( gf ,uv,metric);
-// double du[2] = {n->u - p->u,n->v - p->v};
-// double ldu = sqrt(DSQR(du[0])+DSQR(du[1]));
-// du[0]/=ldu;
-// du[1]/=ldu;
-// double fact = 1./sqrt (metric[0] * du[0] * du[0] +
-// 2 * metric[1] * du[0] * du[1] +
-// metric[2] * du[1] * du[1]);
- double fact = 1.0;
-
+ double fact = 1.0;
sTot += fact;
U += n->u * fact;
V += n->v * fact;
@@ -1366,30 +1454,35 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
double tot_length = 0;
double LC = 0;
- std::list<BDS_Edge*>::iterator eit = p->edges.begin();
- while(eit != p->edges.end()) {
- if((*eit)->numfaces() == 1) return false;
- BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
- const double l_e = (*eit)->length();
- U += op->u * l_e;
- V += op->v * l_e;
- tot_length += l_e;
- LC += op->lc();
- ++eit;
- }
-
- U /= tot_length;
- V /= tot_length;
- LC /= p->edges.size();
std::list<BDS_Face*> ts;
p->getTriangles(ts);
std::list<BDS_Face*>::iterator it = ts.begin();
std::list<BDS_Face*>::iterator ite = ts.end();
+
while(it != ite) {
BDS_Face *t = *it;
- if(!test_move_point_parametric_triangle(p, U, V, t))
+ BDS_Point *n[4];
+ t->getNodes(n);
+ for (int i = 0; i<t->numEdges();i++){
+ U += n[i]->u;
+ V += n[i]->v;
+ LC += n[i]->lc();
+ tot_length += 1;
+ }
+ ++it;
+ }
+ U /= tot_length;
+ V /= tot_length;
+ LC /= p->edges.size();
+
+ it = ts.begin();
+ while(it != ite) {
+ BDS_Face *t = *it;
+ if(!test_move_point_parametric_triangle(p, U, V, t)){
+ printf("coucou %g %g -> %g %g\n", p->u, p->v,U,V);
return false;
+ }
++it;
}
@@ -1400,7 +1493,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
- eit = p->edges.begin();
+ std::list<BDS_Edge*>::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
(*eit)->update();
++eit;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 1622b46c3d..4ffc15521c 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -254,9 +254,13 @@ public:
}
inline BDS_Face * otherFace(const BDS_Face *f) const
{
- if(numfaces()!=2) throw;
+ if(numfaces()!=2) {
+ printf("otherFace wrong, ony %d faces attached to edge %d %d\n",numfaces(),p1->iD,p2->iD);
+ throw;
+ }
if(f == _faces[0]) return _faces[1];
if(f == _faces[1]) return _faces[0];
+ printf("otherFace wrong : the edge does not belong to the face \n",numfaces(),p1->iD,p2->iD);
throw;
}
inline void del(BDS_Face *t)
@@ -299,6 +303,29 @@ public:
BDS_Edge *e1, *e2, *e3, *e4;
BDS_GeomEntity *g;
inline int numEdges () const { return e4 ? 4 : 3; }
+ inline BDS_Edge *oppositeEdge (BDS_Point *p){
+ if (e4){
+ printf("oppositeEdge to point %d cannot be applied to a quad\n",p->iD);
+ throw;
+ }
+ if (e1->p1 != p && e1->p2 != p)return e1;
+ if (e2->p1 != p && e2->p2 != p)return e2;
+ if (e3->p1 != p && e3->p2 != p)return e3;
+ printf("point %d does not belong to this triangle\n",p->iD);
+ throw;
+ }
+ inline BDS_Point *oppositeVertex (BDS_Edge *e){
+ if (e4){
+ printf("oppositeVertex to edge %d %d cannot be applied to a quad\n",e->p1->iD,e->p2->iD);
+ throw;
+ }
+
+ if (e == e1)return e2->commonvertex(e3);
+ if (e == e2)return e1->commonvertex(e3);
+ if (e == e3)return e1->commonvertex(e2);
+ printf("edge %d %d does not belong to this triangle\n",e->p1->iD,e->p2->iD);
+ throw;
+ }
inline void getNodes(BDS_Point *n[4]) const
{
if (!e4){
@@ -459,6 +486,7 @@ public:
// 2D operators
BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0,
std::set<EdgeToRecover> *not_recovered=0);
+ BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2);
bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest);
bool collapse_edge_parametric(BDS_Edge*, BDS_Point*);
void snap_point(BDS_Point*, BDS_Mesh *geom = 0);
@@ -478,4 +506,8 @@ public:
void outputScalarField(std::list<BDS_Face*> t, const char *fn, int param, GFace *gf=0);
void recur_tag(BDS_Face *t, BDS_GeomEntity *g);
+int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
+ double x3, double y3, double x4, double y4,
+ double x[2]);
+
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 83b8771738..70eb301412 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.29 2008-05-06 21:11:47 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.30 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -622,11 +622,20 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std::vector<MTetrahedron*> tetrahedra2;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
- std::vector<MVertex*> ve;
+ std::vector<MVertex*> ve,vf;
getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
- tetrahedra2.push_back
- (new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5]));
+ if (nPts == 1){
+ tetrahedra2.push_back
+ (new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5]));
+ }
+ else{
+ getFaceVertices(gr, t, vf, faceVertices, linear, nPts);
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ tetrahedra2.push_back
+ (new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
+ ve, nPts + 1));
+ }
delete t;
}
gr->tetrahedra = tetrahedra2;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 97582f8029..2277626b21 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.224 2008-06-07 17:20:48 geuzaine Exp $
+# $Id: Makefile,v 1.225 2008-06-10 08:37:34 remacle Exp $
#
# Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
#
@@ -44,6 +44,7 @@ SRC = Generator.cpp \
meshGFaceBDS.cpp \
meshGFaceDelaunayInsertion.cpp \
meshGFaceOptimize.cpp \
+ meshGFaceQuadrilateralize.cpp \
meshGRegion.cpp \
meshGRegionDelaunayInsertion.cpp \
meshGRegionTransfinite.cpp \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 116a55ba72..a4cb15bd91 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.135 2008-06-05 11:52:50 samtech Exp $
+// $Id: meshGFace.cpp,v 1.136 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -23,6 +23,7 @@
#include "meshGFace.h"
#include "meshGFaceBDS.h"
#include "meshGFaceDelaunayInsertion.h"
+#include "meshGFaceQuadrilateralize.h"
#include "meshGFaceOptimize.h"
#include "DivideAndConquer.h"
#include "BackgroundMesh.h"
@@ -102,8 +103,8 @@ void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
for (int i = 0; i < N; i++){
MVertex *v1 = itr->ge->lines[i]->getVertex(0);
MVertex *v2 = itr->ge->lines[i]->getVertex(1);
- if ((v1->getNum() == p1 && v2->getNum() == p2) ||
- (v1->getNum() == p2 && v2->getNum() == p1)){
+ if ((v1->getIndex() == p1 && v2->getIndex() == p2) ||
+ (v1->getIndex() == p2 && v2->getIndex() == p1)){
double t1;
double lc1 = -1;
if (v1->onWhat() == g1) t1 = bb.low();
@@ -242,23 +243,23 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
if(pass_ == 1){
- e2r->insert(EdgeToRecover(vstart->getNum(), vend->getNum(), ge));
+ e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
return true;
}
else{
- BDS_Point *pstart = m->find_point(vstart->getNum());
- BDS_Point *pend = m->find_point(vend->getNum());
+ BDS_Point *pstart = m->find_point(vstart->getIndex());
+ BDS_Point *pend = m->find_point(vend->getIndex());
if(!pstart->g){
- m->add_geom (vstart->getNum(), 0);
- BDS_GeomEntity *g0 = m->get_geom(vstart->getNum(), 0);
+ m->add_geom (vstart->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
pstart->g = g0;
}
if(!pend->g){
- m->add_geom(vend->getNum(), 0);
- BDS_GeomEntity *g0 = m->get_geom(vend->getNum(), 0);
+ m->add_geom(vend->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
pend->g = g0;
}
- BDS_Edge * e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
+ BDS_Edge * e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("The unrecoverable edge is on model edge %d", ge->tag());
@@ -272,15 +273,15 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->mesh_vertices.begin());
- if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getNum(), vend->getNum(), ge));
+ if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
- BDS_Point *pstart = m->find_point(vstart->getNum());
+ BDS_Point *pstart = m->find_point(vstart->getIndex());
if(!pstart->g){
- m->add_geom (vstart->getNum(), 0);
- BDS_GeomEntity *g0 = m->get_geom(vstart->getNum(), 0);
+ m->add_geom (vstart->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
pstart->g = g0;
}
- e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
+ e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("The unrecoverable edge is on model edge %d", ge->tag());
@@ -291,9 +292,9 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
for (unsigned int i = 1; i < ge->mesh_vertices.size(); i++){
vstart = ge->mesh_vertices[i - 1];
vend = ge->mesh_vertices[i];
- if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getNum(), vend->getNum(), ge));
+ if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
- e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
+ e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
@@ -305,9 +306,9 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
}
vstart = vend;
vend = *(ge->getEndVertex()->mesh_vertices.begin());
- if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getNum(), vend->getNum(), ge));
+ if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
- e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
+ e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e)e->g = g;
else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
@@ -315,10 +316,10 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
// ge->mesh_vertices.size(), ge->mesh_vertices.size());
// return false;
}
- BDS_Point *pend = m->find_point(vend->getNum());
+ BDS_Point *pend = m->find_point(vend->getIndex());
if(!pend->g){
- m->add_geom (vend->getNum(), 0);
- BDS_GeomEntity *g0 = m->get_geom(vend->getNum(), 0);
+ m->add_geom (vend->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
pend->g = g0;
}
}
@@ -402,8 +403,8 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
}
U_[count] = param.x();
V_[count] = param.y();
- (*itv)->setNum(count);
- numbered_vertices[(*itv)->getNum()] = *itv;
+ (*itv)->setIndex(count);
+ numbered_vertices[(*itv)->getIndex()] = *itv;
bbox += SPoint3(param.x(), param.y(), 0);
count++;
++itv;
@@ -468,7 +469,7 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
for(int i = 0; i < doc.numPoints; i++){
MVertex *here = (MVertex *)doc.points[i].data;
- int num = here->getNum();
+ int num = here->getIndex();
double U, V;
if(num < 0){ // fake bbox points
U = bb[-1 - num]->x();
@@ -486,7 +487,7 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
MVertex *V1 = (MVertex*)doc.points[doc.triangles[i].a].data;
MVertex *V2 = (MVertex*)doc.points[doc.triangles[i].b].data;
MVertex *V3 = (MVertex*)doc.points[doc.triangles[i].c].data;
- m->add_triangle(V1->getNum(), V2->getNum(), V3->getNum());
+ m->add_triangle(V1->getIndex(), V2->getIndex(), V3->getIndex());
}
@@ -597,7 +598,7 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
Msg::Debug("Computing mesh size field at mesh vertices", edgesToRecover.size());
for(int i = 0; i < doc.numPoints; i++){
MVertex *here = (MVertex *)doc.points[i].data;
- int num = here->getNum();
+ int num = here->getIndex();
GEntity *ge = here->onWhat();
BDS_Point *pp = m->find_point(num);
if(ge->dim() == 0){
@@ -1337,6 +1338,8 @@ void meshGFace::operator() (GFace *gf)
Msg::Error("Impossible to mesh face %d", gf->tag());
}
+ // gmshQMorph(gf);
+
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
}
diff --git a/Mesh/meshGFaceBDS.h b/Mesh/meshGFaceBDS.h
index dd2af73680..994658785d 100644
--- a/Mesh/meshGFaceBDS.h
+++ b/Mesh/meshGFaceBDS.h
@@ -41,4 +41,5 @@ void gmshDelaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap);
void gmshCollapseSmallEdges(GModel &gm);
BDS_Mesh *gmsh2BDS(std::list<GFace*> &l);
double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2);
+void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q);
#endif
diff --git a/Mesh/meshGFaceQuadrilateralize.cpp b/Mesh/meshGFaceQuadrilateralize.cpp
new file mode 100644
index 0000000000..ac1df83780
--- /dev/null
+++ b/Mesh/meshGFaceQuadrilateralize.cpp
@@ -0,0 +1,586 @@
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "meshGFaceQuadrilateralize.h"
+#include "Message.h"
+#include "Numeric.h"
+#include "meshGFaceDelaunayInsertion.h"
+#include "meshGFaceOptimize.h"
+#include "meshGFaceBDS.h"
+#include "BDS.h"
+#include "SVector3.h"
+
+class gmshEdgeFront {
+public:
+ typedef std::set<BDS_Edge *>::const_iterator eiter;
+private:
+ BDS_Mesh *m;
+ GFace *gf;
+ void getFrontEdges (BDS_Point *p, eiter & it1, eiter & it2) const;
+ void getFrontEdges (BDS_Point *p, std::vector<eiter> & f) const;
+public:
+ std::set<BDS_Edge *> edges;
+ std::set<BDS_Edge*> stat[5];
+ eiter begin() {return edges.begin();}
+ eiter end() {return edges.end();}
+ gmshEdgeFront (BDS_Mesh *_m, GFace *_gf) : m(_m), gf(_gf){
+ }
+ // initiate edges in the front i.e.
+ // take all edges that have one neighbor
+ // and all edges that have a quad and a triangle
+ // neighbor
+ void initiate ();
+ // compute normal vector to an edge that points
+ // inside the domain
+ SVector3 normal (BDS_Edge *e) const;
+ // compute the state of a front edge
+ // 0 both vertices have edge angles < 60 deg
+ // 1 e->p1 have and edge angle > 60 deg
+ // 2 e->p2 have and edge angle > 60 deg
+ // 3 both vertices have edge angles > 60 deg
+ int computeStatus (BDS_Edge *e) const;
+ inline bool inFront (BDS_Edge *e) const{return getStatus(e) != -1;}
+ inline int getStatus (BDS_Edge *e) const{
+ for (int i=0;i<5;i++){
+ if (stat[i].find(e) != stat[i].end())return i;
+ }
+ return -1;
+ }
+ // get one edge of the front that is tagged "tag"
+ // and do whatever has to be done to remove it from
+ // the front/tag. return false if the front is not empty.
+ bool emptyFront (int tag);
+ // update the status of the edge
+ void updateStatus (BDS_Edge *e);
+ // form a quad now
+ bool formQuad (BDS_Edge *e, BDS_Edge *left, BDS_Edge *right);
+ // delete the cavity delimitated by 4 edges
+ void emptyCavity (BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left, BDS_Edge *right);
+ // delete an edge from the front
+ void deleteFromFront(BDS_Edge *e){
+ edges.erase(e);
+ for (int i=0;i<5;i++){
+ std::set<BDS_Edge*>::iterator it = stat[i].find(e);
+ if (it !=stat[i].end()){
+ stat[i].erase(it);
+ return;
+ }
+ }
+ }
+ // taking a point from the front, find the optimal edge
+ BDS_Edge *findOptimalEdge(BDS_Point *p,BDS_Point *avoid);
+};
+
+void gmshEdgeFront::updateStatus (BDS_Edge *e){
+ for (int i=0;i<5;i++){
+ std::set<BDS_Edge*>::iterator it = stat[i].find(e);
+ if (it !=stat[i].end()){
+ stat[i].erase(it);
+ stat[computeStatus(e)].insert(e);
+ return;
+ }
+ }
+ throw;
+}
+
+SVector3 norm_edge (BDS_Point *p1, BDS_Point *p2){
+ SVector3 n (p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
+ n.normalize();
+ return n;
+}
+
+void recur_empty_cavity (BDS_Face *f,
+ BDS_Edge *be[4],
+ BDS_Point *bv[4],
+ std::set<BDS_Face*> & faces,
+ std::set<BDS_Edge*> & edges,
+ std::set<BDS_Point*> & vertices){
+ if (faces.find(f) != faces.end())return;
+ faces.insert(f);
+ BDS_Edge *ee[3] = {f->e1,f->e2,f->e3};
+ for (int i=0;i<3;i++){
+ BDS_Edge *e = ee[i];
+ if (e != be[0] && e != be[1] && e != be[2] && e != be[3])
+ {
+ edges.insert(e);
+ BDS_Face *of = e->otherFace(f);
+ recur_empty_cavity (of,be,bv,faces,edges,vertices);
+ }
+ }
+}
+
+
+void gmshEdgeFront::emptyCavity (BDS_Edge *bottom, BDS_Edge *top, BDS_Edge *left, BDS_Edge *right){
+ // not optimal for now, will be improved further away
+ BDS_Face *f=0;
+ if (bottom->faces(0) && bottom->faces(0)->numEdges() == 3)f=bottom->faces(0);
+ else if (bottom->faces(1))f = bottom->faces(1);
+
+ std::set<BDS_Face*> m_faces;
+ std::set<BDS_Edge*> m_edges;
+ std::set<BDS_Point*> m_vertices;
+ BDS_Edge *be[4] = {bottom,top,left,right};
+ BDS_Point *bv[4] = {bottom->commonvertex(left),
+ left->commonvertex(top),
+ top->commonvertex(right),
+ right->commonvertex(bottom)};
+
+
+ recur_empty_cavity (f,be,bv,m_faces,m_edges,m_vertices);
+ // printf("%d edges %d faces\n",m_edges.size(),m_faces.size());
+ for (std::set<BDS_Face*> :: iterator it = m_faces.begin() ; it != m_faces.end() ; ++it)m->del_face(*it);
+ for (std::set<BDS_Edge*> :: iterator it = m_edges.begin() ; it != m_edges.end() ; ++it)m->del_edge(*it);
+}
+
+
+SVector3 gmshEdgeFront::normal (BDS_Edge*e) const{
+ BDS_Face *t = e->faces(0);
+ if (e->numfaces() == 2 && e->faces(1)->numEdges() == 3)t=e->faces(1);
+
+ // printf("%d faces e->faces(0)->numEdges() = %d\n",e->numfaces(),e->faces(0)->numEdges());
+
+ /*
+ points p1, p2 and p3
+ p3 does NOT belong to the edge
+ p = p1 (1-u-v) + p2 u + p3 v;
+ J = [p2-p1 p3-p1 (p2-p1)x(p3-p1)]
+ N = - J^-1 * {0,1,0}
+ n = N/|N|
+ */
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+ BDS_Point *p3;
+ if (e == t->e1)
+ p3 = t->e2->commonvertex(t->e3);
+ else if (e == t->e2)
+ p3 = t->e1->commonvertex(t->e3);
+ else if (e == t->e3)
+ p3 = t->e2->commonvertex(t->e1);
+ else
+ throw;
+
+ SVector3 t1 (p2->X-p1->X,p2->Y-p1->Y,p2->Z-p1->Z);
+ SVector3 t2 (p3->X-p1->X,p3->Y-p1->Y,p3->Z-p1->Z);
+ SVector3 t3 = crossprod(t1,t2);
+ double m[3][3] = {{t1.x(),t2.x(),t3.x()},
+ {t1.y(),t2.y(),t3.y()},
+ {t1.z(),t2.z(),t3.z()}};
+ double im[3][3];
+ inv3x3(m,im);
+ SVector3 n (im[1][0],im[1][1],im[1][2]);
+ n.normalize();
+ return n;
+}
+
+void gmshEdgeFront::getFrontEdges (BDS_Point *p, std::vector<eiter> & fe) const{
+ for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+ eiter it = edges.find(*itp);
+ if ( it != edges.end())
+ fe.push_back(it);
+ }
+}
+
+
+void gmshEdgeFront::getFrontEdges (BDS_Point *p, eiter & it1, eiter & it2) const{
+ int count = 0;
+ for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+ if (count == 0){
+ it1 = edges.find(*itp);
+ if ( it1 != edges.end()) count++;
+ }
+ else if (count == 1){
+ it2 = edges.find(*itp);
+ if ( it2 != edges.end()) return;
+ }
+ }
+ printf("point %d is in the front but has only %d edges\n",p->iD,count);
+
+ throw;
+}
+
+void gmshEdgeFront::initiate () {
+ edges.clear();
+ for (int i=0;i<5;i++)stat[i].clear();
+ std::list<BDS_Edge*>::iterator it = m->edges.begin();
+ while (it!= m->edges.end()){
+ if (((*it)->numfaces() == 1 && (*it)->faces(0)->e4 == 0) ||
+ ((*it)->numfaces() == 2 && (*it)->numTriangles() == 1)) {
+ edges.insert(*it);
+ }
+ ++it;
+ }
+ for (eiter it2 = begin(); it2 !=end() ; ++it2){
+ int status = computeStatus(*it2);
+ stat[status].insert(*it2);
+ }
+}
+
+
+double angle3Points ( BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 ){
+ SVector3 a(p1->X-p2->X,p1->Y-p2->Y,p1->Z-p2->Z);
+ SVector3 b(p3->X-p2->X,p3->Y-p2->Y,p3->Z-p2->Z);
+ SVector3 c = crossprod(a,b);
+ double sinA = c.norm();
+ double cosA = dot(a,b);
+ // printf("%d %d %d -> %g %g\n",p1->iD,p2->iD,p3->iD,cosA,sinA);
+ return atan2 (sinA,cosA);
+}
+
+int gmshEdgeFront::computeStatus (BDS_Edge *e) const {
+ eiter it11, it12, it21,it22;
+ getFrontEdges (e->p1, it11, it12);
+ getFrontEdges (e->p2, it21, it22);
+
+ BDS_Edge *e1 = (*it11) == e ? *it12 : *it11;
+ BDS_Edge *e2 = (*it21) == e ? *it22 : *it21;
+
+ double angle1 = angle3Points ( (*it11)->othervertex(e->p1), e->p1, (*it12)->othervertex(e->p1));
+ double angle2 = angle3Points ( (*it21)->othervertex(e->p2), e->p2, (*it22)->othervertex(e->p2));
+
+ SVector3 n1 = normal(e);
+ SVector3 n2 = norm_edge (e->p1,e1->othervertex(e->p1));
+ SVector3 n3 = norm_edge (e->p2,e2->othervertex(e->p2));
+ if (dot(n1,n2) < 0)angle1 = M_PI;
+ if (dot(n1,n3) < 0)angle2 = M_PI;
+
+ // printf("edge %d %d angles %g (%d %d) %g (%d %d)\n",e->p1->iD,e->p2->iD,angle1*180/M_PI,e1->p1->iD,e1->p2->iD,angle2*180/M_PI,
+ // e2->p1->iD,e2->p2->iD);
+
+ const double angleLimit = 3*M_PI/4.;
+
+ if (angle1 > angleLimit && angle2 > angleLimit)return 0;
+ if (angle1 > angleLimit && angle2 < angleLimit)return 1;
+ if (angle1 < angleLimit && angle2 > angleLimit)return 2;
+ return 3;
+}
+
+
+bool gmshEdgeFront::formQuad (BDS_Edge *e,
+ BDS_Edge *left,
+ BDS_Edge *right){
+
+ printf("e (%d,%d), l(%d,%d), r(%d,%d)\n",
+ e->p1->iD,e->p2->iD,
+ left->p1->iD,left->p2->iD,
+ right->p1->iD,right->p2->iD);
+
+
+ // outputScalarField(m->triangles, "deb_before.pos", 0);
+
+
+ std::vector<BDS_Point*> toUpdate;
+
+ BDS_Point *pleft = left->othervertex(e->p1);
+ BDS_Point *pright = right->othervertex(e->p2);
+
+ // recover edge pleft->pright
+ BDS_Edge *top = m->find_edge(pleft,pright);
+
+ /*
+ We first have to ensure that, if left or right
+ are closing the front, then even number of edges should remain in
+ both left and right parts of the front
+ */
+
+ if (!top) {
+ // printf("recover\n");
+ // top = m->recover_edge_fast (pleft,pright);
+ // if (!top)
+ top = m->recover_edge (pleft->iD,pright->iD);
+ // printf("recover done %p\n",top);
+ if (!top)return false;
+ }
+ // delete internal triangles
+ emptyCavity (e,top,left,right);
+ // add the quad
+ m->add_quadrangle(e,left,top,right);
+ // update the edge status
+ // bottom edge leaves the front
+ deleteFromFront(e);
+ // top edge becomes part of the front
+
+ /* printf("(%d,%d),(%d,%d),(%d,%d),(%d,%d)\n",
+ e->p1->iD,e->p2->iD,
+ left->p1->iD,left->p2->iD,
+ top->p1->iD,top->p2->iD,
+ right->p1->iD,right->p2->iD);
+ */
+ // outputScalarField(m->triangles, "deb.pos", 0);
+
+
+ // if left edge was in the front, then it leaves the front
+ // because it has either 2 neighboring quads or one quad
+ // and the void
+ // if the left edge was not in the front, then it becomes
+ // part of it.
+ // printf("coucou1\n");
+ if (inFront(left))deleteFromFront(left);
+ else edges.insert(left);
+
+ // printf("coucou2\n");
+ if (inFront(right))deleteFromFront(right);
+ else edges.insert(right);
+
+ // printf("coucou3\n");
+ if (inFront(top))deleteFromFront(top);
+ else edges.insert(top);
+
+ toUpdate.push_back(e->p1);
+ toUpdate.push_back(e->p2);
+ toUpdate.push_back(pleft);
+ toUpdate.push_back(pright);
+
+ for (int i=0;i<toUpdate.size();i++){
+ toUpdate[i]->config_modified = true;
+ bool done = m->smooth_point_parametric(toUpdate[i], gf);
+ // printf("smooth done %d (g %d)\n",done,toUpdate[i]->g->classif_degree);
+ }
+
+ for (int i=0;i<toUpdate.size();i++){
+ BDS_Point *p = toUpdate[i];
+ for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+ if (inFront(*itp)){
+ updateStatus(*itp);
+ }
+ }
+ }
+ return true;
+}
+
+BDS_Edge *gmshEdgeFront::findOptimalEdge(BDS_Point *p, BDS_Point *avoid){
+
+ eiter it1,it2;
+ getFrontEdges (p, it1, it2);
+ // compute bissector of front edges, this is the optimal direction
+ SVector3 n1 = normal(*it1);
+ SVector3 n2 = normal(*it2);
+ SVector3 n = n1+n2;
+ n.normalize();
+
+ // printf("POINT %g %g %g N %g %g %g\n",p->X,p->Y,p->Z,n.x(),n.y(),n.z());
+
+ double lowerBound = cos (M_PI/6.0);
+ BDS_Edge *found = 0;
+ for (std::list<BDS_Edge*>::iterator itp = p->edges.begin(); itp != p->edges.end() ; ++ itp){
+ // the edge is not in the front and is not bounded by quads only
+ if (*it1 != *itp && *it2 != *itp && (*itp)->numTriangles()){
+ BDS_Point *q = (*itp)->othervertex(p);
+ SVector3 d (q->X-p->X,q->Y-p->Y,q->Z-p->Z);
+ d.normalize();
+ double COS = dot(n,d);
+ if (COS > lowerBound && q != avoid){
+ lowerBound = COS;
+ found = *itp;
+ }
+ }
+ }
+ if (found) return found;
+ else{
+
+ std::list<BDS_Face*> ts;
+ double x[2];
+ const double L = 0.5*sqrt(3.0)*((*it2)->length() * (*it1)->length()) ;
+ p->getTriangles(ts);
+ std::list<BDS_Face*>::iterator it = ts.begin();
+ std::list<BDS_Face*>::iterator ite = ts.end();
+ while(it != ite) {
+ BDS_Face *t = *it;
+ if (!t->e4){
+ BDS_Edge *e = t->oppositeEdge(p);
+ if (e->numfaces() == 2){
+ BDS_Face *f = e->otherFace(t);
+ if (!f->e4){
+ BDS_Point *target = f->oppositeVertex(e);
+ // ONLY WORKS IN 2D for now !!!!!!!!!!!!!!!!!!!
+ Intersect_Edges_2d ( e->p1->X,e->p1->Y,
+ e->p2->X,e->p2->Y,
+ p->X,p->Y,
+ p->X+n.x(),p->Y + n.y(),x);
+ if ( x[0] >= 0 && x[0] <= 1){
+ SVector3 d (target->X-p->X,target->Y-p->Y,target->Z-p->Z);
+ d.normalize();
+ double COS = dot(n,d);
+ double L2 = sqrt ((target->X - p->X) *(target->X - p->X) +
+ (target->X - p->Y) *(target->X - p->Y) +
+ (target->X - p->Z) *(target->X - p->Z) );
+
+ // swapping the edge alllow to find an edgge that has the right direction and
+ // right size
+ if (COS > cos (M_PI/6.0) && L2 < L){
+ m->swap_edge( e, BDS_SwapEdgeTestQuality(false,false));
+ BDS_Edge *newE = m->find_edge(p,target);
+ // printf("swapping -> %p\n",newE);
+ return newE;
+ }
+ // split the edge
+ else{
+ BDS_Point *mid;
+ mid = m->add_point(++m->MAXPOINTNUMBER,(1.-x[0])*e->p1->u + x[0]*e->p2->u,
+ (1.-x[0])*e->p1->v + x[0]*e->p2->v,gf);
+ mid->lc() = 0.5 * (p->lc() + target->lc());
+ mid->g = e->p1->g;
+ m->split_edge(e, mid);
+ BDS_Edge *newE = m->find_edge(p,mid);
+ // printf("splitting -> %p %p\n",newE,e->p1->g);
+ // m->cleanup();
+ return newE;
+ }
+ }
+ }
+ }
+ }
+ ++it;
+ }
+ }
+ printf("zarbi\n");
+ return 0;
+}
+
+
+
+bool gmshEdgeFront::emptyFront (int tag){
+ // front edges tagged "tag" is empty
+ if (stat[tag].size() == 0)return true;
+ BDS_Edge *e = *(stat[tag].begin());
+ BDS_Edge *left,*right;
+ eiter it1,it2;
+
+ std::vector<eiter> fe1,fe2;
+
+ printf("front edges %d %d tag %d\n",e->p1->iD,e->p2->iD,tag);
+
+ switch(tag){
+ // both left and right neighboring edges are
+ // sufficiently angled in order to allow to
+ // form the quad
+ case 3:
+ {
+ getFrontEdges (e->p1, it1, it2);
+ if (*it1 == e)left=*it2;
+ else left = *it1;
+ getFrontEdges (e->p2, it1, it2);
+ if (*it1 == e)right=*it2;
+ else right = *it1;
+ // printf("case 3\n");
+ }
+ break;
+ // right edge is angled with current edge
+ // we therefore find a new front edge in the
+ // mesh that allows to move to tag 3
+ case 2:
+ getFrontEdges (e->p1, it1, it2);
+ if (*it1 == e)left=*it2;
+ else left = *it1;
+ // printf("case 2 left edge %p\n",left);
+ right = findOptimalEdge(e->p2,left->othervertex(e->p1));
+ if (right)getFrontEdges (right->othervertex(e->p2),fe2);
+ break;
+ // left edge is angled with current edge
+ // we therefore find a new front edge in the
+ // mesh that allows to move to tag 3
+ case 1:
+ getFrontEdges (e->p2, it1, it2);
+ if (*it2 == e)right=*it1;
+ else right = *it2;
+ // printf("case 1 right edge %p %p\n",e,right);
+ left = findOptimalEdge( e->p1,right->othervertex(e->p2));
+ if (left)getFrontEdges (left->othervertex(e->p1),fe1);
+ break;
+ // no neighboring edge is angled with current edge
+ // we therefore find a new front edge in the
+ // mesh that allows to move to tag 3
+ case 0:
+ left = findOptimalEdge( e->p1,0);
+ if (left)right= findOptimalEdge( e->p2,left->othervertex(e->p1));
+ if(right)getFrontEdges (right->othervertex(e->p2),fe2);
+ if(left)getFrontEdges (left->othervertex(e->p1),fe1);
+ // printf("strange\n");
+ break;
+ default:
+ throw;
+ }
+
+ if ( fe1.size() || fe2.size() || !left || !right || !formQuad (e, left, right)){
+ // printf("front cloeses : algo has to be finished\n");
+ stat[tag].erase(stat[tag].begin());
+ stat[4].insert(e);
+ }
+ return false;
+}
+
+static int numQuads ( BDS_Mesh *m)
+{
+ int N = 0;
+ for (std::list<BDS_Face*> ::iterator it = m->triangles.begin() ; it != m->triangles.end() ; ++it){
+ if ((*it)->e4) N++;
+ }
+ return N;
+}
+
+int gmshQMorph (GFace *gf) {
+
+ // assert first that there exist a triangulation of
+ // the face
+ if (!gf->triangles.size()){
+ Msg::Warning("Cannot Quadrilaterize a face that has not been triangulated");
+ return -1;
+ }
+
+ // create data structures for mesh manipulations
+ std::list<GFace*> l; l.push_back(gf);
+ BDS_Mesh *pm = gmsh2BDS(l);
+
+ // create the front
+ gmshEdgeFront front (pm,gf);
+ front.initiate();
+
+ int ITER = 1;
+
+ // empty the front for front edges tagged 3, 2, 1 then 0
+ int _numQuads = 0;
+ while (1){
+ if (front.emptyFront ( 3 )){
+ if (front.emptyFront ( 2 )){
+ if (front.emptyFront ( 1 )){
+ if (front.emptyFront ( 0 )){
+ int ns;
+ smoothVertexPass(gf,*pm,ns,false);
+ printf("nex row iter %6d->>>\n",ITER);
+ front.initiate();
+ int _numQuadsNew = numQuads (pm);
+ if (front.edges.size() == 0 || _numQuads == _numQuadsNew)
+ break;
+ _numQuads = _numQuadsNew;
+ }
+ }
+ }
+ }
+ ITER++;
+ if (1 || ITER%100 == 0){
+ char name[256];
+ sprintf(name,"qmorph-face%d-iter%d.pos",gf->tag(),ITER);
+ outputScalarField(pm->triangles, name, 0);
+ }
+ // if (ITER == 1123)break;
+ }
+ // delete the BDS
+ delete pm;
+
+}
+
diff --git a/Mesh/meshGFaceQuadrilateralize.h b/Mesh/meshGFaceQuadrilateralize.h
new file mode 100644
index 0000000000..636adb9081
--- /dev/null
+++ b/Mesh/meshGFaceQuadrilateralize.h
@@ -0,0 +1,25 @@
+#ifndef _GMSH_QUADRILATERALIZE_
+#define _GMSH_QUADRILATERALIZE_
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+class GFace;
+int gmshQMorph (GFace *gf);
+
+#endif
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index fa995cd482..44b84b07fd 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.48 2008-06-05 11:52:50 samtech Exp $
+// $Id: meshGRegion.cpp,v 1.49 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -118,6 +118,12 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
std::vector<MVertex*> &numberedV)
{
+ // Improvement has to be done here :
+ // netgen splits some of the existing edges of the
+ // mesh. If those edges are classified on some
+ // model faces, new points SHOULD be classified
+ // on the model face and get the right set of parametric coordinates.
+
for(int i = numberedV.size(); i < out.numberofpoints; i++){
MVertex *v = new MVertex(out.pointlist[i * 3 + 0],
out.pointlist[i * 3 + 1],
@@ -158,13 +164,46 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[1] = numberedV[out.trifacelist[i * 3 + 1] - 1];
v[2] = numberedV[out.trifacelist[i * 3 + 2] - 1];
GFace *gf = gr->model()->getFaceByTag(out.trifacemarkerlist[i]);
+
+ double guess[2] = {0,0};
+ int Count = 0;
+ for(int j = 0; j < 3; j++){
+ if(v[j]->onWhat()->dim() == 2){
+ v[j]->getParameter(0,guess[0]);
+ v[j]->getParameter(1,guess[1]);
+ Count++;
+ }
+ }
+ 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
(std::find(v[j]->onWhat()->mesh_vertices.begin(),
v[j]->onWhat()->mesh_vertices.end(),
v[j]));
- MVertex *v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
+ MVertex *v1b;
+ if (CTX.mesh.order > 1){
+ // PARAMETRIC COORDINATES SHOULD BE SET for the vertex !!!!!!!!!!!!!
+ // This is not 100 % safe yet, so we reserve that operation for high order
+ // meshes.
+ GPoint gp = gf->closestPoint (SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
+ // printf("ah que coucou\n");
+ // const double U(0),V(0);
+ // MVertex *v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf,U,V);
+ Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher",gf->tag());
+ Msg::Debug("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());
+ // To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
+ v[j]->x() = gp.x();
+ v[j]->y() = gp.y();
+ v[j]->z() = gp.z();
+ v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf,gp.u(),gp.v());
+ }
+ else{
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z());
+ }
+
gf->mesh_vertices.push_back(v1b);
numberedV[out.trifacelist[i * 3 + j] - 1] = v1b;
delete v[j];
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index d340f2b483..cba41e37c8 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -1,4 +1,4 @@
-// $Id: FunctionSpace.cpp,v 1.5 2008-06-05 18:25:18 geuzaine Exp $
+// $Id: FunctionSpace.cpp,v 1.6 2008-06-10 08:37:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -585,6 +585,26 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
F.monomials = generatePascalTriangle(5);
F.points = gmshGeneratePointsTriangle(5, false);
break;
+ case MSH_TET_4 :
+ F.monomials = generatePascalTetrahedron(1);
+ F.points = gmshGeneratePointsTetrahedron(5, false);
+ break;
+ case MSH_TET_10 :
+ F.monomials = generatePascalTetrahedron(2);
+ F.points = gmshGeneratePointsTetrahedron(2, false);
+ break;
+ case MSH_TET_20 :
+ F.monomials = generatePascalTetrahedron(3);
+ F.points = gmshGeneratePointsTetrahedron(3, false);
+ break;
+ case MSH_TET_35 :
+ F.monomials = generatePascalTetrahedron(4);
+ F.points = gmshGeneratePointsTetrahedron(4, false);
+ break;
+ case MSH_TET_56 :
+ F.monomials = generatePascalTetrahedron(5);
+ F.points = gmshGeneratePointsTetrahedron(5, false);
+ break;
default :
throw;
}
--
GitLab