diff --git a/Geo/GenericEdge.cpp b/Geo/GenericEdge.cpp
index e191af49bb50b2bd9c07662ee99ccb36c006ca2c..595460e06f043977641b5ec66728fdb028dbd637 100644
--- a/Geo/GenericEdge.cpp
+++ b/Geo/GenericEdge.cpp
@@ -1,7 +1,9 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
#include <limits>
#include "GmshConfig.h"
@@ -11,7 +13,6 @@
#include "GenericFace.h"
#include "Context.h"
-
GenericEdge::ptrfunction_int_double_refvector GenericEdge::EdgeEvalXYZFromT = NULL;
GenericEdge::ptrfunction_int_refdouble_refdouble GenericEdge::EdgeEvalParBounds = NULL;
GenericEdge::ptrfunction_int_refstring GenericEdge::EdgeGeomType = NULL;
@@ -22,39 +23,39 @@ GenericEdge::ptrfunction_int_refvector_refdouble_refvector_refbool GenericEdge::
GenericEdge::ptrfunction_int_refbool GenericEdge::EdgeIs3D = NULL;
GenericEdge::ptrfunction_int_int_double_int_refvector GenericEdge::EdgeReparamOnFace = NULL;
-
-GenericEdge::GenericEdge(GModel *m, int num, int _native_id, GVertex *v1, GVertex *v2, bool _isseam):GEdge(m, num, v1, v2), id(_native_id),is_seam(_isseam){
+GenericEdge::GenericEdge(GModel *m, int num, int _native_id, GVertex *v1, GVertex *v2, bool _isseam)
+ : GEdge(m, num, v1, v2), id(_native_id),is_seam(_isseam)
+{
if ((!EdgeEvalParBounds)||(!EdgeEvalXYZFromT)) Msg::Error("GenericEdge::ERROR: Callback not set");
bool ok = EdgeEvalParBounds(id,s0,s1);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalParBounds ! " );
}
-
-GenericEdge::~GenericEdge(){
+GenericEdge::~GenericEdge()
+{
}
-
-Range<double> GenericEdge::parBounds(int i) const{
+Range<double> GenericEdge::parBounds(int i) const
+{
return Range<double>(s0, s1);
}
-
-SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const{
- vector<double> res(2,0.);
+SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const
+{
+ std::vector<double> res(2,0.);
if (!EdgeReparamOnFace) Msg::Error("GenericEdge::ERROR: Callback EdgeReparamOnFace not set");
bool ok = EdgeReparamOnFace(id,face->getNativeInt(),par, dir,res);
if (!ok){
- cout << "GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id " << getNativeInt() << endl;
Msg::Error("GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id %d",getNativeInt() );
}
return SPoint2(res[0],res[1]);;
}
-
-GPoint GenericEdge::closestPoint(const SPoint3 &qp, double ¶m) const{
- vector<double> queryPoint(3,0.);
+GPoint GenericEdge::closestPoint(const SPoint3 &qp, double ¶m) const
+{
+ std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i];
- vector<double> res(3,0.);
+ std::vector<double> res(3,0.);
bool is_on_edge;
if (!EdgeClosestPoint) Msg::Error("GenericEdge::ERROR: Callback EdgeClosestPoint not set");
bool ok = EdgeClosestPoint(id,queryPoint,param,res,is_on_edge);
@@ -63,33 +64,32 @@ GPoint GenericEdge::closestPoint(const SPoint3 &qp, double ¶m) const{
return GPoint(res[0], res[1], res[2], this, param);
}
-
-bool GenericEdge::isSeam(const GFace *face) const{
-// return false;
+bool GenericEdge::isSeam(const GFace *face) const
+{
return is_seam;
}
-
-GPoint GenericEdge::point(double par) const{
- vector<double> res(3,0.);
+GPoint GenericEdge::point(double par) const
+{
+ std::vector<double> res(3,0.);
if (!EdgeEvalXYZFromT) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalXYZFromT not set");
bool ok = EdgeEvalXYZFromT(id,par,res);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalXYZFromT ! " );
return GPoint(res[0], res[1], res[2], this, par);
}
-
-SVector3 GenericEdge::firstDer(double par) const{
- vector<double> res(3,0.);
+SVector3 GenericEdge::firstDer(double par) const
+{
+ std::vector<double> res(3,0.);
if (!EdgeEvalFirstDer) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalFirstDer not set");
bool ok = EdgeEvalFirstDer(id,par,res);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalFirstDer ! " );
return SVector3(res[0],res[1],res[2]);
}
-
-GEntity::GeomType GenericEdge::geomType() const{
- string s;
+GEntity::GeomType GenericEdge::geomType() const
+{
+ std::string s;
if (!EdgeGeomType) Msg::Error("GenericEdge::ERROR: Callback EdgeGeomType not set");
bool ok = EdgeGeomType(id,s);
if (!ok){
@@ -121,8 +121,8 @@ GEntity::GeomType GenericEdge::geomType() const{
return Unknown;
}
-
-double GenericEdge::curvature(double par) const{
+double GenericEdge::curvature(double par) const
+{
double res;
if (!EdgeEvalCurvature) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalCurvature not set");
bool ok = EdgeEvalCurvature(id,par,res);
@@ -130,8 +130,8 @@ double GenericEdge::curvature(double par) const{
return res;
}
-
-bool GenericEdge::is3D() const{
+bool GenericEdge::is3D() const
+{
bool res;
if (!EdgeIs3D) Msg::Error("GenericEdge::ERROR: Callback EdgeIs3D not set");
bool ok = EdgeIs3D(id,res);
@@ -139,8 +139,8 @@ bool GenericEdge::is3D() const{
return res;
}
-
-bool GenericEdge::degenerate(int) const{
+bool GenericEdge::degenerate(int) const
+{
bool res=false;
if (!EdgeDegenerated) Msg::Error("GenericEdge::ERROR: Callback EdgeDegenerated not set");
bool ok = EdgeDegenerated(id,res);
@@ -148,7 +148,6 @@ bool GenericEdge::degenerate(int) const{
return res;
}
-
int GenericEdge::minimumDrawSegments() const
{
if(geomType() == Line)
@@ -157,63 +156,55 @@ int GenericEdge::minimumDrawSegments() const
return CTX::instance()->geom.numSubEdges * GEdge::minimumDrawSegments();
}
-
-
-
-
-
-
-LinearSeamEdge::LinearSeamEdge(GModel *m, int num, GVertex *v1, GVertex *v2):GEdge(m, num, v1, v2){
+LinearSeamEdge::LinearSeamEdge(GModel *m, int num, GVertex *v1, GVertex *v2):GEdge(m, num, v1, v2)
+{
s0=0.;
s1=v1->xyz().distance(v2->xyz());
first_der = SVector3(v1->xyz(),v2->xyz());
first_der.normalize();
}
-
-LinearSeamEdge::~LinearSeamEdge(){
+LinearSeamEdge::~LinearSeamEdge()
+{
}
-
-Range<double> LinearSeamEdge::parBounds(int i) const{
+Range<double> LinearSeamEdge::parBounds(int i) const
+{
return Range<double>(s0, s1);
}
-
-GPoint LinearSeamEdge::point(double par) const{
+GPoint LinearSeamEdge::point(double par) const
+{
SVector3 res = v0->xyz() + par*first_der;
return GPoint(res[0], res[1], res[2], this, par);
}
-
-SVector3 LinearSeamEdge::firstDer(double par) const{
+SVector3 LinearSeamEdge::firstDer(double par) const
+{
return first_der;
}
-
-GEntity::GeomType LinearSeamEdge::geomType() const{
+GEntity::GeomType LinearSeamEdge::geomType() const
+{
return Line;
}
-
-double LinearSeamEdge::curvature(double par) const{
+double LinearSeamEdge::curvature(double par) const
+{
return 0.;
}
-
-bool LinearSeamEdge::is3D() const{
+bool LinearSeamEdge::is3D() const
+{
return false;
}
-
-bool LinearSeamEdge::degenerate(int) const{
+bool LinearSeamEdge::degenerate(int) const
+{
return false;
}
-
GPoint LinearSeamEdge::closestPoint(const SPoint3 &q, double &t) const
{
return GEdge::closestPoint(q,t);
}
-
-
diff --git a/Geo/GenericEdge.h b/Geo/GenericEdge.h
index e2ffd71b8ec95cffc7d81b7c8f9b6a2022067abf..cdd5a759da0192017844ae086313cee77a08069f 100644
--- a/Geo/GenericEdge.h
+++ b/Geo/GenericEdge.h
@@ -1,7 +1,9 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_EDGE_H
#define _GENERIC_EDGE_H
@@ -13,109 +15,97 @@
#include "Range.h"
#include <vector>
-using namespace std;
-
class GenericFace;
-/*The set of Generic Entities is a generic interface to any other modeler.
- Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
+/* The set of Generic Entities is a generic interface to any other modeler.
+ Callbacks (function pointers) are given, sending requests, enquiries, to the
+ native modeler. */
class GenericEdge : public GEdge {
- protected:
- double s0, s1;
- int id;
- const bool is_seam;
- public:
- // callbacks typedef
- typedef bool (*ptrfunction_int_double_refvector)(int, double, vector<double>&);
- typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&);
- typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&);
- typedef bool (*ptrfunction_int_refstring)(int, string&);
- typedef bool (*ptrfunction_int_refbool)(int, bool&);
- typedef bool (*ptrfunction_int_refvector_refdouble_refvector_refbool)(int, const vector<double> &, double &, vector<double>&, bool &);
- typedef bool (*ptrfunction_int_int_double_int_refvector)(int, int, double, int, vector<double> &);
-
- GenericEdge(GModel *model, int num,int _native_id, GVertex *v1, GVertex *v2, bool _isseam=false);
- virtual ~GenericEdge();
-
- virtual Range<double> parBounds(int i) const;
- virtual GeomType geomType() const;
- virtual bool degenerate(int) const;
- virtual GPoint point(double p) const;
- virtual SVector3 firstDer(double par) const;
- virtual double curvature (double par) const;
- virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
- virtual GPoint closestPoint(const SPoint3 &queryPoint, double ¶m) const;
-
- ModelType getNativeType() const { return GenericModel; }
- virtual int getNativeInt()const{return id;};
-
- virtual int minimumDrawSegments () const;// for output
-
- bool is3D() const;
- bool isSeam(const GFace *) const;
-
- // sets the callbacks, to be given by the user
- static void setEdgeEvalXYZFromT(ptrfunction_int_double_refvector fct){EdgeEvalXYZFromT = fct;};
- static void setEdgeEvalParBounds(ptrfunction_int_refdouble_refdouble fct){EdgeEvalParBounds = fct;};
- static void setEdgeGeomType(ptrfunction_int_refstring fct){EdgeGeomType = fct;};
- static void setEdgeDegenerated(ptrfunction_int_refbool fct){EdgeDegenerated = fct;};
- static void setEdgeEvalFirstDer(ptrfunction_int_double_refvector fct){EdgeEvalFirstDer = fct;};
- static void setEdgeEvalCurvature(ptrfunction_int_double_refdouble fct){EdgeEvalCurvature = fct;};
- static void setEdgeClosestPoint(ptrfunction_int_refvector_refdouble_refvector_refbool fct){EdgeClosestPoint = fct;};
- static void setEdgeIs3D(ptrfunction_int_refbool fct){EdgeIs3D = fct;};
- static void setEdgeReparamOnFace(ptrfunction_int_int_double_int_refvector fct){EdgeReparamOnFace = fct;};
-
- private:
- // the callbacks:
- // --------------
- static ptrfunction_int_double_refvector EdgeEvalXYZFromT;
- static ptrfunction_int_refdouble_refdouble EdgeEvalParBounds;
- static ptrfunction_int_refstring EdgeGeomType;
- static ptrfunction_int_refbool EdgeDegenerated;
- static ptrfunction_int_double_refvector EdgeEvalFirstDer;
- static ptrfunction_int_double_refdouble EdgeEvalCurvature;
- // the first vector is a query point xyz, fills the second vector with closest point
- // on edge using orthogonal projection. Fills double param with parametric coordinate of end point projection.
- static ptrfunction_int_refvector_refdouble_refvector_refbool EdgeClosestPoint;
- static ptrfunction_int_refbool EdgeIs3D;
- static ptrfunction_int_int_double_int_refvector EdgeReparamOnFace;
+protected:
+ double s0, s1;
+ int id;
+ const bool is_seam;
+public:
+ // callbacks typedef
+ typedef bool (*ptrfunction_int_double_refvector)(int, double, std::vector<double>&);
+ typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&);
+ typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&);
+ typedef bool (*ptrfunction_int_refstring)(int, std::string&);
+ typedef bool (*ptrfunction_int_refbool)(int, bool&);
+ typedef bool (*ptrfunction_int_refvector_refdouble_refvector_refbool)(int, const std::vector<double> &, double &, std::vector<double>&, bool &);
+ typedef bool (*ptrfunction_int_int_double_int_refvector)(int, int, double, int, std::vector<double> &);
+
+ GenericEdge(GModel *model, int num,int _native_id, GVertex *v1, GVertex *v2, bool _isseam=false);
+ virtual ~GenericEdge();
+
+ virtual Range<double> parBounds(int i) const;
+ virtual GeomType geomType() const;
+ virtual bool degenerate(int) const;
+ virtual GPoint point(double p) const;
+ virtual SVector3 firstDer(double par) const;
+ virtual double curvature (double par) const;
+ virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
+ virtual GPoint closestPoint(const SPoint3 &queryPoint, double ¶m) const;
+
+ ModelType getNativeType() const { return GenericModel; }
+ virtual int getNativeInt()const{return id;};
+
+ virtual int minimumDrawSegments () const;// for output
+
+ bool is3D() const;
+ bool isSeam(const GFace *) const;
+
+ // sets the callbacks, to be given by the user
+ static void setEdgeEvalXYZFromT(ptrfunction_int_double_refvector fct){EdgeEvalXYZFromT = fct;};
+ static void setEdgeEvalParBounds(ptrfunction_int_refdouble_refdouble fct){EdgeEvalParBounds = fct;};
+ static void setEdgeGeomType(ptrfunction_int_refstring fct){EdgeGeomType = fct;};
+ static void setEdgeDegenerated(ptrfunction_int_refbool fct){EdgeDegenerated = fct;};
+ static void setEdgeEvalFirstDer(ptrfunction_int_double_refvector fct){EdgeEvalFirstDer = fct;};
+ static void setEdgeEvalCurvature(ptrfunction_int_double_refdouble fct){EdgeEvalCurvature = fct;};
+ static void setEdgeClosestPoint(ptrfunction_int_refvector_refdouble_refvector_refbool fct){EdgeClosestPoint = fct;};
+ static void setEdgeIs3D(ptrfunction_int_refbool fct){EdgeIs3D = fct;};
+ static void setEdgeReparamOnFace(ptrfunction_int_int_double_int_refvector fct){EdgeReparamOnFace = fct;};
+
+private:
+ // the callbacks:
+ // --------------
+ static ptrfunction_int_double_refvector EdgeEvalXYZFromT;
+ static ptrfunction_int_refdouble_refdouble EdgeEvalParBounds;
+ static ptrfunction_int_refstring EdgeGeomType;
+ static ptrfunction_int_refbool EdgeDegenerated;
+ static ptrfunction_int_double_refvector EdgeEvalFirstDer;
+ static ptrfunction_int_double_refdouble EdgeEvalCurvature;
+ // the first vector is a query point xyz, fills the second vector with closest point
+ // on edge using orthogonal projection. Fills double param with parametric coordinate of end point projection.
+ static ptrfunction_int_refvector_refdouble_refvector_refbool EdgeClosestPoint;
+ static ptrfunction_int_refbool EdgeIs3D;
+ static ptrfunction_int_int_double_int_refvector EdgeReparamOnFace;
};
-
-
-
-
-
-
-
-
-
-
class LinearSeamEdge : public GEdge {
- protected:
- double s0, s1;
- SVector3 first_der;
- public:
- LinearSeamEdge(GModel *model, int num, GVertex *v1, GVertex *v2);
- virtual ~LinearSeamEdge();
-
- virtual Range<double> parBounds(int i) const;
- virtual GeomType geomType() const;
- virtual bool degenerate(int) const;
- virtual GPoint point(double p) const;
- virtual SVector3 firstDer(double par) const;
- virtual double curvature (double par) const;
- virtual bool isSeam(const GFace *face) const { return true; }
- bool is3D() const;
- virtual GPoint closestPoint(const SPoint3 &queryPoint, double ¶m) const;
-
- ModelType getNativeType() const { return GenericModel; }
+protected:
+ double s0, s1;
+ SVector3 first_der;
+public:
+ LinearSeamEdge(GModel *model, int num, GVertex *v1, GVertex *v2);
+ virtual ~LinearSeamEdge();
+
+ virtual Range<double> parBounds(int i) const;
+ virtual GeomType geomType() const;
+ virtual bool degenerate(int) const;
+ virtual GPoint point(double p) const;
+ virtual SVector3 firstDer(double par) const;
+ virtual double curvature (double par) const;
+ virtual bool isSeam(const GFace *face) const { return true; }
+ bool is3D() const;
+ virtual GPoint closestPoint(const SPoint3 &queryPoint, double ¶m) const;
+
+ ModelType getNativeType() const { return GenericModel; }
};
-
#endif
diff --git a/Geo/GenericFace.cpp b/Geo/GenericFace.cpp
index 59df180c893d35e51ed0dc9726e259ec48ef504a..b2e89ba0a411a81df8800d85b76229d859d92d94 100644
--- a/Geo/GenericFace.cpp
+++ b/Geo/GenericFace.cpp
@@ -1,9 +1,11 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
-#include "GmshConfig.h"
+#include <math.h>
#include "GmshMessage.h"
#include "GModel.h"
#include "GEdgeLoop.h"
@@ -12,9 +14,6 @@
#include "GenericFace.h"
#include "Numeric.h"
#include "Context.h"
-#include <math.h>
-
-
GenericFace::ptrfunction_int_refstring GenericFace::FaceGeomType = NULL;
GenericFace::ptrfunction_int_refvector_refvector GenericFace::FaceUVFromXYZ = NULL;
@@ -29,31 +28,31 @@ GenericFace::ptrfunction_int_refvector_refvector_refvector GenericFace::FaceFirs
GenericFace::ptrfunction_int_refvector_refvector_refvector_refvector GenericFace::FaceSecondDer = NULL;
GenericFace::ptrfunction_int_refbool_refbool_refdouble_refdouble GenericFace::FacePeriodicInfo = NULL;
-
-GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_native_id){
+GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_native_id)
+{
if (!FaceParBounds) Msg::Fatal("Genericface::ERROR: Callback FaceParBounds not set");
bool ok = FaceParBounds(id,0,umin,umax);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalParBounds ! " );
ok = FaceParBounds(id,1,vmin,vmax);
-
+
computePeriodicity();
}
-
-GenericFace::~GenericFace(){
+GenericFace::~GenericFace()
+{
}
-
-Range<double> GenericFace::parBounds(int i) const{
+Range<double> GenericFace::parBounds(int i) const
+{
if(i == 0) return Range<double>(umin, umax);
return Range<double>(vmin, vmax);
}
-
-SVector3 GenericFace::normal(const SPoint2 ¶m) const{
- vector<double> res(3,0.);
- vector<double> par(2,0.);
+SVector3 GenericFace::normal(const SPoint2 ¶m) const
+{
+ std::vector<double> res(3,0.);
+ std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceEvalNormal) Msg::Fatal("Genericface::ERROR: Callback FaceEvalNormal not set");
bool ok = FaceEvalNormal(id,par,res);
@@ -61,12 +60,12 @@ SVector3 GenericFace::normal(const SPoint2 ¶m) const{
return SVector3(res[0],res[1],res[2]);
}
-
-Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 ¶m) const{
+Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 ¶m) const
+{
if (!FaceFirstDer) Msg::Fatal("Genericface::ERROR: Callback FaceFirstDer not set");
- vector<double> deru(3,0.);
- vector<double> derv(3,0.);
- vector<double> par(2,0.);
+ std::vector<double> deru(3,0.);
+ std::vector<double> derv(3,0.);
+ std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i];
bool ok = FaceFirstDer(id,par,deru,derv);
if (!ok) Msg::Error("GenericFace::ERROR from FaceFirstDer ! " );
@@ -74,12 +73,12 @@ Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 ¶m) const{
SVector3(derv[0],derv[1],derv[2]));
}
-
-void GenericFace::secondDer(const SPoint2 ¶m,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{
- vector<double> deruu(3,0.);
- vector<double> dervv(3,0.);
- vector<double> deruv(3,0.);
- vector<double> par(2,0.);
+void GenericFace::secondDer(const SPoint2 ¶m,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ std::vector<double> deruu(3,0.);
+ std::vector<double> dervv(3,0.);
+ std::vector<double> deruv(3,0.);
+ std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceSecondDer) Msg::Fatal("Genericface::ERROR: Callback FaceSecondDer not set");
bool ok = FaceSecondDer(id,par,deruu,dervv,deruv);
@@ -90,12 +89,12 @@ void GenericFace::secondDer(const SPoint2 ¶m,SVector3 *dudu, SVector3 *dvdv,
return;
}
-
-GPoint GenericFace::point(double par1, double par2) const{
- vector<double> uv(2,0.);
+GPoint GenericFace::point(double par1, double par2) const
+{
+ std::vector<double> uv(2,0.);
uv[0] = par1;
uv[1] = par2;
- vector<double> xyz(3,0.);
+ std::vector<double> xyz(3,0.);
double pp[2] = {par1, par2};
if (!FaceXYZFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceXYZFromUV not set");
bool ok = FaceXYZFromUV(id,uv,xyz);
@@ -109,10 +108,11 @@ GPoint GenericFace::point(double par1, double par2) const{
}
-GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const{
- vector<double> uvres(2,0.);
- vector<double> xyzres(3,0.);
- vector<double> queryPoint(3,0.);
+GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const
+{
+ std::vector<double> uvres(2,0.);
+ std::vector<double> xyzres(3,0.);
+ std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i];
if (!FaceClosestPoint) Msg::Fatal("Genericface::ERROR: Callback FaceClosestPoint not set");
bool ok = FaceClosestPoint(id,queryPoint,xyzres,uvres);// orthogonal projection
@@ -121,10 +121,6 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]
// check if the projected point lies on the face...
bool on_the_face;
- // ok = FaceContainsPointFromXYZ(id,xyzres,on_the_face);// check if the projected point lies on the face...
- // if (!ok) cout << "ERROR from FaceContainsPointFromXYZ ! " << endl;
- // if (!on_the_face) cout << "GenericFace::closestPoint::WARNING (using XYZ) !!!! The returned point does not lies on the face ! " << endl;
-
if (!FaceContainsPointFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromUV not set");
ok = FaceContainsPointFromUV(id,uvres,on_the_face);
if (!ok) Msg::Error("GenericFace::ERROR from FaceContainsPointFromUV ! " );
@@ -133,16 +129,16 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]
return GPoint(xyzres[0], xyzres[1], xyzres[2], this, pp);
}
-
-SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{
- vector<double> uvres(2,0.);
- vector<double> xyzres(3,0.);
- vector<double> queryPoint(3,0.);
+SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
+{
+ std::vector<double> uvres(2,0.);
+ std::vector<double> xyzres(3,0.);
+ std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i];
bool ok=true;
if (onSurface){
if (!FaceUVFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceUVFromXYZ not set");
- ok = FaceUVFromXYZ(id,queryPoint,uvres);// assuming point is on surface
+ ok = FaceUVFromXYZ(id,queryPoint,uvres);// assuming point is on surface
if (!ok) Msg::Error("GenericFace::ERROR from FaceUVFromXYZ ! " );
}
if ((!onSurface)||(!ok)){// if not on surface
@@ -158,9 +154,9 @@ SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{
return SPoint2(uvres[0],uvres[1]);
}
-
-GEntity::GeomType GenericFace::geomType() const{
- string s;
+GEntity::GeomType GenericFace::geomType() const
+{
+ std::string s;
if (!FaceGeomType) Msg::Fatal("Genericface::ERROR: Callback FaceGeomType not set");
bool ok = FaceGeomType(id,s);
if (!ok){
@@ -188,12 +184,12 @@ GEntity::GeomType GenericFace::geomType() const{
return Unknown;
}
-
-double GenericFace::curvatureMax(const SPoint2 ¶m) const{
- vector<double> dirMax(3,0.);
- vector<double> dirMin(3,0.);
+double GenericFace::curvatureMax(const SPoint2 ¶m) const
+{
+ std::vector<double> dirMax(3,0.);
+ std::vector<double> dirMin(3,0.);
double curvMax,curvMin;
- vector<double> par(2,0.);
+ std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set");
bool ok = FaceCurvatures(id,par,dirMax,dirMin,curvMax,curvMin);
@@ -201,12 +197,13 @@ double GenericFace::curvatureMax(const SPoint2 ¶m) const{
return std::max(fabs(curvMax), fabs(curvMin));
}
-
-double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,double *curvMax,double *curvMin) const{
- vector<double> param(2,0.);
+double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,
+ double *curvMax,double *curvMin) const
+{
+ std::vector<double> param(2,0.);
for (int i=0;i<3;i++) param[i] = _param[i];
- vector<double> dirMax(3,0.);
- vector<double> dirMin(3,0.);
+ std::vector<double> dirMax(3,0.);
+ std::vector<double> dirMin(3,0.);
if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set");
bool ok = FaceCurvatures(id,param,dirMax,dirMin,*curvMax,*curvMin);
@@ -217,10 +214,10 @@ double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3
return *curvMax;
}
-
-bool GenericFace::containsPoint(const SPoint3 &pt) const{
+bool GenericFace::containsPoint(const SPoint3 &pt) const
+{
bool res;
- vector<double> queryPoint(3,0.);
+ std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = pt[i];
if (!FaceContainsPointFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromXYZ not set");
bool ok = FaceContainsPointFromXYZ(id,queryPoint,res);
@@ -228,59 +225,42 @@ bool GenericFace::containsPoint(const SPoint3 &pt) const{
return res;
}
-
-void GenericFace::computePeriodicity(){
+void GenericFace::computePeriodicity()
+{
if (!FacePeriodicInfo) Msg::Fatal("Genericface::ERROR: Callback FacePeriodicInfo not set");
bool ok = FacePeriodicInfo(id,_periodic[0], _periodic[1],_period[0],_period[1]);
if (!ok) Msg::Error("GenericFace::ERROR from FacePeriodicInfo ! " );
-
-
}
-
-void GenericFace::createLoops(){
- const bool debug = false;
-
- if (debug){
- cout << "Initial loops:";
- for (std::list<GEdge *>::iterator it = l_edges.begin();it!=l_edges.end();it++){
- cout << (*it)->tag() << " ";
- }
- cout << endl;
- }
-
-
+void GenericFace::createLoops()
+{
edgeLoops.clear();
l_edges.clear();
l_dirs.clear();
- if (debug) cout << "Creating loops, face " << tag() << " (" << getNativeInt() << ")" << endl;
- for (set<int>::iterator it_loop = loopsnumber.begin();it_loop!=loopsnumber.end();it_loop++){// for each loop
- if (debug) cout << "Loop #" << *it_loop << " made of edges ";
- pair<multimap<int, pair<GEdge*,int> >::iterator, multimap<int, pair<GEdge*,int> >::iterator> range = bnd.equal_range(*it_loop);
- list<GEdge*> l_wire;
- for (multimap<int, pair<GEdge*,int> >::iterator it = range.first;it!=range.second;it++){// for all edges
- if (debug) cout << it->second.first->tag() << " ";
- l_wire.push_back(it->second.first);
+ for (std::set<int>::iterator it_loop = loopsnumber.begin();
+ it_loop!=loopsnumber.end();it_loop++){// for each loop
+ std::pair<std::multimap<int, std::pair<GEdge*,int> >::iterator,
+ std::multimap<int, std::pair<GEdge*,int> >::iterator>
+ range = bnd.equal_range(*it_loop);
+ std::list<GEdge*> l_wire;
+ for (std::multimap<int, std::pair<GEdge*,int> >::iterator it = range.first;
+ it!=range.second;it++){// for all edges
+ l_wire.push_back(it->second.first);
}
- if (debug) cout << endl << "-> GEdgeLoop made of edges (sign) ";
GEdgeLoop el(l_wire);
for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){
l_edges.push_back(it->ge);
l_dirs.push_back(it->_sign);
- if (debug) cout << it->ge->tag() << "(" << ((it->_sign<0.) ? "-" : "+") << ") " ;
if (el.count() == 2){
- it->ge->meshAttributes.minimumMeshSegments =
+ it->ge->meshAttributes.minimumMeshSegments =
std::max(it->ge->meshAttributes.minimumMeshSegments,2);
}
if (el.count() == 1){
- it->ge->meshAttributes.minimumMeshSegments =
+ it->ge->meshAttributes.minimumMeshSegments =
std::max(it->ge->meshAttributes.minimumMeshSegments,3);
}
}
- if (debug) cout << endl;
- edgeLoops.push_back(el);
+ edgeLoops.push_back(el);
}
}
-
-
diff --git a/Geo/GenericFace.h b/Geo/GenericFace.h
index 395b982f208fcc6ca2c6c4019af69651d78f4dba..3aec81819aa76208f3ba5db6f7474950ff1e8220 100644
--- a/Geo/GenericFace.h
+++ b/Geo/GenericFace.h
@@ -1,7 +1,9 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_FACE_H_
#define _GENERIC_FACE_H_
@@ -14,90 +16,89 @@
#include "Range.h"
#include <vector>
-using namespace std;
-
-/*The set of Generic Entities is a generic interface to any other modeler.
- Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
+/* The set of Generic Entities is a generic interface to any other modeler.
+ Callbacks (function pointers) are given, sending requests, enquiries, to the
+ native modeler. */
class GenericFace : public GFace {
- protected:
- int id;
- double umin, umax, vmin, vmax;// face uv bounds
- bool _periodic[2];// is periodic in u, v
- double _period[2];
-
- public:
- // callbacks typedef
- typedef bool (*ptrfunction_int_refstring)(int, string&);
- typedef bool (*ptrfunction_int_refbool_refbool_refdouble_refdouble)(int, bool&, bool&, double&, double&);
- typedef bool (*ptrfunction_int_refvector_refvector)(const int , const vector<double> &, vector<double> &);
- typedef bool (*ptrfunction_int_refvector_refvector_refvector)(const int , const vector<double> &, vector<double> &, vector<double> &);
- typedef bool (*ptrfunction_int_refvector_refbool)(const int , const vector<double> &, bool &);
- typedef bool (*ptrfunction_int_int_refdouble_refdouble)(const int, const int, double &, double &);
- typedef bool (*ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble)(const int, const vector<double> &, vector<double> &, vector<double> &, double &, double &);
- typedef bool (*ptrfunction_int_refvector_refvector_refvector_refvector)(const int , const vector<double> &, vector<double> &, vector<double> &, vector<double> &);
-
- GenericFace(GModel *m, int num, int _native_id);
- virtual ~GenericFace();
-
- Range<double> parBounds(int i) const;
- virtual GPoint point(double par1, double par2) const;
- virtual GPoint closestPoint(const SPoint3 & queryPoint,
- const double initialGuess[2]) const;
- virtual bool containsPoint(const SPoint3 &pt) const;
- virtual SVector3 normal(const SPoint2 ¶m) const;
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
- virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
- virtual GEntity::GeomType geomType() const;
- virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
- virtual double curvatureMax(const SPoint2 ¶m) const;
- virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
- double *curvMax, double *curvMin) const;
- virtual bool periodic(int dir) const { return false;}//_periodic[dir]; }// TODO ?
- virtual double period(int dir) const{ return _period[dir]; };
-
- ModelType getNativeType() const { return GenericModel; }
- virtual int getNativeInt()const{return id;};
-
- void addBndInfo(int loop, GEdge *ptr, int sign){
- bnd.insert(make_pair(loop, make_pair(ptr,sign)));
- l_dirs.push_back(sign);
- l_edges.push_back(ptr);
- loopsnumber.insert(loop);
- ptr->addFace(this);
- };
- void createLoops();
-
- void computePeriodicity();
-
- // sets callbacks
- static void setFaceGeomType(ptrfunction_int_refstring fct){FaceGeomType = fct;};
- static void setFaceUVFromXYZ(ptrfunction_int_refvector_refvector fct){FaceUVFromXYZ = fct;};
- static void setFaceClosestPoint(ptrfunction_int_refvector_refvector_refvector fct){FaceClosestPoint = fct;};
- static void setFaceContainsPointFromXYZ(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromXYZ = fct;};
- static void setFaceContainsPointFromUV(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromUV = fct;};
- static void setFaceXYZFromUV(ptrfunction_int_refvector_refvector fct){FaceXYZFromUV = fct;};
- static void setFaceParBounds(ptrfunction_int_int_refdouble_refdouble fct){FaceParBounds = fct;};
- static void setFaceCurvatures(ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble fct){FaceCurvatures = fct;};
- static void setFaceEvalNormal(ptrfunction_int_refvector_refvector fct){FaceEvalNormal = fct;};
- static void setFaceFirstDer(ptrfunction_int_refvector_refvector_refvector fct){FaceFirstDer = fct;};
- static void setFaceSecondDer(ptrfunction_int_refvector_refvector_refvector_refvector fct){FaceSecondDer = fct;};
- static void setFacePeriodicInfo(ptrfunction_int_refbool_refbool_refdouble_refdouble fct){FacePeriodicInfo = fct;};
-
- private:
- multimap<int, pair<GEdge*,int> > bnd;
- set<int> loopsnumber;
-
- // the callbacks:
- static ptrfunction_int_refstring FaceGeomType;
- static ptrfunction_int_refvector_refvector FaceUVFromXYZ,FaceXYZFromUV,FaceEvalNormal;
- static ptrfunction_int_refvector_refvector_refvector FaceClosestPoint,FaceFirstDer;
- static ptrfunction_int_refvector_refbool FaceContainsPointFromXYZ,FaceContainsPointFromUV;
- static ptrfunction_int_int_refdouble_refdouble FaceParBounds;
- static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures;
- static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer;
- static ptrfunction_int_refbool_refbool_refdouble_refdouble FacePeriodicInfo;
-
+protected:
+ int id;
+ double umin, umax, vmin, vmax;// face uv bounds
+ bool _periodic[2];// is periodic in u, v
+ double _period[2];
+
+public:
+ // callbacks typedef
+ typedef bool (*ptrfunction_int_refstring)(int, std::string&);
+ typedef bool (*ptrfunction_int_refbool_refbool_refdouble_refdouble)(int, bool&, bool&, double&, double&);
+ typedef bool (*ptrfunction_int_refvector_refvector)(const int , const std::vector<double> &, std::vector<double> &);
+ typedef bool (*ptrfunction_int_refvector_refvector_refvector)(const int , const std::vector<double> &, std::vector<double> &, std::vector<double> &);
+ typedef bool (*ptrfunction_int_refvector_refbool)(const int , const std::vector<double> &, bool &);
+ typedef bool (*ptrfunction_int_int_refdouble_refdouble)(const int, const int, double &, double &);
+ typedef bool (*ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble)(const int, const std::vector<double> &, std::vector<double> &, std::vector<double> &, double &, double &);
+ typedef bool (*ptrfunction_int_refvector_refvector_refvector_refvector)(const int , const std::vector<double> &, std::vector<double> &, std::vector<double> &, std::vector<double> &);
+
+ GenericFace(GModel *m, int num, int _native_id);
+ virtual ~GenericFace();
+
+ Range<double> parBounds(int i) const;
+ virtual GPoint point(double par1, double par2) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint,
+ const double initialGuess[2]) const;
+ virtual bool containsPoint(const SPoint3 &pt) const;
+ virtual SVector3 normal(const SPoint2 ¶m) const;
+ virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
+ virtual GEntity::GeomType geomType() const;
+ virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
+ virtual double curvatureMax(const SPoint2 ¶m) const;
+ virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const;
+ virtual bool periodic(int dir) const { return false;}//_periodic[dir]; }// TODO ?
+ virtual double period(int dir) const{ return _period[dir]; };
+
+ ModelType getNativeType() const { return GenericModel; }
+ virtual int getNativeInt()const{return id;};
+
+ void addBndInfo(int loop, GEdge *ptr, int sign)
+ {
+ bnd.insert(std::make_pair(loop, std::make_pair(ptr,sign)));
+ l_dirs.push_back(sign);
+ l_edges.push_back(ptr);
+ loopsnumber.insert(loop);
+ ptr->addFace(this);
+ };
+ void createLoops();
+
+ void computePeriodicity();
+
+ // sets callbacks
+ static void setFaceGeomType(ptrfunction_int_refstring fct){FaceGeomType = fct;};
+ static void setFaceUVFromXYZ(ptrfunction_int_refvector_refvector fct){FaceUVFromXYZ = fct;};
+ static void setFaceClosestPoint(ptrfunction_int_refvector_refvector_refvector fct){FaceClosestPoint = fct;};
+ static void setFaceContainsPointFromXYZ(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromXYZ = fct;};
+ static void setFaceContainsPointFromUV(ptrfunction_int_refvector_refbool fct){FaceContainsPointFromUV = fct;};
+ static void setFaceXYZFromUV(ptrfunction_int_refvector_refvector fct){FaceXYZFromUV = fct;};
+ static void setFaceParBounds(ptrfunction_int_int_refdouble_refdouble fct){FaceParBounds = fct;};
+ static void setFaceCurvatures(ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble fct){FaceCurvatures = fct;};
+ static void setFaceEvalNormal(ptrfunction_int_refvector_refvector fct){FaceEvalNormal = fct;};
+ static void setFaceFirstDer(ptrfunction_int_refvector_refvector_refvector fct){FaceFirstDer = fct;};
+ static void setFaceSecondDer(ptrfunction_int_refvector_refvector_refvector_refvector fct){FaceSecondDer = fct;};
+ static void setFacePeriodicInfo(ptrfunction_int_refbool_refbool_refdouble_refdouble fct){FacePeriodicInfo = fct;};
+
+private:
+ std::multimap<int, std::pair<GEdge*,int> > bnd;
+ std::set<int> loopsnumber;
+
+ // the callbacks:
+ static ptrfunction_int_refstring FaceGeomType;
+ static ptrfunction_int_refvector_refvector FaceUVFromXYZ,FaceXYZFromUV,FaceEvalNormal;
+ static ptrfunction_int_refvector_refvector_refvector FaceClosestPoint,FaceFirstDer;
+ static ptrfunction_int_refvector_refbool FaceContainsPointFromXYZ,FaceContainsPointFromUV;
+ static ptrfunction_int_int_refdouble_refdouble FaceParBounds;
+ static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures;
+ static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer;
+ static ptrfunction_int_refbool_refbool_refdouble_refdouble FacePeriodicInfo;
};
#endif
diff --git a/Geo/GenericRegion.cpp b/Geo/GenericRegion.cpp
index 4031c11d9bd4affec9c9143075fad4c29c6b7585..e8325d2f98108c902ac207292e12a243686f4c0d 100644
--- a/Geo/GenericRegion.cpp
+++ b/Geo/GenericRegion.cpp
@@ -1,27 +1,25 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
-#include "GmshConfig.h"
-#include "GmshMessage.h"
#include "GModel.h"
#include "GenericVertex.h"
#include "GenericEdge.h"
#include "GenericFace.h"
#include "GenericRegion.h"
-
-GenericRegion::GenericRegion(GModel *m, int num, int _native_id):GRegion(m, num), id(_native_id)
+GenericRegion::GenericRegion(GModel *m, int num, int _native_id)
+ : GRegion(m, num), id(_native_id)
{
}
-
GenericRegion::~GenericRegion()
{
}
-
GEntity::GeomType GenericRegion::geomType() const
{
return Unknown;
diff --git a/Geo/GenericRegion.h b/Geo/GenericRegion.h
index 74850cc8b2e7cb3535f48a167b283391adc97e48..1ae038d3feb3a21eb6294131188f0ca2ef77a65d 100644
--- a/Geo/GenericRegion.h
+++ b/Geo/GenericRegion.h
@@ -1,7 +1,9 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_REGION_H_
#define _GENERIC_REGION_H_
@@ -9,29 +11,31 @@
#include "GmshConfig.h"
#include "GRegion.h"
-/*The set of Generic Entities is a generic interface to any other modeler.
- Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
+/* The set of Generic Entities is a generic interface to any other modeler.
+ Callbacks (function pointers) are given, sending requests, enquiries, to the
+ native modeler. */
class GenericRegion : public GRegion {
- public:
- GenericRegion(GModel *m, int num, int _native_id);
- virtual ~GenericRegion();
-
- virtual GeomType geomType() const;
-
- ModelType getNativeType() const { return GenericModel; }
- virtual int getNativeInt()const{return id;};
-
- // TODO: When using GRegion->l_dirs and l_faces, what is the convention for l_dirs ? For now, assuming positive value for normals pointing inside the region.
- void addFace(GenericFace *ptr, int sign){
- l_dirs.push_back(sign);
- l_faces.push_back(ptr);
- ptr->addRegion(this);
- };
-
-
- private:
- int id;
+public:
+ GenericRegion(GModel *m, int num, int _native_id);
+ virtual ~GenericRegion();
+
+ virtual GeomType geomType() const;
+
+ ModelType getNativeType() const { return GenericModel; }
+ virtual int getNativeInt()const{return id;};
+
+ // TODO: When using GRegion->l_dirs and l_faces, what is the convention for
+ // l_dirs ? For now, assuming positive value for normals pointing inside the
+ // region.
+ void addFace(GenericFace *ptr, int sign){
+ l_dirs.push_back(sign);
+ l_faces.push_back(ptr);
+ ptr->addRegion(this);
+ };
+
+private:
+ int id;
};
#endif
diff --git a/Geo/GenericVertex.cpp b/Geo/GenericVertex.cpp
index 49500ce09a5479b46df9f3576c880bdd6d3ee019..03a4320338d57074da2236ffa68247de576398e8 100644
--- a/Geo/GenericVertex.cpp
+++ b/Geo/GenericVertex.cpp
@@ -1,27 +1,28 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
-#include "GmshConfig.h"
+#include <algorithm>
#include "GModel.h"
#include "MVertex.h"
#include "MPoint.h"
#include "GenericVertex.h"
#include "GenericEdge.h"
#include "GenericFace.h"
-#include<algorithm>
-
GenericVertex::ptrfunction_int_vector GenericVertex::VertexXYZ = NULL;
GenericVertex::ptrfunction_int_doubleptr_voidptr GenericVertex::VertexMeshSize = NULL;
-
-GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num), id(_native_id){
+GenericVertex::GenericVertex(GModel *m, int num, int _native_id)
+ : GVertex(m, num), id(_native_id)
+{
if (!VertexXYZ)
Msg::Fatal("GenericVertex::ERROR: Callback not set");
- vector<double> vec(3,0.);
+ std::vector<double> vec(3,0.);
bool ok = VertexXYZ(id,vec);
if (!ok) Msg::Error("GenericVertex::ERROR from callback VertexXYZ ");
_x=vec[0];
@@ -29,30 +30,26 @@ GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num)
_z=vec[2];
}
-
-GenericVertex::GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec):GVertex(m, num), id(_native_id){
+GenericVertex::GenericVertex(GModel *m, int num, int _native_id, const std::vector<double> &vec)
+ : GVertex(m, num), id(_native_id)
+{
if (!VertexXYZ)
Msg::Fatal("GenericVertex::ERROR: Callback not set");
-
_x=vec[0];
_y=vec[1];
_z=vec[2];
}
-
GenericVertex::~GenericVertex(){
-}
+}
SPoint2 GenericVertex::reparamOnFace(const GFace *gf, int dir) const
{
-
SPoint3 pt(_x,_y,_z);
return gf->parFromPoint(pt,true);
-
}
-
void GenericVertex::setPosition(GPoint &p)
{
_x = p.x();
diff --git a/Geo/GenericVertex.h b/Geo/GenericVertex.h
index 46f5c60f4343d1fd840311f4acdeaf9446f153e3..0bee21fdb419e0d4b2bf77eaf77e12bd47988450 100644
--- a/Geo/GenericVertex.h
+++ b/Geo/GenericVertex.h
@@ -1,7 +1,9 @@
-// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+//
+// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_VERTEX_H_
#define _GENERIC_VERTEX_H_
@@ -11,64 +13,61 @@
#include "GVertex.h"
#include "Context.h"
#include <vector>
-#include <iostream>
-
-using namespace std;
-/*The set of Generic Entities is a generic interface to any other modeler.
- Callbacks (function pointers) are given, sending requests, enquiries, to the native modeler. */
+/* The set of Generic Entities is a generic interface to any other modeler.
+ Callbacks (function pointers) are given, sending requests, enquiries, to the
+ native modeler. */
class GenericVertex : public GVertex {
- protected:
- int id;
- double _x, _y, _z;
- public:
- // callbacks typedef
- typedef bool (*ptrfunction_int_vector)(int, vector<double>&);
- typedef bool (*ptrfunction_int_doubleptr_voidptr)(int, double*, void*);
+protected:
+ int id;
+ double _x, _y, _z;
+public:
+ // callbacks typedef
+ typedef bool (*ptrfunction_int_vector)(int, std::vector<double>&);
+ typedef bool (*ptrfunction_int_doubleptr_voidptr)(int, double*, void*);
- GenericVertex(GModel *m, int num, int _native_id);
- GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec);
- virtual ~GenericVertex();
+ GenericVertex(GModel *m, int num, int _native_id);
+ GenericVertex(GModel *m, int num, int _native_id, const std::vector<double> &vec);
+ virtual ~GenericVertex();
- virtual GPoint point() const { return GPoint(x(), y(), z()); }
- virtual double x() const { return _x; }
- virtual double y() const { return _y; }
- virtual double z() const { return _z; }
+ virtual GPoint point() const { return GPoint(x(), y(), z()); }
+ virtual double x() const { return _x; }
+ virtual double y() const { return _y; }
+ virtual double z() const { return _z; }
- virtual void setPosition(GPoint &p);
- virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
+ virtual void setPosition(GPoint &p);
+ virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
- ModelType getNativeType() const { return GenericModel; }
- virtual int getNativeInt()const{return id;};
+ ModelType getNativeType() const { return GenericModel; }
+ virtual int getNativeInt()const{return id;};
- // sets the callbacks
- static void setVertexXYZ(ptrfunction_int_vector fct){VertexXYZ = fct;};
- static void setVertexMeshSize(ptrfunction_int_doubleptr_voidptr fct){VertexMeshSize = fct;};
+ // sets the callbacks
+ static void setVertexXYZ(ptrfunction_int_vector fct){VertexXYZ = fct;};
+ static void setVertexMeshSize(ptrfunction_int_doubleptr_voidptr fct){VertexMeshSize = fct;};
- // meshing-related methods:
- virtual void setPrescribedMeshSizeAtVertex(double l) {
- cout << "GenericVertex:: setPrescribedMeshSizeAtVertex !!!, aborting" << endl;
- throw;
+ // meshing-related methods:
+ virtual void setPrescribedMeshSizeAtVertex(double l)
+ {
+ Msg::Error("GenericVertex::setPrescribedMeshSizeAtVertex");
+ }
+ virtual inline double prescribedMeshSizeAtVertex() const
+ {
+ double size;
+ void *chose = NULL;
+ if (!VertexMeshSize(id,&size,chose)){
+ Msg::Error("GenericVertex::ERROR from callback VertexMeshSize");
+ return CTX::instance()->lc;
}
- virtual inline double prescribedMeshSizeAtVertex() const {
- double size;
- void *chose = NULL;
- if (!VertexMeshSize(id,&size,chose)){
- cout<< "GenericVertex::ERROR from callback VertexMeshSize "<< endl;
- Msg::Warning("GenericVertex::ERROR from callback VertexMeshSize ");
- cout << "Check for ring edges !!!" << endl;
- return CTX::instance()->lc;
- }
- return size;
- }
-
- private:
- // the callbacks:
- // --------------
- // fills vector xyz for vertex of int id
- static ptrfunction_int_vector VertexXYZ;
- static ptrfunction_int_doubleptr_voidptr VertexMeshSize;
+ return size;
+ }
+
+private:
+ // the callbacks:
+ // --------------
+ // fills vector xyz for vertex of int id
+ static ptrfunction_int_vector VertexXYZ;
+ static ptrfunction_int_doubleptr_voidptr VertexMeshSize;
};
#endif
diff --git a/Mesh/BackgroundMesh3D.cpp b/Mesh/BackgroundMesh3D.cpp
index c4f4e9f73d5534a82fe7864707b280370132e2b0..d2026669a89b8611493376b9be1ef7bdc94a481d 100644
--- a/Mesh/BackgroundMesh3D.cpp
+++ b/Mesh/BackgroundMesh3D.cpp
@@ -1,23 +1,16 @@
-#include "BackgroundMesh3D.h"
#include <fstream>
#include <algorithm>
#include <iostream>
#include <string>
#include <numeric>
-
-#include <time.h>
-
+#include "BackgroundMesh3D.h"
#include "MElement.h"
#include "GFace.h"
#include "GRegion.h"
#include "BackgroundMeshManager.h"
#include "BackgroundMesh2D.h"
-
#include "MElementOctree.h"
-
-//#include "pointInsertion.h"
-
#include "MTetrahedron.h"
#include "OS.h"
@@ -28,12 +21,13 @@
#include "linearSystemFull.h"
#endif
+using namespace std;
-int signof(int i){
+static int signof(int i)
+{
return ((i < 0) ? -1 : 1);
}
-
// TODO: virer les trucs "vertextorank", mettre cette classe-ci :
//class evalDiffusivity3DFunction : public simpleFunction<double>{
@@ -47,18 +41,18 @@ int signof(int i){
// const double threshold;
//};
-
-backgroundMesh3D::backgroundMesh3D(GRegion *_gr):BGMBase(3,_gr),debug(false),verbose(false){
+backgroundMesh3D::backgroundMesh3D(GRegion *_gr)
+ : BGMBase(3,_gr), debug(false), verbose(false)
+{
computeSizeField();
}
-
-backgroundMesh3D::~backgroundMesh3D(){
+backgroundMesh3D::~backgroundMesh3D()
+{
}
-
-void backgroundMesh3D::computeSizeField(){
-
+void backgroundMesh3D::computeSizeField()
+{
cout << "backgroundMesh3D::computeSizeField() " << endl;
// fills dirichlet BC
@@ -78,7 +72,7 @@ void backgroundMesh3D::computeSizeField(){
SPoint2 p;
reparamMeshVertexOnFace(v, face, p);
sizeField[v] = bgm2d->size(p.x(),p.y());
-// cout << "sets sizefield for vertex " << v << endl;
+ // cout << "sets sizefield for vertex " << v << endl;
}
}
}
@@ -87,11 +81,13 @@ void backgroundMesh3D::computeSizeField(){
simpleFunction<double> ONE(1.0);
propagateValues(sizeField,ONE);
-// cout << "backgroundMesh3D::size of sizeField: " << sizeField.size() << endl;
+ // cout << "backgroundMesh3D::size of sizeField: " << sizeField.size() << endl;
}
-
-void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane){
+void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet,
+ simpleFunction<double> &eval_diffusivity,
+ bool in_parametric_plane)
+{
// same as Size_field::solve()
GRegion *gr = dynamic_cast<GRegion*>(gf);
@@ -161,42 +157,44 @@ void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
#endif
}
-
-GPoint backgroundMesh3D::get_GPoint_from_MVertex(const MVertex *v)const{
+GPoint backgroundMesh3D::get_GPoint_from_MVertex(const MVertex *v)const
+{
return GPoint(v->x(),v->y(),v->z(),dynamic_cast<GRegion*>(gf));
}
-
-MVertex* backgroundMesh3D::get_nearest_neighbor(const MVertex *v){
+MVertex* backgroundMesh3D::get_nearest_neighbor(const MVertex *v)
+{
double distance;
return get_nearest_neighbor(v->x(),v->y(),v->z(), distance);
}
-
-MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz){
+MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz)
+{
double distance;
return get_nearest_neighbor(xyz, distance);
}
-
-MVertex* backgroundMesh3D::get_nearest_neighbor(double x, double y, double z){
+MVertex* backgroundMesh3D::get_nearest_neighbor(double x, double y, double z)
+{
double distance;
return get_nearest_neighbor(x,y,z,distance);
}
-
-MVertex* backgroundMesh3D::get_nearest_neighbor(double x, double y, double z, double &distance ){
+MVertex* backgroundMesh3D::get_nearest_neighbor(double x, double y, double z,
+ double &distance )
+{
double xyz[3] = {x,y,z};
return get_nearest_neighbor(xyz, distance);
}
-
-MElementOctree* backgroundMesh3D::getOctree(){
+MElementOctree* backgroundMesh3D::getOctree()
+{
if(!octree){
GRegion *gr = dynamic_cast<GRegion*>(gf);
Msg::Debug("Rebuilding BackgroundMesh element octree");
std::vector<MElement*> copy;// !!!
- for (std::vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();it!=gr->tetrahedra.end();it++){
+ for (std::vector<MTetrahedron*>::iterator it = gr->tetrahedra.begin();
+ it!=gr->tetrahedra.end();it++){
copy.push_back(*it);
}
octree = new MElementOctree(copy);
@@ -204,8 +202,8 @@ MElementOctree* backgroundMesh3D::getOctree(){
return octree;
}
-
-MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & distance ){
+MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & distance )
+{
// using the octree instead of ANN, faster.
MElement *elem = const_cast<MElement*>(findElement(xyz[0],xyz[1],xyz[2]));
if (!elem) return NULL;
@@ -220,7 +218,6 @@ MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & dist
vector<double>::iterator itmax = std::max_element(distances.begin(),distances.end());
return candidates[std::distance(distances.begin(),itmax)];
-
// map<double,MVertex*> distances;
// SPoint3 p(xyz[0],xyz[1],xyz[2]);
// for (int i=0;i<elem->getNumVertices();i++){
@@ -232,13 +229,10 @@ MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & dist
// return it->second;
}
-
vector<montripletbis> frameFieldBackgroundMesh3D::permutation = vector<montripletbis>();
-
-frameFieldBackgroundMesh3D::frameFieldBackgroundMesh3D(GRegion *_gr):backgroundMesh3D(_gr){
-
-
+frameFieldBackgroundMesh3D::frameFieldBackgroundMesh3D(GRegion *_gr):backgroundMesh3D(_gr)
+{
// readValue("param.dat","SMOOTHCF",smooth_the_crossfield);
smooth_the_crossfield = true;
@@ -271,16 +265,16 @@ frameFieldBackgroundMesh3D::frameFieldBackgroundMesh3D(GRegion *_gr):backgroundM
}
}
-
-frameFieldBackgroundMesh3D::~frameFieldBackgroundMesh3D(){
+frameFieldBackgroundMesh3D::~frameFieldBackgroundMesh3D()
+{
#if defined(HAVE_ANN)
if(annTreeBnd) delete annTreeBnd;
if(dataPtsBnd) annDeallocPts(dataPtsBnd);
#endif
}
-
-void frameFieldBackgroundMesh3D::initiate_ANN_research(){
+void frameFieldBackgroundMesh3D::initiate_ANN_research()
+{
#ifdef HAVE_ANN
// ANN research for 2D !!!
int maxPts = listOfBndVertices.size();
@@ -293,13 +287,13 @@ void frameFieldBackgroundMesh3D::initiate_ANN_research(){
dataPtsBnd[i][k] = (v->point())[k];
++i;
}
- annTreeBnd=new ANNkd_tree(dataPtsBnd, maxPts , 3);
+ annTreeBnd=new ANNkd_tree(dataPtsBnd, maxPts , 3);
#endif
return;
}
-
-void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){
+void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries()
+{
// this is used when the crossfield is not smoothed !!!
// otherwise, the smoothness is computed on the fly while smoothing the cf !
@@ -311,7 +305,8 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){
double mean_angle=0.;
vector<double> vectorial_smoothness(3);
- for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++){// for all vertices
+ for (vert2elemtype::iterator it_vertex=vert2elem.begin();
+ it_vertex!=vert2elem.end();it_vertex++){// for all vertices
themap.clear();
neighbors.clear();
current = it_vertex->first;
@@ -328,13 +323,15 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){
STensor3 &ref = itcurrent->second;
multimap<double,MVertex*>::iterator it_neighbors_begin = themap.begin();
- crossFieldSmoothness[current] = compare_to_neighbors(current->point(), ref, it_neighbors_begin, themap.end(), mean_axis,mean_angle,vectorial_smoothness);
+ crossFieldSmoothness[current] = compare_to_neighbors
+ (current->point(), ref, it_neighbors_begin, themap.end(), mean_axis,
+ mean_angle,vectorial_smoothness);
// crossFieldVectorialSmoothness[current] = vectorial_smoothness;
}
}
-
-void frameFieldBackgroundMesh3D::computeCrossField(){
+void frameFieldBackgroundMesh3D::computeCrossField()
+{
// TODO: mettre des parametres façon gmsh ???
// reading parameters
int Niter = 3;
@@ -347,7 +344,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
smoothness_threshold=0.85;
double percentage_threshold = 0.98;
-
vector<double> scal_prods;
vector<double> signs;
vector<double> angles;
@@ -360,7 +356,8 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
map<MVertex*,bool> vertex_is_still;
map<MVertex*,double> vertex_movement;
MVertex *current;
- for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++){
+ for (vert2elemtype::iterator it_vertex=vert2elem.begin();
+ it_vertex!=vert2elem.end();it_vertex++){
current = it_vertex->first;
if (current->onWhat()->dim()<=2)
vertex_is_still[current] = true;
@@ -371,15 +368,14 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
// OLD - NEW COMPARISON
map<MVertex*,double> vertex_to_rank;
- for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++){// for all vertices
+ for (vert2elemtype::iterator it_vertex=vert2elem.begin();
+ it_vertex!=vert2elem.end();it_vertex++){// for all vertices
//vertex_to_rank[it_vertex->first] = 0.;
vertex_to_rank[it_vertex->first] = 1.;
rank.insert(make_pair(0.,it_vertex->first));
}
// END OLD - NEW COMPARISON
-
-
double percentage_not_done=0.;
double norme = 10.*norme_threshold;
double norme_inf=10.*norme_threshold;
@@ -389,34 +385,44 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
vector<double> vectorial_smoothness(3);
// vector<int> nb_local_iterations;
-
// main smoothing loop
- while (((norme_inf>norme_threshold) && (percentage_not_done<percentage_threshold)) && (iIter<Niter)){// for maximum Niter iterations, or until convergence
+ while (((norme_inf>norme_threshold) &&
+ (percentage_not_done<percentage_threshold)) &&
+ (iIter<Niter)){// for maximum Niter iterations, or until convergence
// nb_local_iterations.clear();
angles.clear();
norme_inf=0.;
int counter_done=0;
int counter_not_done=0;
- while (rank.size()){// for all vertices, i.e. all cavities
- // for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++)// for all vertices, i.e. all cavities
- //MVertex *current = it_vertex->first;
+ while (rank.size()){// for all vertices, i.e. all cavities
+ // for (vert2elemtype::iterator
+ //it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++)//
+ //for all vertices, i.e. all cavities MVertex *current = it_vertex->first;
current = (rank.begin())->second;
rank.erase((rank.begin()));
-
// smooth 3D vertices only
if (current->onWhat()->dim()<=2){
- if (verbose) cout << "-------------------- discard point " << current->getNum() << " on dim " << current->onWhat()->dim() << " coord (" << current->x() << "," << current->y() << "," << current->z()<< ")" << endl;
+ if (verbose) cout << "-------------------- discard point "
+ << current->getNum() << " on dim "
+ << current->onWhat()->dim() << " coord ("
+ << current->x() << "," << current->y() << ","
+ << current->z()<< ")" << endl;
continue;
}
TensorStorageType::iterator itcurrent = crossField.find(current);
STensor3 &ref = itcurrent->second;
- if (verbose) cout << "-------------------- working on point " << current->getNum() << " on dim " << current->onWhat()->dim() << " coord (" << current->x() << "," << current->y() << "," << current->z()<< ")" << endl;
+ if (verbose) cout << "-------------------- working on point "
+ << current->getNum() << " on dim "
+ << current->onWhat()->dim() << " coord ("
+ << current->x() << "," << current->y() << ","
+ << current->z()<< ")" << endl;
pair<graphtype::iterator, graphtype::iterator> range = graph.equal_range(current);
graphtype::iterator itgraph = range.first;
- bool all_neighbors_still=true;// if nothing has changed since previous iteration: nothing to do !
+ bool all_neighbors_still=true;// if nothing has changed since previous
+ // iteration: nothing to do !
for (;itgraph!=range.second;itgraph++){// for all neighbors
if (vertex_is_still[itgraph->second.second]==false){
all_neighbors_still = false;
@@ -424,23 +430,24 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
}
}
- if (all_neighbors_still && vertex_is_still[current]){// neighbors didn't move, and current didn't move either -> nothing to do !
+ // neighbors didn't move, and current didn't move either -> nothing to do !
+ if (all_neighbors_still && vertex_is_still[current]){
vertex_movement[current] = 0.;
rank_temp.insert(make_pair(0.,current));
counter_not_done++;
continue;
}
-
-
// OLD - NEW COMPARISON
multimap<double,MVertex*> neighbors_to_trust;
itgraph = range.first;
for (;itgraph!=range.second;itgraph++){// for all neighbors
- neighbors_to_trust.insert(make_pair(vertex_to_rank[itgraph->second.second],itgraph->second.second));
- // if (vertex_to_rank[itgraph->second.second] <= 1. /180.*M_PI){
- // neighbors_to_trust.insert(make_pair(vertex_to_rank[itgraph->second.second],itgraph->second.second));
- // }
+ neighbors_to_trust.insert(make_pair(vertex_to_rank[itgraph->second.second],
+ itgraph->second.second));
+ // if (vertex_to_rank[itgraph->second.second] <= 1. /180.*M_PI){
+ // neighbors_to_trust.insert(make_pair(vertex_to_rank[itgraph->second.second],
+ // itgraph->second.second));
+ // }
}
if (!neighbors_to_trust.size()){
rank_temp.insert(make_pair(M_PI,current));
@@ -449,11 +456,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
}
// END OLD - NEW COMPARISON
-
-
-
-
-
counter_done++;
double angle_to_go;;
@@ -463,26 +465,27 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
for (;Nlocaliter<20;Nlocaliter++){// iterations, convergence of the local cavity...
multimap<double,MVertex*>::iterator it_neighbors_to_trust = neighbors_to_trust.begin();
- crossFieldSmoothness[current] = compare_to_neighbors(current->point(), ref, it_neighbors_to_trust, neighbors_to_trust.end(), mean_axis,mean_angle,vectorial_smoothness);
- // crossFieldVectorialSmoothness[current] = vectorial_smoothness;
+ crossFieldSmoothness[current] = compare_to_neighbors
+ (current->point(), ref, it_neighbors_to_trust, neighbors_to_trust.end(),
+ mean_axis,mean_angle,vectorial_smoothness);
+ // crossFieldVectorialSmoothness[current] = vectorial_smoothness;
// rotating directions to align closest vectors...
angle_to_go = mean_angle;
ref = apply_rotation(mean_axis,angle_to_go,ref);
- // cout << " iter " << Nlocaliter << ": mean_angle = " << mean_angle << endl;
+ // cout << " iter " << Nlocaliter << ": mean_angle = " << mean_angle <<
+ // endl;
if (fabs(mean_angle)<localangle_threshold/3.) break;
}
- // nb_local_iterations.push_back(Nlocaliter+1);
+ // nb_local_iterations.push_back(Nlocaliter+1);
if (verbose) cout << "iIter " << iIter << ", Nlocaliter = " << Nlocaliter << endl;
-
- // computing the total rotation
+ // computing the total rotation
//get_min_rotation_matrix(ref_original, ref, mean_angle, mean_axis,localangle_threshold,true);
get_min_rotation_matrix(ref_original, ref, mean_angle, mean_axis,localangle_threshold);//,true);
norme_inf = max(norme_inf,mean_angle);
// cout << " #local_iter=" << Nlocaliter << " final mean_angle = " << mean_angle << " threshold=" << localangle_threshold << " condition :" << (mean_angle <= localangle_threshold) << endl;
-
angles.push_back(mean_angle);
vertex_is_still[current] = (mean_angle <= localangle_threshold);
vertex_movement[current] = mean_angle;
@@ -502,8 +505,13 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
norme = std::accumulate(angles.begin(),angles.end(),0.);
if (angles.size()) norme = norme/M_PI*180./angles.size();
percentage_not_done = ((double)counter_not_done)/(counter_done+counter_not_done);
- cout << " --- iter " << iIter << " NormeInf=" << norme_inf << " mean Angle=" << norme << " counter_not_done=" << counter_not_done << " NotDone (%)=" << (percentage_not_done*100.) << " %" << endl;
- // cout << "Average number of local iterations: " << ((double)std::accumulate(nb_local_iterations.begin(),nb_local_iterations.end(),0))/nb_local_iterations.size() << endl;
+ cout << " --- iter " << iIter << " NormeInf=" << norme_inf
+ << " mean Angle=" << norme << " counter_not_done="
+ << counter_not_done << " NotDone (%)="
+ << (percentage_not_done*100.) << " %" << endl;
+ // cout << "Average number of local iterations: "
+ // << ((double)std::accumulate(nb_local_iterations.begin(),
+ // nb_local_iterations.end(),0))/nb_local_iterations.size() << endl;
//if (debug){
double temp = log10(Niter);
@@ -515,10 +523,9 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
iIter++;
}// end Niter iterations
-
-
// also computes smoothness for boundary points
- for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++){
+ for (vert2elemtype::iterator it_vertex=vert2elem.begin();
+ it_vertex!=vert2elem.end();it_vertex++){
current = it_vertex->first;
if (current->onWhat()->dim()<=2){
TensorStorageType::iterator itcurrent = crossField.find(current);
@@ -526,26 +533,23 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
pair<graphtype::iterator, graphtype::iterator> range = graph.equal_range(current);
graphtype::iterator itgraph = range.first;
-
multimap<double,MVertex*> neighbors_to_trust;
itgraph = range.first;
for (;itgraph!=range.second;itgraph++){// for all neighbors
neighbors_to_trust.insert(make_pair(0.,itgraph->second.second));
}
- crossFieldSmoothness[current] = compare_to_neighbors(current->point(), ref, neighbors_to_trust.begin(), neighbors_to_trust.end(), mean_axis,mean_angle,vectorial_smoothness);
+ crossFieldSmoothness[current] = compare_to_neighbors
+ (current->point(), ref, neighbors_to_trust.begin(),
+ neighbors_to_trust.end(), mean_axis,mean_angle,vectorial_smoothness);
- //crossFieldSmoothness[current] = compare_to_neighbors(current->point(), ref, itgraph, range.second, mean_axis,mean_angle,vectorial_smoothness);
+ //crossFieldSmoothness[current] = compare_to_neighbors
+ // (current->point(), ref, itgraph, range.second, mean_axis,
+ // mean_angle,vectorial_smoothness);
// crossFieldVectorialSmoothness[current] = vectorial_smoothness;
}
}
-
-
-
- return;
-
}
-
//STensor3 frameFieldBackgroundMesh3D::get_random_cross()const{
// SVector3 vec1(rand()%10000/9999. , rand()%10000/9999. , rand()%10000/9999. );
// vec1.normalize();
@@ -563,19 +567,21 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
// return random_cross;
//}
-
-void frameFieldBackgroundMesh3D::initiate_crossfield(){
+void frameFieldBackgroundMesh3D::initiate_crossfield()
+{
crossField.clear();
MVertex *v;
// first, for boundaries :
GRegion *gr = dynamic_cast<GRegion*>(gf);
- list<GFace*> faces = gr->faces();
- // here, not using the gm2D since we are interested by the new 2D vertices, not the old (now erased) ones... alternative would be to reset the 2DBGM...
+ list<GFace*> faces = gr->faces();
+ // here, not using the gm2D since we are interested by the new 2D vertices,
+ // not the old (now erased) ones... alternative would be to reset the 2DBGM...
for (list<GFace*>::iterator it=faces.begin();it!=faces.end();it++){// for all GFace
GFace *face = *it;
- frameFieldBackgroundMesh2D *bgm2d = dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face));
- // initializing the vertices on the faces
+ frameFieldBackgroundMesh2D *bgm2d =
+ dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face));
+ // initializing the vertices on the faces
for (unsigned int i=0;i<face->getNumMeshElements();i++){// for all elements
MElement *e = face->getMeshElement(i);
if (e->getDim()!=2) continue;// of dim=2
@@ -594,7 +600,7 @@ void frameFieldBackgroundMesh3D::initiate_crossfield(){
}
}
- // then, for volume :
+ // then, for volume :
for (unsigned int i=0;i<gr->getNumMeshElements();i++){// for all elements
MElement *e = gr->getMeshElement(i);
if (e->getDim()!=3) continue;
@@ -613,14 +619,15 @@ void frameFieldBackgroundMesh3D::initiate_crossfield(){
}
}
-
-MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v){
+MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v)
+{
double distance;
return get_nearest_neighbor_on_boundary(v,distance);
}
-
-MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v, double &distance){
+MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary
+ (MVertex* v, double &distance)
+{
#ifdef HAVE_ANN
ANNpoint q = annAllocPt(3);
for (int k=0; k< 3; ++k)
@@ -642,33 +649,35 @@ MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v
#endif
}
-
-double frameFieldBackgroundMesh3D::get_smoothness(double x, double y, double z){
+double frameFieldBackgroundMesh3D::get_smoothness(double x, double y, double z)
+{
return get_field_value(x,y,z,crossFieldSmoothness);
-};
-
+};
-double frameFieldBackgroundMesh3D::get_approximate_smoothness(double x, double y, double z){
+double frameFieldBackgroundMesh3D::get_approximate_smoothness(double x, double y, double z)
+{
return crossFieldSmoothness[get_nearest_neighbor(x,y,z)];
}
-
-double frameFieldBackgroundMesh3D::get_approximate_smoothness(MVertex *v){
+double frameFieldBackgroundMesh3D::get_approximate_smoothness(MVertex *v)
+{
return get_approximate_smoothness(v->x(),v->y(),v->z());
}
-
-double frameFieldBackgroundMesh3D::get_smoothness(MVertex *v){
+double frameFieldBackgroundMesh3D::get_smoothness(MVertex *v)
+{
return get_nodal_value(v,crossFieldSmoothness);
};
-
-void frameFieldBackgroundMesh3D::eval_approximate_crossfield(double x, double y, double z, STensor3 &cf){
+void frameFieldBackgroundMesh3D::eval_approximate_crossfield(double x, double y,
+ double z, STensor3 &cf)
+{
cf = crossField[get_nearest_neighbor(x,y,z)];
};
-
-void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert, STensor3 &cf){
+void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert,
+ STensor3 &cf)
+{
// check if vertex is ours...
TensorStorageType::const_iterator itfind = crossField.find(const_cast<MVertex*>(vert));
if (itfind != crossField.end())
@@ -678,31 +687,35 @@ void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert
};
-void frameFieldBackgroundMesh3D::get_vectorial_smoothness(double x, double y, double z, SVector3 &cf){
+void frameFieldBackgroundMesh3D::get_vectorial_smoothness(double x, double y, double z,
+ SVector3 &cf)
+{
vector<double>res = get_field_value(x,y,z,crossFieldVectorialSmoothness);
for (int i=0;i<3;i++) cf(i)=res[i];
};
-
-void frameFieldBackgroundMesh3D::get_vectorial_smoothness(const MVertex *vert, SVector3 &cf){
+void frameFieldBackgroundMesh3D::get_vectorial_smoothness(const MVertex *vert, SVector3 &cf)
+{
vector<double> res = get_nodal_value(vert,crossFieldVectorialSmoothness);
for (int i=0;i<3;i++) cf(i)=res[i];
};
-
-double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, double x, double y, double z){
+double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, double x,
+ double y, double z)
+{
vector<double>res = get_field_value(x,y,z,crossFieldVectorialSmoothness);
return res[idir];
};
-
-double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, const MVertex *vert){
+double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir,
+ const MVertex *vert)
+{
vector<double> res = get_nodal_value(vert,crossFieldVectorialSmoothness);
return res[idir];
};
-
-void frameFieldBackgroundMesh3D::build_vertex_to_element_table(){
+void frameFieldBackgroundMesh3D::build_vertex_to_element_table()
+{
GRegion *gr = dynamic_cast<GRegion*>(gf);
MElement *e;
MVertex *v;
@@ -726,28 +739,29 @@ void frameFieldBackgroundMesh3D::build_vertex_to_element_table(){
}
}
-
-const MElement* backgroundMesh3D::getElement(unsigned int i)const{
+const MElement* backgroundMesh3D::getElement(unsigned int i)const
+{
GRegion *gr = dynamic_cast<GRegion*>(gf);
return gr->getMeshElement(i);
}
-
-unsigned int backgroundMesh3D::getNumMeshElements()const{
+unsigned int backgroundMesh3D::getNumMeshElements()const
+{
GRegion *gr = dynamic_cast<GRegion*>(gf);
return gr->getNumMeshElements();
}
-
-// this function fills the multimap "graph": vertex to direct neighbors, or indirect neighbors, depending on the "max_recursion_level".
-void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level){
-
+// this function fills the multimap "graph": vertex to direct neighbors, or
+// indirect neighbors, depending on the "max_recursion_level".
+void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level)
+{
set<MVertex*> visited,start;
set<MElement*> visited_elements;
multimap<int,MVertex*> proximity;
//int counter=0;
- for (vert2elemtype::iterator it_vertex=vert2elem.begin();it_vertex!=vert2elem.end();it_vertex++){// for all vertices
+ for (vert2elemtype::iterator it_vertex=vert2elem.begin();
+ it_vertex!=vert2elem.end();it_vertex++){// for all vertices
MVertex *current_vertex = it_vertex->first;
visited.clear();
visited_elements.clear();
@@ -781,16 +795,20 @@ void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level)
}
}
-
-void frameFieldBackgroundMesh3D::get_recursive_neighbors(set<MVertex*> &start, set<MVertex*> &visited, set<MElement*> &visited_elements, multimap<int,MVertex*> &proximity, int max_level, int level){
+void frameFieldBackgroundMesh3D::get_recursive_neighbors
+ (set<MVertex*> &start, set<MVertex*> &visited, set<MElement*> &visited_elements,
+ multimap<int,MVertex*> &proximity, int max_level, int level)
+{
level++;
if (level>max_level) return;
set<MVertex*> new_vertices;
- for (set<MVertex*>::iterator it_start=start.begin();it_start!=start.end();it_start++){// for all initial vertices
+ for (set<MVertex*>::iterator it_start=start.begin();
+ it_start!=start.end();it_start++){// for all initial vertices
MVertex *current = *it_start;
- // cout << "get_recursive_neighbors : on vertex " << current->getNum() << " (" << current << ")" << endl;
+ // cout << "get_recursive_neighbors : on vertex " << current->getNum()
+ // << " (" << current << ")" << endl;
vert2elemtype::iterator itfind = vert2elem.find(current);
if (itfind==vert2elem.end()) continue;
set<MElement*>::iterator it_elem = itfind->second.begin();
@@ -815,12 +833,15 @@ void frameFieldBackgroundMesh3D::get_recursive_neighbors(set<MVertex*> &start, s
}
}
- get_recursive_neighbors(new_vertices, visited, visited_elements, proximity, max_level, level);
-
+ get_recursive_neighbors(new_vertices, visited, visited_elements,
+ proximity, max_level, level);
}
-
-double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor3 &ref, multimap<double,MVertex*>::iterator itbegin, multimap<double,MVertex*>::iterator itend, SVector3 &mean_axis,double &mean_angle, vector<double> &vectorial_smoothness){
+double frameFieldBackgroundMesh3D::compare_to_neighbors
+ (SPoint3 current, STensor3 &ref, multimap<double,MVertex*>::iterator itbegin,
+ multimap<double,MVertex*>::iterator itend, SVector3 &mean_axis,
+ double &mean_angle, vector<double> &vectorial_smoothness)
+{
for (int i=0;i<3;i++) mean_axis(i)=0.;
multimap<double,MVertex*>::iterator it = itbegin;
SVector3 rotation_axis;
@@ -835,14 +856,15 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor
vector<double> ponderations_vec;
- vector<double> alternative;
+ vector<double> alternative;
for (;it!=itend;it++){// for all trusted neighbors
//neighbor = it->second.second;
neighbor = it->second;
//ponderations_vec.push_back(1.);
ponderations_vec.push_back((fabs(it->first) >= smoothness_threshold) ? 1. : 1.e-3);
- //ponderations_vec.push_back(((fabs(it->first) >= smoothness_threshold) ? 1. : 1.e-3) / (current.distance(neighbor->point())));
+ //ponderations_vec.push_back(((fabs(it->first) >= smoothness_threshold) ?
+ //1. : 1.e-3) / (current.distance(neighbor->point())));
itneighbor = crossField.find(neighbor);
STensor3 &neibcross = itneighbor->second;
@@ -882,9 +904,6 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor
// double smoothness_scalar = (*std::min_element(vectorial_smoothness.begin(),vectorial_smoothness.end()));
double smoothness_scalar = 1. - (std::accumulate(alternative.begin(),alternative.end(),0.)/alternative.size())/M_PI*4.;// APPROXIMATELY between 0 (not smooth) and 1 (smooth), (sometimes <0, always > 1).
-
-
-
vector<double>::iterator itan=all_angle.begin();
vector<double>::iterator itpond=ponderations_vec.begin();
@@ -898,13 +917,12 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor
mean_angle = mean_axis.norm()/pond_total;
mean_axis.normalize();
-
-
return smoothness_scalar;
}
-
-STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const double &angle, const STensor3 &thecross){
+STensor3 frameFieldBackgroundMesh3D::apply_rotation
+ (const SVector3 &axis, const double &angle, const STensor3 &thecross)
+{
double rotmat[3][3];
STensor3 res2;
get_rotation_matrix(angle,axis,&rotmat[0][0]);
@@ -919,7 +937,6 @@ STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const
// }
//
-
for (int i=0;i<3;i++)
for (int j=0;j<3;j++)
for (int k=0;k<3;k++)
@@ -985,9 +1002,11 @@ STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const
}
-// TODO: ce genre de fct... mettre dans une classe FrameField ? Ou bien juste un set de fcts à part ? Parce que ça devient bizarre d'avoir ça dans un BGM ???
-void frameFieldBackgroundMesh3D::get_rotation_matrix(const double &angle_to_go, const SVector3 &rotvec, double* rotmat){
-
+// TODO: ce genre de fct... mettre dans une classe FrameField ? Ou bien juste un
+// set de fcts à part ? Parce que ça devient bizarre d'avoir ça dans un BGM ???
+void frameFieldBackgroundMesh3D::get_rotation_matrix
+ (const double &angle_to_go, const SVector3 &rotvec, double* rotmat)
+{
double c = cos(angle_to_go);
double s = sin(angle_to_go);
double temp01 = rotvec[0]*rotvec[1]*(1.-c);
@@ -1007,9 +1026,10 @@ void frameFieldBackgroundMesh3D::get_rotation_matrix(const double &angle_to_go,
rotmat[8] = square2 +(1.-square2)*c;
}
-
-void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, double threshold, bool debugflag){
-
+void frameFieldBackgroundMesh3D::get_min_rotation_matrix
+ (const STensor3 &reference, const STensor3 &thecross, double &minimum_angle,
+ SVector3 &rotation_axis, double threshold, bool debugflag)
+{
minimum_angle = M_PI/2.;
pair<int,int> p;
@@ -1021,7 +1041,9 @@ void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &referen
get_rotation_angle_and_axis(reference,thecross,a,v,i_rotation_perm,permutation[iperm]);
if (debugflag){
if (fabs(a)<M_PI/2.){
- cout << " temp parameters: angle=" << a*180./M_PI << "pair=(" << iperm << "," << i_rotation_perm << ") axis=(" << v(0) << "," << v(1) << "," << v(2) << ")" << endl;
+ cout << " temp parameters: angle=" << a*180./M_PI
+ << "pair=(" << iperm << "," << i_rotation_perm << ") axis=("
+ << v(0) << "," << v(1) << "," << v(2) << ")" << endl;
}
else
cout << " temp parameters: angle=" << a*180./M_PI << endl;
@@ -1032,23 +1054,26 @@ void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &referen
rotation_axis = v;
}
-
if (minimum_angle<threshold) break;
}
}
// cout << "pair=(" << p.first << "," << p.second << ")" << endl;
if (debugflag){
- cout << " ---> MIN parameters: angle=" << minimum_angle*180./M_PI << " axis=(" << rotation_axis(0) << "," << rotation_axis(1) << "," << rotation_axis(2) << ")" << endl;
+ cout << " ---> MIN parameters: angle=" << minimum_angle*180./M_PI
+ << " axis=(" << rotation_axis(0) << "," << rotation_axis(1)
+ << "," << rotation_axis(2) << ")" << endl;
}
return;
}
-// compute the angle and axis of rotation between two (unit-orthogonal) cross fields
-// using crossfield vectors permutations defined by modulo and trip
-// The rotation matrix between cross fields C and M is R = C M^T
-void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, int modulo, montripletbis &trip){
-
+// compute the angle and axis of rotation between two (unit-orthogonal) cross
+// fields using crossfield vectors permutations defined by modulo and trip The
+// rotation matrix between cross fields C and M is R = C M^T
+void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis
+ (const STensor3 &reference, const STensor3 &thecross, double &minimum_angle,
+ SVector3 &rotation_axis, int modulo, montripletbis &trip)
+{
//double MT[3][3],C[3][3],R[3][3];
double C[3][3],R[3][3];
@@ -1069,7 +1094,7 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &ref
// ////////////////////
// computing the trace for the angle...
- // MT is transpose of reference !!!
+ // MT is transpose of reference !!!
// R[0][0] = C[0][0]*MT[0][0] + C[0][1]*MT[1][0] + C[0][2]*MT[2][0];
// R[1][1] = C[1][0]*MT[0][1] + C[1][1]*MT[1][1] + C[1][2]*MT[2][1];
// R[2][2] = C[2][0]*MT[0][2] + C[2][1]*MT[1][2] + C[2][2]*MT[2][2];
@@ -1085,7 +1110,7 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &ref
if (fabs(minimum_angle)>M_PI/2.) return;// no need to continue...
- // computing rotation axis
+ // computing rotation axis
// computing the remaining terms, not on diagonal
if (fabs(minimum_angle)<1.e-6){
rotation_axis(0)=0.;
@@ -1108,8 +1133,8 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &ref
return;
}
-
-void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filename){
+void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filename)
+{
FILE *f = Fopen (filename.c_str(),"w");
fprintf(f,"View \"Background Mesh\"{\n");
@@ -1126,12 +1151,11 @@ void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filenam
eval_approximate_crossfield(v,cf);
for (int i=0;i<3;i++){
double vs = get_vectorial_smoothness(i,v);
- fprintf(f,"VP(%g,%g,%g){%g,%g,%g};\n",v->x(),v->y(),v->z(),cf(0,i)*vs,cf(1,i)*vs,cf(2,i)*vs);
+ fprintf(f,"VP(%g,%g,%g){%g,%g,%g};\n", v->x(), v->y(), v->z(),
+ cf(0,i)*vs,cf(1,i)*vs,cf(2,i)*vs);
}
}
}
fprintf(f,"};\n");
fclose(f);
}
-
-
diff --git a/Mesh/BackgroundMesh3D.h b/Mesh/BackgroundMesh3D.h
index 51226b62c381714c9ab2c54e4a4101b873927089..344aa8a7c850df2e20c8f36f9ad5d0d2a9a28f7f 100644
--- a/Mesh/BackgroundMesh3D.h
+++ b/Mesh/BackgroundMesh3D.h
@@ -14,14 +14,7 @@
#include <iostream>
#include <iomanip>
#include <numeric>
-
-
-
-//#include "SVector3.h"
-//#include "STensor3.h"
#include "BGMBase.h"
-
-
#include "pointInsertion.h"
#if defined(HAVE_ANN)
@@ -29,145 +22,140 @@
class ANNkd_tree;
#endif
+// TODO: do we really need to build the graph vertex2elem-> vertex 2 neighbors ?
+// This was useful for the smoothness computation, if we want to take into
+// account vertices at a given topological distance... but using distance=1, it
+// seems to work... ?!
-
-// TODO: do we really need to build the graph vertex2elem-> vertex 2 neighbors ? This was useful for the smoothness computation, if we want to take into account vertices at a given topological distance... but using distance=1, it seems to work... ?!
+class montripletbis{
+public:
+ ~montripletbis(){};
+ montripletbis(int a,int b, int c):vec(3){
+ vec[0] = a;
+ vec[1] = b;
+ vec[2] = c;
+ };
+ inline int operator()(int i)const{return vec[i];}
+private:
+ std::vector<int> vec;
+};
-using namespace std;
+// difference with BackgroundMesh2D: no copy of components, working directly on the vertices and elements of GRegion
+class backgroundMesh3D : public BGMBase {
+protected:
+ virtual void computeSizeField();
+ virtual void propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane = false);
+ virtual GPoint get_GPoint_from_MVertex(const MVertex *) const;
+ virtual const MElement* getElement(unsigned int i)const;
+ virtual unsigned int getNumMeshElements()const;
+ const bool debug;
+ const bool verbose;
-class montripletbis{
- public:
- ~montripletbis(){};
- montripletbis(int a,int b, int c):vec(3){
- vec[0] = a;
- vec[1] = b;
- vec[2] = c;
- };
- inline int operator()(int i)const{return vec[i];}
- private:
- vector<int> vec;
-};
+public:
+ backgroundMesh3D(GRegion *_gr);
+ virtual ~backgroundMesh3D();
-// difference with BackgroundMesh2D: no copy of components, working directly on the vertices and elements of GRegion
+ virtual MElementOctree* getOctree();
-class backgroundMesh3D : public BGMBase {
- protected:
- virtual void computeSizeField();
- virtual void propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane = false);
- virtual GPoint get_GPoint_from_MVertex(const MVertex *) const;
-
- virtual const MElement* getElement(unsigned int i)const;
- virtual unsigned int getNumMeshElements()const;
-
- const bool debug;
- const bool verbose;
-
- public:
-
- backgroundMesh3D(GRegion *_gr);
- virtual ~backgroundMesh3D();
-
- virtual MElementOctree* getOctree();
-
- virtual MVertex* get_nearest_neighbor(const double* xyz, double & distance );
- virtual MVertex* get_nearest_neighbor(const double* xyz);
- virtual MVertex* get_nearest_neighbor(const MVertex *v);
- virtual MVertex* get_nearest_neighbor(double x, double y, double z);
- virtual MVertex* get_nearest_neighbor(double x, double y, double z, double &distance );
+ virtual MVertex* get_nearest_neighbor(const double* xyz, double & distance );
+ virtual MVertex* get_nearest_neighbor(const double* xyz);
+ virtual MVertex* get_nearest_neighbor(const MVertex *v);
+ virtual MVertex* get_nearest_neighbor(double x, double y, double z);
+ virtual MVertex* get_nearest_neighbor(double x, double y, double z, double &distance );
};
class frameFieldBackgroundMesh3D : public backgroundMesh3D{
- public:
-// typedef tr1::unordered_map<hash_key_ptr, set<MElement*> > vert2elemtype;
-// typedef tr1::unordered_map<MElement*, set<MVertex*> > elem2verttype;
- typedef std::map<hash_key_ptr, set<MElement*> > vert2elemtype;
- typedef std::map<MElement*, set<MVertex*> > elem2verttype;
- typedef multimap<MVertex*, pair<int,MVertex*> > graphtype;
-
- protected:
- bool smooth_the_crossfield;
-
- void initiate_crossfield();
- void initiate_ANN_research();
- void computeCrossField();
- void computeSmoothnessOnlyFromBoundaries();
-
- void build_vertex_to_element_table();
-
- // fills the multimap "graph": vertex to direct neighbors, or indirect neighbors, depending on the "max_recursion_level".
- void build_neighbors(const int &max_recursion_level);
-
- // function called only by "build_neighbors", recursively recovering vertices neighbors.
- void get_recursive_neighbors(set<MVertex*> &start, set<MVertex*> &visited, set<MElement*> &visited_elements, multimap<int,MVertex*> &proximity, int max_level, int level=0);
-
-
- double compare_to_neighbors(SPoint3 current, STensor3 &ref, multimap<double,MVertex*>::iterator itbegin, multimap<double,MVertex*>::iterator itend, SVector3 &mean_axis,double &mean_angle, vector<double> &vectorial_smoothness);
- STensor3 apply_rotation(const SVector3 &axis, const double &angle, const STensor3 &thecross);
- void get_rotation_matrix(const double &angle_to_go, const SVector3 &rotvec, double* rotmat);
- void get_min_rotation_matrix(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, double threshold=-1., bool debugflag=false);
- void get_rotation_angle_and_axis(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, int modulo, montripletbis &trip);
-
-
-
- TensorStorageType crossField;
- DoubleStorageType crossFieldSmoothness;
- VectorStorageType crossFieldVectorialSmoothness;
- vert2elemtype vert2elem;
- elem2verttype elem2vert;
-
- set<MVertex*> listOfBndVertices;
- graphtype graph;
- double smoothness_threshold;
-
- static vector<montripletbis> permutation;
+public:
+ // typedef tr1::unordered_map<hash_key_ptr, std::set<MElement*> > vert2elemtype;
+ // typedef tr1::unordered_map<MElement*, std::set<MVertex*> > elem2verttype;
+ typedef std::map<hash_key_ptr, std::set<MElement*> > vert2elemtype;
+ typedef std::map<MElement*, std::set<MVertex*> > elem2verttype;
+ typedef std::multimap<MVertex*, std::pair<int,MVertex*> > graphtype;
+
+protected:
+ bool smooth_the_crossfield;
+
+ void initiate_crossfield();
+ void initiate_ANN_research();
+ void computeCrossField();
+ void computeSmoothnessOnlyFromBoundaries();
+
+ void build_vertex_to_element_table();
+
+ // fills the multimap "graph": vertex to direct neighbors, or indirect neighbors, depending on the "max_recursion_level".
+ void build_neighbors(const int &max_recursion_level);
+
+ // function called only by "build_neighbors", recursively recovering vertices neighbors.
+ void get_recursive_neighbors(std::set<MVertex*> &start, std::set<MVertex*> &visited, std::set<MElement*> &visited_elements, std::multimap<int,MVertex*> &proximity, int max_level, int level=0);
+
+
+ double compare_to_neighbors(SPoint3 current, STensor3 &ref, std::multimap<double,MVertex*>::iterator itbegin, std::multimap<double,MVertex*>::iterator itend, SVector3 &mean_axis,double &mean_angle, std::vector<double> &vectorial_smoothness);
+ STensor3 apply_rotation(const SVector3 &axis, const double &angle, const STensor3 &thecross);
+ void get_rotation_matrix(const double &angle_to_go, const SVector3 &rotvec, double* rotmat);
+ void get_min_rotation_matrix(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, double threshold=-1., bool debugflag=false);
+ void get_rotation_angle_and_axis(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, int modulo, montripletbis &trip);
+
+
+
+ TensorStorageType crossField;
+ DoubleStorageType crossFieldSmoothness;
+ VectorStorageType crossFieldVectorialSmoothness;
+ vert2elemtype vert2elem;
+ elem2verttype elem2vert;
+
+ std::set<MVertex*> listOfBndVertices;
+ graphtype graph;
+ double smoothness_threshold;
+
+ static std::vector<montripletbis> permutation;
#if defined(HAVE_ANN)
- ANNkd_tree* annTree;
- ANNpointArray dataPts;
- ANNkd_tree* annTreeBnd;
- ANNpointArray dataPtsBnd;
+ ANNkd_tree* annTree;
+ ANNpointArray dataPts;
+ ANNkd_tree* annTreeBnd;
+ ANNpointArray dataPtsBnd;
#endif
- public:
- frameFieldBackgroundMesh3D(GRegion *_gf);
- virtual ~frameFieldBackgroundMesh3D();
-
- MVertex* get_nearest_neighbor_on_boundary(MVertex* v);
- MVertex* get_nearest_neighbor_on_boundary(MVertex* v, double &distance);
-
-// bool findvertex(const MVertex *v)const{
-// map<MVertex*,double>::const_iterator it = sizeField.find(const_cast<MVertex*>(v));
-// if (it==sizeField.end()) return false;
-// return true;
-// };
-
- double get_smoothness(double x, double y, double z);
- double get_smoothness(MVertex *v);
-
- double get_approximate_smoothness(double x, double y, double z);
- double get_approximate_smoothness(MVertex *v);
-
- void eval_approximate_crossfield(double x, double y, double z, STensor3 &cf);
- void eval_approximate_crossfield(const MVertex *vert, STensor3 &cf);
-
- void get_vectorial_smoothness(double x, double y, double z, SVector3 &cf);
- void get_vectorial_smoothness(const MVertex *vert, SVector3 &cf);
- double get_vectorial_smoothness(const int idir, double x, double y, double z);
- double get_vectorial_smoothness(const int idir, const MVertex *vert);
-
- inline void exportCrossField(const string &filename){export_tensor_as_vectors(filename,crossField);};
- inline void exportSmoothness(const string &filename) const{export_scalar(filename,crossFieldSmoothness);};
- void exportVectorialSmoothness(const string &filename);
-
-// private:
-// STensor3 get_random_cross()const;
+public:
+ frameFieldBackgroundMesh3D(GRegion *_gf);
+ virtual ~frameFieldBackgroundMesh3D();
+
+ MVertex* get_nearest_neighbor_on_boundary(MVertex* v);
+ MVertex* get_nearest_neighbor_on_boundary(MVertex* v, double &distance);
+
+ // bool findvertex(const MVertex *v)const{
+ // std::map<MVertex*,double>::const_iterator it = sizeField.find(const_cast<MVertex*>(v));
+ // if (it==sizeField.end()) return false;
+ // return true;
+ // };
+
+ double get_smoothness(double x, double y, double z);
+ double get_smoothness(MVertex *v);
+
+ double get_approximate_smoothness(double x, double y, double z);
+ double get_approximate_smoothness(MVertex *v);
+
+ void eval_approximate_crossfield(double x, double y, double z, STensor3 &cf);
+ void eval_approximate_crossfield(const MVertex *vert, STensor3 &cf);
+
+ void get_vectorial_smoothness(double x, double y, double z, SVector3 &cf);
+ void get_vectorial_smoothness(const MVertex *vert, SVector3 &cf);
+ double get_vectorial_smoothness(const int idir, double x, double y, double z);
+ double get_vectorial_smoothness(const int idir, const MVertex *vert);
+
+ inline void exportCrossField(const std::string &filename){export_tensor_as_vectors(filename,crossField);};
+ inline void exportSmoothness(const std::string &filename) const{export_scalar(filename,crossFieldSmoothness);};
+ void exportVectorialSmoothness(const std::string &filename);
+
+ // private:
+ // STensor3 get_random_cross()const;
};
diff --git a/Mesh/pointInsertionRTreeTools.h b/Mesh/pointInsertionRTreeTools.h
index 8ac501c870e6664205c6ae20588606aec44e4ea1..42945a89ee9496482d89547b55d46ab0ac4656a5 100644
--- a/Mesh/pointInsertionRTreeTools.h
+++ b/Mesh/pointInsertionRTreeTools.h
@@ -206,7 +206,7 @@ public:
virtual bool empty(){ return points.empty(); };
protected:
- set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
+ std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
};
class listOfPointsVectorialSmoothness : public listOfPointsScalarSmoothness{
@@ -218,7 +218,7 @@ public:
};
virtual SVector3 get_first_direction(){ return (*points.begin())->direction; };
protected:
- set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
+ std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
};
class listOfPointsFifo : public listOfPoints{
@@ -248,7 +248,7 @@ public:
virtual bool empty(){ return points.empty(); };
protected:
- queue<smoothness_vertex_pair*> points;
+ std::queue<smoothness_vertex_pair*> points;
};
#endif