diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index b181ec6effaaf2c760458cadd4cdce4fb1d64b1a..e9d281fc03504d83d3f6c0d76232a1047d967a91 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -45,8 +45,8 @@ int GmshInitialize(int argc, char **argv)
GMSH_PluginManager::instance()->registerDefaultPlugins();
#endif
- // Check for buggy obsolete GSL versions
- check_gsl();
+ // Initialize numeric library (gsl, robust predicates)
+ Init_Numeric();
return 1;
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index ede3847246d09b85496378505c22fc9944f6955a..445a327db835e7f2c872359664efd1d1066dfa55 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -12,9 +12,6 @@
#include "meshGFaceDelaunayInsertion.h"
#include "qualityMeasures.h"
-bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t);
-bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t);
-
void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf)
{
FILE *f = fopen(iii, "w");
@@ -85,8 +82,8 @@ BDS_Vector::BDS_Vector(const BDS_Point &p2, const BDS_Point &p1)
{
}
-void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
- double c[3])
+static void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
+ double c[3])
{
double a[3] = {p1->X - p2->X, p1->Y - p2->Y, p1->Z - p2->Z};
double b[3] = {p1->X - p3->X, p1->Y - p3->Y, p1->Z - p3->Z};
@@ -95,8 +92,8 @@ void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
c[0] = a[1] * b[2] - a[2] * b[1];
}
-void vector_triangle_parametric(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
- double &c)
+static void vector_triangle_parametric(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
+ double &c)
{
double a[2] = {p1->u - p2->u, p1->v - p2->v};
double b[2] = {p1->u - p3->u, p1->v - p3->v};
@@ -110,14 +107,14 @@ void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
norme(c);
}
-double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+static double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
{
double c[3];
vector_triangle(p1, p2, p3, c);
return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
}
-double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
+static double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
{
double c;
vector_triangle_parametric(p1, p2, p3, c);
@@ -967,6 +964,91 @@ int BDS_Edge::numTriangles() const
return NT;
}
+// use robust predicates for not allowing to revert a triangle by
+// moving one of its vertices
+
+static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t)
+{
+ BDS_Point *pts[4];
+ t->getNodes(pts);
+
+ double pa[2] = {pts[0]->u, pts[0]->v};
+ double pb[2] = {pts[1]->u, pts[1]->v};
+ double pc[2] = {pts[2]->u, pts[2]->v};
+ double pd[2] = {pts[3]->u, pts[3]->v};
+
+ double ori_init1 = gmsh::orient2d(pa, pb, pc);
+ double ori_init2 = gmsh::orient2d(pc, pd, pa);
+
+ if(p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
+ }
+ else if(p == pts[1]){
+ pb[0] = u;
+ pb[1] = v;
+ }
+ else if(p == pts[2]){
+ pc[0] = u;
+ pc[1] = v;
+ }
+ else if(p == pts[3]){
+ pd[0] = u;
+ pd[1] = v;
+ }
+ else{
+ Msg::Error("Something wrong in move_point_parametric_quad");
+ return false;
+ }
+
+ double ori_final1 = gmsh::orient2d(pa, pb, pc);
+ double ori_final2 = gmsh::orient2d(pc, pd, pa);
+ // allow to move a point when a triangle was flat
+ return ori_init1*ori_final1 > 0 && ori_init2*ori_final2 > 0 ;
+}
+
+static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t)
+{
+ if (t->e4)
+ return test_move_point_parametric_quad(p,u,v,t);
+ BDS_Point *pts[4];
+ t->getNodes(pts);
+
+ double pa[2] = {pts[0]->u, pts[0]->v};
+ double pb[2] = {pts[1]->u, pts[1]->v};
+ double pc[2] = {pts[2]->u, pts[2]->v};
+
+ double a[2] = {pb[0] - pa[0], pb[1] - pa[1]};
+ double b[2] = {pc[0] - pa[0], pc[1] - pa[1]};
+
+ double area_init = fabs(a[0] * b[1] - a[1] * b[0]);
+
+ if(area_init == 0.0) return true;
+
+ double ori_init = gmsh::orient2d(pa, pb, pc);
+
+ if(p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
+ }
+ else if(p == pts[1]){
+ pb[0] = u;
+ pb[1] = v;
+ }
+ else if(p == pts[2]){
+ pc[0] = u;
+ pc[1] = v;
+ }
+ else
+ return false;
+
+ double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
+ if(area_final < 0.1 * area_init) return false;
+ double ori_final = gmsh::orient2d(pa, pb, pc);
+ // allow to move a point when a triangle was flat
+ return ori_init*ori_final > 0;
+}
+
bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
{
if(e->numfaces() != 2)
@@ -1069,186 +1151,6 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
return true;
}
-// use robust predicates for not allowing to revert a triangle by
-// moving one of its vertices
-
-bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t)
-{
- if (t->e4)
- return test_move_point_parametric_quad(p,u,v,t);
- BDS_Point *pts[4];
- t->getNodes(pts);
-
- double pa[2] = {pts[0]->u, pts[0]->v};
- double pb[2] = {pts[1]->u, pts[1]->v};
- double pc[2] = {pts[2]->u, pts[2]->v};
-
- double a[2] = {pb[0] - pa[0], pb[1] - pa[1]};
- double b[2] = {pc[0] - pa[0], pc[1] - pa[1]};
-
- double area_init = fabs(a[0] * b[1] - a[1] * b[0]);
-
- if(area_init == 0.0) return true;
-
- double ori_init = gmsh::orient2d(pa, pb, pc);
-
- if(p == pts[0]){
- pa[0] = u;
- pa[1] = v;
- }
- else if(p == pts[1]){
- pb[0] = u;
- pb[1] = v;
- }
- else if(p == pts[2]){
- pc[0] = u;
- pc[1] = v;
- }
- else
- return false;
-
- double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
- if(area_final < 0.1 * area_init) return false;
- double ori_final = gmsh::orient2d(pa, pb, pc);
- // allow to move a point when a triangle was flat
- return ori_init*ori_final > 0;
-}
-
-bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t)
-{
- BDS_Point *pts[4];
- t->getNodes(pts);
-
- double pa[2] = {pts[0]->u, pts[0]->v};
- double pb[2] = {pts[1]->u, pts[1]->v};
- double pc[2] = {pts[2]->u, pts[2]->v};
- double pd[2] = {pts[3]->u, pts[3]->v};
-
- double ori_init1 = gmsh::orient2d(pa, pb, pc);
- double ori_init2 = gmsh::orient2d(pc, pd, pa);
-
- if(p == pts[0]){
- pa[0] = u;
- pa[1] = v;
- }
- else if(p == pts[1]){
- pb[0] = u;
- pb[1] = v;
- }
- else if(p == pts[2]){
- pc[0] = u;
- pc[1] = v;
- }
- else if(p == pts[3]){
- pd[0] = u;
- pd[1] = v;
- }
- else{
- Msg::Error("Something wrong in move_point_parametric_quad");
- return false;
- }
-
- double ori_final1 = gmsh::orient2d(pa, pb, pc);
- double ori_final2 = gmsh::orient2d(pc, pd, pa);
- // allow to move a point when a triangle was flat
- return ori_init1*ori_final1 > 0 && ori_init2*ori_final2 > 0 ;
-}
-
-// d^2_i = (x^2_i - x)^T M (x_i - x)
-// = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2
-
-struct smoothVertexData{
- BDS_Point *p;
- GFace *gf;
- double scalu, scalv;
- std::list<BDS_Face*> ts;
-};
-
-double smoothing_objective_function(double U, double V, BDS_Point *v,
- std::list<BDS_Face*> &ts, double su, double sv,
- GFace *gf)
-{
- GPoint gp = gf->point(U * su, V * sv);
-
- const double oldX = v->X;
- const double oldY = v->Y;
- const double oldZ = v->Z;
- v->X = gp.x();
- v->Y = gp.y();
- v->Z = gp.z();
-
- std::list<BDS_Face*>::iterator it = ts.begin();
- std::list<BDS_Face*>::iterator ite = ts.end();
- double qMin = 1.;
- while(it != ite) {
- BDS_Face *t = *it;
- qMin = std::min(qmTriangle(*it, QMTRI_RHO), qMin);
- ++it;
- }
- v->X = oldX;
- v->Y = oldY;
- v->Z = oldZ;
- return -qMin;
-}
-
-void deriv_smoothing_objective_function(double U, double V,
- double &F, double &dFdU, double &dFdV,
- void *data)
-{
- smoothVertexData *svd = (smoothVertexData*)data;
- BDS_Point *v = svd->p;
- const double LARGE = 1.e5;
- const double SMALL = 1./LARGE;
- F = smoothing_objective_function(U, V, v, svd->ts,
- svd->scalu, svd->scalv, svd->gf);
- double F_U = smoothing_objective_function(U + SMALL, V, v, svd->ts,
- svd->scalu, svd->scalv, svd->gf);
- double F_V = smoothing_objective_function(U, V + SMALL, v, svd->ts,
- svd->scalu, svd->scalv, svd->gf);
- dFdU = (F_U - F) * LARGE;
- dFdV = (F_V - F) * LARGE;
-}
-
-double smooth_obj(double U, double V, void *data)
-{
- smoothVertexData *svd = (smoothVertexData*)data;
- return smoothing_objective_function(U, V, svd->p, svd->ts,
- svd->scalu, svd->scalv, svd->gf);
-}
-
-void optimize_vertex_position(GFace *GF, BDS_Point *data, double su, double sv)
-{
- if(data->g && data->g->classif_degree <= 1) return;
- smoothVertexData vd;
- vd.p = data;
- vd.scalu = su;
- vd.scalv = sv;
- vd.gf = GF;
- data->getTriangles(vd.ts);
- double U = data->u, V = data->v, val;
-
- val = smooth_obj(U, V, &vd);
- if(val < -.90) return;
-
- minimize_2(smooth_obj, deriv_smoothing_objective_function, &vd, 5, U,V,val);
- std::list<BDS_Face*>::iterator it = vd.ts.begin();
- std::list<BDS_Face*>::iterator ite = vd.ts.end();
- while(it != ite) {
- BDS_Face *t = *it;
- if(!test_move_point_parametric_triangle(data, U, V, t)){
- return;
- }
- ++it;
- }
-
- data->u = U;
- data->v = V;
- GPoint gp = GF->point(U * su, V * sv);
- data->X = gp.x();
- data->Y = gp.y();
- data->Z = gp.z();
-}
-
bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
{
if(!p->config_modified) return false;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index ca16a8c848155b5951688aa4aa3148e24eef11b9..e8ba1572ac7ee5448b7888cabe5c14e58167ece2 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -28,17 +28,6 @@ class GFace;
class GEdge;
class GVertex;
-void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
-void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
-double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
-double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
-void optimize_vertex_position(GFace *GF, BDS_Point *data, double su, double sv);
-void swap_config(BDS_Edge *e,
- BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
- BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
- BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
- BDS_Point **p41, BDS_Point **p42, BDS_Point **p43);
-
class BDS_GeomEntity
{
public:
@@ -468,6 +457,12 @@ class BDS_Mesh
void cleanup();
};
+void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
+void swap_config(BDS_Edge *e,
+ BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
+ BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
+ BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
+ BDS_Point **p41, BDS_Point **p42, BDS_Point **p43);
void outputScalarField(std::list<BDS_Face*> t, const char *fn, int param, GFace *gf=0);
void recur_tag(BDS_Face *t, BDS_GeomEntity *g);
int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 5aa4e049fc03bfdde81374001c94609f992806aa..e5ea84bcb15879546b83da5e2edc36422a615bad 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -380,10 +380,9 @@ void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
double F = -smooth_obj_HighOrderN(uv, &vdN);
if (F < eps){
double val;
- minimize_N(2 * vdN.v.size(),
- smooth_obj_HighOrderN,
+ minimize_N(2 * vdN.v.size(), smooth_obj_HighOrderN,
deriv_smoothing_objective_function_HighOrderN,
- &vdN, 1, uv,val);
+ &vdN, 1, uv, val);
double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
printf("%12.5E %12.5E\n", F, Fafter);
if (F < Fafter){
@@ -399,20 +398,18 @@ void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
}
}
-double angle3Points ( MVertex *p1, MVertex *p2, MVertex *p3 )
+double angle3Points(MVertex *p1, MVertex *p2, MVertex *p3)
{
- SVector3 a(p1->x()-p2->x(),p1->y()-p2->y(),p1->z()-p2->z());
- SVector3 b(p3->x()-p2->x(),p3->y()-p2->y(),p3->z()-p2->z());
- SVector3 c = crossprod(a,b);
+ SVector3 a(p1->x() - p2->x(), p1->y() - p2->y(), p1->z() - p2->z());
+ SVector3 b(p3->x() - p2->x(), p3->y() - p2->y(), p3->z() - p2->z());
+ SVector3 c = crossprod(a, b);
double sinA = c.norm();
- double cosA = dot(a,b);
+ double cosA = dot(a, b);
// printf("%d %d %d -> %g %g\n",p1->iD,p2->iD,p3->iD,cosA,sinA);
- return atan2 (sinA,cosA);
+ return atan2 (sinA, cosA);
}
-/*
-A curvilinear edge smooth and swap
-*/
+// A curvilinear edge smooth and swap
typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
@@ -517,47 +514,7 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
e2t[p].push_back(t);
}
}
- /*
- v2t_cont :: iterator it = adjv.begin();
- while (it != adjv.end()){
- MVertex *ver= it->first;
- GEntity *ge = ver->onWhat();
- if (ge->dim() == 2){
- double initu,initv;
- ver->getParameter(0, initu);
- ver->getParameter(1, initv);
-
- smoothVertexDataHON vdN;
- vdN.ts = it->second;
- for (int i=0;i<vdN.ts.size();i++){
- MTriangle *t = vdN.ts[i];
- }
- vdN.v = e;
- vdN.gf = gf;
-
- double val;
- double F = -smooth_obj_HighOrder(initu,initv, &vd);
- if (F < .2){
- 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){
- success = true;
- ver->setParameter(0,initu);
- ver->setParameter(1,initv);
- GPoint gp = gf->point(initu,initv);
- ver->x() = gp.x();
- ver->y() = gp.y();
- ver->z() = gp.z();
- }
- }
- }
- ++it;
- }
- */
- bool success = false;
-
for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
std::pair<MVertex*, MVertex*> edge = it->first;
std::vector<MVertex*> e;
@@ -584,7 +541,7 @@ bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
}
}
- return success;
+ return true;
}
static void parametricCoordinates(MVertex *v, GFace *gf, double &uu, double &vv)
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 2dddd409d2367acc1586bb20eb3aee40ab95d757..9f2623bf81decec21248eca4f12e18226aacd137 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -142,7 +142,6 @@ void BuildSwapPattern4(SwapPattern *sc)
sc->trianguls = trgul ;
}
-
void BuildSwapPattern5(SwapPattern *sc)
{
static int trgl[][3] =
@@ -202,10 +201,10 @@ void BuildSwapPattern7(SwapPattern *sc)
sc->trianguls = trgul ;
}
-bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
- MTet4 *tet,
- int iLocalEdge,
- const gmshQualityMeasure4Tet &cr)
+bool gmshEdgeSwap(std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr)
{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
@@ -753,28 +752,29 @@ double smoothing_objective_function_3D(double X, double Y, double Z,
return -qMin;
}
-void deriv_smoothing_objective_function_3D(double X, double Y, double Z,
- double &F,
- double &dFdX, double &dFdY, double &dFdZ,
- void *data)
+void deriv_smoothing_objective_function_3D(double *XYZ, double *dF,
+ double &F, void *data)
{
smoothVertexData3D *svd = (smoothVertexData3D*)data;
MVertex *v = svd->v;
- const double LARGE = svd->LC*1.e5;
- const double SMALL = 1./LARGE;
- F = smoothing_objective_function_3D(X, Y, Z, v, svd->ts);
- double F_X = smoothing_objective_function_3D(X + SMALL, Y, Z, v, svd->ts);
- double F_Y = smoothing_objective_function_3D(X, Y + SMALL, Z, v, svd->ts);
- double F_Z = smoothing_objective_function_3D(X, Y, Z + SMALL, v, svd->ts);
- dFdX = (F_X - F) * LARGE;
- dFdY = (F_Y - F) * LARGE;
- dFdZ = (F_Z - F) * LARGE;
+ const double LARGE = svd->LC * 1.e5;
+ const double SMALL = 1. / LARGE;
+ F = smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2], v, svd->ts);
+ double F_X = smoothing_objective_function_3D
+ (XYZ[0] + SMALL, XYZ[1], XYZ[2], v, svd->ts);
+ double F_Y = smoothing_objective_function_3D
+ (XYZ[0], XYZ[1] + SMALL, XYZ[2], v, svd->ts);
+ double F_Z = smoothing_objective_function_3D
+ (XYZ[0], XYZ[1], XYZ[2] + SMALL, v, svd->ts);
+ dF[0] = (F_X - F) * LARGE;
+ dF[1] = (F_Y - F) * LARGE;
+ dF[2] = (F_Z - F) * LARGE;
}
-double smooth_obj_3D(double X, double Y, double Z, void *data)
+double smooth_obj_3D(double *XYZ, void *data)
{
smoothVertexData3D *svd = (smoothVertexData3D*)data;
- return smoothing_objective_function_3D(X, Y, Z, svd->v, svd->ts);
+ return smoothing_objective_function_3D(XYZ[0], XYZ[1], XYZ[2], svd->v, svd->ts);
}
bool gmshSmoothVertexOptimize(MTet4 *t,
@@ -789,13 +789,11 @@ bool gmshSmoothVertexOptimize(MTet4 *t,
vd.LC = 1.0; // WRONG
gmshBuildVertexCavity_recur(t, t->tet()->getVertex(iVertex), vd.ts);
- double xopti = vd.v->x();
- double yopti = vd.v->y();
- double zopti = vd.v->z();
+ double xyzopti[3] = {vd.v->x(), vd.v->y(), vd.v->z()};
double val;
- minimize_3(smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
- xopti, yopti, zopti, val);
+ minimize_N(3, smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
+ xyzopti, val);
double vTot = 0;
@@ -810,9 +808,9 @@ bool gmshSmoothVertexOptimize(MTet4 *t,
double y = t->tet()->getVertex(iVertex)->y();
double z = t->tet()->getVertex(iVertex)->z();
- t->tet()->getVertex(iVertex)->x() = xopti;
- t->tet()->getVertex(iVertex)->y() = yopti;
- t->tet()->getVertex(iVertex)->z() = zopti;
+ t->tet()->getVertex(iVertex)->x() = xyzopti[0];
+ t->tet()->getVertex(iVertex)->y() = xyzopti[1];
+ t->tet()->getVertex(iVertex)->z() = xyzopti[2];
double newQuals[2000];
if(vd.ts.size() >= 2000){
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index dc72f506634934aa03a6f2dfe4f41bcc4467a320..5a359830021e7fe5319e9b409676dc92264ffebb 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -12,26 +12,22 @@
#if defined(HAVE_BLAS)
extern "C" {
- void dgemm_(const char *transa, const char *transb,
- int *l, int *n, int *m, double *alpha,
- const void *a, int *lda, void *b, int *ldb,
- double *beta, void *c, int *ldc);
- void dgemv_(const char *trans, int *m, int *n, double *alpha,
- void *a, int *lda, void *x, int *incx,
- double *beta, void *y, int *incy);
+ void dgemm_(const char *transa, const char *transb, int *m, int *n, int *k,
+ double *alpha, double *a, int *lda, double *b, int *ldb,
+ double *beta, double *c, int *ldc);
+ void dgemv_(const char *trans, int *m, int *n, double *alpha, double *a,
+ int *lda, double *x, int *incx, double *beta, double *y, int *incy);
}
#endif
#if defined(HAVE_LAPACK)
extern "C" {
- void dgesv_(const int *N, const int *nrhs, double *A, const int *lda,
- int *ipiv, double *b, const int *ldb, int *info);
- void dgetrf_(const int *M, const int *N, double *A, const int *lda,
- int *ipiv, int *info);
- void dgesvd_(const char* jobu, const char *jobvt, const int *M, const int *N,
- double *A, const int *lda, double *S, double* U, const int *ldu,
- double *VT, const int *ldvt, double *work, const int *lwork,
- const int *info);
+ void dgesv_(int *N, int *nrhs, double *A, int *lda, int *ipiv,
+ double *b, int *ldb, int *info);
+ void dgetrf_(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
+ void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
+ double *A, int *lda, double *S, double* U, int *ldu,
+ double *VT, int *ldvt, double *work, int *lwork, int *info);
}
#endif
@@ -132,11 +128,12 @@ class Gmsh_Matrix
}
inline void mult(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c)
{
-#if 0 // defined(HAVE_BLAS)
- void dgemm_(const char *transa, const char *transb,
- int *l, int *n, int *m, double *alpha,
- const void *a, int *lda, void *b, int *ldb,
- double *beta, void *c, int *ldc);
+#if defined(HAVE_BLAS)
+ int M = c.size1(), N = c.size2(), K = _c;
+ int LDA = _r, LDB = b.size1(), LDC = c.size1();
+ double alpha = 1., beta = 0.;
+ dgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB,
+ &beta, c._data, &LDC);
#else
c.scale(0.);
for(int i = 0; i < _r; i++)
@@ -145,19 +142,19 @@ class Gmsh_Matrix
c._data[i + _r * j] += (*this)(i, k) * b(k, j);
#endif
}
- inline void blas_dgemm(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c,
- const SCALAR alpha=1., const SCALAR beta=1.)
+ inline void blas_dgemm(Gmsh_Matrix<SCALAR> &a, Gmsh_Matrix<SCALAR> &b,
+ SCALAR alpha=1., SCALAR beta=1.)
{
-#if 0 // defined(HAVE_BLAS)
- void dgemm_(const char *transa, const char *transb,
- int *l, int *n, int *m, double *alpha,
- const void *a, int *lda, void *b, int *ldb,
- double *beta, void *c, int *ldc);
+#if defined(HAVE_BLAS)
+ int M = size1(), N = size2(), K = a.size2();
+ int LDA = a.size1(), LDB = b.size1(), LDC = size1();
+ dgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB,
+ &beta, _data, &LDC);
#else
- Gmsh_Matrix<SCALAR> temp(b.size1(), c.size2());
- temp.mult(b, c);
- scale(beta);
+ Gmsh_Matrix<SCALAR> temp(a.size1(), b.size2());
+ a.mult(b, temp);
temp.scale(alpha);
+ scale(beta);
add(temp);
#endif
}
@@ -184,10 +181,11 @@ class Gmsh_Matrix
}
inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &y)
{
-#if 0 //defined(HAVE_BLAS)
- void dgemv_(const char *trans, int *m, int *n, double *alpha,
- void *a, int *lda, void *x, int *incx,
- double *beta, void *y, int *incy);
+#if defined(HAVE_BLAS)
+ int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
+ double alpha = 1., beta = 0.;
+ dgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
+ &beta, y._data, &INCY);
#else
y.scale(0.);
for(int i = 0; i < _r; i++)
@@ -362,10 +360,10 @@ class GSL_Matrix
{
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., _data, b._data, 0., c._data);
}
- inline void blas_dgemm(const GSL_Matrix &x, GSL_Matrix &b,
- const double alpha=1., const double beta=1.)
+ inline void blas_dgemm(GSL_Matrix &a, GSL_Matrix &b,
+ double alpha=1., double beta=1.)
{
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, x._data, b._data, beta, _data);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, a._data, b._data, beta, _data);
}
inline void set_all(const double &m)
{
@@ -384,9 +382,9 @@ class GSL_Matrix
{
gsl_matrix_add(_data, m._data);
}
- inline void mult(const GSL_Vector &x, GSL_Vector &b)
+ inline void mult(const GSL_Vector &x, GSL_Vector &y)
{
- gsl_blas_dgemv(CblasNoTrans, 1., _data, x._data, 0., b._data);
+ gsl_blas_dgemv(CblasNoTrans, 1., _data, x._data, 0., y._data);
}
inline GSL_Matrix transpose()
{
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 35b4256978912cfdfdfb54727ddc0d2ebcaa02f0..b6f63a77f10823e3cd49280fe46809feca377c80 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -26,7 +26,7 @@ void my_gsl_msg(const char *reason, const char *file, int line,
Msg::Error("GSL: %s (%s, line %d)", reason, file, line);
}
-int check_gsl()
+int Init_Numeric()
{
// set new error handler
gsl_set_error_handler(&my_gsl_msg);
@@ -143,7 +143,7 @@ double brent(double ax, double bx, double cx,
b = gsl_min_fminimizer_x_minimum(s);
if(status != GSL_SUCCESS)
- Msg::Error("minimization did not converge: f(%g) = %g", b, fn1(b, NULL));
+ Msg::Error("Minimization did not converge: f(%g) = %g", b, fn1(b, NULL));
*xmin = b;
gsl_min_fminimizer_free(s);
@@ -235,194 +235,29 @@ void newt(double x[], int n, int *check,
gsl_vector_free(xx);
}
-static double (*f_stat) (double, double, void *data);
-static void (*df_stat) (double, double, double &, double &, double &, void *data);
-static double (*f_stat3) (double, double, double, void *data);
-static void (*df_stat3) (double, double, double, double &, double &,
- double &, double &,void *data);
static double (*f_statN) (double*, void *data);
-static void (*df_statN) (double*, double *, double &,void *data);
-static double fobj(const gsl_vector * x, void * data)
-{
- double u, v;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- return f_stat(u, v, data);
-}
-
-static double fobj3(const gsl_vector * x, void * data)
-{
- double u, v,w;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- return f_stat3(u, v, w, data);
-}
+static void (*df_statN) (double*, double*, double &, void *data);
-static double fobjN (const gsl_vector * x, void * data)
+static double fobjN (const gsl_vector *x, void *data)
{
return f_statN(x->data, data);
}
-static void dfobj (const gsl_vector * x, void * params, gsl_vector * g)
-{
- double u, v, f, duf, dvf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- df_stat(u, v, f, duf, dvf, params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
-}
-
-static void dfobj3 (const gsl_vector * x, void * params, gsl_vector * g)
-{
- double u, v, w,f, duf, dvf,dwf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- df_stat3(u, v, w, f, duf, dvf, dwf, params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
- gsl_vector_set(g, 2, dwf);
-}
-
-static void dfobjN(const gsl_vector * x, void * params, gsl_vector * g)
+static void dfobjN(const gsl_vector *x, void *params, gsl_vector *g)
{
double f;
- df_statN(x->data, g->data,f,params);
+ df_statN(x->data, g->data, f, params);
}
-static void fdfobj(const gsl_vector * x, void * params, double * f, gsl_vector * g)
+static void fdfobjN(const gsl_vector *x, void *params, double *f, gsl_vector *g)
{
- double u, v, duf, dvf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- df_stat(u, v, *f, duf, dvf, params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
+ df_statN(x->data, g->data, *f, params);
}
-static void fdfobj3(const gsl_vector * x, void * params, double * f, gsl_vector * g)
-{
- double u, v, w,duf, dvf,dwf;
- u = gsl_vector_get(x, 0);
- v = gsl_vector_get(x, 1);
- w = gsl_vector_get(x, 2);
- df_stat3(u, v, w, *f, duf, dvf, dwf, params);
- gsl_vector_set(g, 0, duf);
- gsl_vector_set(g, 1, dvf);
- gsl_vector_set(g, 2, dwf);
-}
-
-static void fdfobjN(const gsl_vector * x, void * params, double * f, gsl_vector * g)
-{
- df_statN(x->data,g->data,*f,params);
-}
-
-void minimize_2(double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &res)
-{
- f_stat = f;
- df_stat = df;
-
- int iter = 0;
- int status;
-
- const gsl_multimin_fdfminimizer_type *T;
- gsl_multimin_fdfminimizer *s;
- gsl_vector *x;
- gsl_multimin_function_fdf my_func;
-
- my_func.f = &fobj;
- my_func.df = &dfobj;
- my_func.fdf = &fdfobj;
- my_func.n = 2;
- my_func.params = data;
-
- x = gsl_vector_alloc (2);
- gsl_vector_set (x, 0, U);
- gsl_vector_set (x, 1, V);
-
- T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc (T, 2);
-
- gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
-
- do{
- iter++;
- status = gsl_multimin_fdfminimizer_iterate(s);
- if(status)
- break;
- status = gsl_multimin_test_gradient(s->gradient, 1e-3);
- }
- while (status == GSL_CONTINUE && iter < niter);
-
- U = gsl_vector_get(s->x, 0);
- V = gsl_vector_get(s->x, 1);
- res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
-}
-
-void minimize_3(double (*f) (double, double, double,void *data),
- void (*df) (double , double , double ,
- double &, double &, double &, double &,
- void *data) ,
- void *data,int niter,
- double &U, double &V, double &W, double &res)
-{
- f_stat3 = f;
- df_stat3 = df;
-
- int iter = 0;
- int status;
-
- const gsl_multimin_fdfminimizer_type *T;
- gsl_multimin_fdfminimizer *s;
- gsl_vector *x;
- gsl_multimin_function_fdf my_func;
-
- my_func.f = &fobj3;
- my_func.df = &dfobj3;
- my_func.fdf = &fdfobj3;
- my_func.n = 3;
- my_func.params = data;
-
- x = gsl_vector_alloc (3);
- gsl_vector_set(x, 0, U);
- gsl_vector_set(x, 1, V);
- gsl_vector_set(x, 2, W);
-
- T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc(T, 3);
-
- gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
-
- do{
- iter++;
- status = gsl_multimin_fdfminimizer_iterate (s);
- if (status)
- break;
- status = gsl_multimin_test_gradient (s->gradient, 1e-3);
- }
- while (status == GSL_CONTINUE && iter < niter);
-
- U = gsl_vector_get (s->x, 0);
- V = gsl_vector_get (s->x, 1);
- W = gsl_vector_get (s->x, 2);
- res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
-}
-
-void minimize_N(int N,
- double (*f) (double*,void *data),
- void (*df) (double*,double*,double &,void *data) ,
- void *data,int niter,
- double *U, double &res)
+void minimize_N(int N, double (*f)(double*, void *data),
+ void (*df)(double*, double*, double &, void *data),
+ void *data, int niter, double *U, double &res)
{
f_statN = f;
df_statN = df;
@@ -441,38 +276,37 @@ void minimize_N(int N,
my_func.n = N;
my_func.params = data;
- x = gsl_vector_alloc (N);
- for (int i=0;i<N;i++)
- gsl_vector_set (x, i, U[i]);
+ x = gsl_vector_alloc(N);
+ for (int i = 0; i < N; i++)
+ gsl_vector_set(x, i, U[i]);
T = gsl_multimin_fdfminimizer_conjugate_fr;
- s = gsl_multimin_fdfminimizer_alloc (T, N);
+ s = gsl_multimin_fdfminimizer_alloc(T, N);
- gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
+ gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
do{
iter++;
- status = gsl_multimin_fdfminimizer_iterate (s);
+ status = gsl_multimin_fdfminimizer_iterate(s);
if (status)
break;
- status = gsl_multimin_test_gradient (s->gradient, 1e-3);
+ status = gsl_multimin_test_gradient(s->gradient, 1e-3);
}
while (status == GSL_CONTINUE && iter < niter);
- for (int i=0;i<N;i++)
- U[i] = gsl_vector_get (s->x, i);
+ for (int i = 0; i < N; i++)
+ U[i] = gsl_vector_get(s->x, i);
res = s->f;
- gsl_multimin_fdfminimizer_free (s);
- gsl_vector_free (x);
+ gsl_multimin_fdfminimizer_free(s);
+ gsl_vector_free(x);
}
#else
-int check_gsl()
+int Init_Numeric()
{
- // initilize robust geometric predicates
gmsh::exactinit() ;
return 1;
}
@@ -489,23 +323,6 @@ void newt(double x[], int n, int *check,
Msg::Error("Gmsh must be compiled with GSL support for newt");
}
-void minimize_2(double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &res)
-{
- Msg::Error("Gmsh must be compiled with GSL support for minimize_2");
-}
-
-void minimize_3(double (*f) (double, double, double, void *data),
- void (*df) (double, double, double , double &, double &,
- double &, double &, void *data) ,
- void *data,int niter,
- double &U, double &V, double &W, double &res)
-{
- Msg::Error("Gmsh must be compiled with GSL support for minimize_3");
-}
-
void minimize_N(int N, double (*f) (double*, void *data),
void (*df) (double*, double*, double &, void *data) ,
void *data,int niter,
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 4dc536ba36853c274213c25ed8d5be8ac83736c9..7d3b29e374305f4eb4795decc701d0f21242bb7d 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -8,20 +8,11 @@
#include "NumericEmbedded.h"
-int check_gsl();
+int Init_Numeric();
double brent(double ax, double bx, double cx,
double (*f) (double), double tol, double *xmin);
void newt(double x[], int n, int *check,
void (*func) (int, double[], double[]));
-void minimize_2(double (*f) (double, double, void *data),
- void (*df) (double, double, double &, double &, double &, void *data),
- void *data,int niter,
- double &U, double &V, double &res);
-void minimize_3(double (*f) (double, double, double, void *data),
- void (*df) (double, double, double , double &, double &,
- double &, double &, void *data),
- void *data,int niter,
- double &U, double &V, double &W, double &res);
void minimize_N(int N, double (*f) (double*, void *data),
void (*df) (double*, double*, double &, void *data),
void *data,int niter,
diff --git a/configure b/configure
index 6b661cfdcce59e7edc2e31e6f2fd104e7356a323..7735f0be9b289ccaa119cb2429a3ad659fc2674b 100755
--- a/configure
+++ b/configure
@@ -576,7 +576,7 @@ PACKAGE_VERSION=
PACKAGE_STRING=
PACKAGE_BUGREPORT=
-ac_unique_file="Parser/Gmsh.y"
+ac_unique_file="Geo/GModel.h"
# Factoring default headers for most tests.
ac_includes_default="\
#include <stdio.h>
@@ -5254,8 +5254,6 @@ fi
echo "${ECHO_T}$ac_cv_lib_fftw3_main" >&6; }
if test $ac_cv_lib_fftw3_main = yes; then
FFTW3="yes"
-else
- FFTW3="no"
fi
if test "x${FFTW3}" != "xyes"; then
diff --git a/configure.in b/configure.in
index 3f7eb290f4ffcc5b60e19ed8c04ad2104f6456c5..8e7e30ffe543b9424ba219b6a342c609946b39af 100644
--- a/configure.in
+++ b/configure.in
@@ -5,12 +5,8 @@ dnl bugs and problems to <gmsh@geuz.org>.
dnl Process this file with autoconf to produce the configure script.
-dnl Note: each CHECK_LIB adds the found library to LIBS, which serves
-dnl implicitly when checking for the next one. So one should check
-dnl in reverse order! (start with -lm, ...)
-
dnl Check that this is the gmsh source tree
-AC_INIT(Parser/Gmsh.y)
+AC_INIT(Geo/GModel.h)
dnl Parse '--with' command-line options
AC_ARG_WITH(fltk-prefix,
@@ -679,7 +675,7 @@ if test "x$enable_fm" != "xno"; then
if test "x${FFTW3_PREFIX}" != "x"; then
LDFLAGS="-L${FFTW3_PREFIX}/lib ${LDFLAGS}"
fi
- AC_CHECK_LIB(fftw3,main,FFTW3="yes",FFTW3="no")
+ AC_CHECK_LIB(fftw3,main,FFTW3="yes")
if test "x${FFTW3}" != "xyes"; then
FM=no
AC_MSG_WARN([Could not find FFTW3: disabling FourierModel.])