diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 9ae7a5e18e14544ee61f88f1034245d0f9a88c17..2b6f73db4a1655db9c467dc5471d552ac641e4a0 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.21 2006-11-20 12:44:09 remacle Exp $
+// $Id: GFace.cpp,v 1.22 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,7 +23,7 @@
#include "GFace.h"
#include "GEdge.h"
#include "Message.h"
-#include "Utils.h"
+#include "Numeric.h"
#if defined(HAVE_GSL)
#include <gsl/gsl_vector.h>
@@ -232,7 +232,7 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
double ex[3], t1[3], t2[3];
// check coherence of results for non-plane surfaces
- if(geomType() != GEntity::Plane) {
+ if(geomType() != GEntity::Plane && geomType() != GEntity::DiscreteSurface) {
double res2[3], c[3], cosc, sinc, angplan;
double eps = 1.e-3;
@@ -352,6 +352,13 @@ void GFace::getMeanPlaneData(double VX[3], double VY[3],
z = meanPlane.z;
}
+void GFace::getMeanPlaneData(double plan[3][3]) const
+{
+ for(int i = 0; i < 3; i++)
+ for(int j = 0; j < 3; j++)
+ plan[i][j] = meanPlane.plan[i][j];
+}
+
double GFace::curvature (const SPoint2 ¶m) const
{
if (geomType() == Plane)return 0;
diff --git a/Geo/GFace.h b/Geo/GFace.h
index b25055dd9e08b2afb70739b70b49b89154422ba4..6666e8c091f99953a2bf41fdc61874229bd9c48d 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -134,6 +134,7 @@ class GFace : public GEntity
// Get the mean plane info
void getMeanPlaneData(double VX[3], double VY[3],
double &x, double &y, double &z) const;
+ void getMeanPlaneData(double plan[3][3]) const;
// a crude graphical representation using a "cross" defined by pairs
// of start/end points
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index d5b2f8bfe8161e4f8f1b7a0c8665bb812d53b2db..0730b9e563aa143e4515b7e3b2a5b6120bdae878 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.60 2006-11-25 17:07:45 geuzaine Exp $
+// $Id: Geo.cpp,v 1.61 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,7 +23,6 @@
#include "Numeric.h"
#include "Geo.h"
#include "GeoInterpolation.h"
-#include "Utils.h"
#include "Context.h"
extern Mesh *THEM;
@@ -3058,3 +3057,21 @@ bool ProjectPointOnSurface(Surface * s, Vertex * p, double *u, double *v)
return true;
}
+void direction(Vertex * v1, Vertex * v2, double d[3])
+{
+ d[0] = v2->Pos.X - v1->Pos.X;
+ d[1] = v2->Pos.Y - v1->Pos.Y;
+ d[2] = v2->Pos.Z - v1->Pos.Z;
+}
+
+void Projette(Vertex * v, double mat[3][3])
+{
+ double X, Y, Z;
+
+ X = v->Pos.X * mat[0][0] + v->Pos.Y * mat[0][1] + v->Pos.Z * mat[0][2];
+ Y = v->Pos.X * mat[1][0] + v->Pos.Y * mat[1][1] + v->Pos.Z * mat[1][2];
+ Z = v->Pos.X * mat[2][0] + v->Pos.Y * mat[2][1] + v->Pos.Z * mat[2][2];
+ v->Pos.X = X;
+ v->Pos.Y = Y;
+ v->Pos.Z = Z;
+}
diff --git a/Geo/Geo.h b/Geo/Geo.h
index c46a8a4d5d7b82ebf4dfb7aa71eee21b3f9afbd7..c1095956c7db625919fc15ce7b96ce996fe65b79 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -331,4 +331,7 @@ int recognize_seg(int typ, List_T *liste, int *seg);
int recognize_loop(List_T *liste, int *loop);
int recognize_surfloop(List_T *liste, int *loop);
+void direction(Vertex * v1, Vertex * v2, double d[3]);
+void Projette(Vertex * v, double mat[3][3]);
+
#endif
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index abfb02b088a82fe20c5720c06c1667eb139e9983..ab53ed1f2a2f41b992678792f8b348da4a276830 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.2 2006-11-25 17:07:45 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.3 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -24,7 +24,6 @@
#include "GeoInterpolation.h"
#include "GeoStringInterface.h"
#include "GeoUtils.h"
-#include "Utils.h"
#include "Numeric.h"
extern Mesh *THEM;
diff --git a/Geo/Makefile b/Geo/Makefile
index ee7d667839eefdc582dd247db19f0fbcc217c507..d7cf8727db41c1d4b9f22eeaddd000475bb3f587 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.115 2006-11-25 16:52:43 geuzaine Exp $
+# $Id: Makefile,v 1.116 2006-11-25 18:03:49 geuzaine Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -96,8 +96,7 @@ GFace.o: GFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h MFace.h \
../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
ExtrudeParams.h GFace.h GEdgeLoop.h Pair.h GRegion.h \
- ../Common/SmoothNormals.h ../Common/Message.h ../Mesh/Utils.h \
- ../Geo/Geo.h ../DataStr/Tree.h ../DataStr/avl.h ../Geo/ExtrudeParams.h
+ ../Common/SmoothNormals.h ../Common/Message.h
GRegion.o: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h SPoint2.h \
GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h MFace.h \
@@ -117,7 +116,7 @@ gmshFace.o: gmshFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
ExtrudeParams.h GFace.h GEdgeLoop.h Pair.h GRegion.h \
../Common/SmoothNormals.h gmshVertex.h Geo.h ../DataStr/Tree.h \
../DataStr/avl.h gmshEdge.h gmshFace.h GeoInterpolation.h \
- ../Mesh/Utils.h ../Geo/Geo.h ../Common/Message.h
+ ../Common/Message.h
gmshRegion.o: gmshRegion.cpp GModel.h GVertex.h GEntity.h Range.h \
SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h \
SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
@@ -144,8 +143,7 @@ OCCFace.o: OCCFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
ExtrudeParams.h GFace.h GEdgeLoop.h Pair.h GRegion.h \
../Common/SmoothNormals.h OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h \
- ../Common/Message.h ../Mesh/Utils.h ../Geo/Geo.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../Geo/ExtrudeParams.h
+ ../Common/Message.h
OCCRegion.o: OCCRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h SPoint2.h \
GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h MFace.h \
@@ -212,7 +210,7 @@ Geo.o: Geo.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
../DataStr/List.h ../DataStr/Tree.h ../Numeric/Numeric.h Geo.h \
../Common/GmshDefines.h ExtrudeParams.h GeoInterpolation.h \
- ../Mesh/Utils.h ../Geo/Geo.h ../Common/Context.h
+ ../Common/Context.h
GeoStringInterface.o: GeoStringInterface.cpp ../Common/Gmsh.h \
../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
@@ -226,8 +224,7 @@ GeoInterpolation.o: GeoInterpolation.cpp ../Common/Gmsh.h \
../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
../DataStr/Tree.h Geo.h ../Common/GmshDefines.h ExtrudeParams.h \
- GeoInterpolation.h GeoStringInterface.h GeoUtils.h ../Mesh/Utils.h \
- ../Geo/Geo.h ../Numeric/Numeric.h
+ GeoInterpolation.h GeoStringInterface.h GeoUtils.h ../Numeric/Numeric.h
GeoUtils.o: GeoUtils.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 1af8003970b8d914362861cbecca9369d563b378..a80802a30564544a0cf6b9297bd955169924e619 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.13 2006-11-25 17:07:45 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.14 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -25,7 +25,7 @@
#include "OCCEdge.h"
#include "OCCFace.h"
#include "Message.h"
-#include "Utils.h"
+#include "Numeric.h"
#if defined(HAVE_OCC)
#include "Geom_CylindricalSurface.hxx"
@@ -213,10 +213,7 @@ int OCCFace::containsPoint(const SPoint3 &pt) const
pl.Coefficients ( n[0], n[1], n[2], c );
norme(n);
double angle = 0.;
- Vertex v,v1,v2;
- v.Pos.X = pt.x();
- v.Pos.Y = pt.y();
- v.Pos.Z = pt.z();
+ double v[3] = {pt.x(), pt.y(), pt.z()};
for(std::list<GEdge*>::const_iterator it = l_edges.begin(); it != l_edges.end(); it++) {
GEdge *c = *it;
int N=10;
@@ -226,13 +223,9 @@ int OCCFace::containsPoint(const SPoint3 &pt) const
double u2 = (double)(j + 1) / (double)N;
GPoint pp1 = c->point(range.low() + u1 * (range.high() - range.low() ));
GPoint pp2 = c->point(range.low() + u2 * (range.high() - range.low() ));
- v1.Pos.X = pp1.x();
- v1.Pos.Y = pp1.y();
- v1.Pos.Z = pp1.z();
- v2.Pos.X = pp2.x();
- v2.Pos.Y = pp2.y();
- v2.Pos.Z = pp2.z();
- angle += angle_plan(&v, &v1, &v2, n);
+ double v1[3] = {pp1.x(), pp1.y(), pp1.z()};
+ double v2[3] = {pp2.x(), pp2.y(), pp2.z()};
+ angle += angle_plan(v, v1, v2, n);
}
}
// we're inside if angle equals 2 * pi
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 29e7047dc12e4004a1a4a5c7471e843cc97cfb5e..adc962d18e5d889a90036087c686c98d7551a5cb 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.22 2006-11-25 16:52:43 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.23 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -25,7 +25,7 @@
#include "gmshFace.h"
#include "Geo.h"
#include "GeoInterpolation.h"
-#include "Utils.h"
+#include "Numeric.h"
#include "Message.h"
extern Mesh *THEM;
@@ -232,10 +232,7 @@ int gmshFace::containsPoint(const SPoint3 &pt) const
// for the (possibly wrong) orientation at the end
double n[3] = {meanPlane.a, meanPlane.b, meanPlane.c};
double angle = 0.;
- Vertex v;
- v.Pos.X = pt.x();
- v.Pos.Y = pt.y();
- v.Pos.Z = pt.z();
+ double v[3] = {pt.x(), pt.y(), pt.z()};
for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
Curve *c;
List_Read(s->Generatrices, i, &c);
@@ -245,7 +242,9 @@ int gmshFace::containsPoint(const SPoint3 &pt) const
double u2 = (double)(j + 1) / (double)N;
Vertex p1 = InterpolateCurve(c, u1, 0);
Vertex p2 = InterpolateCurve(c, u2, 0);
- angle += angle_plan(&v, &p1, &p2, n);
+ double v1[3] = {p1.Pos.X, p1.Pos.Y, p1.Pos.Z};
+ double v2[3] = {p2.Pos.X, p2.Pos.Y, p2.Pos.Z};
+ angle += angle_plan(v, v1, v2, n);
}
}
// we're inside if angle equals 2 * pi
diff --git a/Mesh/Makefile b/Mesh/Makefile
index d08c771ada79bf06a73ecb3d74f84cc714626f67..a5c8636c2ec9a45035a1c10726ae33e80d359097 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.140 2006-11-25 16:52:44 geuzaine Exp $
+# $Id: Makefile,v 1.141 2006-11-25 18:03:49 geuzaine Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -34,7 +34,6 @@ SRC = DivideAndConquer.cpp \
BackgroundMesh.cpp \
BDS.cpp \
Generator.cpp \
- Utils.cpp \
meshGEdge.cpp \
meshGFace.cpp \
meshGFaceTransfinite.cpp \
@@ -120,12 +119,6 @@ Generator.o: Generator.cpp ../Common/Gmsh.h ../Common/Message.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \
../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h BackgroundMesh.h \
SecondOrder.h
-Utils.o: Utils.cpp ../Common/Gmsh.h ../Common/Message.h \
- ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
- ../Geo/Geo.h ../Common/GmshDefines.h ../Geo/ExtrudeParams.h \
- ../Geo/GeoInterpolation.h ../Geo/Geo.h Utils.h ../Numeric/Numeric.h \
- ../Common/Context.h
meshGEdge.o: meshGEdge.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -139,8 +132,7 @@ meshGEdge.o: meshGEdge.cpp ../Common/Gmsh.h ../Common/Message.h \
../Geo/SVector3.h ../Numeric/Numeric.h ../Common/Context.h \
../Geo/ExtrudeParams.h ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h \
- ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h Utils.h \
- ../Geo/Geo.h ../Geo/ExtrudeParams.h BackgroundMesh.h
+ ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h BackgroundMesh.h
meshGFace.o: meshGFace.cpp meshGFace.h DivideAndConquer.h \
BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
@@ -153,11 +145,10 @@ meshGFace.o: meshGFace.cpp meshGFace.h DivideAndConquer.h \
../Common/Context.h ../DataStr/List.h ../Geo/ExtrudeParams.h \
../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/Pair.h ../Geo/ExtrudeParams.h Utils.h ../Geo/Geo.h \
- ../DataStr/Tree.h ../DataStr/avl.h ../Geo/ExtrudeParams.h \
- ../Common/Message.h BDS.h ../Common/Views.h ../Common/ColorTable.h \
- ../Common/VertexArray.h ../Common/SmoothNormals.h \
- ../Common/AdaptiveViews.h ../Common/GmshMatrix.h
+ ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Common/Message.h BDS.h \
+ ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
+ ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
+ ../Common/GmshMatrix.h
meshGFaceTransfinite.o: meshGFaceTransfinite.cpp meshGFace.h \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
deleted file mode 100644
index 6ee7c7a85a3a92bd4513886d01bb5de83c06a039..0000000000000000000000000000000000000000
--- a/Mesh/Utils.cpp
+++ /dev/null
@@ -1,724 +0,0 @@
-// $Id: Utils.cpp,v 1.35 2006-11-25 16:52:44 geuzaine Exp $
-//
-// Copyright (C) 1997-2006 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>.
-//
-// Contributor(s):
-// Pierre Badel
-//
-
-#include "Gmsh.h"
-#include "Geo.h"
-#include "GeoInterpolation.h"
-#include "Utils.h"
-#include "Numeric.h"
-#include "Context.h"
-
-#if defined(HAVE_GSL)
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
-#else
-#define NRANSI
-#include "nrutil.h"
-void dsvdcmp(double **a, int m, int n, double w[], double **v);
-#endif
-
-extern Context_T CTX;
-
-void direction(Vertex * v1, Vertex * v2, double d[3])
-{
- d[0] = v2->Pos.X - v1->Pos.X;
- d[1] = v2->Pos.Y - v1->Pos.Y;
- d[2] = v2->Pos.Z - v1->Pos.Z;
-}
-
-void Projette(Vertex * v, double mat[3][3])
-{
- double X, Y, Z;
-
- X = v->Pos.X * mat[0][0] + v->Pos.Y * mat[0][1] + v->Pos.Z * mat[0][2];
- Y = v->Pos.X * mat[1][0] + v->Pos.Y * mat[1][1] + v->Pos.Z * mat[1][2];
- Z = v->Pos.X * mat[2][0] + v->Pos.Y * mat[2][1] + v->Pos.Z * mat[2][2];
- v->Pos.X = X;
- v->Pos.Y = Y;
- v->Pos.Z = Z;
-}
-
-/*
- Le concept d'un plan moyen calcule au sens des moidres carres n'est
- pas le bon pour les surfaces non-planes : imagine un quart de cercle
- extrude d'une faible hauteur. Le plan moyen sera dans le plan du
- cercle! En attendant mieux, il y a donc un test de coherence
- supplementaire pour les surfaces non-planes. */
-
-void MeanPlane(List_T * points, Surface * s)
-{
- int i, j, min, ndata, na;
- double res[4], ex[3], t1[3], t2[3], svd[3];
- Vertex *v;
- double xm = 0., ym = 0., zm = 0.;
-
- ndata = List_Nbr(points);
- na = 3;
- for(i = 0; i < ndata; i++) {
- List_Read(points, i, &v);
- xm += v->Pos.X;
- ym += v->Pos.Y;
- zm += v->Pos.Z;
- }
- xm /= (double)ndata;
- ym /= (double)ndata;
- zm /= (double)ndata;
-
-#if defined(HAVE_GSL)
- gsl_matrix *U = gsl_matrix_alloc(ndata, na);
- gsl_matrix *V = gsl_matrix_alloc(na, na);
- gsl_vector *W = gsl_vector_alloc(na);
- gsl_vector *TMPVEC = gsl_vector_alloc(na);
- for(i = 0; i < ndata; i++) {
- List_Read(points, i, &v);
- gsl_matrix_set(U, i, 0, v->Pos.X - xm);
- gsl_matrix_set(U, i, 1, v->Pos.Y - ym);
- gsl_matrix_set(U, i, 2, v->Pos.Z - zm);
- }
- gsl_linalg_SV_decomp(U, V, W, TMPVEC);
- svd[0] = gsl_vector_get(W, 0);
- svd[1] = gsl_vector_get(W, 1);
- svd[2] = gsl_vector_get(W, 2);
- if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
- min = 0;
- else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
- min = 1;
- else
- min = 2;
- res[0] = gsl_matrix_get(V, 0, min);
- res[1] = gsl_matrix_get(V, 1, min);
- res[2] = gsl_matrix_get(V, 2, min);
- norme(res);
- gsl_matrix_free(U);
- gsl_matrix_free(V);
- gsl_vector_free(W);
- gsl_vector_free(TMPVEC);
-#else
- double **U = dmatrix(1, ndata, 1, na);
- double **V = dmatrix(1, na, 1, na);
- double *W = dvector(1, na);
- for(i = 0; i < ndata; i++) {
- List_Read(points, i, &v);
- U[i + 1][1] = v->Pos.X - xm;
- U[i + 1][2] = v->Pos.Y - ym;
- U[i + 1][3] = v->Pos.Z - zm;
- }
- dsvdcmp(U, ndata, na, W, V);
- if(fabs(W[1]) < fabs(W[2]) && fabs(W[1]) < fabs(W[3]))
- min = 1;
- else if(fabs(W[2]) < fabs(W[1]) && fabs(W[2]) < fabs(W[3]))
- min = 2;
- else
- min = 3;
- svd[0] = W[1];
- svd[1] = W[2];
- svd[2] = W[3];
- res[0] = V[1][min];
- res[1] = V[2][min];
- res[2] = V[3][min];
- norme(res);
- free_dmatrix(U, 1, ndata, 1, na);
- free_dmatrix(V, 1, na, 1, na);
- free_dvector(W, 1, na);
-#endif
-
- // check coherence of results for non-plane surfaces
- if(s->Typ != MSH_SURF_PLAN) {
- double res2[3], c[3], cosc, sinc, angplan;
- double eps = 1.e-3;
- Vertex v1, v2, v3;
- v1 = InterpolateSurface(s, 0.5, 0.5, 0, 0);
- v2 = InterpolateSurface(s, 0.5 + eps, 0.5, 0, 0);
- v3 = InterpolateSurface(s, 0.5, 0.5 + eps, 0, 0);
- t1[0] = v2.Pos.X - v1.Pos.X;
- t1[1] = v2.Pos.Y - v1.Pos.Y;
- t1[2] = v2.Pos.Z - v1.Pos.Z;
- t2[0] = v3.Pos.X - v1.Pos.X;
- t2[1] = v3.Pos.Y - v1.Pos.Y;
- t2[2] = v3.Pos.Z - v1.Pos.Z;
- norme(t1);
- norme(t2);
- // prodve(t1, t2, res2);
- // Warning: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?)
- prodve(t2, t1, res2);
- norme(res2);
- prodve(res, res2, c);
- prosca(res, res2, &cosc);
- sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
- angplan = myatan2(sinc, cosc);
- angplan = angle_02pi(angplan) * 180. / Pi;
- if((angplan > 70 && angplan < 110) || (angplan > 260 && angplan < 280)) {
- Msg(INFO, "SVD failed (angle=%g): using rough algo...", angplan);
- res[0] = res2[0];
- res[1] = res2[1];
- res[2] = res2[2];
- goto end;
- }
- }
-
- ex[0] = ex[1] = ex[2] = 0.0;
- if(res[0] == 0.)
- ex[0] = 1.0;
- else if(res[1] == 0.)
- ex[1] = 1.0;
- else
- ex[2] = 1.0;
-
- prodve(res, ex, t1);
- norme(t1);
- prodve(t1, res, t2);
- norme(t2);
-
-end:
- res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
-
- for(i = 0; i < 3; i++)
- s->plan[0][i] = t1[i];
- for(i = 0; i < 3; i++)
- s->plan[1][i] = t2[i];
- for(i = 0; i < 3; i++)
- s->plan[2][i] = res[i];
-
- s->a = res[0];
- s->b = res[1];
- s->c = res[2];
- s->d = res[3];
-
- Msg(DEBUG1, "Surface: %d", s->Num);
- Msg(DEBUG2, "SVD : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min);
- Msg(DEBUG2, "Plane : (%g x + %g y + %g z = %g)", s->a, s->b, s->c, s->d);
- Msg(DEBUG2, "Normal : (%g , %g , %g )", s->a, s->b, s->c);
- Msg(DEBUG3, "t1 : (%g , %g , %g )", t1[0], t1[1], t1[2]);
- Msg(DEBUG3, "t2 : (%g , %g , %g )", t2[0], t2[1], t2[2]);
-
- for(i = 0; i < 3; i++) {
- for(j = 0; j < 3; j++) {
- s->invplan[i][j] = s->plan[j][i];
- }
- }
-
- //check coherence for plane surfaces
- if(s->Typ == MSH_SURF_PLAN) {
- Curve *c;
- for(i = 0; i < List_Nbr(s->Generatrices); i++) {
- List_Read(s->Generatrices, i, &c);
- if(c->Num > 0) {
- List_Read(s->Generatrices, i, &c);
- for(j = 0; j < List_Nbr(c->Control_Points); j++) {
- List_Read(c->Control_Points, j, &v);
- double d =
- s->a * v->Pos.X + s->b * v->Pos.Y + s->c * v->Pos.Z - s->d;
- if(fabs(d) > CTX.lc * 1.e-3) {
- Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!",
- s->Num, s->a, s->b, s->c, s->d);
- Msg(GERROR3, "Control point %d = (%g,%g,%g), val=%g",
- v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, d);
- return;
- }
- }
- }
- }
- }
-
-}
-
-
-#define Precision 1.e-8
-#define MaxIter 25
-#define NumInitGuess 11
-
-void find_bestuv(Surface * s, double X, double Y,
- double *U, double *V, double *Z, int N)
-{
- double d, mina, min = 0., minu = 0., minv = 0., Unew, Vnew;
- static int i, j;
- Vertex P;
-
- d = 1. / (double)N;
-
- for(i = 0; i <= N; i++) {
- for(j = 0; j <= N; j++) {
- Unew = ((double)i) * d;
- Vnew = ((double)j) * d;
- P = InterpolateSurface(s, Unew, Vnew, 0, 0);
- if(!i && !j) {
- min = myhypot(X - P.Pos.X, Y - P.Pos.Y);
- minu = Unew;
- minv = Vnew;
- *Z = P.Pos.Z;
- }
- else {
- if((mina = myhypot(X - P.Pos.X, Y - P.Pos.Y)) < min) {
- min = mina;
- minu = Unew;
- minv = Vnew;
- *Z = P.Pos.Z;
- }
- }
- }
- }
- *U = minu;
- *V = minv;
-}
-
-void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
-{
- int i, j, k, n = 3;
- double TT[3][3];
-
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- II[i - 1][j - 1] = 0.0;
- TT[i - 1][j - 1] = 0.0;
- }
- }
-
-#if defined(HAVE_GSL)
- gsl_matrix *M = gsl_matrix_alloc(3, 3);
- gsl_matrix *V = gsl_matrix_alloc(3, 3);
- gsl_vector *W = gsl_vector_alloc(3);
- gsl_vector *TMPVEC = gsl_vector_alloc(3);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- gsl_matrix_set(M, i - 1, j - 1, MM[i - 1][j - 1]);
- }
- }
- gsl_linalg_SV_decomp(M, V, W, TMPVEC);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- double ww = gsl_vector_get(W, i - 1);
- if(fabs(ww) > 1.e-16) { //singular value precision
- TT[i - 1][j - 1] += gsl_matrix_get(M, j - 1, i - 1) / ww;
- }
- }
- }
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- for(k = 1; k <= n; k++) {
- II[i - 1][j - 1] +=
- gsl_matrix_get(V, i - 1, k - 1) * TT[k - 1][j - 1];
- }
- }
- }
- gsl_matrix_free(M);
- gsl_matrix_free(V);
- gsl_vector_free(W);
- gsl_vector_free(TMPVEC);
-#else
- double **M = dmatrix(1, 3, 1, 3);
- double **V = dmatrix(1, 3, 1, 3);
- double *W = dvector(1, 3);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- M[i][j] = MM[i - 1][j - 1];
- }
- }
- dsvdcmp(M, n, n, W, V);
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- if(fabs(W[i]) > 1.e-16) { //singular value precision
- TT[i - 1][j - 1] += M[j][i] / W[i];
- }
- }
- }
- for(i = 1; i <= n; i++) {
- for(j = 1; j <= n; j++) {
- for(k = 1; k <= n; k++) {
- II[i - 1][j - 1] += V[i][k] * TT[k - 1][j - 1];
- }
- }
- }
- free_dmatrix(M, 1, n, 1, n);
- free_dmatrix(V, 1, n, 1, n);
- free_dvector(W, 1, n);
-#endif
-}
-
-
-/*
- X = X ( u, v )
- Y = Y ( u, v )
- Z = Z ( u, v )
-
-*/
-
-
-void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
- double relax)
-{
- double Unew = 0., Vnew = 0., err,err2;
- int iter;
- Vertex D_u, D_v, P;
- double mat[3][3], jac[3][3];
- double umin, umax, vmin, vmax;
- double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1,0,1};
-
- if(s->Typ == MSH_SURF_NURBS) {
- umin = s->ku[0];
- umax = s->ku[s->OrderU + s->Nu];
- vmin = s->kv[0];
- vmax = s->kv[s->OrderV + s->Nv];
- }
- else {
- umin = vmin = 0.0-1.e-8;
- umax = vmax = 1.0+1.e-8;
- }
-
- for(int i = 0; i < NumInitGuess; i++){
- for(int j = 0; j < NumInitGuess; j++){
-
- *U = init[i];
- *V = init[j];
- err = 1.0;
- iter = 1;
-
- while(err > Precision && iter < MaxIter) {
- P = InterpolateSurface(s, *U, *V, 0, 0);
- D_u = InterpolateSurface(s, *U, *V, 1, 1);
- D_v = InterpolateSurface(s, *U, *V, 1, 2);
-
- mat[0][0] = D_u.Pos.X;
- mat[0][1] = D_u.Pos.Y;
- mat[0][2] = D_u.Pos.Z;
- mat[1][0] = D_v.Pos.X;
- mat[1][1] = D_v.Pos.Y;
- mat[1][2] = D_v.Pos.Z;
- mat[2][0] = 0.;
- mat[2][1] = 0.;
- mat[2][2] = 0.;
- invert_singular_matrix3x3(mat, jac);
-
- Unew = *U + relax *
- (jac[0][0] * (X - P.Pos.X) + jac[1][0] * (Y - P.Pos.Y) +
- jac[2][0] * (Z - P.Pos.Z));
- Vnew = *V + relax *
- (jac[0][1] * (X - P.Pos.X) + jac[1][1] * (Y - P.Pos.Y) +
- jac[2][1] * (Z - P.Pos.Z));
-
- err = DSQR(Unew - *U) + DSQR(Vnew - *V);
- // A BETTER TEST !! (JFR/AUG 2006)
- err2 = DSQR(X - P.Pos.X) + DSQR(Y - P.Pos.Y) + DSQR(Z - P.Pos.Z);
-
- iter++;
- *U = Unew;
- *V = Vnew;
- }
-
-
-
- if(iter < MaxIter && err <= Precision && err2 <= 1.e-5 &&
- Unew <= umax && Vnew <= vmax &&
- Unew >= umin && Vnew >= vmin){
- if (err2 > Precision)
- Msg(WARNING,"converged for i=%d j=%d (err=%g iter=%d) BUT err2 = %g", i, j, err, iter,err2);
- return;
- }
- }
- }
-
- if(relax < 1.e-6)
- Msg(GERROR, "Could not converge: surface mesh will be wrong");
- else {
- Msg(INFO, "Relaxation factor = %g", 0.75 * relax);
- XYZtoUV(s, X, Y, Z, U, V, 0.75 * relax);
- }
-
-}
-
-void XYtoUV(Surface * s, double *X, double *Y,
- double *U, double *V, double *Z, double relax)
-{
- double det, Unew = 0., Vnew = 0., err, mat[2][2], jac[2][2];
- int iter;
- Vertex D_u, D_v, P;
- double umin, umax, vmin, vmax;
- double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1};
-
- if(s->Typ == MSH_SURF_NURBS) {
- umin = s->ku[0];
- umax = s->ku[s->OrderU + s->Nu];
- vmin = s->kv[0];
- vmax = s->kv[s->OrderV + s->Nv];
- }
- else {
- umin = vmin = 0.0;
- umax = vmax = 1.0;
- }
-
- for(int i = 0; i < NumInitGuess; i++){
- for(int j = 0; j < NumInitGuess; j++){
-
- *U = init[i];
- *V = init[j];
- err = 1.0;
- iter = 1;
-
- while(err > Precision && iter < MaxIter) {
- P = InterpolateSurface(s, *U, *V, 0, 0);
- D_u = InterpolateSurface(s, *U, *V, 1, 1);
- D_v = InterpolateSurface(s, *U, *V, 1, 2);
- mat[0][0] = D_u.Pos.X;
- mat[0][1] = D_u.Pos.Y;
- mat[1][0] = D_v.Pos.X;
- mat[1][1] = D_v.Pos.Y;
- det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
-
- if(det == 0.0) {
- iter = MaxIter;
- break;
- }
-
- jac[0][0] = mat[1][1] / det;
- jac[0][1] = -mat[0][1] / det;
- jac[1][0] = -mat[1][0] / det;
- jac[1][1] = mat[0][0] / det;
-
- Unew = *U + relax * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
- Vnew = *V + relax * (jac[0][1] * (*X - P.Pos.X) + jac[1][1] * (*Y - P.Pos.Y));
-
- err = DSQR(Unew - *U) + DSQR(Vnew - *V);
-
- iter++;
- *U = Unew;
- *V = Vnew;
- }
-
- *Z = P.Pos.Z;
-
- if(iter < MaxIter && err <= Precision &&
- Unew <= umax && Vnew <= vmax &&
- Unew >= umin && Vnew >= vmin){
- //printf("converged for i=%d j=%d (err=%g iter=%d)\n", i, j, err, iter);
- return;
- }
- }
- }
-
- if(relax < 1.e-6) {
- Msg(GERROR, "Could not converge: surface mesh will probably be wrong");
- if(Unew > umax)
- *U = umax;
- if(Vnew > vmax)
- *V = vmax;
- if(Unew < umin)
- *U = umin;
- if(Vnew < vmin)
- *V = vmin;
- find_bestuv(s, *X, *Y, U, V, Z, 30);
- P = InterpolateSurface(s, *U, *V, 0, 0);
- *X = P.Pos.X;
- *Y = P.Pos.Y;
- *Z = P.Pos.Z;
- }
- else {
- Msg(INFO, "Relaxation factor = %g", 0.75 * relax);
- XYtoUV(s, X, Y, U, V, Z, 0.75 * relax);
- }
-
-}
-
-int Oriente(List_T * cu, double n[3])
-{
- int N, i, a, b, c;
- double cosa, sina, sum, v[3], w[3], u[3];
- Vertex *ver[3];
-
- N = List_Nbr(cu);
-
- if(N < 3){
- Msg(GERROR, "Unable to orient contour with less than 3 vertices");
- n[0] = 0.;
- n[1] = 0.;
- n[2] = 1.;
- return 0;
- }
-
- sum = 0.0;
- for(i = 0; i < N; i++) {
- if(i == N - 1) {
- a = N - 1;
- b = 1;
- c = 2;
- }
- else if(i == N - 2) {
- a = N - 2;
- b = N - 1;
- c = 1;
- }
- else {
- a = i;
- b = i + 1;
- c = i + 2;
- }
- List_Read(cu, a, &ver[0]);
- List_Read(cu, b, &ver[1]);
- List_Read(cu, c, &ver[2]);
-
- u[0] = ver[1]->Pos.X - ver[0]->Pos.X;
- u[1] = ver[1]->Pos.Y - ver[0]->Pos.Y;
- u[2] = ver[1]->Pos.Z - ver[0]->Pos.Z;
-
- v[0] = ver[2]->Pos.X - ver[1]->Pos.X;
- v[1] = ver[2]->Pos.Y - ver[1]->Pos.Y;
- v[2] = ver[2]->Pos.Z - ver[1]->Pos.Z;
- norme(u);
- norme(v);
- prodve(u, v, w);
- prosca(w, n, &sina);
- prosca(u, v, &cosa);
- sum += myatan2(sina, cosa);
- }
-
- if(sum < 0)
- return (1);
- else
- return (0);
-}
-
-double angle_3p(Vertex * V, Vertex * P1, Vertex * P2)
-{
- double PA[3], PB[3], angplan;
- double cosc, sinc, c[3];
-
- PA[0] = P1->Pos.X - V->Pos.X;
- PA[1] = P1->Pos.Y - V->Pos.Y;
- PA[2] = P1->Pos.Z - V->Pos.Z;
-
- PB[0] = P2->Pos.X - V->Pos.X;
- PB[1] = P2->Pos.Y - V->Pos.Y;
- PB[2] = P2->Pos.Z - V->Pos.Z;
-
- norme(PA);
- norme(PB);
-
- prodve(PA, PB, c);
-
- prosca(PA, PB, &cosc);
- sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
-
- angplan = myatan2(sinc, cosc);
-
- return angplan;
-}
-
-double angle_plan(Vertex * V, Vertex * P1, Vertex * P2, double n[3])
-{
- double PA[3], PB[3], angplan;
- double cosc, sinc, c[3];
-
- PA[0] = P1->Pos.X - V->Pos.X;
- PA[1] = P1->Pos.Y - V->Pos.Y;
- PA[2] = P1->Pos.Z - V->Pos.Z;
-
- PB[0] = P2->Pos.X - V->Pos.X;
- PB[1] = P2->Pos.Y - V->Pos.Y;
- PB[2] = P2->Pos.Z - V->Pos.Z;
-
- norme(PA);
- norme(PB);
-
- prodve(PA, PB, c);
-
- prosca(PA, PB, &cosc);
- prosca(c, n, &sinc);
- angplan = myatan2(sinc, cosc);
-
- return angplan;
-}
-
-double angle_3pts(Vertex * a, Vertex * b, Vertex * c)
-{
- double L, prosca, angle;
-
- L = myhypot((a->Pos.X - b->Pos.X), (a->Pos.Y - b->Pos.Y)) *
- myhypot((b->Pos.X - c->Pos.X), (b->Pos.Y - c->Pos.Y));
-
- prosca = ((a->Pos.X - b->Pos.X) * (c->Pos.X - b->Pos.X) +
- (a->Pos.Y - b->Pos.Y) * (c->Pos.Y - b->Pos.Y)) / L;
-
- angle = myacos(prosca) * 180. / Pi;
- return (angle);
-}
-
-double trapeze(IntPoint * P1, IntPoint * P2)
-{
- return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
-}
-
-
-void RecursiveIntegration(IntPoint * from, IntPoint * to,
- double (*f) (double X), List_T * pPoints,
- double Prec, int *depth)
-{
- IntPoint P, p1;
- double err, val1, val2, val3;
-
- (*depth)++;
-
- P.t = 0.5 * (from->t + to->t);
- P.lc = f(P.t);
-
- val1 = trapeze(from, to);
- val2 = trapeze(from, &P);
- val3 = trapeze(&P, to);
-
- err = fabs(val1 - val2 - val3);
- // Msg(INFO,"Int %22.15 E %22.15 E %22.15 E\n", val1,val2,val3);
- if(((err < Prec) && (*depth > 1)) || (*depth > 25)) {
- List_Read(pPoints, List_Nbr(pPoints) - 1, &p1);
- P.p = p1.p + val2;
- List_Add(pPoints, &P);
-
- List_Read(pPoints, List_Nbr(pPoints) - 1, &p1);
- to->p = p1.p + val3;
- List_Add(pPoints, to);
- }
- else {
- RecursiveIntegration(from, &P, f, pPoints, Prec, depth);
- RecursiveIntegration(&P, to, f, pPoints, Prec, depth);
- }
- (*depth)--;
-}
-
-double Integration(double t1, double t2, double (*f) (double X),
- List_T * pPoints, double Prec)
-{
- int depth;
- IntPoint from, to;
-
- depth = 0;
-
- from.t = t1;
- from.lc = f(from.t);
- from.p = 0.0;
- List_Add(pPoints, &from);
-
- to.t = t2;
- to.lc = f(to.t);
- RecursiveIntegration(&from, &to, f, pPoints, Prec, &depth);
-
- List_Read(pPoints, List_Nbr(pPoints) - 1, &to);
- return (to.p);
-}
diff --git a/Mesh/Utils.h b/Mesh/Utils.h
deleted file mode 100644
index ed7a83441514c49b1f5e9eee555c135bf95b0815..0000000000000000000000000000000000000000
--- a/Mesh/Utils.h
+++ /dev/null
@@ -1,52 +0,0 @@
-#ifndef _UTILS_H_
-#define _UTILS_H_
-
-// Copyright (C) 1997-2006 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 "Geo.h"
-#include "List.h"
-
-void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
-void direction (Vertex * v1, Vertex * v2, double d[3]);
-void Projette (Vertex * v, double mat[3][3]);
-void MeanPlane(List_T *point, Surface *s);
-void find_bestuv (Surface * s, double X, double Y,
- double *U, double *V, double *Z, int N);
-void XYtoUV (Surface * s, double *X, double *Y,
- double *U, double *V, double *Z, double relax);
-void XYZtoUV (Surface *s, double X, double Y, double Z,
- double *U, double *V, double relax);
-int Oriente (List_T * cu, double n[3]);
-double angle_3p (Vertex * V, Vertex * P1, Vertex * P2);
-double angle_plan (Vertex * V, Vertex * P1, Vertex * P2, double n[3]);
-double angle_3pts (Vertex * a, Vertex * b, Vertex * c);
-
-typedef struct{
- int Num;
- double t, lc, p;
-}IntPoint;
-
-double trapeze (IntPoint * P1, IntPoint * P2);
-void RecursiveIntegration (IntPoint * from, IntPoint * to, double (*f) (double X),
- List_T * pPoints, double Prec, int *depth);
-double Integration (double t1, double t2, double (*f) (double X),
- List_T * pPoints, double Prec);
-
-#endif
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 06473a926b1ef71e780e73b34cb56f22a78025a5..c8e70dc1b1e888b48cbb3cf97ea00b030bbc7c86 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.18 2006-11-25 16:52:44 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.19 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,7 +23,6 @@
#include "meshGEdge.h"
#include "GEdge.h"
#include "GFace.h"
-#include "Utils.h"
#include "BackgroundMesh.h"
#include "Context.h"
#include "Message.h"
@@ -128,6 +127,71 @@ double F_One_bis(double t)
return norm(der);
}
+typedef struct{
+ int Num;
+ double t, lc, p;
+}IntPoint;
+
+double trapeze(IntPoint * P1, IntPoint * P2)
+{
+ return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
+}
+
+void RecursiveIntegration(IntPoint * from, IntPoint * to,
+ double (*f) (double X), List_T * pPoints,
+ double Prec, int *depth)
+{
+ IntPoint P, p1;
+ double err, val1, val2, val3;
+
+ (*depth)++;
+
+ P.t = 0.5 * (from->t + to->t);
+ P.lc = f(P.t);
+
+ val1 = trapeze(from, to);
+ val2 = trapeze(from, &P);
+ val3 = trapeze(&P, to);
+
+ err = fabs(val1 - val2 - val3);
+ // Msg(INFO,"Int %22.15 E %22.15 E %22.15 E\n", val1,val2,val3);
+ if(((err < Prec) && (*depth > 1)) || (*depth > 25)) {
+ List_Read(pPoints, List_Nbr(pPoints) - 1, &p1);
+ P.p = p1.p + val2;
+ List_Add(pPoints, &P);
+
+ List_Read(pPoints, List_Nbr(pPoints) - 1, &p1);
+ to->p = p1.p + val3;
+ List_Add(pPoints, to);
+ }
+ else {
+ RecursiveIntegration(from, &P, f, pPoints, Prec, depth);
+ RecursiveIntegration(&P, to, f, pPoints, Prec, depth);
+ }
+ (*depth)--;
+}
+
+double Integration(double t1, double t2, double (*f) (double X),
+ List_T * pPoints, double Prec)
+{
+ int depth;
+ IntPoint from, to;
+
+ depth = 0;
+
+ from.t = t1;
+ from.lc = f(from.t);
+ from.p = 0.0;
+ List_Add(pPoints, &from);
+
+ to.t = t2;
+ to.lc = f(to.t);
+ RecursiveIntegration(&from, &to, f, pPoints, Prec, &depth);
+
+ List_Read(pPoints, List_Nbr(pPoints) - 1, &to);
+ return (to.p);
+}
+
void deMeshGEdge :: operator() (GEdge *ge)
{
for (unsigned int i=0;i<ge->mesh_vertices.size();i++) delete ge->mesh_vertices[i];
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a923e9f487765cf406cd5189d8558936c5355105..4a5aec99d5238748acfc5c3dbd259b11c07ef81f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.27 2006-11-25 17:07:45 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.28 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -28,7 +28,6 @@
#include "MVertex.h"
#include "MElement.h"
#include "Context.h"
-#include "Utils.h"
#include "GPoint.h"
#include "Message.h"
#include "Numeric.h"
@@ -87,8 +86,6 @@ int Orientation (std::vector<MVertex*> &cu)
return (0);
}
-
-
class fromCartesianToParametric
{
GFace *gf;
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 9e685eeea25577745f293d4bdbf324cd3e7aa3f3..43dc210da40dcbefe24221b2ad1f580c2e42a569 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.28 2006-09-11 17:58:19 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.29 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -34,6 +34,8 @@
#include <gsl/gsl_version.h>
#include <gsl/gsl_errno.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
void new_handler(const char *reason, const char *file, int line,
int gsl_errno)
@@ -59,21 +61,20 @@ int check_gsl()
// set new error handler
gsl_set_error_handler(&new_handler);
-
// initilize robust geometric predicates
gmsh::exactinit() ;
-
return 1;
}
#else
+#define NRANSI
+#include "nrutil.h"
+void dsvdcmp(double **a, int m, int n, double w[], double **v);
int check_gsl()
{
-
// initilize robust geometric predicates
gmsh::exactinit() ;
-
return 1;
}
@@ -280,6 +281,79 @@ double inv3x3(double mat[3][3], double inv[3][3])
return det;
}
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
+{
+ int i, j, k, n = 3;
+ double TT[3][3];
+
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ II[i - 1][j - 1] = 0.0;
+ TT[i - 1][j - 1] = 0.0;
+ }
+ }
+
+#if defined(HAVE_GSL)
+ gsl_matrix *M = gsl_matrix_alloc(3, 3);
+ gsl_matrix *V = gsl_matrix_alloc(3, 3);
+ gsl_vector *W = gsl_vector_alloc(3);
+ gsl_vector *TMPVEC = gsl_vector_alloc(3);
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ gsl_matrix_set(M, i - 1, j - 1, MM[i - 1][j - 1]);
+ }
+ }
+ gsl_linalg_SV_decomp(M, V, W, TMPVEC);
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ double ww = gsl_vector_get(W, i - 1);
+ if(fabs(ww) > 1.e-16) { //singular value precision
+ TT[i - 1][j - 1] += gsl_matrix_get(M, j - 1, i - 1) / ww;
+ }
+ }
+ }
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ for(k = 1; k <= n; k++) {
+ II[i - 1][j - 1] +=
+ gsl_matrix_get(V, i - 1, k - 1) * TT[k - 1][j - 1];
+ }
+ }
+ }
+ gsl_matrix_free(M);
+ gsl_matrix_free(V);
+ gsl_vector_free(W);
+ gsl_vector_free(TMPVEC);
+#else
+ double **M = dmatrix(1, 3, 1, 3);
+ double **V = dmatrix(1, 3, 1, 3);
+ double *W = dvector(1, 3);
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ M[i][j] = MM[i - 1][j - 1];
+ }
+ }
+ dsvdcmp(M, n, n, W, V);
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ if(fabs(W[i]) > 1.e-16) { //singular value precision
+ TT[i - 1][j - 1] += M[j][i] / W[i];
+ }
+ }
+ }
+ for(i = 1; i <= n; i++) {
+ for(j = 1; j <= n; j++) {
+ for(k = 1; k <= n; k++) {
+ II[i - 1][j - 1] += V[i][k] * TT[k - 1][j - 1];
+ }
+ }
+ }
+ free_dmatrix(M, 1, n, 1, n);
+ free_dmatrix(V, 1, n, 1, n);
+ free_dvector(W, 1, n);
+#endif
+}
+
double angle_02pi(double A3)
{
double DP = 2 * Pi;
@@ -292,6 +366,31 @@ double angle_02pi(double A3)
return A3;
}
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3])
+{
+ double PA[3], PB[3], angplan;
+ double cosc, sinc, c[3];
+
+ PA[0] = P1[0] - V[0];
+ PA[1] = P1[1] - V[1];
+ PA[2] = P1[2] - V[2];
+
+ PB[0] = P2[0] - V[0];
+ PB[1] = P2[1] - V[1];
+ PB[2] = P2[2] - V[2];
+
+ norme(PA);
+ norme(PB);
+
+ prodve(PA, PB, c);
+
+ prosca(PA, PB, &cosc);
+ prosca(c, n, &sinc);
+ angplan = myatan2(sinc, cosc);
+
+ return angplan;
+}
+
double triangle_area(double p0[3], double p1[3], double p2[3])
{
double a[3], b[3], c[3];
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 514fdd7d631f559d245bb3a6e65d580f8396403e..27309d5125d2c502977a9b35a11d03db6a7ff881 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -74,7 +74,9 @@ double det3x3(double mat[3][3]);
double trace3x3(double mat[3][3]);
double trace2 (double mat[3][3]);
double inv3x3(double mat[3][3], double inv[3][3]);
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
double angle_02pi(double A3);
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
double triangle_area(double p0[3], double p1[3], double p2[3]);
char float2char(float f);
float char2float(char c);
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 68c06fcfca47dfa70a4ff013a24621f88d723fe4..537ed2d427b6bff244c0330fb4e5e88bf4dda62f 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.115 2006-11-25 16:52:53 geuzaine Exp $
+# $Id: Makefile,v 1.116 2006-11-25 18:03:49 geuzaine Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -156,8 +156,18 @@ Triangulate.o: Triangulate.cpp ../Common/Gmsh.h ../Common/Message.h \
Plugin.h ../Common/Options.h ../Common/Views.h ../Common/ColorTable.h \
../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
../Common/AdaptiveViews.h ../Common/GmshMatrix.h Triangulate.h \
- ../Common/Context.h ../Geo/Geo.h ../Common/GmshDefines.h \
- ../Geo/ExtrudeParams.h ../Mesh/Utils.h
+ ../Common/Context.h ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/gmshFace.h \
+ ../Geo/Geo.h ../Common/GmshDefines.h ../Geo/ExtrudeParams.h \
+ ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+ ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+ ../Geo/GEntity.h ../Geo/MVertex.h ../Geo/GPoint.h ../Geo/SPoint2.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+ ../Geo/MElement.h ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/ExtrudeParams.h ../Geo/MElement.h \
+ ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
+ ../Geo/gmshVertex.h ../Geo/Geo.h ../Geo/GVertex.h ../Geo/Range.h
Warp.o: Warp.cpp Plugin.h ../Common/Options.h ../Common/Message.h \
../Common/Views.h ../Common/ColorTable.h ../DataStr/List.h \
../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
diff --git a/Plugin/Triangulate.cpp b/Plugin/Triangulate.cpp
index 79a2241b251e8430f34e70276679cd884e88c4dc..0016526292ceee062bb5dd36ebd18f6011d41c84 100644
--- a/Plugin/Triangulate.cpp
+++ b/Plugin/Triangulate.cpp
@@ -1,4 +1,4 @@
-// $Id: Triangulate.cpp,v 1.31 2006-11-25 16:52:53 geuzaine Exp $
+// $Id: Triangulate.cpp,v 1.32 2006-11-25 18:03:49 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -19,16 +19,15 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+#include <vector>
#include "Gmsh.h"
#include "Plugin.h"
#include "Triangulate.h"
-#include "List.h"
-#include "Tree.h"
#include "Views.h"
#include "Context.h"
#include "Malloc.h"
-#include "Geo.h"
-#include "Utils.h"
+#include "MVertex.h"
+#include "gmshFace.h"
extern Context_T CTX;
@@ -86,8 +85,8 @@ void GMSH_TriangulatePlugin::catchErrorMessage(char *errorMessage) const
#if !defined(HAVE_TRIANGLE)
-void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
- int nbTimeStep, int nbComp)
+static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
+ int nbTimeStep, int nbComp)
{
Msg(GERROR, "Triangle is not compiled in this version of Gmsh");
}
@@ -103,36 +102,42 @@ extern "C"
#include "triangle.h"
}
-void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
- int nbTimeStep, int nbComp)
+static void Project(MVertex *v, double mat[3][3])
+{
+ double X = v->x() * mat[0][0] + v->y() * mat[0][1] + v->z() * mat[0][2];
+ double Y = v->x() * mat[1][0] + v->y() * mat[1][1] + v->z() * mat[1][2];
+ double Z = v->x() * mat[2][0] + v->y() * mat[2][1] + v->z() * mat[2][2];
+ v->x() = X;
+ v->y() = Y;
+ v->z() = Z;
+}
+
+static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
+ int nbTimeStep, int nbComp)
{
if(nbIn < 3)
return;
- List_T *points = List_Create(nbIn, 1, sizeof(Vertex *));
+ std::vector<MVertex*> points;
+
int nb = List_Nbr(inList) / nbIn;
int j = 0;
- for(int i = 0; i < List_Nbr(inList); i += nb) {
- Vertex *v = Create_Vertex(j++,
- *(double *)List_Pointer_Fast(inList, i),
- *(double *)List_Pointer_Fast(inList, i + 1),
- *(double *)List_Pointer_Fast(inList, i + 2), 1., 0.);
- List_Add(points, &v);
- }
-
- Surface *s = Create_Surface(1, MSH_SURF_PLAN);
- MeanPlane(points, s);
-
- for(int i = 0; i < List_Nbr(points); i++) {
- Vertex *v;
- List_Read(points, i, &v);
- Projette(v, s->plan);
- }
-
- Free_Surface(&s, NULL);
+ for(int i = 0; i < List_Nbr(inList); i += nb)
+ points.push_back(new MVertex(*(double *)List_Pointer_Fast(inList, i),
+ *(double *)List_Pointer_Fast(inList, i + 1),
+ *(double *)List_Pointer_Fast(inList, i + 2),
+ 0, j++));
+
+ gmshFace *s = new gmshFace(0, 1);
+ s->computeMeanPlane(points);
+ double plan[3][3];
+ s->getMeanPlaneData(plan);
+ for(unsigned int i = 0; i < points.size(); i++)
+ Project(points[i], plan);
+ delete s;
struct triangulateio in;
- in.numberofpoints = List_Nbr(points);
+ in.numberofpoints = points.size();
in.pointlist = (REAL *) Malloc(in.numberofpoints * 2 * sizeof(REAL));
in.numberofpointattributes = 0;
in.pointattributelist = NULL;
@@ -146,20 +151,14 @@ void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
in.holelist = NULL;
j = 0;
- for(int i = 0; i < List_Nbr(points); i++) {
- Vertex *v;
- List_Read(points, i, &v);
- in.pointlist[j] = v->Pos.X;
- in.pointlist[j + 1] = v->Pos.Y;
+ for(int i = 0; i < points.size(); i++) {
+ in.pointlist[j] = points[i]->x();
+ in.pointlist[j + 1] = points[i]->y();
j += 2;
}
- for(int i = 0; i < List_Nbr(points); i++) {
- Vertex *v;
- List_Read(points, i, &v);
- Free_Vertex(&v, NULL);
- }
- List_Delete(points);
+ for(int i = 0; i < points.size(); i++)
+ delete points[i];
struct triangulateio out;
out.pointlist = NULL;