diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 0c599353e20186a97932da349bf295c676311ec4..f80268ca45fb35ff48eeea97070d7a0bf0462955 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.38 2008-02-23 15:30:07 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.39 2008-03-01 01:32:02 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -297,7 +297,6 @@ int GModel::readMSH(const std::string &name)
if(!fgets(str, sizeof(str), fp)) return 0;
int format, size;
if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) return 0;
-
if(format){
binary = true;
Msg(INFO, "Mesh is in binary format");
@@ -340,10 +339,8 @@ int GModel::readMSH(const std::string &name)
int numVertices;
if(sscanf(str, "%d", &numVertices) != 1) return 0;
Msg(INFO, "%d vertices", numVertices);
-
vertexVector.clear();
vertexMap.clear();
-
int progress = (numVertices > 100000) ? numVertices / 25 : 0;
int minVertex = numVertices + 1, maxVertex = -1;
for(int i = 0; i < numVertices; i++) {
@@ -368,7 +365,6 @@ int GModel::readMSH(const std::string &name)
Msg(PROGRESS, "Read %d vertices", i + 1);
}
if(progress) Msg(PROGRESS, "");
-
// If the vertex numbering is dense, tranfer the map into a
// vector to speed up element creation
if((int)vertexMap.size() == numVertices &&
@@ -393,7 +389,6 @@ int GModel::readMSH(const std::string &name)
int numElements;
sscanf(str, "%d", &numElements);
Msg(INFO, "%d elements", numElements);
-
int progress = (numElements > 100000) ? numElements / 25 : 0;
if(!binary){
for(int i = 0; i < numElements; i++) {
@@ -470,13 +465,21 @@ int GModel::readMSH(const std::string &name)
}
else if(!strncmp(&str[1], "NodeData", 8)) {
- // there's some post-processing data to read later on, so cache
- // the vertex indexing data
+
+ // there's some nodal post-processing data to read later on, so
+ // cache the vertex indexing data
if(vertexVector.size())
_vertexVectorCache = vertexVector;
else
_vertexMapCache = vertexMap;
postpro = true;
+
+ }
+ else if(!strncmp(&str[1], "ElementData", 11)) {
+
+ // there's some element post-processing data to read later on
+ postpro = true;
+
}
do {
diff --git a/Graphics/Post.cpp b/Graphics/Post.cpp
index d16e3b8144f2ffcd719d3b0e58b887110e953daa..0b69dc3cb2a9ebeffde99b75adfa2b262d45d22a 100644
--- a/Graphics/Post.cpp
+++ b/Graphics/Post.cpp
@@ -1,4 +1,4 @@
-// $Id: Post.cpp,v 1.153 2008-02-24 16:18:19 geuzaine Exp $
+// $Id: Post.cpp,v 1.154 2008-03-01 01:32:02 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -31,7 +31,6 @@
#include "Context.h"
#include "gl2ps.h"
-
#if defined(HAVE_MATH_EVAL)
#include "matheval.h"
#endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 1d92fdffdbc5d65ac557fe36b81afd93ebad7c9c..a2f6397488e1055fa1f787d40c71838986a2d074 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.102 2008-02-21 12:11:12 geuzaine Exp $
+// $Id: BDS.cpp,v 1.103 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -240,7 +240,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
if(!p1 || !p2) throw;
- Msg(DEBUG, " edge %d %d has to be recovered", num1, num2);
+ Msg(DEBUG, "edge %d %d has to be recovered", num1, num2);
int ix = 0;
int ixMax = 300;
@@ -260,7 +260,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0));
std::set<EdgeToRecover>::iterator itr2 =
e2r->find(EdgeToRecover(num1, num2, 0));
- Msg(DEBUG2, " edge %d %d on model edge %d cannot be recovered because"
+ Msg(DEBUG2, "edge %d %d on model edge %d cannot be recovered because"
" it intersects %d %d on model edge %d", num1, num2, itr2->ge->tag(),
e->p1->iD, e->p2->iD, itr1->ge->tag());
// now throw a class that contains the diagnostic
@@ -274,7 +274,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
}
// if(ix > 300){
-// Msg(WARNING," edge %d %d cannot be recovered after %d iterations, trying again",
+// Msg(WARNING, "edge %d %d cannot be recovered after %d iterations, trying again",
// num1, num2, ix);
// ix = 0;
// }
@@ -285,7 +285,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
if(!eee){
outputScalarField(triangles, "debugp.pos", 1);
outputScalarField(triangles, "debugr.pos", 0);
- Msg(GERROR," edge %d %d cannot be recovered at all, look at debugp.pos "
+ Msg(GERROR, "edge %d %d cannot be recovered at all, look at debugp.pos "
"and debugr.pos", num1, num2);
return 0;
}
@@ -1434,7 +1434,7 @@ recombine_T::recombine_T(const BDS_Edge *_e)
void BDS_Mesh::recombineIntoQuads(const double angle_limit, GFace *gf)
{
- Msg(INFO,"Recombining triangles for surface %d", gf->tag());
+ 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();
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index 66da06851507ad6daee5bbded0f0badfda9b4236..8956f8749d078183f57ef183c4d1e6188102ef70 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -1,4 +1,4 @@
-// $Id: DivideAndConquer.cpp,v 1.15 2008-02-21 19:27:58 geuzaine Exp $
+// $Id: DivideAndConquer.cpp,v 1.16 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -19,21 +19,17 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
-/*
- A L G O R I T H M E D I V I D E A N D C O N Q U E R
-
- le noeud de cette methode est de pouvoir fusionner deux
- triangulations de Delaunay en une seule (routine merge) on procede
- alors recursivement en separant les points en deux groupes puis en
- separant les groupes en 2 ... jusqu'a n'obtenir que 1 2 ou 3 points
- (la triangulation est alors triviale)
-
- Dans le mailleur, on utilise cet algorithme pour construire le
- maillage initial
-
- !!! il faut PERTURBER les points d'une faible valeur aleatoire pour
- eviter d'avoir 3 points alignes ou 4 points cocycliques !!!
-*/
+// Triangulation using a divide and conquer algorithm
+//
+// The main idea is to be able to merge two Delaunay triangulations
+// into a single one (see the 'merge' function). Points are then
+// separated into two groups, then each group into two, ... until
+// having 1, 2 or 3 points (the triangulation is then trivial).
+//
+// This is used to contruct the initial mesh.
+//
+// Warning: point positions must be PERTURBED by a small random
+// value to avoid 3 aligned points or 4 cocyclical points!
#include "Message.h"
#include "Numeric.h"
@@ -102,7 +98,7 @@ PointNumero First(PointNumero x)
return (pPointArray[x].adjacent)->point_num;
}
-int Is_left_of(PointNumero x, PointNumero y, PointNumero check)
+int IsLeftOf(PointNumero x, PointNumero y, PointNumero check)
{
double pa[2] = {(double)pPointArray[x].where.h, (double)pPointArray[x].where.v};
double pb[2] = {(double)pPointArray[y].where.h, (double)pPointArray[y].where.v};
@@ -114,9 +110,9 @@ int Is_left_of(PointNumero x, PointNumero y, PointNumero check)
return result > 0;
}
-int Is_right_of(PointNumero x, PointNumero y, PointNumero check)
+int IsRightOf(PointNumero x, PointNumero y, PointNumero check)
{
- return Is_left_of(y, x, check);
+ return IsLeftOf(y, x, check);
}
Segment LowerCommonTangent(DT vl, DT vr)
@@ -130,12 +126,12 @@ Segment LowerCommonTangent(DT vl, DT vr)
z1 = First(x);
z2 = Predecessor(x, z1);
for(;;) {
- if(Is_right_of(x, y, z)) {
+ if(IsRightOf(x, y, z)) {
temp = z;
z = Successor(z, y);
y = temp;
}
- else if(Is_right_of(x, y, z2)) {
+ else if(IsRightOf(x, y, z2)) {
temp = z2;
z2 = Predecessor(z2, x);
x = temp;
@@ -159,12 +155,12 @@ Segment UpperCommonTangent(DT vl, DT vr)
z1 = First(x);
z2 = Predecessor(y, z);
for(;;) {
- if(Is_left_of(x, y, z2)) {
+ if(IsLeftOf(x, y, z2)) {
temp = z2;
z2 = Predecessor(z2, y);
y = temp;
}
- else if(Is_left_of(x, y, z1)) {
+ else if(IsLeftOf(x, y, z1)) {
temp = z1;
z1 = Successor(z1, x);
x = temp;
@@ -215,7 +211,7 @@ int merge(DT vl, DT vr)
r1 = Predecessor(r, l);
if(r1 == -1)
return 0;
- if(Is_right_of(l, r, r1))
+ if(IsRightOf(l, r, r1))
a = 1;
else {
out = 0;
@@ -231,7 +227,7 @@ int merge(DT vl, DT vr)
if(!Delete(r, r1))
return 0;
r1 = r2;
- if(Is_right_of(l, r, r1))
+ if(IsRightOf(l, r, r1))
out = a = 1;
}
}
@@ -240,7 +236,7 @@ int merge(DT vl, DT vr)
l1 = Successor(l, r);
if(l1 == -1)
return 0;
- if(Is_left_of(r, l, l1))
+ if(IsLeftOf(r, l, l1))
b = 1;
else {
out = 0;
@@ -256,7 +252,7 @@ int merge(DT vl, DT vr)
if(!Delete(l, l1))
return 0;
l1 = l2;
- if(Is_left_of(r, l, l1))
+ if(IsLeftOf(r, l, l1))
out = b = 1;
}
}
@@ -307,7 +303,7 @@ DT recur_trig(PointNumero left, PointNumero right)
Insert(left, right);
Insert(left, left + 1);
Insert(left + 1, right);
- if(Is_right_of(left, right, left + 1)) {
+ if(IsRightOf(left, right, left + 1)) {
FixFirst(left, left + 1);
FixFirst(left + 1, right);
FixFirst(right, left);
@@ -340,18 +336,12 @@ int comparePoints(const void *i, const void *j)
return (x < 0.) ? -1 : 1;
}
-// this fonction builds the delaunay triangulation and the voronoi for
-// a window. All error handling is done here
-int DelaunayAndVoronoi(DocPeek doc)
+// this fonction builds the delaunay triangulation for a window
+int BuildDelaunay(DocRecord *doc)
{
pPointArray = doc->points;
-
- if(doc->numPoints < 2)
- return 1;
-
qsort(doc->points, doc->numPoints, sizeof(PointRecord), comparePoints);
recur_trig(0, doc->numPoints - 1);
-
return 1;
}
@@ -515,19 +505,8 @@ PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n)
return ptr;
}
-void filldel(Delaunay *deladd, int aa, int bb, int cc,
- PointRecord *points)
-{
- deladd->t.a = aa;
- deladd->t.b = bb;
- deladd->t.c = cc;
- deladd->v.voisin1 = NULL;
- deladd->v.voisin2 = NULL;
- deladd->v.voisin3 = NULL;
-}
-
// Convertir les listes d'adjacence en triangles
-int Conversion(DocPeek doc)
+int ConvertDListToTriangles(DocRecord *doc)
{
// on suppose que n >= 3. gPointArray est suppose OK.
@@ -544,10 +523,7 @@ int Conversion(DocPeek doc)
// nombre de triangles que l'on doit obtenir
count2 = 2 * (n - 1) - count2;
- if(doc->delaunay)
- Free(doc->delaunay);
-
- doc->delaunay = (Delaunay *) Malloc(2 * count2 * sizeof(Delaunay));
+ doc->triangles = (Triangle *) Malloc(2 * count2 * sizeof(Triangle));
for(i = 0; i < n; i++) {
// on cree une liste de points connectes au point i (t) + nombre
@@ -563,11 +539,13 @@ int Conversion(DocPeek doc)
for(i = 0; i < n; i++) {
for(j = 0; j < striangle[i].t_length; j++) {
if((striangle[i].t[j] > i) && (striangle[i].t[j + 1] > i) &&
- (Is_right_of(i, striangle[i].t[j], striangle[i].t[j + 1]))) {
+ (IsRightOf(i, striangle[i].t[j], striangle[i].t[j + 1]))) {
aa = i;
bb = striangle[i].t[j];
cc = striangle[i].t[j + 1];
- filldel(&doc->delaunay[count], aa, bb, cc, gPointArray);
+ doc->triangles[count].a = aa;
+ doc->triangles[count].b = bb;
+ doc->triangles[count].c = cc;
count++;
}
}
@@ -581,7 +559,7 @@ int Conversion(DocPeek doc)
}
// Cette routine efface toutes les listes d'adjacence du pPointArray
-void remove_all_dlist(int n, PointRecord *pPointArray)
+void RemoveAllDList(int n, PointRecord *pPointArray)
{
int i;
DListPeek p, temp;
@@ -598,13 +576,23 @@ void remove_all_dlist(int n, PointRecord *pPointArray)
}
}
-void Make_Mesh_With_Points(DocRecord *ptr, PointRecord *Liste, int Numpoints)
+DocRecord::DocRecord(int n)
+ : numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
+{
+ if(numPoints)
+ points = (PointRecord*)Malloc(numPoints * sizeof(PointRecord));
+}
+
+DocRecord::~DocRecord()
+{
+ Free(points);
+ Free(triangles);
+}
+
+void DocRecord::MakeMeshWithPoints()
{
- ptr->numTriangles = 0;
- ptr->points = Liste;
- ptr->numPoints = Numpoints;
- ptr->delaunay = 0;
- DelaunayAndVoronoi(ptr);
- Conversion(ptr);
- remove_all_dlist(ptr->numPoints, ptr->points);
+ if(numPoints < 3) return;
+ BuildDelaunay(this);
+ ConvertDListToTriangles(this);
+ RemoveAllDList(numPoints, points);
}
diff --git a/Mesh/DivideAndConquer.h b/Mesh/DivideAndConquer.h
index 44511c488b255204a465d6dd8003c2ba72642a48..28c2ad380049ec4c20c416c6bf4bbccd6013ba84 100644
--- a/Mesh/DivideAndConquer.h
+++ b/Mesh/DivideAndConquer.h
@@ -20,29 +20,19 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
-typedef struct _POINT PointRecord, *PointPeek;
-typedef struct _DOC DocRecord, *DocPeek;
typedef struct _CDLIST DListRecord, *DListPeek;
-typedef struct _DELAUNAY Delaunay, *delpeek;
typedef int PointNumero;
-struct _DOC{
- PointRecord *points; // points to triangulate
- int numPoints; // number of points
- int numTriangles; // number of triangles
- Delaunay *delaunay; // 2D results
-};
-
typedef struct{
double v;
double h;
}MPoint;
-struct _POINT{
+typedef struct{
MPoint where;
DListPeek adjacent;
void *data;
-};
+}PointRecord;
struct _CDLIST{
PointNumero point_num;
@@ -74,15 +64,15 @@ typedef struct{
PointNumero a, b, c;
}Triangle;
-typedef struct {
- Delaunay *voisin1, *voisin2, *voisin3;
-}Voronoi;
-
-struct _DELAUNAY{
- Triangle t;
- Voronoi v;
+class DocRecord{
+ public:
+ int numPoints; // number of points
+ PointRecord *points; // points to triangulate
+ int numTriangles; // number of triangles
+ Triangle *triangles; // 2D results
+ DocRecord(int n);
+ ~DocRecord();
+ void MakeMeshWithPoints();
};
-void Make_Mesh_With_Points(DocRecord *ptr, PointRecord *Liste, int Numpoints);
-
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 888a297aadf05102a23b17d45ebb97c98defc3d3..562e77dfbc89db2a2eb92ce6f75db821d7e7bace 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.23 2008-02-27 12:39:08 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.24 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -98,7 +98,7 @@ bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double
// N is the total number of points (3 for quadratic, 4 for cubic ...)
// u[0] = u0;
// u[N-1] = uN;
- // dist ( ge(u[i]), ge(u[i+1]) ) = dist ( ge(u[i+1]), ge(u[i+2]) ) , i = 0,...,N-2
+ // dist(ge(u[i]), ge(u[i+1])) = dist(ge(u[i+1]), ge(u[i+2])), i = 0,...,N-2
// initialize as equidistant in parameter space
u[0] = u0;
@@ -891,7 +891,7 @@ void deriv_smoothing_objective_function_HighOrderN(double *uv, double *dF,
}
}
-void optimizeNodeLocations ( GFace *gf, smoothVertexDataHON &vdN , double eps = .2)
+void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
{
if(!vdN.v.size()) return;
double uv[20];
@@ -958,8 +958,8 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
GEntity *ge = ver->onWhat();
if (ge->dim() == 2){
double initu,initv;
- ver->getParameter ( 0,initu);
- ver->getParameter ( 1,initv);
+ ver->getParameter(0, initu);
+ ver->getParameter(1, initv);
smoothVertexDataHON vdN;
vdN.ts = it->second;
@@ -973,7 +973,7 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
double val;
double F = -smooth_obj_HighOrder(initu,initv, &vd);
if (F < .2){
- minimize_2 ( smooth_obj_HighOrder,
+ minimize_2(smooth_obj_HighOrder,
deriv_smoothing_objective_function_HighOrder, &vd, 1, initu,initv,val);
double Fafter = -smooth_obj_HighOrder(initu,initv, &vd);
if (F < Fafter){
@@ -1098,28 +1098,28 @@ bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
MVertex *vert4 = (n4 < n1) ? e4[e4.size() - i - 1] : e4[i];
MVertex *vert2 = (n2 < n3) ? e2[i] : e2[e2.size() - i - 1];
double U1,V1,U2,V2,U3,V3,U4,V4,U,V,nU1,nV1,nU2,nV2,nU3,nV3,nU4,nV4;
- parametricCoordinates(vert , gf,U,V);
- parametricCoordinates(vert1, gf,U1,V1);
- parametricCoordinates(vert2, gf,U2,V2);
- parametricCoordinates(vert3, gf,U3,V3);
- parametricCoordinates(vert4, gf,U4,V4);
- parametricCoordinates(n1, gf,nU1,nV1);
- parametricCoordinates(n2, gf,nU2,nV2);
- parametricCoordinates(n3, gf,nU3,nV3);
- parametricCoordinates(n4, gf,nU4,nV4);
+ parametricCoordinates(vert , gf, U, V);
+ parametricCoordinates(vert1, gf, U1, V1);
+ parametricCoordinates(vert2, gf, U2, V2);
+ parametricCoordinates(vert3, gf, U3, V3);
+ parametricCoordinates(vert4, gf, U4, V4);
+ parametricCoordinates(n1, gf, nU1, nV1);
+ parametricCoordinates(n2, gf, nU2, nV2);
+ parametricCoordinates(n3, gf, nU3, nV3);
+ parametricCoordinates(n4, gf, nU4, nV4);
- Unew[k][i] = U + relax * ( (1.-u) * U4 + u * U2 +
- (1.-v) * U1 + v * U3 -
- ( (1.-u)*(1.-v) * nU1
- + u * (1.-v) * nU2
- + u*v*nU3
- + (1.-u) * v * nU4) - U);
- Vnew[k][i] = V + relax * ( (1.-u) * V4 + u * V2 +
- (1.-v) * V1 + v * V3 -
- ( (1.-u)*(1.-v) * nV1
- + u * (1.-v) * nV2
- + u*v*nV3
- + (1.-u) * v * nV4) - V);
+ Unew[k][i] = U + relax * ((1.-u) * U4 + u * U2 +
+ (1.-v) * U1 + v * U3 -
+ ((1.-u)*(1.-v) * nU1
+ + u * (1.-v) * nU2
+ + u * v * nU3
+ + (1.-u) * v * nU4) - U);
+ Vnew[k][i] = V + relax * ((1.-u) * V4 + u * V2 +
+ (1.-v) * V1 + v * V3 -
+ ((1.-u)*(1.-v) * nV1
+ + u * (1.-v) * nV2
+ + u * v * nV3
+ + (1.-u) * v * nV4) - V);
GPoint gp = gf->point(Unew[k][i],Vnew[k][i]);
vert->x() = gp.x();
vert->y() = gp.y();
@@ -1127,8 +1127,8 @@ bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
}
double minJloc = 1.e22;
double maxJloc = -1.e22;
- getMinMaxJac (t1, minJloc, maxJloc);
- getMinMaxJac (t2, minJloc, maxJloc);
+ getMinMaxJac(t1, minJloc, maxJloc);
+ getMinMaxJac(t2, minJloc, maxJloc);
if (minJloc > minJ){
kopt = k;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 709cfb91d4e9a46dc7d5be71034336a2d997c32b..16438dfcfec42f4d63924f6822bbbf262e8bec30 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.122 2008-02-23 16:19:22 remacle Exp $
+// $Id: meshGFace.cpp,v 1.123 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -259,7 +259,7 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
BDS_Edge * e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
if (e) e->g = g;
else {
- // Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag());
+ // Msg(GERROR, "The unrecoverable edge is on model edge %d",ge->tag());
return false;
}
return true;
@@ -281,7 +281,7 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
if (e) e->g = g;
else {
- // Msg(GERROR,"The unrecoverable edge is on model edge %d",ge->tag());
+ // Msg(GERROR, "The unrecoverable edge is on model edge %d",ge->tag());
// return false;
}
}
@@ -294,7 +294,7 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
if (e) e->g = g;
else {
- // Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",
+ // Msg(GERROR, "Unable to recover an edge %g %g && %g %g (%d/%d)",
// vstart->x(), vstart->y(), vend->x(), vend->y(), i,
// ge->mesh_vertices.size());
// return false;
@@ -308,7 +308,7 @@ bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
e = m->recover_edge(vstart->getNum(), vend->getNum(), e2r, not_recovered);
if (e)e->g = g;
else {
- // Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/%d)",
+ // Msg(GERROR, "Unable to recover an edge %g %g && %g %g (%d/%d)",
// vstart->x(), vstart->y(), vend->x(), vend->y(),
// ge->mesh_vertices.size(), ge->mesh_vertices.size());
// return false;
@@ -361,6 +361,11 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
++it;
}
+ if (all_vertices.size() < 3){
+ Msg(WARNING, "Cannot triangulate less than 3 vertices");
+ return false;
+ }
+
double *U_ = new double[all_vertices.size()];
double *V_ = new double[all_vertices.size()];
@@ -407,112 +412,97 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
SVector3 dd (bbox.max(), bbox.min());
double LC2D = norm(dd);
- // Use a divide & conquer type algorithm to create a triangulation
- // We add to the triangulation a box with 4 points that encoses
- // the domain.
-
- /// Fill the DocRecord Data Structure with the points
- DocRecord doc;
-
- doc.points = (PointRecord*)malloc((all_vertices.size() + 4) * sizeof(PointRecord));
- itv = all_vertices.begin();
- int j = 0;
- while(itv != all_vertices.end()){
- MVertex *here = *itv;
- double XX = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
- double YY = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
- doc.points[j].where.h = U_[j] + XX;
- doc.points[j].where.v = V_[j] + YY;
- doc.points[j].adjacent = NULL;
- doc.points[j].data = here;
- j++;
- ++itv;
- }
-
- // Increase the size of the bounding box
- bbox *= 2.5;
- // add 4 points than encloses the domain (use negative number to
- // distinguish thos fake vertices)
- MVertex *bb[4];
- bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, gf, -1);
- bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, gf, -2);
- bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, gf, -3);
- bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, gf, -4);
-
- for(int ip = 0; ip < 4; ip++){
- doc.points[all_vertices.size() + ip].where.h = bb[ip]->x();
- doc.points[all_vertices.size() + ip].where.v = bb[ip]->y();
- doc.points[all_vertices.size() + ip].adjacent = 0;
- doc.points[all_vertices.size() + ip].data = bb[ip];
- }
-
- // if the number of vertices is less or equal to 2, then no elements are generated
- if (all_vertices.size() <= 2){
- free(doc.points);
- free(doc.delaunay);
- for(int ip = 0; ip < 4; ip++) delete bb[ip];
- }
-
- // Use "fast" inhouse recursive algo to generate the triangulation
- // At this stage the triangulation is not what we need
- // -) It does not necessary recover the boundaries
- // -) It contains triangles outside the domain (the first edge
- // loop is the outer one)
- Msg(DEBUG1,"Meshing of the convex hull (%d points)", all_vertices.size());
- Make_Mesh_With_Points(&doc, doc.points, all_vertices.size() + 4);
-
- // Buid a BDS_Mesh structure that is convenient for doing the actual meshing procedure
+ // Buid a BDS_Mesh structure that is convenient for doing the actual
+ // meshing procedure
BDS_Mesh *m = new BDS_Mesh;
m->scalingU = 1;
m->scalingV = 1;
- for(int i = 0; i < doc.numPoints; i++){
- MVertex *here = (MVertex *)doc.points[i].data;
- int num = here->getNum();
- double U, V;
- // This test was missing in 2.0.0 and led to the seemingly random
- // Windows/Mac slowdowns (we were passing random numbers to curve
- // interpolation, and straight line interpol does "while(not in
- // bounds) i--" which would take forever to get back into a
- // reasonnable interval).
- // JFR : the fix was WRONG, I fixed the fix ;-)
- if(num< 0){ // fake bbox points
- U = bb[-1-num]->x();
- V = bb[-1-num]->y();
- }
- else{
- U = U_[num];
- V = V_[num];
+ // Use a divide & conquer type algorithm to create a triangulation.
+ // We add to the triangulation a box with 4 points that encloses the
+ // domain.
+ {
+ DocRecord doc(all_vertices.size() + 4);
+ itv = all_vertices.begin();
+ int j = 0;
+ while(itv != all_vertices.end()){
+ MVertex *here = *itv;
+ double XX = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
+ double YY = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
+ doc.points[j].where.h = U_[j] + XX;
+ doc.points[j].where.v = V_[j] + YY;
+ doc.points[j].adjacent = NULL;
+ doc.points[j].data = here;
+ j++;
+ ++itv;
}
- BDS_Point *pp = m->add_point(num, U, V, gf);
+ // Increase the size of the bounding box
+ bbox *= 2.5;
+ // add 4 points than encloses the domain (use negative number to
+ // distinguish thos fake vertices)
+ MVertex *bb[4];
+ bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, gf, -1);
+ bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, gf, -2);
+ bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, gf, -3);
+ bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, gf, -4);
- GEntity *ge = here->onWhat();
- if(ge->dim() == 0){
- pp->lcBGM() = BGM_MeshSize(ge,0,0,here->x(),here->y(),here->z());
+ for(int ip = 0; ip < 4; ip++){
+ doc.points[all_vertices.size() + ip].where.h = bb[ip]->x();
+ doc.points[all_vertices.size() + ip].where.v = bb[ip]->y();
+ doc.points[all_vertices.size() + ip].adjacent = 0;
+ doc.points[all_vertices.size() + ip].data = bb[ip];
}
- else if(ge->dim() == 1){
- double u;
- here->getParameter(0,u);
- pp->lcBGM() = BGM_MeshSize(ge,u,0,here->x(),here->y(),here->z());
+
+ // Use "fast" inhouse recursive algo to generate the triangulation
+ // At this stage the triangulation is not what we need
+ // -) It does not necessary recover the boundaries
+ // -) It contains triangles outside the domain (the first edge
+ // loop is the outer one)
+ Msg(DEBUG1, "Meshing of the convex hull (%d points)", all_vertices.size());
+ doc.MakeMeshWithPoints();
+
+ for(int i = 0; i < doc.numPoints; i++){
+ MVertex *here = (MVertex *)doc.points[i].data;
+ int num = here->getNum();
+ double U, V;
+ if(num < 0){ // fake bbox points
+ U = bb[-1 - num]->x();
+ V = bb[-1 - num]->y();
+ }
+ else{
+ U = U_[num];
+ V = V_[num];
+ }
+
+ BDS_Point *pp = m->add_point(num, U, V, gf);
+
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 0){
+ pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
+ }
+ else if(ge->dim() == 1){
+ double u;
+ here->getParameter(0,u);
+ pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
+ }
+ else
+ pp->lcBGM() = 1.e22;
+
+ pp->lc() = pp->lcBGM();
}
- else
- pp->lcBGM() = 1.e22;
- pp->lc() = pp->lcBGM();
- }
-
- Msg(DEBUG1, "Meshing of the convex hull (%d points) done", all_vertices.size());
+ Msg(DEBUG1, "Meshing of the convex hull (%d points) done", all_vertices.size());
+
+ for(int i = 0; i < doc.numTriangles; i++) {
+ 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());
+ }
- for(int i = 0; i < doc.numTriangles; i++) {
- MVertex *V1 = (MVertex*)doc.points[doc.delaunay[i].t.a].data;
- MVertex *V2 = (MVertex*)doc.points[doc.delaunay[i].t.b].data;
- MVertex *V3 = (MVertex*)doc.points[doc.delaunay[i].t.c].data;
- m->add_triangle(V1->getNum(), V2->getNum(), V3->getNum());
+ for(int ip = 0; ip < 4; ip++) delete bb[ip];
}
- free(doc.points);
- free(doc.delaunay);
- for(int ip = 0; ip < 4; ip++) delete bb[ip];
// Recover the boundary edges and compute characteristic lenghts
// using mesh edge spacing
@@ -544,7 +534,7 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
++it;
}
- Msg(DEBUG1,"Recovering %d mesh Edges", edgesToRecover.size());
+ Msg(DEBUG1, "Recovering %d mesh Edges", edgesToRecover.size());
// effectively recover the medge
it = edges.begin();
@@ -566,8 +556,9 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
int p1 = itr->p1;
int p2 = itr->p2;
int tag = itr->ge->tag();
- Msg(WARNING,"MEdge %d %d in GEdge %d",p1,p2,tag);
+ Msg(WARNING, "MEdge %d %d in GEdge %d",p1,p2,tag);
}
+ // delete the mesh
delete m;
delete [] U_;
delete [] V_;
@@ -576,10 +567,10 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
return false;
}
if(RECUR_ITER > 0)
- Msg(WARNING,":-) Gmsh was able to recover all edges after %d ITERATIONS",
+ Msg(WARNING, ":-) Gmsh was able to recover all edges after %d ITERATIONS",
RECUR_ITER);
- // Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
+ // Msg(INFO, "Boundary Edges recovered for surface %d",gf->tag());
// Look for an edge that is on the boundary for which one of the two
// neighbors has a negative number node. The other triangle is
// inside the domain and, because all edges were recovered,
@@ -644,9 +635,9 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
if (debug){
char name[245];
- sprintf(name,"surface%d-recovered-real.pos", gf->tag());
+ sprintf(name, "surface%d-recovered-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0);
- sprintf(name,"surface%d-recovered-param.pos", gf->tag());
+ sprintf(name, "surface%d-recovered-param.pos", gf->tag());
outputScalarField(m->triangles, name, 1);
}
@@ -716,15 +707,14 @@ bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
}
else if (debug){
char name[256];
- sprintf(name,"real%d.pos",gf->tag());
- outputScalarField(m->triangles, name,0,gf);
- sprintf(name,"param%d.pos",gf->tag());
+ sprintf(name, "real%d.pos", gf->tag());
+ outputScalarField(m->triangles, name, 0, gf);
+ sprintf(name, "param%d.pos", gf->tag());
outputScalarField(m->triangles, name,1);
}
// delete the mesh
delete m;
-
delete [] U_;
delete [] V_;
@@ -843,7 +833,7 @@ bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
if (seam && seam_the_first){
coords = ((*it)._sign == 1) ? mesh1d_seam : mesh1d_seam_reversed;
found = (*it);
- Msg(INFO,"This test case would have failed in Previous Gmsh Version ;-)");
+ Msg(INFO, "This test case would have failed in Previous Gmsh Version ;-)");
}
else{
coords = ((*it)._sign == 1) ? mesh1d : mesh1d_reversed;
@@ -1030,64 +1020,64 @@ bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- // build a point record structure to create the initial mesh
- DocRecord doc;
- doc.points = (PointRecord*)malloc((nbPointsTotal + 4) * sizeof(PointRecord));
- int count = 0;
- for (unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
- std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
- for (unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
- BDS_Point *pp = edgeLoop_BDS[j];
- const double U = pp->u;
- const double V = pp->v;
- double XX = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
- double YY = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
- doc.points[count].where.h = U + XX;
- doc.points[count].where.v = V + YY;
- doc.points[count].adjacent = NULL;
- doc.points[count].data = pp;
- count++;
+ // Use a divide & conquer type algorithm to create a triangulation.
+ // We add to the triangulation a box with 4 points that encloses the
+ // domain.
+ {
+ DocRecord doc(nbPointsTotal + 4);
+ int count = 0;
+ for (unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
+ std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
+ for (unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
+ BDS_Point *pp = edgeLoop_BDS[j];
+ const double U = pp->u;
+ const double V = pp->v;
+ double XX = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
+ double YY = CTX.mesh.rand_factor * LC2D * (double)rand() / (double)RAND_MAX;
+ doc.points[count].where.h = U + XX;
+ doc.points[count].where.v = V + YY;
+ doc.points[count].adjacent = NULL;
+ doc.points[count].data = pp;
+ count++;
+ }
}
- }
- // Increase the size of the bounding box, add 4 points that encloses
- // the domain, use negative number to distinguish those fake
- // vertices
- bbox *= 3.5;
- MVertex *bb[4];
- bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, 0, -1);
- bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, 0, -2);
- bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, 0, -3);
- bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, 0, -4);
-
- for ( int ip = 0; ip < 4; ip++){
- BDS_Point *pp = m->add_point(-ip - 1, bb[ip]->x(), bb[ip]->y(), gf);
- m->add_geom(gf->tag(), 2);
- BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
- pp->g = g;
- doc.points[nbPointsTotal+ip].where.h = bb[ip]->x();
- doc.points[nbPointsTotal+ip].where.v = bb[ip]->y();
- doc.points[nbPointsTotal+ip].adjacent = 0;
- doc.points[nbPointsTotal+ip].data = pp;
- }
- for (int ip = 0; ip < 4; ip++) delete bb[ip];
-
- // Use "fast" inhouse recursive algo to generate the triangulation
- // At this stage the triangulation is not what we need
- // -) It does not necessary recover the boundaries
- // -) It contains triangles outside the domain (the first edge
- // loop is the outer one)
- Msg(DEBUG1, "Meshing of the convex hull (%d points)", nbPointsTotal);
- Make_Mesh_With_Points(&doc, doc.points,nbPointsTotal + 4);
- for(int i = 0; i < doc.numTriangles; i++){
- BDS_Point *p1 = (BDS_Point*)doc.points[doc.delaunay[i].t.a].data;
- BDS_Point *p2 = (BDS_Point*)doc.points[doc.delaunay[i].t.b].data;
- BDS_Point *p3 = (BDS_Point*)doc.points[doc.delaunay[i].t.c].data;
- m->add_triangle(p1->iD, p2->iD, p3->iD);
+ // Increase the size of the bounding box, add 4 points that enclose
+ // the domain, use negative number to distinguish those fake
+ // vertices
+ bbox *= 3.5;
+ MVertex *bb[4];
+ bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, 0, -1);
+ bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, 0, -2);
+ bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, 0, -3);
+ bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, 0, -4);
+ for (int ip = 0; ip < 4; ip++){
+ BDS_Point *pp = m->add_point(-ip - 1, bb[ip]->x(), bb[ip]->y(), gf);
+ m->add_geom(gf->tag(), 2);
+ BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
+ pp->g = g;
+ doc.points[nbPointsTotal+ip].where.h = bb[ip]->x();
+ doc.points[nbPointsTotal+ip].where.v = bb[ip]->y();
+ doc.points[nbPointsTotal+ip].adjacent = 0;
+ doc.points[nbPointsTotal+ip].data = pp;
+ }
+ for (int ip = 0; ip < 4; ip++) delete bb[ip];
+
+ // Use "fast" inhouse recursive algo to generate the triangulation
+ // At this stage the triangulation is not what we need
+ // -) It does not necessary recover the boundaries
+ // -) It contains triangles outside the domain (the first edge
+ // loop is the outer one)
+ Msg(DEBUG1, "Meshing of the convex hull (%d points)", nbPointsTotal);
+ doc.MakeMeshWithPoints();
+
+ for(int i = 0; i < doc.numTriangles; i++){
+ BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data;
+ BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data;
+ BDS_Point *p3 = (BDS_Point*)doc.points[doc.triangles[i].c].data;
+ m->add_triangle(p1->iD, p2->iD, p3->iD);
+ }
}
- // Free stuff
- free(doc.points);
- free(doc.delaunay);
// Recover the boundary edges and compute characteristic lenghts
// using mesh edge spacing
@@ -1117,7 +1107,7 @@ bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- // Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
+ // Msg(INFO, "Boundary Edges recovered for surface %d",gf->tag());
// Look for an edge that is on the boundary for which one of the two
// neighbors has a negative number node. The other triangle is
// inside the domain and, because all edges were recovered,
@@ -1244,7 +1234,6 @@ bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- // delete the mesh
if (debug){
char name[245];
sprintf(name, "surface%d-final-real.pos", gf->tag());
@@ -1257,8 +1246,10 @@ bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
insertVerticesInFace(gf, m);
laplaceSmoothing(gf);
}
-
+
+ // delete the mesh
delete m;
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
gf->meshStatistics.best_element_shape,
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index a96a3c866062523303a163b0bd9cdc722fb90830..6695c2e635d22b70960326f936240a0f64a8ab94 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.8 2008-02-21 13:34:40 geuzaine Exp $
+// $Id: meshGFaceBDS.cpp,v 1.9 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -611,7 +611,7 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
t_sm += t6 - t5;
m.cleanup();
IT++;
- Msg(DEBUG1," iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
+ Msg(DEBUG1, " iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
"%6d splits, %6d swaps, %6d collapses, %6d moves",
IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps, nb_smooth);
if (nb_split == 0 && nb_collaps == 0) break;
@@ -619,16 +619,16 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6;
- Msg(DEBUG1," ---------------------------------------");
- Msg(DEBUG1," CPU Report ");
- Msg(DEBUG1," ---------------------------------------");
- Msg(DEBUG1," CPU SWAP %12.5E sec (%3.0f %%)", t_sw,100 * t_sw / t_total);
- Msg(DEBUG1," CPU SPLIT %12.5E sec (%3.0f %%) ", t_spl,100 * t_spl / t_total);
- Msg(DEBUG1," CPU COLLAPS %12.5E sec (%3.0f %%) ", t_col,100 * t_col / t_total);
- Msg(DEBUG1," CPU SMOOTH %12.5E sec (%3.0f %%) ", t_sm,100 * t_sm / t_total);
- Msg(DEBUG1," ---------------------------------------");
- Msg(DEBUG1," CPU TOTAL %12.5E sec ",t_total);
- Msg(DEBUG1," ---------------------------------------");
+ Msg(DEBUG1, " ---------------------------------------");
+ Msg(DEBUG1, " CPU Report ");
+ Msg(DEBUG1, " ---------------------------------------");
+ Msg(DEBUG1, " CPU SWAP %12.5E sec (%3.0f %%)", t_sw,100 * t_sw / t_total);
+ Msg(DEBUG1, " CPU SPLIT %12.5E sec (%3.0f %%) ", t_spl,100 * t_spl / t_total);
+ Msg(DEBUG1, " CPU COLLAPS %12.5E sec (%3.0f %%) ", t_col,100 * t_col / t_total);
+ Msg(DEBUG1, " CPU SMOOTH %12.5E sec (%3.0f %%) ", t_sm,100 * t_sm / t_total);
+ Msg(DEBUG1, " ---------------------------------------");
+ Msg(DEBUG1, " CPU TOTAL %12.5E sec ",t_total);
+ Msg(DEBUG1, " ---------------------------------------");
}
void allowAppearanceofEdge (BDS_Point *p1, BDS_Point *p2)
diff --git a/Plugin/Triangulate.cpp b/Plugin/Triangulate.cpp
index 7083b9c2d4cfe7822957092304a4affe0d940c29..9e79310cedde9df626b75860594ad49d13acaeac 100644
--- a/Plugin/Triangulate.cpp
+++ b/Plugin/Triangulate.cpp
@@ -1,4 +1,4 @@
-// $Id: Triangulate.cpp,v 1.45 2008-02-20 09:20:47 geuzaine Exp $
+// $Id: Triangulate.cpp,v 1.46 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -117,8 +117,7 @@ static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
double lc = norm(SVector3(bbox.max(), bbox.min()));
// build a point record structure for the divide and conquer algorithm
- DocRecord doc;
- doc.points = (PointRecord*)malloc(points.size() * sizeof(PointRecord));
+ DocRecord doc(points.size());
for(unsigned int i = 0; i < points.size(); i++){
double XX = CTX.mesh.rand_factor * lc * (double)rand() / (double)RAND_MAX;
double YY = CTX.mesh.rand_factor * lc * (double)rand() / (double)RAND_MAX;
@@ -130,13 +129,13 @@ static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
}
// triangulate
- Make_Mesh_With_Points(&doc, doc.points, nbIn);
+ doc.MakeMeshWithPoints();
// create output (using unperturbed data)
for(int i = 0; i < doc.numTriangles; i++){
- double *pa = (double*)doc.points[doc.delaunay[i].t.a].data;
- double *pb = (double*)doc.points[doc.delaunay[i].t.b].data;
- double *pc = (double*)doc.points[doc.delaunay[i].t.c].data;
+ double *pa = (double*)doc.points[doc.triangles[i].a].data;
+ double *pb = (double*)doc.points[doc.triangles[i].b].data;
+ double *pc = (double*)doc.points[doc.triangles[i].c].data;
for(int j = 0; j < 3; j++) {
List_Add(outList, pa + j);
List_Add(outList, pb + j);
@@ -149,9 +148,6 @@ static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
}
(*nbOut)++;
}
-
- free(doc.points);
- free(doc.delaunay);
}
PView *GMSH_TriangulatePlugin::execute(PView *v)
diff --git a/Post/PView.cpp b/Post/PView.cpp
index a83789a938c7c6d3775dedf87f784fed533a2980..2392cc74b43614d0c797081f2155dbc18ba054b6 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -1,4 +1,4 @@
-// $Id: PView.cpp,v 1.17 2008-02-20 09:24:41 geuzaine Exp $
+// $Id: PView.cpp,v 1.18 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -300,22 +300,56 @@ bool PView::readMSH(std::string filename, int fileIndex)
return false;
}
- Msg(INFO, "Reading post-processing data from MSH file...");
-
- // FIXME: to be implemented!
- int index = 0;
- PViewDataGModel *d = new PViewDataGModel(GModel::current());
- if(!d->readMSH(fp)){
- Msg(GERROR, "Could not read data in msh file");
- delete d;
- return false;
- }
- else{
- d->setFileName(filename);
- d->setFileIndex(index);
- new PView(d);
+ char str[256];
+ double version;
+ int format, size, index = -1;
+
+ while(1) {
+
+ do {
+ if(!fgets(str, 256, fp) || feof(fp))
+ break;
+ } while(str[0] != '$');
+
+ if(feof(fp))
+ break;
+
+ if(!strncmp(&str[1], "NodeData", 8)) {
+ index++;
+ if(fileIndex < 0 || fileIndex == index){
+ // either get existing viewData, or create new one
+ // check timestep storage limit, do magic if exceeded
+ PView *p = 0;
+ PViewDataGModel *d;
+ // read view name and timestep
+ // if append
+ // p = getView
+ // d = p->getData()
+ // else
+ d = new PViewDataGModel(GModel::current());
+ if(!d->readMSH(fp)){
+ Msg(GERROR, "Could not read data in msh file");
+ delete d;
+ return false;
+ }
+ else{
+ d->setFileName(filename);
+ d->setFileIndex(index);
+ if(!p) new PView(d);
+ }
+ }
+ }
+
+ do {
+ if(!fgets(str, 256, fp) || feof(fp)){
+ Msg(GERROR, "Prematured end of file");
+ break;
+ }
+ } while(str[0] != '$');
+
}
-
+
+ fclose(fp);
return true;
}
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 1e8d5e6e7577e366b818e3ec17095db6bef4e3a0..48c71a2ada5228abedd0084dcc1b46accc375c9d 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModelIO.cpp,v 1.1 2008-02-24 19:59:03 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.2 2008-03-01 01:32:03 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -27,30 +27,14 @@
#include "PViewDataGModel.h"
/*
- if(!fgets(str, sizeof(str), fp)) return 0;
- // name = str[1] + remove final "
- int timeStep, numData, numComponents;
- double time;
- if(_vertexVector.empty() && _vertexMap.empty()){
- Msg(GERROR, "Mesh vertex information missing: impossible to load dataset");
- return false;
- }
-
- if(fscanf(fp, "%d %lf %d %d", &timeStep, &time, &numData, &numComponents) != 4)
- return 0;
- Msg(INFO, "%d node data", numData);
-
- //std::map<int, int> nodeNumber, nodeIndex
- PViewDataGModel *p = getPViewDataGModel(name)
- if(p){ // add data to existing view
- if(!p.count(timeStep)){
- // we don't have any data for this time step
- p[timeStep] = new nodeData(numNodes);
- }
- data = p[timeStep];
- if(num
- data.scalar.indices.append();
- data.scalar.values.append();
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ // name = str[1] + remove final "
+ int timeStep, numData, numComponents;
+ double time;
+ if(fscanf(fp, "%d %lf %d %d", &timeStep, &time, &numData, &numComponents) != 4)
+ return 0;
+ Msg(INFO, "%d node data", numData);
+ fill in data
*/
bool PViewDataGModel::readMSH(FILE *fp)
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 9d5e0c608fc972652d03d31cf3635b7131345cff..89718988623f6204c4de0e50a6402df20a55c8a1 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,6 +1,7 @@
-$Id: VERSIONS,v 1.398 2008-02-27 12:39:28 geuzaine Exp $
+$Id: VERSIONS,v 1.399 2008-03-01 01:32:03 geuzaine Exp $
-2.1.1 (Feb , 2008): small bug fixes (higher order elements)
+2.1.1 (Feb , 2008): small bug fixes (second order meshes, combine
+views, divide and conquer crash).
2.1.0 (Feb 23, 2008): new post-processing database; complete rewrite
of post-processing drawing code; improved surface mesh algorithms;