Skip to content
Snippets Groups Projects
Commit 7bad641b authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

cleanup

parent f63c2cb5
Branches
Tags
No related merge requests found
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Paul-Emile Bernard
#include <limits> #include <limits>
#include "GmshConfig.h" #include "GmshConfig.h"
...@@ -11,7 +13,6 @@ ...@@ -11,7 +13,6 @@
#include "GenericFace.h" #include "GenericFace.h"
#include "Context.h" #include "Context.h"
GenericEdge::ptrfunction_int_double_refvector GenericEdge::EdgeEvalXYZFromT = NULL; GenericEdge::ptrfunction_int_double_refvector GenericEdge::EdgeEvalXYZFromT = NULL;
GenericEdge::ptrfunction_int_refdouble_refdouble GenericEdge::EdgeEvalParBounds = NULL; GenericEdge::ptrfunction_int_refdouble_refdouble GenericEdge::EdgeEvalParBounds = NULL;
GenericEdge::ptrfunction_int_refstring GenericEdge::EdgeGeomType = NULL; GenericEdge::ptrfunction_int_refstring GenericEdge::EdgeGeomType = NULL;
...@@ -22,39 +23,39 @@ GenericEdge::ptrfunction_int_refvector_refdouble_refvector_refbool GenericEdge:: ...@@ -22,39 +23,39 @@ GenericEdge::ptrfunction_int_refvector_refdouble_refvector_refbool GenericEdge::
GenericEdge::ptrfunction_int_refbool GenericEdge::EdgeIs3D = NULL; GenericEdge::ptrfunction_int_refbool GenericEdge::EdgeIs3D = NULL;
GenericEdge::ptrfunction_int_int_double_int_refvector GenericEdge::EdgeReparamOnFace = 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)
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){ : GEdge(m, num, v1, v2), id(_native_id),is_seam(_isseam)
{
if ((!EdgeEvalParBounds)||(!EdgeEvalXYZFromT)) Msg::Error("GenericEdge::ERROR: Callback not set"); if ((!EdgeEvalParBounds)||(!EdgeEvalXYZFromT)) Msg::Error("GenericEdge::ERROR: Callback not set");
bool ok = EdgeEvalParBounds(id,s0,s1); bool ok = EdgeEvalParBounds(id,s0,s1);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalParBounds ! " ); 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); return Range<double>(s0, s1);
} }
SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const
SPoint2 GenericEdge::reparamOnFace(const GFace *face, double par, int dir) const{ {
vector<double> res(2,0.); std::vector<double> res(2,0.);
if (!EdgeReparamOnFace) Msg::Error("GenericEdge::ERROR: Callback EdgeReparamOnFace not set"); if (!EdgeReparamOnFace) Msg::Error("GenericEdge::ERROR: Callback EdgeReparamOnFace not set");
bool ok = EdgeReparamOnFace(id,face->getNativeInt(),par, dir,res); bool ok = EdgeReparamOnFace(id,face->getNativeInt(),par, dir,res);
if (!ok){ if (!ok){
cout << "GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id " << getNativeInt() << endl;
Msg::Error("GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id %d",getNativeInt() ); Msg::Error("GenericEdge::ERROR from EdgeReparamOnFace ! Edge Native id %d",getNativeInt() );
} }
return SPoint2(res[0],res[1]);; return SPoint2(res[0],res[1]);;
} }
GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const
GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const{ {
vector<double> queryPoint(3,0.); std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i]; 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; bool is_on_edge;
if (!EdgeClosestPoint) Msg::Error("GenericEdge::ERROR: Callback EdgeClosestPoint not set"); if (!EdgeClosestPoint) Msg::Error("GenericEdge::ERROR: Callback EdgeClosestPoint not set");
bool ok = EdgeClosestPoint(id,queryPoint,param,res,is_on_edge); bool ok = EdgeClosestPoint(id,queryPoint,param,res,is_on_edge);
...@@ -63,33 +64,32 @@ GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const{ ...@@ -63,33 +64,32 @@ GPoint GenericEdge::closestPoint(const SPoint3 &qp, double &param) const{
return GPoint(res[0], res[1], res[2], this, param); return GPoint(res[0], res[1], res[2], this, param);
} }
bool GenericEdge::isSeam(const GFace *face) const
bool GenericEdge::isSeam(const GFace *face) const{ {
// return false;
return is_seam; return is_seam;
} }
GPoint GenericEdge::point(double par) const
GPoint GenericEdge::point(double par) const{ {
vector<double> res(3,0.); std::vector<double> res(3,0.);
if (!EdgeEvalXYZFromT) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalXYZFromT not set"); if (!EdgeEvalXYZFromT) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalXYZFromT not set");
bool ok = EdgeEvalXYZFromT(id,par,res); bool ok = EdgeEvalXYZFromT(id,par,res);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalXYZFromT ! " ); if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalXYZFromT ! " );
return GPoint(res[0], res[1], res[2], this, par); return GPoint(res[0], res[1], res[2], this, par);
} }
SVector3 GenericEdge::firstDer(double par) const
SVector3 GenericEdge::firstDer(double par) const{ {
vector<double> res(3,0.); std::vector<double> res(3,0.);
if (!EdgeEvalFirstDer) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalFirstDer not set"); if (!EdgeEvalFirstDer) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalFirstDer not set");
bool ok = EdgeEvalFirstDer(id,par,res); bool ok = EdgeEvalFirstDer(id,par,res);
if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalFirstDer ! " ); if (!ok) Msg::Error("GenericEdge::ERROR from EdgeEvalFirstDer ! " );
return SVector3(res[0],res[1],res[2]); return SVector3(res[0],res[1],res[2]);
} }
GEntity::GeomType GenericEdge::geomType() const
GEntity::GeomType GenericEdge::geomType() const{ {
string s; std::string s;
if (!EdgeGeomType) Msg::Error("GenericEdge::ERROR: Callback EdgeGeomType not set"); if (!EdgeGeomType) Msg::Error("GenericEdge::ERROR: Callback EdgeGeomType not set");
bool ok = EdgeGeomType(id,s); bool ok = EdgeGeomType(id,s);
if (!ok){ if (!ok){
...@@ -121,8 +121,8 @@ GEntity::GeomType GenericEdge::geomType() const{ ...@@ -121,8 +121,8 @@ GEntity::GeomType GenericEdge::geomType() const{
return Unknown; return Unknown;
} }
double GenericEdge::curvature(double par) const
double GenericEdge::curvature(double par) const{ {
double res; double res;
if (!EdgeEvalCurvature) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalCurvature not set"); if (!EdgeEvalCurvature) Msg::Error("GenericEdge::ERROR: Callback EdgeEvalCurvature not set");
bool ok = EdgeEvalCurvature(id,par,res); bool ok = EdgeEvalCurvature(id,par,res);
...@@ -130,8 +130,8 @@ double GenericEdge::curvature(double par) const{ ...@@ -130,8 +130,8 @@ double GenericEdge::curvature(double par) const{
return res; return res;
} }
bool GenericEdge::is3D() const
bool GenericEdge::is3D() const{ {
bool res; bool res;
if (!EdgeIs3D) Msg::Error("GenericEdge::ERROR: Callback EdgeIs3D not set"); if (!EdgeIs3D) Msg::Error("GenericEdge::ERROR: Callback EdgeIs3D not set");
bool ok = EdgeIs3D(id,res); bool ok = EdgeIs3D(id,res);
...@@ -139,8 +139,8 @@ bool GenericEdge::is3D() const{ ...@@ -139,8 +139,8 @@ bool GenericEdge::is3D() const{
return res; return res;
} }
bool GenericEdge::degenerate(int) const
bool GenericEdge::degenerate(int) const{ {
bool res=false; bool res=false;
if (!EdgeDegenerated) Msg::Error("GenericEdge::ERROR: Callback EdgeDegenerated not set"); if (!EdgeDegenerated) Msg::Error("GenericEdge::ERROR: Callback EdgeDegenerated not set");
bool ok = EdgeDegenerated(id,res); bool ok = EdgeDegenerated(id,res);
...@@ -148,7 +148,6 @@ bool GenericEdge::degenerate(int) const{ ...@@ -148,7 +148,6 @@ bool GenericEdge::degenerate(int) const{
return res; return res;
} }
int GenericEdge::minimumDrawSegments() const int GenericEdge::minimumDrawSegments() const
{ {
if(geomType() == Line) if(geomType() == Line)
...@@ -157,63 +156,55 @@ int GenericEdge::minimumDrawSegments() const ...@@ -157,63 +156,55 @@ int GenericEdge::minimumDrawSegments() const
return CTX::instance()->geom.numSubEdges * GEdge::minimumDrawSegments(); 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.; s0=0.;
s1=v1->xyz().distance(v2->xyz()); s1=v1->xyz().distance(v2->xyz());
first_der = SVector3(v1->xyz(),v2->xyz()); first_der = SVector3(v1->xyz(),v2->xyz());
first_der.normalize(); 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); return Range<double>(s0, s1);
} }
GPoint LinearSeamEdge::point(double par) const
GPoint LinearSeamEdge::point(double par) const{ {
SVector3 res = v0->xyz() + par*first_der; SVector3 res = v0->xyz() + par*first_der;
return GPoint(res[0], res[1], res[2], this, par); 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; return first_der;
} }
GEntity::GeomType LinearSeamEdge::geomType() const
GEntity::GeomType LinearSeamEdge::geomType() const{ {
return Line; return Line;
} }
double LinearSeamEdge::curvature(double par) const
double LinearSeamEdge::curvature(double par) const{ {
return 0.; return 0.;
} }
bool LinearSeamEdge::is3D() const
bool LinearSeamEdge::is3D() const{ {
return false; return false;
} }
bool LinearSeamEdge::degenerate(int) const
bool LinearSeamEdge::degenerate(int) const{ {
return false; return false;
} }
GPoint LinearSeamEdge::closestPoint(const SPoint3 &q, double &t) const GPoint LinearSeamEdge::closestPoint(const SPoint3 &q, double &t) const
{ {
return GEdge::closestPoint(q,t); return GEdge::closestPoint(q,t);
} }
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_EDGE_H #ifndef _GENERIC_EDGE_H
#define _GENERIC_EDGE_H #define _GENERIC_EDGE_H
...@@ -13,12 +15,11 @@ ...@@ -13,12 +15,11 @@
#include "Range.h" #include "Range.h"
#include <vector> #include <vector>
using namespace std;
class GenericFace; class GenericFace;
/* The set of Generic Entities is a generic interface to any other 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. */ Callbacks (function pointers) are given, sending requests, enquiries, to the
native modeler. */
class GenericEdge : public GEdge { class GenericEdge : public GEdge {
protected: protected:
...@@ -27,13 +28,13 @@ class GenericEdge : public GEdge { ...@@ -27,13 +28,13 @@ class GenericEdge : public GEdge {
const bool is_seam; const bool is_seam;
public: public:
// callbacks typedef // callbacks typedef
typedef bool (*ptrfunction_int_double_refvector)(int, double, vector<double>&); typedef bool (*ptrfunction_int_double_refvector)(int, double, std::vector<double>&);
typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&); typedef bool (*ptrfunction_int_refdouble_refdouble)(int, double&, double&);
typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&); typedef bool (*ptrfunction_int_double_refdouble)(int, double, double&);
typedef bool (*ptrfunction_int_refstring)(int, string&); typedef bool (*ptrfunction_int_refstring)(int, std::string&);
typedef bool (*ptrfunction_int_refbool)(int, bool&); 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_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, vector<double> &); 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); GenericEdge(GModel *model, int num,int _native_id, GVertex *v1, GVertex *v2, bool _isseam=false);
virtual ~GenericEdge(); virtual ~GenericEdge();
...@@ -84,16 +85,6 @@ class GenericEdge : public GEdge { ...@@ -84,16 +85,6 @@ class GenericEdge : public GEdge {
}; };
class LinearSeamEdge : public GEdge { class LinearSeamEdge : public GEdge {
protected: protected:
double s0, s1; double s0, s1;
...@@ -117,5 +108,4 @@ class LinearSeamEdge : public GEdge { ...@@ -117,5 +108,4 @@ class LinearSeamEdge : public GEdge {
}; };
#endif #endif
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // 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 "GmshMessage.h"
#include "GModel.h" #include "GModel.h"
#include "GEdgeLoop.h" #include "GEdgeLoop.h"
...@@ -12,9 +14,6 @@ ...@@ -12,9 +14,6 @@
#include "GenericFace.h" #include "GenericFace.h"
#include "Numeric.h" #include "Numeric.h"
#include "Context.h" #include "Context.h"
#include <math.h>
GenericFace::ptrfunction_int_refstring GenericFace::FaceGeomType = NULL; GenericFace::ptrfunction_int_refstring GenericFace::FaceGeomType = NULL;
GenericFace::ptrfunction_int_refvector_refvector GenericFace::FaceUVFromXYZ = NULL; GenericFace::ptrfunction_int_refvector_refvector GenericFace::FaceUVFromXYZ = NULL;
...@@ -29,8 +28,8 @@ GenericFace::ptrfunction_int_refvector_refvector_refvector GenericFace::FaceFirs ...@@ -29,8 +28,8 @@ GenericFace::ptrfunction_int_refvector_refvector_refvector GenericFace::FaceFirs
GenericFace::ptrfunction_int_refvector_refvector_refvector_refvector GenericFace::FaceSecondDer = NULL; GenericFace::ptrfunction_int_refvector_refvector_refvector_refvector GenericFace::FaceSecondDer = NULL;
GenericFace::ptrfunction_int_refbool_refbool_refdouble_refdouble GenericFace::FacePeriodicInfo = 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"); if (!FaceParBounds) Msg::Fatal("Genericface::ERROR: Callback FaceParBounds not set");
bool ok = FaceParBounds(id,0,umin,umax); bool ok = FaceParBounds(id,0,umin,umax);
...@@ -40,20 +39,20 @@ GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_ ...@@ -40,20 +39,20 @@ GenericFace::GenericFace(GModel *m, int num, int _native_id):GFace(m, num), id(_
computePeriodicity(); 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); if(i == 0) return Range<double>(umin, umax);
return Range<double>(vmin, vmax); return Range<double>(vmin, vmax);
} }
SVector3 GenericFace::normal(const SPoint2 &param) const
SVector3 GenericFace::normal(const SPoint2 &param) const{ {
vector<double> res(3,0.); std::vector<double> res(3,0.);
vector<double> par(2,0.); std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i]; for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceEvalNormal) Msg::Fatal("Genericface::ERROR: Callback FaceEvalNormal not set"); if (!FaceEvalNormal) Msg::Fatal("Genericface::ERROR: Callback FaceEvalNormal not set");
bool ok = FaceEvalNormal(id,par,res); bool ok = FaceEvalNormal(id,par,res);
...@@ -61,12 +60,12 @@ SVector3 GenericFace::normal(const SPoint2 &param) const{ ...@@ -61,12 +60,12 @@ SVector3 GenericFace::normal(const SPoint2 &param) const{
return SVector3(res[0],res[1],res[2]); return SVector3(res[0],res[1],res[2]);
} }
Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const
Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const{ {
if (!FaceFirstDer) Msg::Fatal("Genericface::ERROR: Callback FaceFirstDer not set"); if (!FaceFirstDer) Msg::Fatal("Genericface::ERROR: Callback FaceFirstDer not set");
vector<double> deru(3,0.); std::vector<double> deru(3,0.);
vector<double> derv(3,0.); std::vector<double> derv(3,0.);
vector<double> par(2,0.); std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i]; for (int i=0;i<3;i++) par[i] = param[i];
bool ok = FaceFirstDer(id,par,deru,derv); bool ok = FaceFirstDer(id,par,deru,derv);
if (!ok) Msg::Error("GenericFace::ERROR from FaceFirstDer ! " ); if (!ok) Msg::Error("GenericFace::ERROR from FaceFirstDer ! " );
...@@ -74,12 +73,12 @@ Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const{ ...@@ -74,12 +73,12 @@ Pair<SVector3,SVector3> GenericFace::firstDer(const SPoint2 &param) const{
SVector3(derv[0],derv[1],derv[2])); SVector3(derv[0],derv[1],derv[2]));
} }
void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{ {
vector<double> deruu(3,0.); std::vector<double> deruu(3,0.);
vector<double> dervv(3,0.); std::vector<double> dervv(3,0.);
vector<double> deruv(3,0.); std::vector<double> deruv(3,0.);
vector<double> par(2,0.); std::vector<double> par(2,0.);
for (int i=0;i<3;i++) par[i] = param[i]; for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceSecondDer) Msg::Fatal("Genericface::ERROR: Callback FaceSecondDer not set"); if (!FaceSecondDer) Msg::Fatal("Genericface::ERROR: Callback FaceSecondDer not set");
bool ok = FaceSecondDer(id,par,deruu,dervv,deruv); bool ok = FaceSecondDer(id,par,deruu,dervv,deruv);
...@@ -90,12 +89,12 @@ void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv, ...@@ -90,12 +89,12 @@ void GenericFace::secondDer(const SPoint2 &param,SVector3 *dudu, SVector3 *dvdv,
return; return;
} }
GPoint GenericFace::point(double par1, double par2) const
GPoint GenericFace::point(double par1, double par2) const{ {
vector<double> uv(2,0.); std::vector<double> uv(2,0.);
uv[0] = par1; uv[0] = par1;
uv[1] = par2; uv[1] = par2;
vector<double> xyz(3,0.); std::vector<double> xyz(3,0.);
double pp[2] = {par1, par2}; double pp[2] = {par1, par2};
if (!FaceXYZFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceXYZFromUV not set"); if (!FaceXYZFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceXYZFromUV not set");
bool ok = FaceXYZFromUV(id,uv,xyz); bool ok = FaceXYZFromUV(id,uv,xyz);
...@@ -109,10 +108,11 @@ GPoint GenericFace::point(double par1, double par2) const{ ...@@ -109,10 +108,11 @@ GPoint GenericFace::point(double par1, double par2) const{
} }
GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const{ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const
vector<double> uvres(2,0.); {
vector<double> xyzres(3,0.); std::vector<double> uvres(2,0.);
vector<double> queryPoint(3,0.); std::vector<double> xyzres(3,0.);
std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i]; for (int i=0;i<3;i++) queryPoint[i] = qp[i];
if (!FaceClosestPoint) Msg::Fatal("Genericface::ERROR: Callback FaceClosestPoint not set"); if (!FaceClosestPoint) Msg::Fatal("Genericface::ERROR: Callback FaceClosestPoint not set");
bool ok = FaceClosestPoint(id,queryPoint,xyzres,uvres);// orthogonal projection bool ok = FaceClosestPoint(id,queryPoint,xyzres,uvres);// orthogonal projection
...@@ -121,10 +121,6 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2] ...@@ -121,10 +121,6 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]
// check if the projected point lies on the face... // check if the projected point lies on the face...
bool 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"); if (!FaceContainsPointFromUV) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromUV not set");
ok = FaceContainsPointFromUV(id,uvres,on_the_face); ok = FaceContainsPointFromUV(id,uvres,on_the_face);
if (!ok) Msg::Error("GenericFace::ERROR from FaceContainsPointFromUV ! " ); if (!ok) Msg::Error("GenericFace::ERROR from FaceContainsPointFromUV ! " );
...@@ -133,11 +129,11 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2] ...@@ -133,11 +129,11 @@ GPoint GenericFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]
return GPoint(xyzres[0], xyzres[1], xyzres[2], this, pp); return GPoint(xyzres[0], xyzres[1], xyzres[2], this, pp);
} }
SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{ {
vector<double> uvres(2,0.); std::vector<double> uvres(2,0.);
vector<double> xyzres(3,0.); std::vector<double> xyzres(3,0.);
vector<double> queryPoint(3,0.); std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = qp[i]; for (int i=0;i<3;i++) queryPoint[i] = qp[i];
bool ok=true; bool ok=true;
if (onSurface){ if (onSurface){
...@@ -158,9 +154,9 @@ SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{ ...@@ -158,9 +154,9 @@ SPoint2 GenericFace::parFromPoint(const SPoint3 &qp, bool onSurface) const{
return SPoint2(uvres[0],uvres[1]); return SPoint2(uvres[0],uvres[1]);
} }
GEntity::GeomType GenericFace::geomType() const
GEntity::GeomType GenericFace::geomType() const{ {
string s; std::string s;
if (!FaceGeomType) Msg::Fatal("Genericface::ERROR: Callback FaceGeomType not set"); if (!FaceGeomType) Msg::Fatal("Genericface::ERROR: Callback FaceGeomType not set");
bool ok = FaceGeomType(id,s); bool ok = FaceGeomType(id,s);
if (!ok){ if (!ok){
...@@ -188,12 +184,12 @@ GEntity::GeomType GenericFace::geomType() const{ ...@@ -188,12 +184,12 @@ GEntity::GeomType GenericFace::geomType() const{
return Unknown; return Unknown;
} }
double GenericFace::curvatureMax(const SPoint2 &param) const
double GenericFace::curvatureMax(const SPoint2 &param) const{ {
vector<double> dirMax(3,0.); std::vector<double> dirMax(3,0.);
vector<double> dirMin(3,0.); std::vector<double> dirMin(3,0.);
double curvMax,curvMin; 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]; for (int i=0;i<3;i++) par[i] = param[i];
if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set"); if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set");
bool ok = FaceCurvatures(id,par,dirMax,dirMin,curvMax,curvMin); bool ok = FaceCurvatures(id,par,dirMax,dirMin,curvMax,curvMin);
...@@ -201,12 +197,13 @@ double GenericFace::curvatureMax(const SPoint2 &param) const{ ...@@ -201,12 +197,13 @@ double GenericFace::curvatureMax(const SPoint2 &param) const{
return std::max(fabs(curvMax), fabs(curvMin)); return std::max(fabs(curvMax), fabs(curvMin));
} }
double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,
double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 *_dirMin,double *curvMax,double *curvMin) const{ double *curvMax,double *curvMin) const
vector<double> param(2,0.); {
std::vector<double> param(2,0.);
for (int i=0;i<3;i++) param[i] = _param[i]; for (int i=0;i<3;i++) param[i] = _param[i];
vector<double> dirMax(3,0.); std::vector<double> dirMax(3,0.);
vector<double> dirMin(3,0.); std::vector<double> dirMin(3,0.);
if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set"); if (!FaceCurvatures) Msg::Fatal("Genericface::ERROR: Callback FaceCurvatures not set");
bool ok = FaceCurvatures(id,param,dirMax,dirMin,*curvMax,*curvMin); bool ok = FaceCurvatures(id,param,dirMax,dirMin,*curvMax,*curvMin);
...@@ -217,10 +214,10 @@ double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3 ...@@ -217,10 +214,10 @@ double GenericFace::curvatures(const SPoint2 &_param,SVector3 *_dirMax,SVector3
return *curvMax; return *curvMax;
} }
bool GenericFace::containsPoint(const SPoint3 &pt) const
bool GenericFace::containsPoint(const SPoint3 &pt) const{ {
bool res; bool res;
vector<double> queryPoint(3,0.); std::vector<double> queryPoint(3,0.);
for (int i=0;i<3;i++) queryPoint[i] = pt[i]; for (int i=0;i<3;i++) queryPoint[i] = pt[i];
if (!FaceContainsPointFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromXYZ not set"); if (!FaceContainsPointFromXYZ) Msg::Fatal("Genericface::ERROR: Callback FaceContainsPointFromXYZ not set");
bool ok = FaceContainsPointFromXYZ(id,queryPoint,res); bool ok = FaceContainsPointFromXYZ(id,queryPoint,res);
...@@ -228,47 +225,33 @@ bool GenericFace::containsPoint(const SPoint3 &pt) const{ ...@@ -228,47 +225,33 @@ bool GenericFace::containsPoint(const SPoint3 &pt) const{
return res; return res;
} }
void GenericFace::computePeriodicity()
void GenericFace::computePeriodicity(){ {
if (!FacePeriodicInfo) Msg::Fatal("Genericface::ERROR: Callback FacePeriodicInfo not set"); if (!FacePeriodicInfo) Msg::Fatal("Genericface::ERROR: Callback FacePeriodicInfo not set");
bool ok = FacePeriodicInfo(id,_periodic[0], _periodic[1],_period[0],_period[1]); bool ok = FacePeriodicInfo(id,_periodic[0], _periodic[1],_period[0],_period[1]);
if (!ok) Msg::Error("GenericFace::ERROR from FacePeriodicInfo ! " ); 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(); edgeLoops.clear();
l_edges.clear(); l_edges.clear();
l_dirs.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 for (std::set<int>::iterator it_loop = loopsnumber.begin();
if (debug) cout << "Loop #" << *it_loop << " made of edges "; it_loop!=loopsnumber.end();it_loop++){// for each loop
pair<multimap<int, pair<GEdge*,int> >::iterator, multimap<int, pair<GEdge*,int> >::iterator> range = bnd.equal_range(*it_loop); std::pair<std::multimap<int, std::pair<GEdge*,int> >::iterator,
list<GEdge*> l_wire; std::multimap<int, std::pair<GEdge*,int> >::iterator>
for (multimap<int, pair<GEdge*,int> >::iterator it = range.first;it!=range.second;it++){// for all edges range = bnd.equal_range(*it_loop);
if (debug) cout << it->second.first->tag() << " "; 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); l_wire.push_back(it->second.first);
} }
if (debug) cout << endl << "-> GEdgeLoop made of edges (sign) ";
GEdgeLoop el(l_wire); GEdgeLoop el(l_wire);
for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){ for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){
l_edges.push_back(it->ge); l_edges.push_back(it->ge);
l_dirs.push_back(it->_sign); l_dirs.push_back(it->_sign);
if (debug) cout << it->ge->tag() << "(" << ((it->_sign<0.) ? "-" : "+") << ") " ;
if (el.count() == 2){ if (el.count() == 2){
it->ge->meshAttributes.minimumMeshSegments = it->ge->meshAttributes.minimumMeshSegments =
std::max(it->ge->meshAttributes.minimumMeshSegments,2); std::max(it->ge->meshAttributes.minimumMeshSegments,2);
...@@ -278,9 +261,6 @@ void GenericFace::createLoops(){ ...@@ -278,9 +261,6 @@ void GenericFace::createLoops(){
std::max(it->ge->meshAttributes.minimumMeshSegments,3); std::max(it->ge->meshAttributes.minimumMeshSegments,3);
} }
} }
if (debug) cout << endl;
edgeLoops.push_back(el); edgeLoops.push_back(el);
} }
} }
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_FACE_H_ #ifndef _GENERIC_FACE_H_
#define _GENERIC_FACE_H_ #define _GENERIC_FACE_H_
...@@ -14,10 +16,9 @@ ...@@ -14,10 +16,9 @@
#include "Range.h" #include "Range.h"
#include <vector> #include <vector>
using namespace std;
/* The set of Generic Entities is a generic interface to any other 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. */ Callbacks (function pointers) are given, sending requests, enquiries, to the
native modeler. */
class GenericFace : public GFace { class GenericFace : public GFace {
protected: protected:
...@@ -28,14 +29,14 @@ class GenericFace : public GFace { ...@@ -28,14 +29,14 @@ class GenericFace : public GFace {
public: public:
// callbacks typedef // callbacks typedef
typedef bool (*ptrfunction_int_refstring)(int, string&); 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_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)(const int , const std::vector<double> &, std::vector<double> &);
typedef bool (*ptrfunction_int_refvector_refvector_refvector)(const int , const vector<double> &, vector<double> &, 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 vector<double> &, bool &); 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_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_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 vector<double> &, vector<double> &, vector<double> &, vector<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); GenericFace(GModel *m, int num, int _native_id);
virtual ~GenericFace(); virtual ~GenericFace();
...@@ -59,8 +60,9 @@ class GenericFace : public GFace { ...@@ -59,8 +60,9 @@ class GenericFace : public GFace {
ModelType getNativeType() const { return GenericModel; } ModelType getNativeType() const { return GenericModel; }
virtual int getNativeInt()const{return id;}; virtual int getNativeInt()const{return id;};
void addBndInfo(int loop, GEdge *ptr, int sign){ void addBndInfo(int loop, GEdge *ptr, int sign)
bnd.insert(make_pair(loop, make_pair(ptr,sign))); {
bnd.insert(std::make_pair(loop, std::make_pair(ptr,sign)));
l_dirs.push_back(sign); l_dirs.push_back(sign);
l_edges.push_back(ptr); l_edges.push_back(ptr);
loopsnumber.insert(loop); loopsnumber.insert(loop);
...@@ -85,8 +87,8 @@ class GenericFace : public GFace { ...@@ -85,8 +87,8 @@ class GenericFace : public GFace {
static void setFacePeriodicInfo(ptrfunction_int_refbool_refbool_refdouble_refdouble fct){FacePeriodicInfo = fct;}; static void setFacePeriodicInfo(ptrfunction_int_refbool_refbool_refdouble_refdouble fct){FacePeriodicInfo = fct;};
private: private:
multimap<int, pair<GEdge*,int> > bnd; std::multimap<int, std::pair<GEdge*,int> > bnd;
set<int> loopsnumber; std::set<int> loopsnumber;
// the callbacks: // the callbacks:
static ptrfunction_int_refstring FaceGeomType; static ptrfunction_int_refstring FaceGeomType;
...@@ -97,7 +99,6 @@ class GenericFace : public GFace { ...@@ -97,7 +99,6 @@ class GenericFace : public GFace {
static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures; static ptrfunction_int_refvector_refvector_refvector_refdouble_refdouble FaceCurvatures;
static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer; static ptrfunction_int_refvector_refvector_refvector_refvector FaceSecondDer;
static ptrfunction_int_refbool_refbool_refdouble_refdouble FacePeriodicInfo; static ptrfunction_int_refbool_refbool_refdouble_refdouble FacePeriodicInfo;
}; };
#endif #endif
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // 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 "GModel.h"
#include "GenericVertex.h" #include "GenericVertex.h"
#include "GenericEdge.h" #include "GenericEdge.h"
#include "GenericFace.h" #include "GenericFace.h"
#include "GenericRegion.h" #include "GenericRegion.h"
GenericRegion::GenericRegion(GModel *m, int num, int _native_id)
GenericRegion::GenericRegion(GModel *m, int num, int _native_id):GRegion(m, num), id(_native_id) : GRegion(m, num), id(_native_id)
{ {
} }
GenericRegion::~GenericRegion() GenericRegion::~GenericRegion()
{ {
} }
GEntity::GeomType GenericRegion::geomType() const GEntity::GeomType GenericRegion::geomType() const
{ {
return Unknown; return Unknown;
......
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_REGION_H_ #ifndef _GENERIC_REGION_H_
#define _GENERIC_REGION_H_ #define _GENERIC_REGION_H_
...@@ -10,7 +12,8 @@ ...@@ -10,7 +12,8 @@
#include "GRegion.h" #include "GRegion.h"
/* The set of Generic Entities is a generic interface to any other 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. */ Callbacks (function pointers) are given, sending requests, enquiries, to the
native modeler. */
class GenericRegion : public GRegion { class GenericRegion : public GRegion {
public: public:
...@@ -22,14 +25,15 @@ class GenericRegion : public GRegion { ...@@ -22,14 +25,15 @@ class GenericRegion : public GRegion {
ModelType getNativeType() const { return GenericModel; } ModelType getNativeType() const { return GenericModel; }
virtual int getNativeInt()const{return id;}; 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. // 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){ void addFace(GenericFace *ptr, int sign){
l_dirs.push_back(sign); l_dirs.push_back(sign);
l_faces.push_back(ptr); l_faces.push_back(ptr);
ptr->addRegion(this); ptr->addRegion(this);
}; };
private: private:
int id; int id;
}; };
......
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // 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 "GModel.h"
#include "MVertex.h" #include "MVertex.h"
#include "MPoint.h" #include "MPoint.h"
#include "GenericVertex.h" #include "GenericVertex.h"
#include "GenericEdge.h" #include "GenericEdge.h"
#include "GenericFace.h" #include "GenericFace.h"
#include<algorithm>
GenericVertex::ptrfunction_int_vector GenericVertex::VertexXYZ = NULL; GenericVertex::ptrfunction_int_vector GenericVertex::VertexXYZ = NULL;
GenericVertex::ptrfunction_int_doubleptr_voidptr GenericVertex::VertexMeshSize = NULL; GenericVertex::ptrfunction_int_doubleptr_voidptr GenericVertex::VertexMeshSize = NULL;
GenericVertex::GenericVertex(GModel *m, int num, int _native_id)
GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num), id(_native_id){ : GVertex(m, num), id(_native_id)
{
if (!VertexXYZ) if (!VertexXYZ)
Msg::Fatal("GenericVertex::ERROR: Callback not set"); Msg::Fatal("GenericVertex::ERROR: Callback not set");
vector<double> vec(3,0.); std::vector<double> vec(3,0.);
bool ok = VertexXYZ(id,vec); bool ok = VertexXYZ(id,vec);
if (!ok) Msg::Error("GenericVertex::ERROR from callback VertexXYZ "); if (!ok) Msg::Error("GenericVertex::ERROR from callback VertexXYZ ");
_x=vec[0]; _x=vec[0];
...@@ -29,30 +30,26 @@ GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num) ...@@ -29,30 +30,26 @@ GenericVertex::GenericVertex(GModel *m, int num, int _native_id):GVertex(m, num)
_z=vec[2]; _z=vec[2];
} }
GenericVertex::GenericVertex(GModel *m, int num, int _native_id, const std::vector<double> &vec)
GenericVertex::GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec):GVertex(m, num), id(_native_id){ : GVertex(m, num), id(_native_id)
{
if (!VertexXYZ) if (!VertexXYZ)
Msg::Fatal("GenericVertex::ERROR: Callback not set"); Msg::Fatal("GenericVertex::ERROR: Callback not set");
_x=vec[0]; _x=vec[0];
_y=vec[1]; _y=vec[1];
_z=vec[2]; _z=vec[2];
} }
GenericVertex::~GenericVertex(){ GenericVertex::~GenericVertex(){
}
}
SPoint2 GenericVertex::reparamOnFace(const GFace *gf, int dir) const SPoint2 GenericVertex::reparamOnFace(const GFace *gf, int dir) const
{ {
SPoint3 pt(_x,_y,_z); SPoint3 pt(_x,_y,_z);
return gf->parFromPoint(pt,true); return gf->parFromPoint(pt,true);
} }
void GenericVertex::setPosition(GPoint &p) void GenericVertex::setPosition(GPoint &p)
{ {
_x = p.x(); _x = p.x();
......
// 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 // See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>. // bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributed by Paul-Emile Bernard
#ifndef _GENERIC_VERTEX_H_ #ifndef _GENERIC_VERTEX_H_
#define _GENERIC_VERTEX_H_ #define _GENERIC_VERTEX_H_
...@@ -11,12 +13,10 @@ ...@@ -11,12 +13,10 @@
#include "GVertex.h" #include "GVertex.h"
#include "Context.h" #include "Context.h"
#include <vector> #include <vector>
#include <iostream>
using namespace std;
/* The set of Generic Entities is a generic interface to any other 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. */ Callbacks (function pointers) are given, sending requests, enquiries, to the
native modeler. */
class GenericVertex : public GVertex { class GenericVertex : public GVertex {
protected: protected:
...@@ -24,11 +24,11 @@ class GenericVertex : public GVertex { ...@@ -24,11 +24,11 @@ class GenericVertex : public GVertex {
double _x, _y, _z; double _x, _y, _z;
public: public:
// callbacks typedef // callbacks typedef
typedef bool (*ptrfunction_int_vector)(int, vector<double>&); typedef bool (*ptrfunction_int_vector)(int, std::vector<double>&);
typedef bool (*ptrfunction_int_doubleptr_voidptr)(int, double*, void*); 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);
GenericVertex(GModel *m, int num, int _native_id, const vector<double> &vec); GenericVertex(GModel *m, int num, int _native_id, const std::vector<double> &vec);
virtual ~GenericVertex(); virtual ~GenericVertex();
virtual GPoint point() const { return GPoint(x(), y(), z()); } virtual GPoint point() const { return GPoint(x(), y(), z()); }
...@@ -47,17 +47,16 @@ class GenericVertex : public GVertex { ...@@ -47,17 +47,16 @@ class GenericVertex : public GVertex {
static void setVertexMeshSize(ptrfunction_int_doubleptr_voidptr fct){VertexMeshSize = fct;}; static void setVertexMeshSize(ptrfunction_int_doubleptr_voidptr fct){VertexMeshSize = fct;};
// meshing-related methods: // meshing-related methods:
virtual void setPrescribedMeshSizeAtVertex(double l) { virtual void setPrescribedMeshSizeAtVertex(double l)
cout << "GenericVertex:: setPrescribedMeshSizeAtVertex !!!, aborting" << endl; {
throw; Msg::Error("GenericVertex::setPrescribedMeshSizeAtVertex");
} }
virtual inline double prescribedMeshSizeAtVertex() const { virtual inline double prescribedMeshSizeAtVertex() const
{
double size; double size;
void *chose = NULL; void *chose = NULL;
if (!VertexMeshSize(id,&size,chose)){ if (!VertexMeshSize(id,&size,chose)){
cout<< "GenericVertex::ERROR from callback VertexMeshSize "<< endl; Msg::Error("GenericVertex::ERROR from callback VertexMeshSize");
Msg::Warning("GenericVertex::ERROR from callback VertexMeshSize ");
cout << "Check for ring edges !!!" << endl;
return CTX::instance()->lc; return CTX::instance()->lc;
} }
return size; return size;
......
#include "BackgroundMesh3D.h"
#include <fstream> #include <fstream>
#include <algorithm> #include <algorithm>
#include <iostream> #include <iostream>
#include <string> #include <string>
#include <numeric> #include <numeric>
#include "BackgroundMesh3D.h"
#include <time.h>
#include "MElement.h" #include "MElement.h"
#include "GFace.h" #include "GFace.h"
#include "GRegion.h" #include "GRegion.h"
#include "BackgroundMeshManager.h" #include "BackgroundMeshManager.h"
#include "BackgroundMesh2D.h" #include "BackgroundMesh2D.h"
#include "MElementOctree.h" #include "MElementOctree.h"
//#include "pointInsertion.h"
#include "MTetrahedron.h" #include "MTetrahedron.h"
#include "OS.h" #include "OS.h"
...@@ -28,12 +21,13 @@ ...@@ -28,12 +21,13 @@
#include "linearSystemFull.h" #include "linearSystemFull.h"
#endif #endif
using namespace std;
int signof(int i){ static int signof(int i)
{
return ((i < 0) ? -1 : 1); return ((i < 0) ? -1 : 1);
} }
// TODO: virer les trucs "vertextorank", mettre cette classe-ci : // TODO: virer les trucs "vertextorank", mettre cette classe-ci :
//class evalDiffusivity3DFunction : public simpleFunction<double>{ //class evalDiffusivity3DFunction : public simpleFunction<double>{
...@@ -47,18 +41,18 @@ int signof(int i){ ...@@ -47,18 +41,18 @@ int signof(int i){
// const double threshold; // const double threshold;
//}; //};
backgroundMesh3D::backgroundMesh3D(GRegion *_gr)
backgroundMesh3D::backgroundMesh3D(GRegion *_gr):BGMBase(3,_gr),debug(false),verbose(false){ : BGMBase(3,_gr), debug(false), verbose(false)
{
computeSizeField(); computeSizeField();
} }
backgroundMesh3D::~backgroundMesh3D()
backgroundMesh3D::~backgroundMesh3D(){ {
} }
void backgroundMesh3D::computeSizeField()
void backgroundMesh3D::computeSizeField(){ {
cout << "backgroundMesh3D::computeSizeField() " << endl; cout << "backgroundMesh3D::computeSizeField() " << endl;
// fills dirichlet BC // fills dirichlet BC
...@@ -90,8 +84,10 @@ void backgroundMesh3D::computeSizeField(){ ...@@ -90,8 +84,10 @@ void backgroundMesh3D::computeSizeField(){
// cout << "backgroundMesh3D::size of sizeField: " << sizeField.size() << endl; // cout << "backgroundMesh3D::size of sizeField: " << sizeField.size() << endl;
} }
void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet,
void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet, simpleFunction<double> &eval_diffusivity, bool in_parametric_plane){ simpleFunction<double> &eval_diffusivity,
bool in_parametric_plane)
{
// same as Size_field::solve() // same as Size_field::solve()
GRegion *gr = dynamic_cast<GRegion*>(gf); GRegion *gr = dynamic_cast<GRegion*>(gf);
...@@ -161,42 +157,44 @@ void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet, simpleFunct ...@@ -161,42 +157,44 @@ void backgroundMesh3D::propagateValues(DoubleStorageType &dirichlet, simpleFunct
#endif #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)); 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; double distance;
return get_nearest_neighbor(v->x(),v->y(),v->z(), 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; double distance;
return get_nearest_neighbor(xyz, 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; double distance;
return get_nearest_neighbor(x,y,z,distance); return get_nearest_neighbor(x,y,z,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 ){ double &distance )
{
double xyz[3] = {x,y,z}; double xyz[3] = {x,y,z};
return get_nearest_neighbor(xyz, distance); return get_nearest_neighbor(xyz, distance);
} }
MElementOctree* backgroundMesh3D::getOctree()
MElementOctree* backgroundMesh3D::getOctree(){ {
if(!octree){ if(!octree){
GRegion *gr = dynamic_cast<GRegion*>(gf); GRegion *gr = dynamic_cast<GRegion*>(gf);
Msg::Debug("Rebuilding BackgroundMesh element octree"); Msg::Debug("Rebuilding BackgroundMesh element octree");
std::vector<MElement*> copy;// !!! 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); copy.push_back(*it);
} }
octree = new MElementOctree(copy); octree = new MElementOctree(copy);
...@@ -204,8 +202,8 @@ MElementOctree* backgroundMesh3D::getOctree(){ ...@@ -204,8 +202,8 @@ MElementOctree* backgroundMesh3D::getOctree(){
return octree; 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. // using the octree instead of ANN, faster.
MElement *elem = const_cast<MElement*>(findElement(xyz[0],xyz[1],xyz[2])); MElement *elem = const_cast<MElement*>(findElement(xyz[0],xyz[1],xyz[2]));
if (!elem) return NULL; if (!elem) return NULL;
...@@ -220,7 +218,6 @@ MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & dist ...@@ -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()); vector<double>::iterator itmax = std::max_element(distances.begin(),distances.end());
return candidates[std::distance(distances.begin(),itmax)]; return candidates[std::distance(distances.begin(),itmax)];
// map<double,MVertex*> distances; // map<double,MVertex*> distances;
// SPoint3 p(xyz[0],xyz[1],xyz[2]); // SPoint3 p(xyz[0],xyz[1],xyz[2]);
// for (int i=0;i<elem->getNumVertices();i++){ // for (int i=0;i<elem->getNumVertices();i++){
...@@ -232,13 +229,10 @@ MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & dist ...@@ -232,13 +229,10 @@ MVertex* backgroundMesh3D::get_nearest_neighbor(const double* xyz, double & dist
// return it->second; // return it->second;
} }
vector<montripletbis> frameFieldBackgroundMesh3D::permutation = vector<montripletbis>(); 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); // readValue("param.dat","SMOOTHCF",smooth_the_crossfield);
smooth_the_crossfield = true; smooth_the_crossfield = true;
...@@ -271,16 +265,16 @@ frameFieldBackgroundMesh3D::frameFieldBackgroundMesh3D(GRegion *_gr):backgroundM ...@@ -271,16 +265,16 @@ frameFieldBackgroundMesh3D::frameFieldBackgroundMesh3D(GRegion *_gr):backgroundM
} }
} }
frameFieldBackgroundMesh3D::~frameFieldBackgroundMesh3D()
frameFieldBackgroundMesh3D::~frameFieldBackgroundMesh3D(){ {
#if defined(HAVE_ANN) #if defined(HAVE_ANN)
if(annTreeBnd) delete annTreeBnd; if(annTreeBnd) delete annTreeBnd;
if(dataPtsBnd) annDeallocPts(dataPtsBnd); if(dataPtsBnd) annDeallocPts(dataPtsBnd);
#endif #endif
} }
void frameFieldBackgroundMesh3D::initiate_ANN_research()
void frameFieldBackgroundMesh3D::initiate_ANN_research(){ {
#ifdef HAVE_ANN #ifdef HAVE_ANN
// ANN research for 2D !!! // ANN research for 2D !!!
int maxPts = listOfBndVertices.size(); int maxPts = listOfBndVertices.size();
...@@ -298,8 +292,8 @@ void frameFieldBackgroundMesh3D::initiate_ANN_research(){ ...@@ -298,8 +292,8 @@ void frameFieldBackgroundMesh3D::initiate_ANN_research(){
return; return;
} }
void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries()
void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){ {
// this is used when the crossfield is not smoothed !!! // this is used when the crossfield is not smoothed !!!
// otherwise, the smoothness is computed on the fly while smoothing the cf ! // otherwise, the smoothness is computed on the fly while smoothing the cf !
...@@ -311,7 +305,8 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){ ...@@ -311,7 +305,8 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){
double mean_angle=0.; double mean_angle=0.;
vector<double> vectorial_smoothness(3); 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(); themap.clear();
neighbors.clear(); neighbors.clear();
current = it_vertex->first; current = it_vertex->first;
...@@ -328,13 +323,15 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){ ...@@ -328,13 +323,15 @@ void frameFieldBackgroundMesh3D::computeSmoothnessOnlyFromBoundaries(){
STensor3 &ref = itcurrent->second; STensor3 &ref = itcurrent->second;
multimap<double,MVertex*>::iterator it_neighbors_begin = themap.begin(); 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; // crossFieldVectorialSmoothness[current] = vectorial_smoothness;
} }
} }
void frameFieldBackgroundMesh3D::computeCrossField()
void frameFieldBackgroundMesh3D::computeCrossField(){ {
// TODO: mettre des parametres façon gmsh ??? // TODO: mettre des parametres façon gmsh ???
// reading parameters // reading parameters
int Niter = 3; int Niter = 3;
...@@ -347,7 +344,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -347,7 +344,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
smoothness_threshold=0.85; smoothness_threshold=0.85;
double percentage_threshold = 0.98; double percentage_threshold = 0.98;
vector<double> scal_prods; vector<double> scal_prods;
vector<double> signs; vector<double> signs;
vector<double> angles; vector<double> angles;
...@@ -360,7 +356,8 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -360,7 +356,8 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
map<MVertex*,bool> vertex_is_still; map<MVertex*,bool> vertex_is_still;
map<MVertex*,double> vertex_movement; map<MVertex*,double> vertex_movement;
MVertex *current; 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; current = it_vertex->first;
if (current->onWhat()->dim()<=2) if (current->onWhat()->dim()<=2)
vertex_is_still[current] = true; vertex_is_still[current] = true;
...@@ -371,15 +368,14 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -371,15 +368,14 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
// OLD - NEW COMPARISON // OLD - NEW COMPARISON
map<MVertex*,double> vertex_to_rank; 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] = 0.;
vertex_to_rank[it_vertex->first] = 1.; vertex_to_rank[it_vertex->first] = 1.;
rank.insert(make_pair(0.,it_vertex->first)); rank.insert(make_pair(0.,it_vertex->first));
} }
// END OLD - NEW COMPARISON // END OLD - NEW COMPARISON
double percentage_not_done=0.; double percentage_not_done=0.;
double norme = 10.*norme_threshold; double norme = 10.*norme_threshold;
double norme_inf=10.*norme_threshold; double norme_inf=10.*norme_threshold;
...@@ -389,34 +385,44 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -389,34 +385,44 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
vector<double> vectorial_smoothness(3); vector<double> vectorial_smoothness(3);
// vector<int> nb_local_iterations; // vector<int> nb_local_iterations;
// main smoothing loop // 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(); // nb_local_iterations.clear();
angles.clear(); angles.clear();
norme_inf=0.; norme_inf=0.;
int counter_done=0; int counter_done=0;
int counter_not_done=0; int counter_not_done=0;
while (rank.size()){// for all vertices, i.e. all cavities 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 // for (vert2elemtype::iterator
//MVertex *current = it_vertex->first; //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; current = (rank.begin())->second;
rank.erase((rank.begin())); rank.erase((rank.begin()));
// smooth 3D vertices only // smooth 3D vertices only
if (current->onWhat()->dim()<=2){ 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; continue;
} }
TensorStorageType::iterator itcurrent = crossField.find(current); TensorStorageType::iterator itcurrent = crossField.find(current);
STensor3 &ref = itcurrent->second; 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); pair<graphtype::iterator, graphtype::iterator> range = graph.equal_range(current);
graphtype::iterator itgraph = range.first; 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 for (;itgraph!=range.second;itgraph++){// for all neighbors
if (vertex_is_still[itgraph->second.second]==false){ if (vertex_is_still[itgraph->second.second]==false){
all_neighbors_still = false; all_neighbors_still = false;
...@@ -424,22 +430,23 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -424,22 +430,23 @@ 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.; vertex_movement[current] = 0.;
rank_temp.insert(make_pair(0.,current)); rank_temp.insert(make_pair(0.,current));
counter_not_done++; counter_not_done++;
continue; continue;
} }
// OLD - NEW COMPARISON // OLD - NEW COMPARISON
multimap<double,MVertex*> neighbors_to_trust; multimap<double,MVertex*> neighbors_to_trust;
itgraph = range.first; itgraph = range.first;
for (;itgraph!=range.second;itgraph++){// for all neighbors for (;itgraph!=range.second;itgraph++){// for all neighbors
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){ // 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 (!neighbors_to_trust.size()){ if (!neighbors_to_trust.size()){
...@@ -449,11 +456,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -449,11 +456,6 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
} }
// END OLD - NEW COMPARISON // END OLD - NEW COMPARISON
counter_done++; counter_done++;
double angle_to_go;; double angle_to_go;;
...@@ -463,26 +465,27 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -463,26 +465,27 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
for (;Nlocaliter<20;Nlocaliter++){// iterations, convergence of the local cavity... for (;Nlocaliter<20;Nlocaliter++){// iterations, convergence of the local cavity...
multimap<double,MVertex*>::iterator it_neighbors_to_trust = neighbors_to_trust.begin(); 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); 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; // crossFieldVectorialSmoothness[current] = vectorial_smoothness;
// rotating directions to align closest vectors... // rotating directions to align closest vectors...
angle_to_go = mean_angle; angle_to_go = mean_angle;
ref = apply_rotation(mean_axis,angle_to_go,ref); 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; 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; 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);
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); 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; // cout << " #local_iter=" << Nlocaliter << " final mean_angle = " << mean_angle << " threshold=" << localangle_threshold << " condition :" << (mean_angle <= localangle_threshold) << endl;
angles.push_back(mean_angle); angles.push_back(mean_angle);
vertex_is_still[current] = (mean_angle <= localangle_threshold); vertex_is_still[current] = (mean_angle <= localangle_threshold);
vertex_movement[current] = mean_angle; vertex_movement[current] = mean_angle;
...@@ -502,8 +505,13 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -502,8 +505,13 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
norme = std::accumulate(angles.begin(),angles.end(),0.); norme = std::accumulate(angles.begin(),angles.end(),0.);
if (angles.size()) norme = norme/M_PI*180./angles.size(); if (angles.size()) norme = norme/M_PI*180./angles.size();
percentage_not_done = ((double)counter_not_done)/(counter_done+counter_not_done); 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 << " --- iter " << iIter << " NormeInf=" << norme_inf
// cout << "Average number of local iterations: " << ((double)std::accumulate(nb_local_iterations.begin(),nb_local_iterations.end(),0))/nb_local_iterations.size() << endl; << " 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){ //if (debug){
double temp = log10(Niter); double temp = log10(Niter);
...@@ -515,10 +523,9 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -515,10 +523,9 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
iIter++; iIter++;
}// end Niter iterations }// end Niter iterations
// also computes smoothness for boundary points // 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; current = it_vertex->first;
if (current->onWhat()->dim()<=2){ if (current->onWhat()->dim()<=2){
TensorStorageType::iterator itcurrent = crossField.find(current); TensorStorageType::iterator itcurrent = crossField.find(current);
...@@ -526,26 +533,23 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -526,26 +533,23 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
pair<graphtype::iterator, graphtype::iterator> range = graph.equal_range(current); pair<graphtype::iterator, graphtype::iterator> range = graph.equal_range(current);
graphtype::iterator itgraph = range.first; graphtype::iterator itgraph = range.first;
multimap<double,MVertex*> neighbors_to_trust; multimap<double,MVertex*> neighbors_to_trust;
itgraph = range.first; itgraph = range.first;
for (;itgraph!=range.second;itgraph++){// for all neighbors for (;itgraph!=range.second;itgraph++){// for all neighbors
neighbors_to_trust.insert(make_pair(0.,itgraph->second.second)); 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; // crossFieldVectorialSmoothness[current] = vectorial_smoothness;
} }
} }
return;
} }
//STensor3 frameFieldBackgroundMesh3D::get_random_cross()const{ //STensor3 frameFieldBackgroundMesh3D::get_random_cross()const{
// SVector3 vec1(rand()%10000/9999. , rand()%10000/9999. , rand()%10000/9999. ); // SVector3 vec1(rand()%10000/9999. , rand()%10000/9999. , rand()%10000/9999. );
// vec1.normalize(); // vec1.normalize();
...@@ -563,18 +567,20 @@ void frameFieldBackgroundMesh3D::computeCrossField(){ ...@@ -563,18 +567,20 @@ void frameFieldBackgroundMesh3D::computeCrossField(){
// return random_cross; // return random_cross;
//} //}
void frameFieldBackgroundMesh3D::initiate_crossfield()
void frameFieldBackgroundMesh3D::initiate_crossfield(){ {
crossField.clear(); crossField.clear();
MVertex *v; MVertex *v;
// first, for boundaries : // first, for boundaries :
GRegion *gr = dynamic_cast<GRegion*>(gf); GRegion *gr = dynamic_cast<GRegion*>(gf);
list<GFace*> faces = gr->faces(); 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... // 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 for (list<GFace*>::iterator it=faces.begin();it!=faces.end();it++){// for all GFace
GFace *face = *it; GFace *face = *it;
frameFieldBackgroundMesh2D *bgm2d = dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face)); frameFieldBackgroundMesh2D *bgm2d =
dynamic_cast<frameFieldBackgroundMesh2D*>(BGMManager::get(face));
// initializing the vertices on the faces // initializing the vertices on the faces
for (unsigned int i=0;i<face->getNumMeshElements();i++){// for all elements for (unsigned int i=0;i<face->getNumMeshElements();i++){// for all elements
MElement *e = face->getMeshElement(i); MElement *e = face->getMeshElement(i);
...@@ -613,14 +619,15 @@ void frameFieldBackgroundMesh3D::initiate_crossfield(){ ...@@ -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; double distance;
return get_nearest_neighbor_on_boundary(v,distance); return get_nearest_neighbor_on_boundary(v,distance);
} }
MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary
MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v, double &distance){ (MVertex* v, double &distance)
{
#ifdef HAVE_ANN #ifdef HAVE_ANN
ANNpoint q = annAllocPt(3); ANNpoint q = annAllocPt(3);
for (int k=0; k< 3; ++k) for (int k=0; k< 3; ++k)
...@@ -642,33 +649,35 @@ MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v ...@@ -642,33 +649,35 @@ MVertex* frameFieldBackgroundMesh3D::get_nearest_neighbor_on_boundary(MVertex* v
#endif #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); 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)]; 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()); 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); return get_nodal_value(v,crossFieldSmoothness);
}; };
void frameFieldBackgroundMesh3D::eval_approximate_crossfield(double x, double y,
void frameFieldBackgroundMesh3D::eval_approximate_crossfield(double x, double y, double z, STensor3 &cf){ double z, STensor3 &cf)
{
cf = crossField[get_nearest_neighbor(x,y,z)]; cf = crossField[get_nearest_neighbor(x,y,z)];
}; };
void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert,
void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert, STensor3 &cf){ STensor3 &cf)
{
// check if vertex is ours... // check if vertex is ours...
TensorStorageType::const_iterator itfind = crossField.find(const_cast<MVertex*>(vert)); TensorStorageType::const_iterator itfind = crossField.find(const_cast<MVertex*>(vert));
if (itfind != crossField.end()) if (itfind != crossField.end())
...@@ -678,31 +687,35 @@ void frameFieldBackgroundMesh3D::eval_approximate_crossfield(const MVertex *vert ...@@ -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); vector<double>res = get_field_value(x,y,z,crossFieldVectorialSmoothness);
for (int i=0;i<3;i++) cf(i)=res[i]; 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); vector<double> res = get_nodal_value(vert,crossFieldVectorialSmoothness);
for (int i=0;i<3;i++) cf(i)=res[i]; for (int i=0;i<3;i++) cf(i)=res[i];
}; };
double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, double x,
double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, double x, double y, double z){ double y, double z)
{
vector<double>res = get_field_value(x,y,z,crossFieldVectorialSmoothness); vector<double>res = get_field_value(x,y,z,crossFieldVectorialSmoothness);
return res[idir]; return res[idir];
}; };
double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir,
double frameFieldBackgroundMesh3D::get_vectorial_smoothness(const int idir, const MVertex *vert){ const MVertex *vert)
{
vector<double> res = get_nodal_value(vert,crossFieldVectorialSmoothness); vector<double> res = get_nodal_value(vert,crossFieldVectorialSmoothness);
return res[idir]; return res[idir];
}; };
void frameFieldBackgroundMesh3D::build_vertex_to_element_table()
void frameFieldBackgroundMesh3D::build_vertex_to_element_table(){ {
GRegion *gr = dynamic_cast<GRegion*>(gf); GRegion *gr = dynamic_cast<GRegion*>(gf);
MElement *e; MElement *e;
MVertex *v; MVertex *v;
...@@ -726,28 +739,29 @@ void frameFieldBackgroundMesh3D::build_vertex_to_element_table(){ ...@@ -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); GRegion *gr = dynamic_cast<GRegion*>(gf);
return gr->getMeshElement(i); return gr->getMeshElement(i);
} }
unsigned int backgroundMesh3D::getNumMeshElements()const
unsigned int backgroundMesh3D::getNumMeshElements()const{ {
GRegion *gr = dynamic_cast<GRegion*>(gf); GRegion *gr = dynamic_cast<GRegion*>(gf);
return gr->getNumMeshElements(); return gr->getNumMeshElements();
} }
// this function fills the multimap "graph": vertex to direct neighbors, or
// this function fills the multimap "graph": vertex to direct neighbors, or indirect neighbors, depending on the "max_recursion_level". // indirect neighbors, depending on the "max_recursion_level".
void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level){ void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level)
{
set<MVertex*> visited,start; set<MVertex*> visited,start;
set<MElement*> visited_elements; set<MElement*> visited_elements;
multimap<int,MVertex*> proximity; multimap<int,MVertex*> proximity;
//int counter=0; //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; MVertex *current_vertex = it_vertex->first;
visited.clear(); visited.clear();
visited_elements.clear(); visited_elements.clear();
...@@ -781,16 +795,20 @@ void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level) ...@@ -781,16 +795,20 @@ void frameFieldBackgroundMesh3D::build_neighbors(const int &max_recursion_level)
} }
} }
void frameFieldBackgroundMesh3D::get_recursive_neighbors
void frameFieldBackgroundMesh3D::get_recursive_neighbors(set<MVertex*> &start, set<MVertex*> &visited, set<MElement*> &visited_elements, multimap<int,MVertex*> &proximity, int max_level, int level){ (set<MVertex*> &start, set<MVertex*> &visited, set<MElement*> &visited_elements,
multimap<int,MVertex*> &proximity, int max_level, int level)
{
level++; level++;
if (level>max_level) return; if (level>max_level) return;
set<MVertex*> new_vertices; 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; 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); vert2elemtype::iterator itfind = vert2elem.find(current);
if (itfind==vert2elem.end()) continue; if (itfind==vert2elem.end()) continue;
set<MElement*>::iterator it_elem = itfind->second.begin(); set<MElement*>::iterator it_elem = itfind->second.begin();
...@@ -815,12 +833,15 @@ void frameFieldBackgroundMesh3D::get_recursive_neighbors(set<MVertex*> &start, s ...@@ -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
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){ (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.; for (int i=0;i<3;i++) mean_axis(i)=0.;
multimap<double,MVertex*>::iterator it = itbegin; multimap<double,MVertex*>::iterator it = itbegin;
SVector3 rotation_axis; SVector3 rotation_axis;
...@@ -842,7 +863,8 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor ...@@ -842,7 +863,8 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor
neighbor = it->second; neighbor = it->second;
//ponderations_vec.push_back(1.); //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);
//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); itneighbor = crossField.find(neighbor);
STensor3 &neibcross = itneighbor->second; STensor3 &neibcross = itneighbor->second;
...@@ -882,9 +904,6 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor ...@@ -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 = (*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). 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 itan=all_angle.begin();
vector<double>::iterator itpond=ponderations_vec.begin(); vector<double>::iterator itpond=ponderations_vec.begin();
...@@ -898,13 +917,12 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor ...@@ -898,13 +917,12 @@ double frameFieldBackgroundMesh3D::compare_to_neighbors(SPoint3 current, STensor
mean_angle = mean_axis.norm()/pond_total; mean_angle = mean_axis.norm()/pond_total;
mean_axis.normalize(); mean_axis.normalize();
return smoothness_scalar; return smoothness_scalar;
} }
STensor3 frameFieldBackgroundMesh3D::apply_rotation
STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const double &angle, const STensor3 &thecross){ (const SVector3 &axis, const double &angle, const STensor3 &thecross)
{
double rotmat[3][3]; double rotmat[3][3];
STensor3 res2; STensor3 res2;
get_rotation_matrix(angle,axis,&rotmat[0][0]); get_rotation_matrix(angle,axis,&rotmat[0][0]);
...@@ -919,7 +937,6 @@ STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const ...@@ -919,7 +937,6 @@ STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const
// } // }
// //
for (int i=0;i<3;i++) for (int i=0;i<3;i++)
for (int j=0;j<3;j++) for (int j=0;j<3;j++)
for (int k=0;k<3;k++) for (int k=0;k<3;k++)
...@@ -985,9 +1002,11 @@ STensor3 frameFieldBackgroundMesh3D::apply_rotation(const SVector3 &axis, const ...@@ -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 ??? // TODO: ce genre de fct... mettre dans une classe FrameField ? Ou bien juste un
void frameFieldBackgroundMesh3D::get_rotation_matrix(const double &angle_to_go, const SVector3 &rotvec, double* rotmat){ // 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 c = cos(angle_to_go);
double s = sin(angle_to_go); double s = sin(angle_to_go);
double temp01 = rotvec[0]*rotvec[1]*(1.-c); double temp01 = rotvec[0]*rotvec[1]*(1.-c);
...@@ -1007,9 +1026,10 @@ void frameFieldBackgroundMesh3D::get_rotation_matrix(const double &angle_to_go, ...@@ -1007,9 +1026,10 @@ void frameFieldBackgroundMesh3D::get_rotation_matrix(const double &angle_to_go,
rotmat[8] = square2 +(1.-square2)*c; rotmat[8] = square2 +(1.-square2)*c;
} }
void frameFieldBackgroundMesh3D::get_min_rotation_matrix
void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, double threshold, bool debugflag){ (const STensor3 &reference, const STensor3 &thecross, double &minimum_angle,
SVector3 &rotation_axis, double threshold, bool debugflag)
{
minimum_angle = M_PI/2.; minimum_angle = M_PI/2.;
pair<int,int> p; pair<int,int> p;
...@@ -1021,7 +1041,9 @@ void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &referen ...@@ -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]); get_rotation_angle_and_axis(reference,thecross,a,v,i_rotation_perm,permutation[iperm]);
if (debugflag){ if (debugflag){
if (fabs(a)<M_PI/2.){ 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 else
cout << " temp parameters: angle=" << a*180./M_PI << endl; cout << " temp parameters: angle=" << a*180./M_PI << endl;
...@@ -1032,23 +1054,26 @@ void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &referen ...@@ -1032,23 +1054,26 @@ void frameFieldBackgroundMesh3D::get_min_rotation_matrix(const STensor3 &referen
rotation_axis = v; rotation_axis = v;
} }
if (minimum_angle<threshold) break; if (minimum_angle<threshold) break;
} }
} }
// cout << "pair=(" << p.first << "," << p.second << ")" << endl; // cout << "pair=(" << p.first << "," << p.second << ")" << endl;
if (debugflag){ 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; return;
} }
// compute the angle and axis of rotation between two (unit-orthogonal) cross fields // compute the angle and axis of rotation between two (unit-orthogonal) cross
// using crossfield vectors permutations defined by modulo and trip // fields using crossfield vectors permutations defined by modulo and trip The
// The rotation matrix between cross fields C and M is R = C M^T // 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){ 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 MT[3][3],C[3][3],R[3][3];
double C[3][3],R[3][3]; double C[3][3],R[3][3];
...@@ -1108,8 +1133,8 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &ref ...@@ -1108,8 +1133,8 @@ void frameFieldBackgroundMesh3D::get_rotation_angle_and_axis(const STensor3 &ref
return; return;
} }
void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filename)
void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filename){ {
FILE *f = Fopen (filename.c_str(),"w"); FILE *f = Fopen (filename.c_str(),"w");
fprintf(f,"View \"Background Mesh\"{\n"); fprintf(f,"View \"Background Mesh\"{\n");
...@@ -1126,12 +1151,11 @@ void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filenam ...@@ -1126,12 +1151,11 @@ void frameFieldBackgroundMesh3D::exportVectorialSmoothness(const string &filenam
eval_approximate_crossfield(v,cf); eval_approximate_crossfield(v,cf);
for (int i=0;i<3;i++){ for (int i=0;i<3;i++){
double vs = get_vectorial_smoothness(i,v); 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"); fprintf(f,"};\n");
fclose(f); fclose(f);
} }
...@@ -14,14 +14,7 @@ ...@@ -14,14 +14,7 @@
#include <iostream> #include <iostream>
#include <iomanip> #include <iomanip>
#include <numeric> #include <numeric>
//#include "SVector3.h"
//#include "STensor3.h"
#include "BGMBase.h" #include "BGMBase.h"
#include "pointInsertion.h" #include "pointInsertion.h"
#if defined(HAVE_ANN) #if defined(HAVE_ANN)
...@@ -29,15 +22,10 @@ ...@@ -29,15 +22,10 @@
class ANNkd_tree; class ANNkd_tree;
#endif #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
// 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... ?! // account vertices at a given topological distance... but using distance=1, it
// seems to work... ?!
using namespace std;
class montripletbis{ class montripletbis{
public: public:
...@@ -49,7 +37,7 @@ class montripletbis{ ...@@ -49,7 +37,7 @@ class montripletbis{
}; };
inline int operator()(int i)const{return vec[i];} inline int operator()(int i)const{return vec[i];}
private: private:
vector<int> vec; std::vector<int> vec;
}; };
...@@ -85,11 +73,11 @@ class backgroundMesh3D : public BGMBase { ...@@ -85,11 +73,11 @@ class backgroundMesh3D : public BGMBase {
class frameFieldBackgroundMesh3D : public backgroundMesh3D{ class frameFieldBackgroundMesh3D : public backgroundMesh3D{
public: public:
// typedef tr1::unordered_map<hash_key_ptr, set<MElement*> > vert2elemtype; // typedef tr1::unordered_map<hash_key_ptr, std::set<MElement*> > vert2elemtype;
// typedef tr1::unordered_map<MElement*, set<MVertex*> > elem2verttype; // typedef tr1::unordered_map<MElement*, std::set<MVertex*> > elem2verttype;
typedef std::map<hash_key_ptr, set<MElement*> > vert2elemtype; typedef std::map<hash_key_ptr, std::set<MElement*> > vert2elemtype;
typedef std::map<MElement*, set<MVertex*> > elem2verttype; typedef std::map<MElement*, std::set<MVertex*> > elem2verttype;
typedef multimap<MVertex*, pair<int,MVertex*> > graphtype; typedef std::multimap<MVertex*, std::pair<int,MVertex*> > graphtype;
protected: protected:
bool smooth_the_crossfield; bool smooth_the_crossfield;
...@@ -105,10 +93,10 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{ ...@@ -105,10 +93,10 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{
void build_neighbors(const int &max_recursion_level); void build_neighbors(const int &max_recursion_level);
// function called only by "build_neighbors", recursively recovering vertices neighbors. // 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); 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, multimap<double,MVertex*>::iterator itbegin, multimap<double,MVertex*>::iterator itend, SVector3 &mean_axis,double &mean_angle, vector<double> &vectorial_smoothness); 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); 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_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_min_rotation_matrix(const STensor3 &reference, const STensor3 &thecross, double &minimum_angle, SVector3 &rotation_axis, double threshold=-1., bool debugflag=false);
...@@ -122,11 +110,11 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{ ...@@ -122,11 +110,11 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{
vert2elemtype vert2elem; vert2elemtype vert2elem;
elem2verttype elem2vert; elem2verttype elem2vert;
set<MVertex*> listOfBndVertices; std::set<MVertex*> listOfBndVertices;
graphtype graph; graphtype graph;
double smoothness_threshold; double smoothness_threshold;
static vector<montripletbis> permutation; static std::vector<montripletbis> permutation;
#if defined(HAVE_ANN) #if defined(HAVE_ANN)
ANNkd_tree* annTree; ANNkd_tree* annTree;
...@@ -143,7 +131,7 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{ ...@@ -143,7 +131,7 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{
MVertex* get_nearest_neighbor_on_boundary(MVertex* v, double &distance); MVertex* get_nearest_neighbor_on_boundary(MVertex* v, double &distance);
// bool findvertex(const MVertex *v)const{ // bool findvertex(const MVertex *v)const{
// map<MVertex*,double>::const_iterator it = sizeField.find(const_cast<MVertex*>(v)); // std::map<MVertex*,double>::const_iterator it = sizeField.find(const_cast<MVertex*>(v));
// if (it==sizeField.end()) return false; // if (it==sizeField.end()) return false;
// return true; // return true;
// }; // };
...@@ -162,9 +150,9 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{ ...@@ -162,9 +150,9 @@ class frameFieldBackgroundMesh3D : public backgroundMesh3D{
double get_vectorial_smoothness(const int idir, double x, double y, double z); double get_vectorial_smoothness(const int idir, double x, double y, double z);
double get_vectorial_smoothness(const int idir, const MVertex *vert); double get_vectorial_smoothness(const int idir, const MVertex *vert);
inline void exportCrossField(const string &filename){export_tensor_as_vectors(filename,crossField);}; inline void exportCrossField(const std::string &filename){export_tensor_as_vectors(filename,crossField);};
inline void exportSmoothness(const string &filename) const{export_scalar(filename,crossFieldSmoothness);}; inline void exportSmoothness(const std::string &filename) const{export_scalar(filename,crossFieldSmoothness);};
void exportVectorialSmoothness(const string &filename); void exportVectorialSmoothness(const std::string &filename);
// private: // private:
// STensor3 get_random_cross()const; // STensor3 get_random_cross()const;
......
...@@ -206,7 +206,7 @@ public: ...@@ -206,7 +206,7 @@ public:
virtual bool empty(){ return points.empty(); }; virtual bool empty(){ return points.empty(); };
protected: protected:
set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points; std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
}; };
class listOfPointsVectorialSmoothness : public listOfPointsScalarSmoothness{ class listOfPointsVectorialSmoothness : public listOfPointsScalarSmoothness{
...@@ -218,7 +218,7 @@ public: ...@@ -218,7 +218,7 @@ public:
}; };
virtual SVector3 get_first_direction(){ return (*points.begin())->direction; }; virtual SVector3 get_first_direction(){ return (*points.begin())->direction; };
protected: protected:
set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points; std::set<smoothness_vertex_pair*, compareSmoothnessVertexPairs> points;
}; };
class listOfPointsFifo : public listOfPoints{ class listOfPointsFifo : public listOfPoints{
...@@ -248,7 +248,7 @@ public: ...@@ -248,7 +248,7 @@ public:
virtual bool empty(){ return points.empty(); }; virtual bool empty(){ return points.empty(); };
protected: protected:
queue<smoothness_vertex_pair*> points; std::queue<smoothness_vertex_pair*> points;
}; };
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment