From e531a4deacd23e89b179d11dc5de11960d32fac1 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 28 Nov 2008 18:42:30 +0000
Subject: [PATCH] fix orientation of recombined quads + removed BDS recombine
---
Fltk/GUI.cpp | 2 +-
Geo/Makefile | 9 +-
Makefile | 2 +-
Mesh/BDS.cpp | 136 ----------------
Mesh/BDS.h | 2 -
Mesh/Makefile | 30 ++--
Mesh/Refine.cpp | 13 +-
Mesh/meshGFace.cpp | 35 +++--
Mesh/meshGFaceDelaunayInsertion.cpp | 234 ++++++++--------------------
Mesh/meshGFaceOptimize.cpp | 141 +++++++----------
Numeric/GmshMatrix.h | 154 +++++++++---------
Numeric/Makefile | 6 +-
doc/gmsh.html | 4 +-
13 files changed, 250 insertions(+), 518 deletions(-)
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index d5c7870885..356d4323d6 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -315,7 +315,7 @@ Context_Item menu_mesh[] = {
#if defined(HAVE_FOURIER_MODEL)
{"Reparameterize", (Fl_Callback *)mesh_parameterize_cb} ,
#endif
- // {"Reclassify", (Fl_Callback *)mesh_classify_cb} ,
+ //{"Reclassify", (Fl_Callback *)mesh_classify_cb} ,
{"Save", (Fl_Callback *)mesh_save_cb} ,
{0}
};
diff --git a/Geo/Makefile b/Geo/Makefile
index d9b895efe3..50527e0ca8 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -101,9 +101,9 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp GFaceCompound.h GFace.h GEntity.h \
GVertex.h SPoint2.h SVector3.h Pair.h ../Numeric/gmshAssembler.h \
../Numeric/gmshLinearSystem.h ../Numeric/gmshLaplace.h \
../Numeric/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
- ../Common/GmshMessage.h ../Common/Gmsh.h ../Common/GmshMessage.h \
- ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
- ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+ ../Common/GmshMessage.h ../Numeric/gmshFunction.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h \
+ ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
@@ -111,7 +111,8 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp GFaceCompound.h GFace.h GEntity.h \
../Numeric/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h ../Common/Octree.h \
../Common/OctreeInternals.h ../Numeric/gmshLinearSystemGmm.h \
- ../Numeric/gmshLinearSystem.h
+ ../Numeric/gmshLinearSystem.h ../Numeric/gmshLinearSystemFull.h \
+ ../Numeric/gmshLinearSystem.h ../Numeric/GmshMatrix.h
GRegion${OBJEXT}: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
diff --git a/Makefile b/Makefile
index a452e27c27..f29451a8b8 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ include variables
GMSH_MAJOR_VERSION = 2
GMSH_MINOR_VERSION = 2
-GMSH_PATCH_VERSION = 6
+GMSH_PATCH_VERSION = 7
GMSH_EXTRA_VERSION =
GMSH_VERSION = ${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}.${GMSH_PATCH_VERSION}${GMSH_EXTRA_VERSION}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 25a8cd271d..ede3847246 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -967,82 +967,6 @@ int BDS_Edge::numTriangles() const
return NT;
}
-// This function does actually the swap without taking into account
-// the feasability of the operation. Those conditions have to be
-// taken into account before doing the edge swap
-bool BDS_Mesh::recombine_edge(BDS_Edge *e)
-{
- /*
- p1
- / | \
- / | \
- op1/ 0 | 1 \op2
- \ | /
- \ | /
- \ p2/
-
- // op1 p1 op2
- // op1 op2 p2
- */
-
- // we test if the edge is deleted
- if(e->deleted)
- return false;
- if(e->numfaces() != 2 || e->numTriangles() != 2)
- return false;
- if(e->g && e->g->classif_degree == 1)
- return false;
-
- BDS_Point *p1 = e->p1;
- BDS_Point *p2 = e->p2;
- BDS_Point *op[2];
- e->oppositeof(op);
-
- BDS_Point *pts1[4];
- e->faces(0)->getNodes(pts1);
-
- BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
- BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
- BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
- BDS_Edge *op2_p2 = find_edge(op[1], p2, e->faces(1));
-
- BDS_GeomEntity *g = 0;
- if(e->faces(0)){
- g = e->faces(0)->g;
- del_face(e->faces(0));
- }
- // not a bug !!!
- if(e->faces(0)) {
- del_face(e->faces(0));
- }
- del_edge(e);
-
- int orientation = 0;
- for(int i = 0; i < 3; i++) {
- if(pts1[i] == p1) {
- if(pts1[(i + 1) % 3] == p2)
- orientation = 1;
- else
- orientation = -1;
- break;
- }
- }
-
- BDS_Face *f;
- if(orientation < 0)
- f = add_quadrangle(p1_op1, op1_p2, op2_p2, p1_op2);
- else
- f = add_quadrangle(p1_op1, p1_op2, op2_p2, op1_p2);
- f->g = g;
-
- p1->config_modified = true;
- p2->config_modified = true;
- op[0]->config_modified = true;
- op[1]->config_modified = true;
-
- return true;
-}
-
bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
{
if(e->numfaces() != 2)
@@ -1491,63 +1415,3 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
return true;
}
-
-struct recombine_T
-{
- const BDS_Edge *e;
- double angle;
- recombine_T(const BDS_Edge *_e);
- bool operator < (const recombine_T &other) const
- {
- return angle < other.angle;
- }
-};
-
-recombine_T::recombine_T(const BDS_Edge *_e)
- : e(_e)
-{
- BDS_Point *op[2];
- BDS_Point *p1 = e->p1;
- BDS_Point *p2 = e->p2;
- e->oppositeof(op);
- BDS_Vector v1(*p1, *op[0]);
- BDS_Vector v2(*op[0], *p2);
- BDS_Vector v3(*p2, *op[1]);
- BDS_Vector v4(*op[1], *p1);
- angle = fabs(90. - v1.angle_deg(v2));
- angle = std::max(fabs(90. - v2.angle_deg(v3)), angle);
- angle = std::max(fabs(90. - v3.angle_deg(v4)), angle);
- angle = std::max(fabs(90. - v4.angle_deg(v1)), angle);
-}
-
-void BDS_Mesh::recombineIntoQuads(const double angle_limit, GFace *gf)
-{
- Msg::Info("Recombining triangles for surface %d", gf->tag());
- for(int i = 0; i < 5; i++){
- std::set<recombine_T> pairs;
- for(std::list<BDS_Edge*>::const_iterator it = edges.begin();
- it != edges.end(); ++it){
- if(!(*it)->deleted && (*it)->numfaces () == 2)
- pairs.insert(recombine_T(*it));
- }
-
- bool rec = false;
- std::set<recombine_T>::iterator itp = pairs.begin();
- while(itp != pairs.end()){
- // recombine if difference between max quad angle and right
- // angle is smaller than tol
- if(itp->angle < gf->meshAttributes.recombineAngle)
- rec |= recombine_edge((BDS_Edge*)itp->e);
- ++itp;
- }
-
- if(!rec) break;
-
- std::set<BDS_Point*, PointLessThan>::iterator itpt = points.begin();
- while(itpt != points.end()){
- smooth_point_parametric(*itpt, gf);
- ++itpt;
- }
- }
-}
-
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index f5f1501cae..cebcaf348f 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -471,10 +471,8 @@ public:
bool split_edge(BDS_Edge *, BDS_Point *);
bool split_face(BDS_Face *, BDS_Point *);
bool edge_constraint(BDS_Point *p1, BDS_Point *p2);
- bool recombine_edge(BDS_Edge *e);
// Global operators
void cleanup();
- void recombineIntoQuads(const double angle, GFace *gf);
};
void outputScalarField(std::list<BDS_Face*> t, const char *fn, int param, GFace *gf=0);
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 1993e4a2f0..84a4c4b085 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -91,7 +91,8 @@ Generator${OBJEXT}: Generator.cpp ../Common/GmshMessage.h ../Numeric/Numeric.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
meshGEdge.h meshGFace.h meshGFaceBDS.h meshGRegion.h BackgroundMesh.h \
- BoundaryLayers.h HighOrder.h ../Post/PView.h ../Post/PViewData.h
+ BoundaryLayers.h HighOrder.h Generator.h ../Post/PView.h \
+ ../Post/PViewData.h
Field${OBJEXT}: Field.cpp ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h Field.h ../Post/PView.h ../Geo/SPoint3.h \
../Geo/GeoInterpolation.h ../Geo/Geo.h ../Common/GmshDefines.h \
@@ -127,8 +128,8 @@ gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp HighOrder.h \
../Numeric/GmshMatrix.h meshGFaceDelaunayInsertion.h \
gmshSmoothHighOrder.h ../Numeric/gmshAssembler.h \
../Numeric/gmshLinearSystem.h ../Numeric/gmshLaplace.h \
- ../Numeric/gmshTermOfFormulation.h ../Common/Gmsh.h \
- ../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
+ ../Numeric/gmshTermOfFormulation.h ../Numeric/gmshFunction.h \
+ ../Common/Gmsh.h ../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
../Numeric/gmshElasticity.h ../Numeric/gmshTermOfFormulation.h \
../Numeric/GmshMatrix.h ../Numeric/gmshLinearSystemGmm.h \
../Numeric/gmshLinearSystem.h ../Common/Context.h ../Geo/CGNSOptions.h \
@@ -227,20 +228,19 @@ meshGFaceBDS${OBJEXT}: meshGFaceBDS.cpp meshGFace.h meshGFaceOptimize.h \
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h \
qualityMeasures.h Field.h ../Common/OS.h
-meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp BDS.h \
- ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
- ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
- ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
- ../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
+meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
BackgroundMesh.h meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h meshGFaceOptimize.h \
- meshGFace.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
- ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
+ ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Numeric/GmshMatrix.h meshGFaceOptimize.h meshGFace.h ../Geo/GFace.h \
+ ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
+ ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+ ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+ ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+ ../Geo/Pair.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h
meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
@@ -252,7 +252,7 @@ meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/Pair.h BackgroundMesh.h
+ ../Geo/Pair.h BackgroundMesh.h Generator.h
meshGFaceQuadrilateralize${OBJEXT}: meshGFaceQuadrilateralize.cpp \
meshGFaceQuadrilateralize.h ../Common/GmshMessage.h \
../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
diff --git a/Mesh/Refine.cpp b/Mesh/Refine.cpp
index 2346b42580..ad3ad79089 100644
--- a/Mesh/Refine.cpp
+++ b/Mesh/Refine.cpp
@@ -44,10 +44,10 @@ static void Subdivide(GEdge *ge)
ge->deleteVertexArrays();
}
-static void Subdivide(GFace *gf,bool splitTrianglesIntoQuads)
+static void Subdivide(GFace *gf, bool splitTrianglesIntoQuads)
{
- if (! splitTrianglesIntoQuads ){
+ if(!splitTrianglesIntoQuads){
std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
@@ -86,14 +86,13 @@ static void Subdivide(GFace *gf,bool splitTrianglesIntoQuads)
MTriangle *t = gf->triangles[i];
if(t->getNumVertices() == 6){
SPoint2 pt, temp;
- bool reparamOK = true;
- for(int k = 0; k<6; k++){
+ for(int k = 0; k < 6; k++){
reparamMeshVertexOnFace(t->getVertex(k), gf, temp);
- pt[0] += temp[0]/6.;
- pt[1] += temp[1]/6.;
+ pt[0] += temp[0] / 6.;
+ pt[1] += temp[1] / 6.;
}
GPoint gp = gf->point(pt);
- MFaceVertex *newv = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pt[0],pt[1]);
+ MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, pt[0], pt[1]);
gf->mesh_vertices.push_back(newv);
quadrangles2.push_back
(new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5)));
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 712fd40eaa..e6a853f009 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -700,9 +700,6 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
gmshOptimizeMeshBDS(gf, *m, 2);
gmshRefineMeshBDS (gf,*m, CTX.mesh.refine_steps, false);
gmshOptimizeMeshBDS(gf, *m, 2);
- if (gf->meshAttributes.recombine){
- m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
- }
}
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
@@ -750,11 +747,12 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
if(AlgoDelaunay2D(gf)){
- gmshBowyerWatson(gf);
- for (int i=0;i<CTX.mesh.nb_smoothing;i++)laplaceSmoothing(gf);
- if (gf->meshAttributes.recombine){
- gmshRecombineIntoQuads(gf);
- }
+ if (CTX.mesh.algo2d == ALGO_2D_FRONTAL)
+ gmshBowyerWatsonFrontal(gf);
+ else
+ gmshBowyerWatson(gf);
+ for (int i = 0; i < CTX.mesh.nb_smoothing; i++)
+ laplaceSmoothing(gf);
}
else if (debug){
char name[256];
@@ -763,12 +761,15 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
sprintf(name, "param%d.pos", gf->tag());
outputScalarField(m->triangles, name,1);
}
-
+
// delete the mesh
delete m;
delete [] U_;
delete [] V_;
+ if (gf->meshAttributes.recombine)
+ gmshRecombineIntoQuads(gf);
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
gf->meshStatistics.best_element_shape,
@@ -1233,9 +1234,6 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
gmshOptimizeMeshBDS(gf, *m, 2);
gmshRefineMeshBDS (gf, *m, -CTX.mesh.refine_steps, false);
gmshOptimizeMeshBDS(gf, *m, 2, &recover_map);
- if (gf->meshAttributes.recombine){
- m->recombineIntoQuads (gf->meshAttributes.recombineAngle, gf);
- }
// compute mesh statistics
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
gf->meshStatistics.longest_edge_length,
@@ -1294,13 +1292,20 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
if (AlgoDelaunay2D(gf)){
- gmshBowyerWatson(gf);
- laplaceSmoothing(gf);
+ if (CTX.mesh.algo2d == ALGO_2D_FRONTAL)
+ gmshBowyerWatsonFrontal(gf);
+ else
+ gmshBowyerWatson(gf);
+ for (int i = 0; i < CTX.mesh.nb_smoothing; i++)
+ laplaceSmoothing(gf);
}
-
+
// delete the mesh
delete m;
+ if (gf->meshAttributes.recombine)
+ gmshRecombineIntoQuads(gf);
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
gf->meshStatistics.best_element_shape,
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index b3315b208a..019b8cdb23 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -3,7 +3,9 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#include "BDS.h"
+#include <set>
+#include <map>
+#include <algorithm>
#include "BackgroundMesh.h"
#include "meshGFaceDelaunayInsertion.h"
#include "meshGFaceOptimize.h"
@@ -11,27 +13,21 @@
#include "GFace.h"
#include "Numeric.h"
#include "GmshMessage.h"
-#include <set>
-#include <map>
-#include <algorithm>
-
-#include "Context.h"
-extern Context_T CTX;
const double LIMIT_ = 0.5 * sqrt(2.0);
-/* This function controls the frontal algorithm */
-static bool isActive ( MTri3 *t , double limit_, int &active){
+// This function controls the frontal algorithm
+static bool isActive(MTri3 *t, double limit_, int &active)
+{
if (t->isDeleted()) return false;
- for (active=0;active<3;active++){
+ for (active = 0; active < 3; active++){
MTri3 *neigh = t->getNeigh(active);
- if (!neigh || neigh->getRadius() < limit_)return true;
- //if (!neigh || neigh->done)return true;
+ if (!neigh || neigh->getRadius() < limit_) return true;
+ //if (!neigh || neigh->done) return true;
}
return false;
}
-
void circumCenterXY(double *p1, double *p2, double *p3, double *res)
{
double d, a1, a2, a3;
@@ -417,7 +413,6 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
int k = 0;
std::list<edgeXface>::iterator it = shell.begin();
-
bool onePointIsTooClose = false;
while (it != shell.end()){
MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
@@ -430,36 +425,37 @@ bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
MTri3 *t4 = new MTri3(t, LL);
- double d1 = sqrt((it->v[0]->x()-v->x())*(it->v[0]->x()-v->x())+
- (it->v[0]->y()-v->y())*(it->v[0]->y()-v->y())+
- (it->v[0]->z()-v->z())*(it->v[0]->z()-v->z()));
- double d2 = sqrt((it->v[1]->x()-v->x())*(it->v[1]->x()-v->x())+
- (it->v[1]->y()-v->y())*(it->v[1]->y()-v->y())+
- (it->v[1]->z()-v->z())*(it->v[1]->z()-v->z()));
- if (d1 < LL*.25 ||d2 < LL*.25)onePointIsTooClose = true;
+ double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
+ (it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
+ (it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
+ double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
+ (it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
+ (it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
+ if (d1 < LL * .25 || d2 < LL * .25) onePointIsTooClose = true;
- // if (t4->getRadius () < LIMIT_ / 2)onePointIsTooClose = true;
+ // if (t4->getRadius () < LIMIT_ / 2) onePointIsTooClose = true;
newTris[k++] = t4;
- // all new triangles are pushed front in order to
- // ba able to destroy them if the cavity is not
- // star shaped around the new vertex.
- new_cavity.push_back (t4);
+ // all new triangles are pushed front in order to be able to
+ // destroy them if the cavity is not star shaped around the new
+ // vertex.
+ new_cavity.push_back(t4);
MTri3 *otherSide = it->t1->getNeigh(it->i1);
if (otherSide)
- new_cavity.push_back (otherSide);
- double ss = fabs(getSurfUV(t4->tri(),Us,Vs));
- if (ss < 1.e-25)ss = 1256172121;
+ new_cavity.push_back(otherSide);
+ double ss = fabs(getSurfUV(t4->tri(), Us, Vs));
+ if (ss < 1.e-25) ss = 1256172121;
newVolume += ss;
++it;
}
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 && !onePointIsTooClose){
- connectTris(new_cavity.begin(),new_cavity.end());
- allTets.insert(newTris,newTris+shell.size());
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 &&
+ !onePointIsTooClose){
+ connectTris(new_cavity.begin(), new_cavity.end());
+ allTets.insert(newTris, newTris + shell.size());
if (activeTets){
- for (std::list<MTri3*>::iterator i = new_cavity.begin();i!=new_cavity.end();++i){
+ for (std::list<MTri3*>::iterator i = new_cavity.begin(); i != new_cavity.end(); ++i){
int active_edge;
- if(isActive(*i,LIMIT_,active_edge) && (*i)->getRadius() > LIMIT_){
+ if(isActive(*i, LIMIT_, active_edge) && (*i)->getRadius() > LIMIT_){
if ((*activeTets).find(*i) == (*activeTets).end())
(*activeTets).insert(*i);
}
@@ -527,8 +523,7 @@ void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
fclose (ff);
}
-
-static void insertAPoint(GFace *gf,
+static void insertAPoint(GFace *gf,
std::set<MTri3*,compareTri3Ptr>::iterator it,
double center[2],
double metric[3],
@@ -573,7 +568,7 @@ static void insertAPoint(GFace *gf,
// printf("the new point is %g %g\n",p.x(),p.y());
MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
- // printf("%g %g -> %g %g %g\n",center[0], center[1],v->x(), v->y(), v->z());
+ // printf("%g %g -> %g %g %g\n",center[0], center[1],v->x(), v->y(), v->z());
v->setNum(Us.size());
double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] +
@@ -616,10 +611,6 @@ static void insertAPoint(GFace *gf,
void gmshBowyerWatson(GFace *gf)
{
- if (CTX.mesh.algo2d == ALGO_2D_FRONTAL){
- gmshBowyerWatsonFrontal(gf);
- return;
- }
std::set<MTri3*,compareTri3Ptr> AllTris;
std::vector<double> vSizes, vSizesBGM, Us, Vs;
@@ -671,7 +662,6 @@ void gmshBowyerWatson(GFace *gf)
We use the approach of Rebay (JCP 1993) that we extend to anisotropic metrics
The point isertion is done on the Voronoi Edges
-
*/
static double length_metric ( const double p[2], const double q[2], const double metric[3])
@@ -705,8 +695,8 @@ static double length_metric ( const double p[2], const double q[2], const double
*/
-
-void gmshBowyerWatsonFrontal(GFace *gf){
+void gmshBowyerWatsonFrontal(GFace *gf)
+{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
std::vector<double> vSizes, vSizesBGM, Us, Vs;
@@ -732,7 +722,8 @@ void gmshBowyerWatsonFrontal(GFace *gf){
ActiveTris.erase(ActiveTris.begin());
// printf("active_tris.size = %d\n",ActiveTris.size());
- if (!worst->isDeleted() && isActive(worst,LIMIT_,active_edge) && worst->getRadius() > LIMIT_){
+ if (!worst->isDeleted() && isActive(worst, LIMIT_, active_edge) &&
+ worst->getRadius() > LIMIT_){
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
@@ -751,43 +742,43 @@ void gmshBowyerWatsonFrontal(GFace *gf){
// compute the middle point of the edge
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
- // printf("the active edge is %d : %g %g -> %g %g\n",active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
- // base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
+ // printf("the active edge is %d : %g %g -> %g %g\n",
+ // active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
+ // base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
double P[2] = {Us[base->getVertex(ip1)->getNum()],
Vs[base->getVertex(ip1)->getNum()]};
double Q[2] = {Us[base->getVertex(ip2)->getNum()],
Vs[base->getVertex(ip2)->getNum()]};
- double midpoint[2] = {0.5*(P[0]+Q[0]),0.5*(P[1]+Q[1])};
+ double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
// now we have the edge center and the center of the circumcircle,
// we try to find a point that would produce a perfect triangle while
// connecting the 2 points of the active edge
- double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
- double norm = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
- dir[0]/=norm;
- dir[1]/=norm;
- const double RATIO = sqrt( dir[0]*dir[0]*metric[0]+
- 2*dir[1]*dir[0]*metric[1]+
- dir[1]*dir[1]*metric[2]);
+ double dir[2] = {center[0] - midpoint[0], center[1] - midpoint[1]};
+ double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
+ dir[0] /= norm;
+ dir[1] /= norm;
+ const double RATIO = sqrt(dir[0] * dir[0] * metric[0] +
+ 2 * dir[1] * dir[0] * metric[1] +
+ dir[1] * dir[1] * metric[2]);
-
- const double p = 0.5*length_metric (P,Q,metric);// / RATIO;
- const double q = length_metric (center,midpoint,metric);
- const double rhoM1 = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
- const double rhoM2 = 0.5 * (vSizesBGM[base->getVertex(ip1)->getNum()] + vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ const double p = 0.5 * length_metric(P, Q, metric); // / RATIO;
+ const double q = length_metric (center, midpoint, metric);
+ const double rhoM1 = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] +
+ vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ const double rhoM2 = 0.5 * (vSizesBGM[base->getVertex(ip1)->getNum()] +
+ vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
- // const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
+ // const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] +
+ // vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
- const double rhoM_hat = std::min(std::max(rhoM,p),(p*p+q*q)/(2*q));
- const double d = (rhoM_hat + sqrt (rhoM_hat*rhoM_hat - p*p))/RATIO;
+ const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
+ const double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
- double newPoint[2] =
- {
- midpoint[0] + d * dir[0],
- midpoint[1] + d * dir[1]
- };
- insertAPoint(gf,AllTris.end(),newPoint,metric,Us,Vs,vSizes,vSizesBGM,AllTris,&ActiveTris,worst);
+ double newPoint[2] = {midpoint[0] + d * dir[0], midpoint[1] + d * dir[1]};
+ insertAPoint(gf, AllTris.end(), newPoint, metric, Us, Vs, vSizes, vSizesBGM,
+ AllTris, &ActiveTris, worst);
}
// if(ITER % 1000== 0){
// char name[245];
@@ -796,107 +787,8 @@ void gmshBowyerWatsonFrontal(GFace *gf){
// }
}
- char name[245];
- sprintf(name,"frontal%d.pos",gf->tag());
- _printTris (name, AllTris, Us,Vs);
- transferDataStructure(gf, AllTris);
-}
-
-
-/*
-void gmshBowyerWatsonFrontal(GFace *gf){
-
- std::set<MTri3*,compareTri3Ptr> AllTris;
- std::vector<double> vSizes, vSizesBGM, Us, Vs;
- buidMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
-
- // delaunise the initial mesh
- int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
- Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
-
-
- int ITER = 0, active_edge;
- while (1) {
- // compute active triangle
- std::set<MTri3*,compareTri3Ptr> ActiveTris;
- std::set<MTri3*,compareTri3Ptr>::reverse_iterator it = AllTris.rbegin();
- for ( ; it!=AllTris.rend();++it){
- if ((*it)->getRadius() < LIMIT_)
- (*it)->done = true;
- else if(isActive(*it,LIMIT_,active_edge))
- ActiveTris.insert(*it);
- }
- if (ActiveTris.size() == 0)break;
-
- // insert points
- while (1){
- if (!ActiveTris.size())break;
- MTri3 *worst = (*ActiveTris.begin());
- ActiveTris.erase(ActiveTris.begin());
- // printf("active_tris.size = %d\n",ActiveTris.size());
-
- if (!worst->isDeleted() && isActive(worst,LIMIT_,active_edge) && worst->getRadius() > LIMIT_){
- if(ITER++ % 5000 == 0)
- Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
- vSizes.size(), worst->getRadius());
- // compute circum center of that guy
- double center[2],uv[2],metric[3],r2;
- MTriangle *base = worst->tri();
- circUV(base, Us, Vs, center, gf);
- double pa[2] = {(Us[base->getVertex(0)->getNum()] +
- Us[base->getVertex(1)->getNum()] +
- Us[base->getVertex(2)->getNum()]) / 3.,
- (Vs[base->getVertex(0)->getNum()] +
- Vs[base->getVertex(1)->getNum()] +
- Vs[base->getVertex(2)->getNum()]) / 3.};
- buildMetric(gf, pa, metric);
- circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
- // compute the middle point of the edge
- int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
- int ip2 = active_edge;
- // printf("the active edge is %d : %g %g -> %g %g\n",active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
- // base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
- double P[2] = {Us[base->getVertex(ip1)->getNum()],
- Vs[base->getVertex(ip1)->getNum()]};
- double Q[2] = {Us[base->getVertex(ip2)->getNum()],
- Vs[base->getVertex(ip2)->getNum()]};
- double midpoint[2] = {0.5*(P[0]+Q[0]),0.5*(P[1]+Q[1])};
-
- // now we have the edge center and the center of the circumcircle,
- // we try to find a point that would produce a perfect triangle while
- // connecting the 2 points of the active edge
-
- double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
- double norm = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
- dir[0]/=norm;
- dir[1]/=norm;
- const double RATIO = sqrt( dir[0]*dir[0]*metric[0]+
- 2*dir[1]*dir[0]*metric[1]+
- dir[1]*dir[1]*metric[2]);
-
-
- const double p = 0.5*length_metric (P,Q,metric);// / RATIO;
- const double q = length_metric (center,midpoint,metric);
- const double rhoM1 = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
- const double rhoM2 = 0.5 * (vSizesBGM[base->getVertex(ip1)->getNum()] + vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
- const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
- // const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
-
- const double rhoM_hat = std::min(std::max(rhoM,p),(p*p+q*q)/(2*q));
- const double d = (rhoM_hat + sqrt (rhoM_hat*rhoM_hat - p*p))/RATIO;
-
- double newPoint[2] =
- {
- midpoint[0] + d * dir[0],
- midpoint[1] + d * dir[1]
- };
- insertAPoint(gf,AllTris.end(),newPoint,metric,Us,Vs,vSizes,vSizesBGM,AllTris,&ActiveTris,worst);
- }
- }
- }
- _printTris ("frontal.pos", AllTris, Us,Vs);
+ // char name[245];
+ // sprintf(name,"frontal%d.pos", gf->tag());
+ // _printTris (name, AllTris, Us,Vs);
transferDataStructure(gf, AllTris);
}
-*/
-
-
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e8bb0bcbc8..978ef56900 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -26,7 +26,7 @@ static void setLcsInit(MTriangle *t, std::map<MVertex*, double> &vSizes)
}
}
-static void setLcs(MTriangle *t, std::map<MVertex*,double> &vSizes)
+static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes)
{
for (int i = 0; i < 3; i++){
for (int j = i + 1; j < 3; j++){
@@ -44,20 +44,6 @@ static void setLcs(MTriangle *t, std::map<MVertex*,double> &vSizes)
}
}
-static void setLcs(MLine *t, std::map<MVertex*,double> &vSizes)
-{
- MVertex *vi = t->getVertex(0);
- MVertex *vj = t->getVertex(1);
- double dx = vi->x()-vj->x();
- double dy = vi->y()-vj->y();
- double dz = vi->z()-vj->z();
- double l = sqrt(dx * dx + dy * dy + dz * dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
- if (iti->second < 0 || iti->second < l) iti->second = l;
- if (itj->second < 0 || itj->second < l) itj->second = l;
-}
-
void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
@@ -70,15 +56,6 @@ void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr
for (unsigned int i = 0;i < gf->triangles.size(); i++)
setLcsInit(gf->triangles[i], vSizesMap);
-// std::list<GEdge*>::iterator it = edges.begin();
-// while (it != edges.end()){
-// GEdge *ge = *it ;
-// for (int i=0;i<ge->lines.size();i++){
-// setLcs(ge->lines[i], vSizesMap);
-// }
-// ++it;
-// }
-
for (unsigned int i = 0;i < gf->triangles.size(); i++)
setLcs(gf->triangles[i], vSizesMap);
@@ -111,7 +88,7 @@ void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr
AllTris.insert(new MTri3(gf->triangles[i], lc));
}
gf->triangles.clear();
- connectTriangles(AllTris );
+ connectTriangles(AllTris);
}
void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris)
@@ -127,9 +104,9 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris)
AllTris.erase(AllTris.begin());
}
}
+
template <class T>
-void buildVertexToElement(std::vector<T*> &eles,
- v2t_cont &adj)
+void buildVertexToElement(std::vector<T*> &eles, v2t_cont &adj)
{
for (unsigned int i = 0; i < eles.size(); i++){
T *t = eles[i];
@@ -148,13 +125,12 @@ void buildVertexToElement(std::vector<T*> &eles,
}
}
-
-void buildVertexToTriangle(std::vector<MTriangle*> &eles, v2t_cont &adj){
+void buildVertexToTriangle(std::vector<MTriangle*> &eles, v2t_cont &adj)
+{
adj.clear();
buildVertexToElement(eles,adj);
}
-
template <class T>
void buildEdgeToElement(std::vector<T*> &triangles, e2t_cont &adj)
{
@@ -187,8 +163,6 @@ void buildEdgeToTriangle(std::vector<MTriangle*> &tris, e2t_cont &adj){
buildEdgeToElement(tris, adj);
}
-
-
void parametricCoordinates(MElement *t, GFace *gf, double u[4], double v[4])
{
for (unsigned int j = 0; j < t->getNumVertices(); j++){
@@ -232,9 +206,9 @@ void laplaceSmoothing(GFace *gf)
// have to test validity !
}
if (fact != 0.0){
- ver->setParameter(0, cu/fact);
- ver->setParameter(1, cv/fact);
- GPoint pt = gf->point(SPoint2(cu/fact, cv/fact));
+ ver->setParameter(0, cu / fact);
+ ver->setParameter(1, cv / fact);
+ GPoint pt = gf->point(SPoint2(cu / fact, cv / fact));
ver->x() = pt.x();
ver->y() = pt.y();
ver->z() = pt.z();
@@ -439,7 +413,7 @@ bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
if (al <= lMax) return false;
MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
- MVertex *v4 = 0 ;
+ MVertex *v4 = 0;
for (int i = 0; i < 3; i++)
if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
v4 = t2->tri()->getVertex(i);
@@ -638,7 +612,7 @@ bool gmshVertexCollapse(const double lMin, MTri3 *t1, GFace *gf,
MVertex *v;
std::vector<MTri3*> cavity;
std::vector<MTri3*> outside;
- std::vector<MVertex*> ring ;
+ std::vector<MVertex*> ring;
if (!gmshBuildVertexCavity(t1, iLocalVertex, &v, cavity, outside, ring)) return false;
@@ -698,7 +672,7 @@ int edgeSwapPass (GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
const std::vector<double> &Us, const std::vector<double> &Vs,
const std::vector<double> &vSizes, const std::vector<double> &vSizesBGM)
{
- typedef std::set<MTri3*, compareTri3Ptr> CONTAINER ;
+ typedef std::set<MTri3*, compareTri3Ptr> CONTAINER;
int nbSwapTot = 0;
std::set<swapquad> configs;
@@ -734,7 +708,7 @@ int edgeSplitPass(double maxLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &allT
std::vector<double> &Us, std::vector<double> &Vs,
std::vector<double> &vSizes, std::vector<double> &vSizesBGM)
{
- typedef std::set<MTri3*, compareTri3Ptr> CONTAINER ;
+ typedef std::set<MTri3*, compareTri3Ptr> CONTAINER;
std::vector<MTri3*> newTris;
int nbSplit = 0;
@@ -766,7 +740,7 @@ int edgeCollapsePass(double minLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &a
std::vector<double> &Us, std::vector<double> &Vs,
std::vector<double> &vSizes, std::vector<double> &vSizesBGM)
{
- typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ typedef std::set<MTri3*,compareTri3Ptr> CONTAINER;
std::vector<MTri3*> newTris;
int nbCollapse = 0;
@@ -795,31 +769,30 @@ int edgeCollapsePass(double minLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &a
return nbCollapse;
}
-extern double angle3Points ( MVertex *p1, MVertex *p2, MVertex *p3 );
+extern double angle3Points(MVertex *p1, MVertex *p2, MVertex *p3);
struct recombine_triangle
{
MElement *t1, *t2;
double angle;
- MVertex *n1,*n2,*n3,*n4;
+ MVertex *n1, *n2, *n3, *n4;
recombine_triangle(const MEdge &me, MElement *_t1, MElement *_t2)
- : t1(_t1),t2(_t2)
+ : t1(_t1), t2(_t2)
{
n1 = me.getVertex(0);
n2 = me.getVertex(1);
- if (t1->getVertex (0) != n1 && t1->getVertex (0) != n2)n3 = t1->getVertex(0);
- else if (t1->getVertex (1) != n1 && t1->getVertex (1) != n2)n3 = t1->getVertex(1);
- else if (t1->getVertex (2) != n1 && t1->getVertex (2) != n2)n3 = t1->getVertex(2);
- if (t2->getVertex (0) != n1 && t2->getVertex (0) != n2)n4 = t2->getVertex(0);
- else if (t2->getVertex (1) != n1 && t2->getVertex (1) != n2)n4 = t2->getVertex(1);
- else if (t2->getVertex (2) != n1 && t2->getVertex (2) != n2)n4 = t2->getVertex(2);
-
- double a1 = 180*angle3Points(n1,n4,n2)/M_PI;
- double a2 = 180*angle3Points(n4,n2,n3)/M_PI;
- double a3 = 180*angle3Points(n2,n3,n1)/M_PI;
- double a4 = 180*angle3Points(n3,n1,n4)/M_PI;
- // printf("%g %g %g %g\n",a1,a2,a3,a4);
+ if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
+ else if (t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
+ else if (t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
+ if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
+ else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
+ else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
+
+ double a1 = 180 * angle3Points(n1, n4, n2) / M_PI;
+ double a2 = 180 * angle3Points(n4, n2, n3) / M_PI;
+ double a3 = 180 * angle3Points(n2, n3, n1) / M_PI;
+ double a4 = 180 * angle3Points(n3, n1, n4) / M_PI;
angle = fabs(90. - a1);
angle = std::max(fabs(90. - a2),angle);
angle = std::max(fabs(90. - a3),angle);
@@ -831,14 +804,10 @@ struct recombine_triangle
}
};
-
-void gmshQuadrilateralizeAlgoBrutal(GFace *gf){
-}
-
-void _gmshRecombineIntoQuads(GFace *gf)
+static void _gmshRecombineIntoQuads(GFace *gf)
{
e2t_cont adj;
- std::set<MElement*> _touched;
+ std::set<MElement*> touched;
std::set<recombine_triangle> pairs;
buildEdgeToElement(gf->triangles, adj);
e2t_cont::iterator it = adj.begin();
@@ -857,40 +826,50 @@ void _gmshRecombineIntoQuads(GFace *gf)
if(itp->angle < gf->meshAttributes.recombineAngle){
MElement *t1 = itp->t1;
MElement *t2 = itp->t2;
- if (_touched.find(t1) == _touched.end() &&
- _touched.find(t2) == _touched.end()){
- _touched.insert(t1);
- _touched.insert(t2);
- MQuadrangle *q = new MQuadrangle (itp->n1,itp->n3,itp->n2,itp->n4);
+ if (touched.find(t1) == touched.end() &&
+ touched.find(t2) == touched.end()){
+ touched.insert(t1);
+ touched.insert(t2);
+
+ int orientation = 0;
+ for(int i = 0; i < 3; i++) {
+ if(t1->getVertex(i) == itp->n1) {
+ if(t1->getVertex((i + 1) % 3) == itp->n2)
+ orientation = 1;
+ else
+ orientation = -1;
+ break;
+ }
+ }
+ MQuadrangle *q;
+ if(orientation < 0)
+ q = new MQuadrangle(itp->n1, itp->n3, itp->n2, itp->n4);
+ else
+ q = new MQuadrangle(itp->n1, itp->n4, itp->n2, itp->n3);
gf->quadrangles.push_back(q);
}
}
++itp;
}
- std::vector<MTriangle*> _newt;
- for ( int i = 0 ; i<gf->triangles.size();i++){
- if (_touched.find(gf->triangles[i]) == _touched.end()){
- _newt.push_back(gf->triangles[i]);
+ std::vector<MTriangle*> triangles2;
+ for (int i = 0; i < gf->triangles.size(); i++){
+ if (touched.find(gf->triangles[i]) == touched.end()){
+ triangles2.push_back(gf->triangles[i]);
}
else {
delete gf->triangles[i];
}
}
- gf->triangles = _newt;
-
+ gf->triangles = triangles2;
}
-void gmshRecombineIntoQuads(GFace *gf){
- _gmshRecombineIntoQuads (gf);
+void gmshRecombineIntoQuads(GFace *gf)
+{
+ _gmshRecombineIntoQuads(gf);
laplaceSmoothing(gf);
- _gmshRecombineIntoQuads (gf);
+ _gmshRecombineIntoQuads(gf);
laplaceSmoothing(gf);
- _gmshRecombineIntoQuads (gf);
+ _gmshRecombineIntoQuads(gf);
laplaceSmoothing(gf);
}
-
-
-
-
-
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index f4e3b37eaf..1fa555f193 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -60,86 +60,80 @@ class Gmsh_Matrix
{
private:
int r, c;
-
- public:
void _back_substitution(int *indx, SCALAR *b)
{
- int i,ii=-1,ip,j;
+ int i, ii = -1, ip, j;
SCALAR sum;
-
- for(i=0; i< c; i++){
- ip=indx[i];
- sum=b[ip];
- b[ip]=b[i];
+ for(i = 0; i < c; i++){
+ ip = indx[i];
+ sum = b[ip];
+ b[ip] = b[i];
if(ii != -1)
- for(j=ii; j <= i-1; j++) sum -= (*this)(i,j)*b[j];
- else if (sum) ii=i;
- b[i]=sum;
+ for(j = ii; j <= i - 1; j++) sum -= (*this)(i, j) * b[j];
+ else if (sum) ii = i;
+ b[i] = sum;
}
- for(i=c-1; i>=0; i--){
- sum=b[i];
- for(j=i+1;j<c;j++) sum -= (*this)(i,j)*b[j];
- b[i]=sum/(*this)(i,i);
+ for(i = c - 1; i >= 0; i--){
+ sum = b[i];
+ for(j = i + 1; j < c; j++) sum -= (*this)(i, j) * b[j];
+ b[i] = sum / (*this)(i, i);
}
}
-
-
- bool _lu_decomposition(int *indx , SCALAR & determinant)
+ bool _lu_decomposition(int *indx , SCALAR &determinant)
{
if (r != c)
Msg::Fatal("Gmsh_Matrix::_lu_decomposition : cannot lu decompose a non-square matrix");
- int i, imax, j,k;
- SCALAR big,dum,sum,temp;
- SCALAR *vv = new SCALAR [c];
-
- determinant=1.0;
- for(i=0; i<c; i++){
- big=0.0;
- for(j=0;j<c; j++)
- if((temp=fabs((*this)(i,j))) > big) big=temp;
- if(big==0.0) {
+ int i, imax, j, k;
+ SCALAR big, dum, sum, temp;
+ SCALAR *vv = new SCALAR[c];
+ determinant = 1.0;
+ for(i = 0; i < c; i++){
+ big = 0.0;
+ for(j = 0; j < c; j++)
+ if((temp=fabs((*this)(i, j))) > big) big = temp;
+ if(big == 0.0) {
return false;
big = 1.e-12;
}
- vv[i]=1.0/big;
+ vv[i] = 1.0 / big;
}
- for(j=0; j<c;j++){
- for(i=0;i<j;i++){
- sum=(*this)(i,j);
- for(k=0;k<i;k++) sum -= (*this)(i,k)*(*this)(k,j);
- (*this)(i,j) = sum;
+ for(j = 0; j < c; j++){
+ for(i = 0; i < j; i++){
+ sum=(*this)(i, j);
+ for(k = 0; k < i; k++) sum -= (*this)(i, k)*(*this)(k, j);
+ (*this)(i, j) = sum;
}
- big=0.0;
- for(i=j; i<c;i++){
- sum=(*this)(i,j);
- for(k=0;k<j;k++)
- sum -= (*this)(i,k)*(*this)(k,j);
- (*this)(i,j)=sum;
- if((dum=vv[i]*fabs(sum)) >= big){
- big=dum;
- imax=i;
+ big = 0.0;
+ for(i = j; i < c; i++){
+ sum = (*this)(i, j);
+ for(k = 0; k < j; k++)
+ sum -= (*this)(i, k) * (*this)(k, j);
+ (*this)(i, j) = sum;
+ if((dum = vv[i] * fabs(sum)) >= big){
+ big = dum;
+ imax = i;
}
}
if(j != imax){
- for(k=0; k <c; k++){
- dum=(*this)(imax,k);
- (*this)(imax,k)=(*this)(j,k);
- (*this)(j,k) = dum;
+ for(k = 0; k < c; k++){
+ dum = (*this)(imax, k);
+ (*this)(imax, k) = (*this)(j, k);
+ (*this)(j, k) = dum;
}
determinant = -(determinant);
- vv[imax]=vv[j];
+ vv[imax] = vv[j];
}
- indx[j]=imax;
- if( (*this)(j,j) == 0.0) (*this)(j,j) = 1.e-20;
- if(j !=c){
- dum=1.0/((*this)(j,j));
- for(i=j+1;i<c;i++) (*this)(i,j) *= dum;
+ indx[j] = imax;
+ if((*this)(j, j) == 0.0) (*this)(j, j) = 1.e-20;
+ if(j != c){
+ dum = 1.0 / ((*this)(j, j));
+ for(i = j + 1; i < c; i++) (*this)(i, j) *= dum;
}
}
delete [] vv;
return true;
}
-
+ public:
inline int size1() const { return r; }
inline int size2() const { return c; }
SCALAR *data;
@@ -157,7 +151,8 @@ class Gmsh_Matrix
Gmsh_Matrix & operator=(const Gmsh_Matrix<SCALAR> &other)
{
if (this != &other){
- r = other.r; c=other.c;
+ r = other.r;
+ c = other.c;
data = new SCALAR[r * c];
memcpy(other);
}
@@ -192,11 +187,11 @@ class Gmsh_Matrix
inline void blas_dgemm(const Gmsh_Matrix<SCALAR> & x, const Gmsh_Matrix<SCALAR> & b,
const SCALAR c_a = 1.0, const SCALAR c_b = 1.0)
{
- Gmsh_Matrix<SCALAR> temp (x.size1(),b.size2());
- temp.mult(x,b);
+ Gmsh_Matrix<SCALAR> temp (x.size1(), b.size2());
+ temp.mult(x, b);
scale(c_b);
temp.scale(c_a);
- add (temp);
+ add(temp);
}
inline void set_all(const SCALAR &m)
{
@@ -204,13 +199,12 @@ class Gmsh_Matrix
}
inline void lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
{
+ int *indx = new int [c];
SCALAR d;
- int i,*indx;
- indx = new int [c];
- if (!_lu_decomposition(indx,d))
+ if (!_lu_decomposition(indx, d))
Msg::Fatal("Singular matrix in Gmsh_Matrix::_lu_decomposition");
- for(i=0; i < c; i++) result[i] = rhs[i];
- _back_substitution(indx,result.data);
+ for(int i = 0; i < c; i++) result[i] = rhs[i];
+ _back_substitution(indx, result.data);
delete [] indx;
}
Gmsh_Matrix cofactor(int i, int j) const
@@ -218,28 +212,26 @@ class Gmsh_Matrix
int ni = size1();
int nj = size2();
Gmsh_Matrix<SCALAR> cof(ni - 1, nj - 1);
- for (int I=0;I<ni;I++){
- for (int J=0;J<nj;J++){
- if (J!=j && I!=i)
- cof (I<i?I:I-1,J<j?J:J-1) = (*this)(I,J);
+ for (int I = 0; I < ni; I++){
+ for (int J = 0; J < nj; J++){
+ if (J != j && I != i)
+ cof (I < i ? I : I - 1, J < j ? J : J - 1) = (*this)(I, J);
}
}
return cof;
}
inline void invert(Gmsh_Matrix& y)
{
- SCALAR d,*col;
- int i,j,*indx;
-
- col = new SCALAR [c];
- indx = new int [c];
- if (!_lu_decomposition(indx,d))
+ SCALAR *col = new SCALAR[c];
+ int *indx = new int[c];
+ SCALAR d;
+ if (!_lu_decomposition(indx, d))
Msg::Fatal("Singular matrix in Gmsh_Matrix::_lu_decomposition");
- for(j=0;j<c;j++){
- for(i=0; i < c; i++) col[i] = 0.0;
- col[j]=1.0;
- _back_substitution(indx,col);
- for(i=0;i<c;i++) y(i,j)=col[i];
+ for(int j = 0; j < c; j++){
+ for(int i = 0; i < c; i++) col[i] = 0.0;
+ col[j] = 1.0;
+ _back_substitution(indx, col);
+ for(int i = 0; i < c; i++) y(i, j) = col[i];
}
delete [] col;
delete [] indx;
@@ -248,10 +240,10 @@ class Gmsh_Matrix
{
Gmsh_Matrix copy = *this;
SCALAR factor = 1.0;
- int *indx = new int [c];
- if (!copy._lu_decomposition(indx, factor))return 0.0;
+ int *indx = new int[c];
+ if (!copy._lu_decomposition(indx, factor)) return 0.0;
SCALAR det = factor;
- for (int i=0;i<c;i++)det *= copy(i,i);
+ for (int i = 0; i < c; i++) det *= copy(i, i);
delete [] indx;
return det;
}
diff --git a/Numeric/Makefile b/Numeric/Makefile
index a864a72669..31bdf450f7 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -72,9 +72,9 @@ gmshTermOfFormulation${OBJEXT}: gmshTermOfFormulation.cpp ../Common/Gmsh.h \
GmshMatrix.h gmshTermOfFormulation.h gmshFunction.h gmshLinearSystem.h \
gmshAssembler.h
gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
- ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Common/Gmsh.h \
- ../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h \
- ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Numeric/GmshMatrix.h ../Common/GmshMessage.h gmshFunction.h \
+ ../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
+ ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
diff --git a/doc/gmsh.html b/doc/gmsh.html
index 84e3dbb306..1a75b5fd36 100644
--- a/doc/gmsh.html
+++ b/doc/gmsh.html
@@ -195,7 +195,9 @@ thumbnail"></a>
<a href="/gmsh/gallery/piston.png">piston</a>,
<a href="/gmsh/gallery/pump.png">pump</a> (EDF R&D).
<li>Native models:
- <a href="/gmsh/gallery/MURUROA.png">ocean</a>,
+ <a href="/gmsh/gallery/MURUROA.png">ocean</a> (here is
+ a <a href="http://braque.mema.ucl.ac.be/~jonathan/gmsh_ocean">screencast</a>
+ explaining how to create such a mesh),
<a href="/gmsh/gallery/magnetron3.gif">magnetron 1</a>,
<a href="/gmsh/gallery/magnetron4.gif">magnetron 2</a> (P. Lefèvre),
<a href="/gmsh/gallery/breaker.gif">circuit breaker</a> (S. K. Choi),
--
GitLab