From 74fc3b00c17204e95ed491c5e977b03877dd0d56 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 1 May 2016 06:43:48 +0000
Subject: [PATCH] fix compile
---
Geo/discreteDiskFace.cpp | 285 ++++++++++++++++++++-------------------
Geo/discreteDiskFace.h | 95 +++++++------
2 files changed, 190 insertions(+), 190 deletions(-)
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index eab809d66f..aab6839092 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -2,13 +2,13 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-
+
#include <queue> // #mark
#include <stdlib.h>
#include "GmshConfig.h"
-
+
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
-
+
#include "GmshMessage.h"
#include "Octree.h"
#include "discreteDiskFace.h"
@@ -26,18 +26,17 @@
#include "convexCombinationTerm.h" // #FIXME
#include "qualityMeasuresJacobian.h" // #temporary?
-
+
static inline void functionShapes(int p, double Xi[2], double* phi)
{
-
switch(p){
-
+
case 1:
phi[0] = 1. - Xi[0] - Xi[1];
phi[1] = Xi[0];
phi[2] = Xi[1];
break;
-
+
case 2:
phi[0] = 1. - 3.*(Xi[0]+Xi[1]) + 2.*(Xi[0]+Xi[1])*(Xi[0]+Xi[1]);
phi[1] = Xi[0]*(2.*Xi[0]-1);
@@ -46,27 +45,24 @@ static inline void functionShapes(int p, double Xi[2], double* phi)
phi[4] = 4.*Xi[0]*Xi[1];
phi[5] = 4.*Xi[1]*(1.-Xi[0]-Xi[1]);
break;
-
+
default:
- Msg::Error("(discreteDiskFace) static inline functionShapes, only first and second order available; order %d requested.",p);
+ Msg::Error("(discreteDiskFace) static inline functionShapes, only first "
+ "and second order available; order %d requested.", p);
break;
}
-
}
-
-
+
static inline void derivativeShapes(int p, double Xi[2], double phi[6][2])
{
-
switch(p){
-
+
case 1:
phi[0][0] = -1. ; phi[0][1] = -1.;
phi[1][0] = 1. ; phi[1][1] = 0.;
phi[2][0] = 0. ; phi[2][1] = 1.;
-
break;
-
+
case 2:
phi[0][0] = -3. + 4.*(Xi[0]+Xi[1]) ; phi[0][1] = -3. + 4.*(Xi[0]+Xi[1]);
phi[1][0] = 4.*Xi[0] - 1. ; phi[1][1] = 0.;
@@ -75,16 +71,16 @@ static inline void derivativeShapes(int p, double Xi[2], double phi[6][2])
phi[4][0] = 4.*Xi[1] ; phi[4][1] = 4.*Xi[0];
phi[5][0] = -4.*Xi[1] ; phi[5][1] = 4. - 4.*Xi[0] - 8.*Xi[1];
break;
-
+
default:
- Msg::Error("(discreteDiskFace) static inline derivativeShapes, only first and second order available; order %d requested.",p);
+ Msg::Error("(discreteDiskFace) static inline derivativeShapes, only first and "
+ "second order available; order %d requested.",p);
break;
}
-
}
-static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double Xi[2]){
-
+static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double Xi[2])
+{
double M[2][2], R[2];
const SPoint3 p0 = my_ddft->p[0];
const SPoint3 p1 = my_ddft->p[1];
@@ -96,18 +92,18 @@ static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double
R[0] = (U[0] - p0.x());
R[1] = (U[1] - p0.y());
sys2x2(M, R, Xi);
-
+
if (my_ddft->tri->getPolynomialOrder() == 2){
-
+
int iter = 1, maxiter = 20;
double error = 1., tol = 1.e-6;
while (error > tol && iter < maxiter){// Newton-Raphson
-
+
double fs[6];// phi_i
functionShapes(2,Xi,fs);
double ds[6][2];// dPhi_i/dXi_j
derivativeShapes(2,Xi,ds);
-
+
double un[2] = {0.,0.};
double jac[2][2] = {{0.,0.},{0.,0.}}; // d(u,v)/d(xi,eta)
for (int i=0; i<6; i++){
@@ -118,7 +114,7 @@ static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double
jac[k][j] += ui[k] * ds[i][j];
}
}
-
+
double inv[2][2];
inv2x2(jac,inv);
double xin[2] = {Xi[0],Xi[1]};
@@ -128,12 +124,12 @@ static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double
xin[i] += inv[i][j] * (U[j] - un[j]);
error += (xin[i] - Xi[i])*(xin[i] - Xi[i]);
}
-
+
error = sqrt(error);
Xi[0] = xin[0];
Xi[1] = xin[1];
iter++;
-
+
} // end Newton-Raphson
if (iter == maxiter){
return false;
@@ -141,18 +137,17 @@ static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double
}// end order 2
return true;
}
-
-
+
// The three following things are mandatory to manipulate octrees (octree in (u;v)).
static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
-{
+{
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
-
+
mmin[0] = std::min(std::min(t->p[0].x(), t->p[1].x()), t->p[2].x());
mmin[1] = std::min(std::min(t->p[0].y(), t->p[1].y()), t->p[2].y());
mmax[0] = std::max(std::max(t->p[0].x(), t->p[1].x()), t->p[2].x());
mmax[1] = std::max(std::max(t->p[0].y(), t->p[1].y()), t->p[2].y());
-
+
if (t->tri->getPolynomialOrder() == 2){
for (int i=3; i<6; i++){
int j = i==5 ? 0 : i-2;
@@ -167,7 +162,7 @@ static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
mmin[2] = -1;
mmax[2] = 1;
}
-
+
static void discreteDiskFaceCentroid(void *a, double*c)
{
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
@@ -175,10 +170,9 @@ static void discreteDiskFaceCentroid(void *a, double*c)
c[1] = (t->p[0].y() + t->p[1].y() + t->p[2].y()) / 3.0;
c[2] = 0.0;
}
-
+
static int discreteDiskFaceInEle(void *a, double*c)// # mark
-{
-
+{
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
double Xi[2];
double U[2] = {c[0],c[1]};
@@ -187,38 +181,38 @@ static int discreteDiskFaceInEle(void *a, double*c)// # mark
if(ok && Xi[0] > -eps && Xi[1] > -eps && 1. - Xi[0] - Xi[1] > -eps)
return 1;
-
+
return 0;
}
-
-static bool orderVertices(const double &tot_length, const std::vector<MVertex*> &l, std::vector<double> &coord)
-{ // called once by constructor ; organize the vertices for the linear system expressing the mapping
+
+static bool orderVertices(const double &tot_length, const std::vector<MVertex*> &l,
+ std::vector<double> &coord)
+{ // called once by constructor ; organize the vertices for the linear system
+ // expressing the mapping
coord.clear();
coord.push_back(0.);
-
+
MVertex* first = l[0];
-
+
for(unsigned int i=1; i < l.size(); i++){
-
+
MVertex* next = l[i];
-
+
const double length = sqrt( (next->x() - first->x()) * (next->x() - first->x()) +
(next->y() - first->y()) * (next->y() - first->y()) +
(next->z() - first->z()) * (next->z() - first->z()) );
-
+
coord.push_back(coord[coord.size()-1] + length / tot_length);
-
+
first = next;
-
+
}
-
return true;
-
}
-
-/*BUILDER*/
-discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation, int p, std::vector<GFace*> *CAD) :
- GFace(gf->model(),123), _parent (gf),_ddft(NULL)
+
+discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
+ int p, std::vector<GFace*> *CAD) :
+ GFace(gf->model(),123), _parent (gf), _ddft(NULL), oct(NULL)
{
initialTriangulation = diskTriangulation;
std::vector<MElement*> mesh = diskTriangulation->tri;
@@ -229,10 +223,10 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
std::map<MVertex*,MVertex*> v2v;// mesh vertex |-> face vertex
std::map<MEdge,MVertex*,Less_Edge> ed2nodes; // edge to interior node(s)
for (unsigned int i=0;i<mesh.size();i++){ // triangle by triangle
-
+
std::vector<MVertex*> vs; // MTriangle vertices
for (unsigned int j=0; j<3; j++){ // loop over vertices AND edges of the current triangle
-
+
MVertex *v = mesh[i]->getVertex(j);// firstly, edge vertices
if (v->onWhat()->dim() == 2) {
std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v);
@@ -241,7 +235,8 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
if (!CAD) vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), 0, 0);
else{
GFace *cad = (*CAD)[i];
- if(cad != v->onWhat())Msg::Fatal("Line %d FILE %s : erroneous cad list",__LINE__,__FILE__);
+ if(cad != v->onWhat())
+ Msg::Fatal("Line %d FILE %s : erroneous cad list",__LINE__,__FILE__);
double pu,pv; v->getParameter(0,pu);v->getParameter(1,pv);
vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), pu, pv);
}
@@ -290,20 +285,21 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
orderVertices(_totLength, _U0, _coords);
-
+
parametrize(false);
buildOct(CAD);
if (!checkOrientationUV()){
- Msg::Info("discreteDIskFace:: parametrization is not one-to-one; fixing the discrete system.");
+ Msg::Info("discreteDIskFace:: parametrization is not one-to-one; fixing "
+ "the discrete system.");
parametrize(true);
buildOct(CAD);
}
-
+
putOnView();
}
-discreteDiskFace::~discreteDiskFace()
+discreteDiskFace::~discreteDiskFace()
{
triangles.clear();
if (_ddft) delete[] _ddft;
@@ -320,7 +316,7 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
const int maxElePerBucket = 15;
oct = Octree_Create(maxElePerBucket, origin, ssize, discreteDiskFaceBB,
discreteDiskFaceCentroid, discreteDiskFaceInEle);
-
+
_ddft = new discreteDiskFaceTriangle[discrete_triangles.size()];
// if (CAD) printf("-------------> %ld %ld\n",CAD->size(),discrete_triangles.size());
@@ -331,16 +327,15 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
my_ddft->p.resize(_N);
for(int io=0; io<_N; io++)
my_ddft->p[io] = coordinates[t->getVertex(io)];
-
+
my_ddft->gf = CAD ? (*CAD)[i] : _parent;
my_ddft->tri = t;
-
+
Octree_Insert(my_ddft, oct);
}
Octree_Arrange(oct);
}
-
bool discreteDiskFace::parametrize(bool one2one) const
{ // #improveme
@@ -351,20 +346,20 @@ bool discreteDiskFace::parametrize(bool one2one) const
linearSystem<double> *lsys_u;
linearSystem<double> *lsys_v;
-
+
lsys_u = new linearSystemCSRTaucs<double>; // taucs :-)
lsys_v = new linearSystemCSRTaucs<double>; // taucs :-)
-
+
dofManager<double> myAssemblerU(lsys_u); // hashing
dofManager<double> myAssemblerV(lsys_v);
-
+
for(size_t i = 0; i < _U0.size(); i++){
MVertex *v = _U0[i];
const double theta = 2 * M_PI * _coords[i];
myAssemblerU.fixVertex(v, 0, 1, cos(theta));
myAssemblerV.fixVertex(v, 0, 1, sin(theta));
}
-
+
for(size_t i = 0; i < discrete_triangles.size(); ++i){
MElement *t = discrete_triangles[i];
for(int j=0; j<t->getNumVertices(); j++){
@@ -372,20 +367,19 @@ bool discreteDiskFace::parametrize(bool one2one) const
myAssemblerV.numberVertex(t->getVertex(j), 0, 1);
}
}
-
-
+
Msg::Debug("Creating term %d dofs numbered %d fixed",
- myAssemblerU.sizeOfR() + myAssemblerV.sizeOfR(), myAssemblerU.sizeOfF() + myAssemblerV.sizeOfF());
-
+ myAssemblerU.sizeOfR() + myAssemblerV.sizeOfR(),
+ myAssemblerU.sizeOfF() + myAssemblerV.sizeOfF());
+
double t1 = Cpu();
-
+
simpleFunction<double> ONE(1.0);
-
if (one2one){
convexLaplaceTerm mappingU(0, 1, &ONE);
convexLaplaceTerm mappingV(0, 1, &ONE);
-
+
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
SElement se(discrete_triangles[i]);
mappingU.addToMatrix(myAssemblerU, &se);
@@ -395,7 +389,7 @@ bool discreteDiskFace::parametrize(bool one2one) const
else{
laplaceTerm mappingU(0, 1, &ONE);
laplaceTerm mappingV(0, 1, &ONE);
-
+
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
SElement se(discrete_triangles[i]);
mappingU.addToMatrix(myAssemblerU, &se);
@@ -408,7 +402,7 @@ bool discreteDiskFace::parametrize(bool one2one) const
lsys_u->systemSolve();
lsys_v->systemSolve();
Msg::Debug("Systems solved");
-
+
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value_U, value_V;
@@ -426,21 +420,20 @@ bool discreteDiskFace::parametrize(bool one2one) const
itf->second[1]= value_V;
}
}
-
+
delete lsys_u;
delete lsys_v;
-
+
return true;
-
}
-
+
void discreteDiskFace::getTriangleUV(const double u,const double v,
discreteDiskFaceTriangle **mt,
double &_xi, double &_eta)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
if (!(*mt)) return;
-
+
double Xi[2];
double U[2] = {u,v};
uv2xi(*mt,U,Xi);
@@ -448,9 +441,8 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
_eta = Xi[1];
}
-
-bool discreteDiskFace::checkOrientationUV(){
-
+bool discreteDiskFace::checkOrientationUV()
+{
discreteDiskFaceTriangle *ct;
if(_order==1){
@@ -474,27 +466,27 @@ bool discreteDiskFace::checkOrientationUV(){
}
return true;
}
- else{
+ else{
double min, max;
std::vector<MVertex*> localVertices;
localVertices.resize(_N);
- for(unsigned int i=0; i<discrete_triangles.size(); i++){
- ct = &_ddft[i];
+ for(unsigned int i=0; i<discrete_triangles.size(); i++){
+ ct = &_ddft[i];
for(int j=0; j<_N; j++)
localVertices[j] = new MVertex(ct->p[j].x(),ct->p[j].y(),0.);
MTriangle6 mt6(localVertices);
jacobianBasedQuality::minMaxJacobianDeterminant(&mt6,min,max);
for(int j=0; j<_N; j++)
- delete localVertices[j];
+ delete localVertices[j];
if (min < 0){
- return false;
+ return false;
break;
}
}
return true;
}
}
-
+
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
@@ -507,7 +499,7 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const
GPoint gp = GPoint(1.e22,1.e22,1.e22,_parent,par);
gp.setNoSuccess();
return gp;
- }
+ }
std::vector<double> eval(_N);
functionShapes(_order,Xi,&eval[0]);
@@ -519,7 +511,7 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const
}
return GPoint(X,Y,Z,_parent,par);
}
-
+
SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
{
if (v->onWhat() == this){
@@ -528,7 +520,7 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
v->getParameter(1,vv);
return SPoint2(uu,vv);
}
-
+
std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
// The 1D mesh has been re-done
@@ -550,80 +542,81 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
}
else if (v->onWhat()->dim()==0)
- Msg::Fatal("discreteDiskFace::parFromVertex vertex classified on a model vertex that is not part of the face");
-
+ Msg::Fatal("discreteDiskFace::parFromVertex vertex classified on a model "
+ "vertex that is not part of the face");
+
return SPoint2(0,0);
}
-
+
SVector3 discreteDiskFace::normal(const SPoint2 ¶m) const
{
return GFace::normal(param);
}
-
+
double discreteDiskFace::curvatureMax(const SPoint2 ¶m) const
{
throw;
return false;
}
-
+
double discreteDiskFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const
{
throw;
return false;
}
-
+
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
{
double xi,eta;
discreteDiskFaceTriangle* ddft;
getTriangleUV(param.x(),param.y(),&ddft,xi,eta);
-
+
MElement* tri = NULL;
if (ddft) tri = ddft->tri;
else {
Msg::Warning("discreteDiskFace::firstDer << triangle not found %g %g",param[0],param[1]);
return Pair<SVector3, SVector3>(SVector3(1,0,0), SVector3(0,1,0));
}
-
+
double Xi[2] = {xi,eta};
double df[6][2];
derivativeShapes(_order,Xi,df);
double dxdxi[3][2] = {{0.,0.},{0.,0.},{0.,0.}};
-
+
double dudxi[2][2] = {{0.,0.},{0.,0.}};
-
+
for (int io=0; io<_N; io++){
-
+
double X = tri->getVertex(io)->x();
double Y = tri->getVertex(io)->y();
double Z = tri->getVertex(io)->z();
-
+
double U = ddft->p[io].x();
double V = ddft->p[io].y();
-
+
dxdxi[0][0] += X*df[io][0];
dxdxi[0][1] += X*df[io][1];
-
+
dxdxi[1][0] += Y*df[io][0];
dxdxi[1][1] += Y*df[io][1];
-
+
dxdxi[2][0] += Z*df[io][0];
dxdxi[2][1] += Z*df[io][1];
-
+
dudxi[0][0] += U*df[io][0];
dudxi[0][1] += U*df[io][1];
-
+
dudxi[1][0] += V*df[io][0];
dudxi[1][1] += V*df[io][1];
-
+
}
-
+
double dxidu[2][2];
inv2x2(dudxi,dxidu);
-
+
double dxdu[3][2];
-
+
for (int i=0;i<3;i++){
for (int j=0;j<2;j++){
dxdu[i][j]=0.;
@@ -632,36 +625,40 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
}
}
}
-
+
return Pair<SVector3, SVector3>(SVector3(dxdu[0][0],dxdu[1][0],dxdu[2][0]),
SVector3(dxdu[0][1],dxdu[1][1],dxdu[2][1]));
}
-
+
void discreteDiskFace::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
-{ // cf Sauvage's thesis
+{ // cf Sauvage's thesis
return;
}
-
void discreteDiskFace::putOnView()
{
char mybuffer [64];
- snprintf(mybuffer,sizeof(mybuffer),"param_u_part%d_order%d.pos",initialTriangulation->idNum,_order);
+ snprintf(mybuffer,sizeof(mybuffer),"param_u_part%d_order%d.pos",
+ initialTriangulation->idNum,_order);
FILE* view_u = Fopen(mybuffer,"w");
- snprintf(mybuffer,sizeof(mybuffer),"param_v_part%d_order%d.pos",initialTriangulation->idNum,_order);
+ snprintf(mybuffer,sizeof(mybuffer),"param_v_part%d_order%d.pos",
+ initialTriangulation->idNum,_order);
FILE* view_v = Fopen(mybuffer,"w");
-
- snprintf(mybuffer,sizeof(mybuffer),"UVx_part%d_order%d.pos",initialTriangulation->idNum,_order);
+
+ snprintf(mybuffer,sizeof(mybuffer),"UVx_part%d_order%d.pos",
+ initialTriangulation->idNum,_order);
FILE* UVx = Fopen(mybuffer,"w");
- snprintf(mybuffer,sizeof(mybuffer),"UVy_part%d_order%d.pos",initialTriangulation->idNum,_order);
+ snprintf(mybuffer,sizeof(mybuffer),"UVy_part%d_order%d.pos",
+ initialTriangulation->idNum,_order);
FILE* UVy = Fopen(mybuffer,"w");
- snprintf(mybuffer,sizeof(mybuffer),"UVz_part%d_order%d.pos",initialTriangulation->idNum,_order);
+ snprintf(mybuffer,sizeof(mybuffer),"UVz_part%d_order%d.pos",
+ initialTriangulation->idNum,_order);
FILE* UVz = Fopen(mybuffer,"w");
if(view_u && view_v && UVx && UVy && UVz){
@@ -676,7 +673,7 @@ void discreteDiskFace::putOnView()
fprintf(view_u,"ST(");
fprintf(view_v,"ST(");
fprintf(UVx,"ST(");
- fprintf(UVy,"ST(");
+ fprintf(UVy,"ST(");
fprintf(UVz,"ST(");
}
else{
@@ -687,14 +684,18 @@ void discreteDiskFace::putOnView()
fprintf(UVz,"ST%d(",_order);
}
for (int j=0; j<_N-1; j++){
- fprintf(view_u,"%g,%g,%g,",my_ddft->tri->getVertex(j)->x(),my_ddft->tri->getVertex(j)->y(),my_ddft->tri->getVertex(j)->z());
- fprintf(view_v,"%g,%g,%g,",my_ddft->tri->getVertex(j)->x(),my_ddft->tri->getVertex(j)->y(),my_ddft->tri->getVertex(j)->z());
+ fprintf(view_u,"%g,%g,%g,",my_ddft->tri->getVertex(j)->x(),
+ my_ddft->tri->getVertex(j)->y(),my_ddft->tri->getVertex(j)->z());
+ fprintf(view_v,"%g,%g,%g,",my_ddft->tri->getVertex(j)->x(),
+ my_ddft->tri->getVertex(j)->y(),my_ddft->tri->getVertex(j)->z());
fprintf(UVx,"%g,%g,%g,",my_ddft->p[j].x(),my_ddft->p[j].y(),0.);
fprintf(UVy,"%g,%g,%g,",my_ddft->p[j].x(),my_ddft->p[j].y(),0.);
fprintf(UVz,"%g,%g,%g,",my_ddft->p[j].x(),my_ddft->p[j].y(),0.);
}
- fprintf(view_u,"%g,%g,%g){",my_ddft->tri->getVertex(_N-1)->x(),my_ddft->tri->getVertex(_N-1)->y(),my_ddft->tri->getVertex(_N-1)->z());
- fprintf(view_v,"%g,%g,%g){",my_ddft->tri->getVertex(_N-1)->x(),my_ddft->tri->getVertex(_N-1)->y(),my_ddft->tri->getVertex(_N-1)->z());
+ fprintf(view_u,"%g,%g,%g){",my_ddft->tri->getVertex(_N-1)->x(),
+ my_ddft->tri->getVertex(_N-1)->y(),my_ddft->tri->getVertex(_N-1)->z());
+ fprintf(view_v,"%g,%g,%g){",my_ddft->tri->getVertex(_N-1)->x(),
+ my_ddft->tri->getVertex(_N-1)->y(),my_ddft->tri->getVertex(_N-1)->z());
fprintf(UVx,"%g,%g,%g){",my_ddft->p[_N-1].x(),my_ddft->p[_N-1].y(),0.);
fprintf(UVy,"%g,%g,%g){",my_ddft->p[_N-1].x(),my_ddft->p[_N-1].y(),0.);
fprintf(UVz,"%g,%g,%g){",my_ddft->p[_N-1].x(),my_ddft->p[_N-1].y(),0.);
@@ -744,8 +745,8 @@ void discreteDiskFace::putOnView()
#endif
*/
}
-
-// useful for mesh generators ----------------------------------------
+
+// useful for mesh generators
// Intersection of a circle and a plane
GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
@@ -765,9 +766,12 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
// the point is situated in the half plane defined
// by direction n1 and point p (exclude the one on the
// other side)
- SVector3 v0(ct->tri->getVertex(0)->x(),ct->tri->getVertex(0)->y(),ct->tri->getVertex(0)->z());
- SVector3 v1(ct->tri->getVertex(1)->x(),ct->tri->getVertex(1)->y(),ct->tri->getVertex(1)->z());
- SVector3 v2(ct->tri->getVertex(2)->x(),ct->tri->getVertex(2)->y(),ct->tri->getVertex(2)->z());
+ SVector3 v0(ct->tri->getVertex(0)->x(),ct->tri->getVertex(0)->y(),
+ ct->tri->getVertex(0)->z());
+ SVector3 v1(ct->tri->getVertex(1)->x(),ct->tri->getVertex(1)->y(),
+ ct->tri->getVertex(1)->z());
+ SVector3 v2(ct->tri->getVertex(2)->x(),ct->tri->getVertex(2)->y(),
+ ct->tri->getVertex(2)->z());
SVector3 t1 = v1 - v0;
SVector3 t2 = v2 - v0;
// let us first compute point q to be the intersection
@@ -799,7 +803,7 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
// they correspond to two points s_1 = q + ta m and s_2 = q + tb m
// we choose the one that corresponds to (s_i - p) . n1 > 0
SVector3 m = crossprod(w,n);
-
+
const double a = dot(m,m);
const double b = 2*dot(m,q-p);
const double c = dot(q-p,q-p) - d*d;
@@ -817,14 +821,14 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
s = s2;
}
else continue;
-
+
// we have now to look if the point is inside the triangle
// s = v1 + u t1 + v t2
// we know the system has a solution because s is in the plane
// defined by v1 and v2 we solve then
// (s - v1) . t1 = u t1^2 + v t2 . t1
// (s - v1) . t2 = u t1 . t2 + v t2^2
-
+
double mat[2][2], b[2],uv[2];
mat[0][0] = dot(t1,t1);
mat[1][1] = dot(t2,t2);
@@ -848,6 +852,5 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
Msg::Debug("ARGG no success intersection circle");
return pp;
}
-
-
+
#endif
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 56fbb1fb36..20124e03a9 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -2,14 +2,14 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-
+
#ifndef _DISCRETE_DISK_FACE_H_
#define _DISCRETE_DISK_FACE_H_
-
+
#include "GmshConfig.h"
-
+
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
-
+
#include <list>
#include <map>
#include "GModel.h"
@@ -23,27 +23,25 @@
#include "PView.h"
#include "robustPredicates.h"
-/*inline utilities*/
-inline int nodeLocalNum(MElement* e, MVertex* v) {// const
+inline int nodeLocalNum(MElement* e, MVertex* v)
+{
for(int i=0; i<e->getNumVertices(); i++)
if (v == e->getVertex(i))
return i;
return -1;
}
-inline int edgeLocalNum(MElement* e, MEdge ed) {// const
+
+inline int edgeLocalNum(MElement* e, MEdge ed)
+{
for(int i=0; i<e->getNumEdges(); i++)
if (ed == e->getEdge(i))
return i;
return -1;
}
-
-
-/*various classes*/
class ANNkd_tree;
class Octree;
class GRbf;
-
class triangulation {
@@ -52,18 +50,20 @@ class triangulation {
// attributes
std::vector<MElement*> tri;// triangles
std::set<MVertex*> vert;// nodes
- std::map<MEdge,std::vector<int>,Less_Edge> ed2tri; // edge to 1 or 2 triangle(s), their num into the vector of MElement*
+ // edge to 1 or 2 triangle(s), their num into the vector of MElement*
+ std::map<MEdge,std::vector<int>,Less_Edge> ed2tri;
std::map<double,std::vector<MVertex*> > bord; //border(s)
std::set<MEdge,Less_Edge> borderEdg; // border edges
GFace *gf;
int idNum; // number of identification, for hashing purposes
-
- //methods
- int genus(){
+
+ int genus()
+ {
return ( -vert.size() + ed2tri.size() - tri.size() + 2 - bord.size() )/2;
}
- void assignVert(){
+ void assignVert()
+ {
for(unsigned int i = 0; i < tri.size(); ++i){
MElement* t = tri[i];
for(int j = 0; j < t->getNumVertices() ; j++){
@@ -74,7 +74,8 @@ class triangulation {
}
}
- void assignEd2tri(){
+ void assignEd2tri()
+ {
for(unsigned int i = 0; i < tri.size(); ++i){
MElement *t = tri[i];
for(int j = 0; j < 3 ; j++){
@@ -84,7 +85,8 @@ class triangulation {
}
}
- void assignBord(){
+ void assignBord()
+ {
for(unsigned int i = 0; i < tri.size(); ++i){
MElement *t = tri[i];
for(int j = 0; j < t->getNumEdges() ; j++){
@@ -109,75 +111,73 @@ class triangulation {
vecver.erase(vecver.begin());
std::map<MVertex*,std::vector<MVertex*> >::iterator im = firstNode2Edge.find(first);
- if (im != firstNode2Edge.end()) Msg::Fatal("Incorrect topology in discreteDiskFace %d", gf->tag());
+ if (im != firstNode2Edge.end())
+ Msg::Error("Incorrect topology in discreteDiskFace %d", gf->tag());
firstNode2Edge[first] = vecver;
firstNode2Edge[first].push_back(last);
}
- while (!firstNode2Edge.empty()) {
+ while (!firstNode2Edge.empty()) {
std::vector<MVertex*> loop;
double length = 0.;
- std::map<MVertex*,std::vector<MVertex*> >::iterator in = firstNode2Edge.begin();
+ std::map<MVertex*,std::vector<MVertex*> >::iterator in = firstNode2Edge.begin();
MVertex* previous = in->first;
while(in != firstNode2Edge.end()) { // it didn't find it
-
-
std::vector<MVertex*> myV = in->second;
for(unsigned int i=0; i<myV.size(); i++){
-
loop.push_back(previous);
-
- MVertex* current = myV[i];
-
+ MVertex* current = myV[i];
length += sqrt( (current->x()-previous->x()) * (current->x()-previous->x()) +
(current->y()-previous->y()) * (current->y()-previous->y()) +
(current->z()-previous->z()) * (current->z()-previous->z()) );
-
+
previous = current;
}
firstNode2Edge.erase(in);
in = firstNode2Edge.find(previous);
}// end while in
- bord.insert(std::make_pair(length,loop)); // it shouldn't be possible to have twice the same length ? actually, it is possible, but quite seldom #fixme ----> multimap ?
+ bord.insert(std::make_pair(length,loop));
+ // it shouldn't be possible to have twice the same length ? actually, it
+ // is possible, but quite seldom #fixme ----> multimap ?
} // end while firstNode2Edge
}// end method
- void assign(){
+ void assign()
+ {
assignVert();
assignEd2tri();
assignBord();
}
- //builder
- triangulation() : gf(0) {}
- triangulation(std::vector<MElement*> input, GFace* gface) : tri(input), gf(gface) {assign();}
+ triangulation() : gf(0) {}
+ triangulation(std::vector<MElement*> input, GFace* gface)
+ : tri(input), gf(gface) { assign(); }
};
-// --------------------
-// triangles in the physical space xyz, with their parametric coordinates
+// triangles in the physical space xyz, with their parametric coordinates
class discreteDiskFaceTriangle {
public:
std::vector<SPoint3> p; // vertices in (u;v) #improveme std::vector instead
//SPoint2 gfp[6]; // CAD model
GFace *gf; // GFace tag
MElement *tri; // mesh triangle in (x;y;z)
- discreteDiskFaceTriangle() : gf(0), tri(0) {}
+ discreteDiskFaceTriangle() : gf(0), tri(0) {}
};
-
-// --------------------
-
+
class discreteDiskFace : public GFace {
GFace *_parent;
void buildOct(std::vector<GFace*> *CAD = NULL) const;
bool parametrize(bool one2one = false) const;
void putOnView();
bool checkOrientationUV();
-
+
public:
- discreteDiskFace(GFace *parent, triangulation* diskTriangulation, int p=1, std::vector<GFace*> *CAD = NULL);// BUILDER
+ discreteDiskFace(GFace *parent, triangulation* diskTriangulation,
+ int p=1, std::vector<GFace*> *CAD = NULL);
virtual ~discreteDiskFace();
- void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v) const;
+ void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt,
+ double &_u, double &_v) const;
GPoint point(double par1, double par2) const;
SPoint2 parFromVertex(MVertex *v) const;
SVector3 normal(const SPoint2&) const;
@@ -191,15 +191,12 @@ class discreteDiskFace : public GFace {
const SVector3 &p, const double &d,
double uv[2]) const;
protected:
- //------------------------------------------------
// a copy of the mesh that should not be destroyed
triangulation* initialTriangulation;
triangulation* geoTriangulation;// parametrized triangulation
- std::vector<MElement*> discrete_triangles;
+ std::vector<MElement*> discrete_triangles;
std::vector<MVertex*> discrete_vertices;
- //------------------------------------------------
- //-----------------------------------------------
int _order;
int _N;// number of dof's for a triangle
double _totLength;
@@ -211,10 +208,10 @@ class discreteDiskFace : public GFace {
mutable v2t_cont adjv;
mutable std::map<MVertex*, Pair<SVector3,SVector3> > firstDerivatives;
mutable discreteDiskFaceTriangle *_ddft;
- mutable Octree *oct = NULL;
+ mutable Octree *oct;
mutable std::vector<double> _coords;
};
-
+
#endif
-
+
#endif
--
GitLab