From cef2e611bd31ad60acdc132f039d376175900883 Mon Sep 17 00:00:00 2001
From: Pierre-Alexandre Beaufort <pierre-alexandre.beaufort@uclouvain.be>
Date: Sat, 30 Apr 2016 13:42:57 +0000
Subject: [PATCH] Partioning of the initial triangulation for reparametrization
(disk-like topology). The mesher cannot be used yet, but the parametrization
does work ! Todo: partioning for large aspect ratio, check if each part is
connected, make discreteEdges (for mesher)
---
Geo/discreteDiskFace.cpp | 565 ++++++++++++++++-----------------------
Geo/discreteDiskFace.h | 182 ++++++++++---
Geo/discreteFace.cpp | 208 +++++++++++++-
Geo/discreteFace.h | 12 +-
4 files changed, 589 insertions(+), 378 deletions(-)
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 7262b86987..eab809d66f 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -1,15 +1,14 @@
-
// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-
+
#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"
@@ -24,20 +23,21 @@
#include "dofManager.h"
#include "laplaceTerm.h"
#include "convexLaplaceTerm.h"
-#include "convexCombinationTerm.h" // #
-#include "BasisFactory.h"
+#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 +46,27 @@ 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);
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,12 +75,12 @@ 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);
break;
}
-
+
}
static inline bool uv2xi(discreteDiskFaceTriangle* my_ddft, double U[2], double Xi[2]){
@@ -96,18 +96,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 +118,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,34 +128,31 @@ 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;
}
- else{
- return true;
- }
}// end order 2
return true;
}
-
-
-// The three things that are mandatory to manipulate octrees (octree in (u;v)).
+
+
+// The three following things are mandatory to manipulate octrees (octree in (u;v)).
static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
-{ // called once by buildOct()
+{
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;
@@ -170,18 +167,18 @@ static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
mmin[2] = -1;
mmax[2] = 1;
}
-
+
static void discreteDiskFaceCentroid(void *a, double*c)
-{ // called once by buildOct()
+{
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
c[0] = (t->p[0].x() + t->p[1].x() + t->p[2].x()) / 3.0;
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
-{ // called once by buildOct()
-
+{
+
discreteDiskFaceTriangle *t = (discreteDiskFaceTriangle *)a;
double Xi[2];
double U[2] = {c[0],c[1]};
@@ -190,49 +187,52 @@ 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
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, std::vector<MTriangle*> &mesh, int p, std::vector<GFace*> *CAD) :
- GFace(gf->model(),123), _parent (gf),_ddft(NULL), oct(NULL)
+discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation, int p, std::vector<GFace*> *CAD) :
+ GFace(gf->model(),123), _parent (gf),_ddft(NULL)
{
+ initialTriangulation = diskTriangulation;
+ std::vector<MElement*> mesh = diskTriangulation->tri;
_order = p;
_N = (p+1)*(p+2)/2;
+ discrete_triangles.resize(mesh.size());
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);
@@ -277,19 +277,25 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, in
}
}// end order == 2
if (_order==1)
- discrete_triangles.push_back(new MTriangle (vs));
+ discrete_triangles[i] = new MTriangle (vs);
else if (_order==2)
- discrete_triangles.push_back(new MTriangle6 (vs));
+ discrete_triangles[i] = new MTriangle6 (vs);
}// end loop over triangles
- buildAllNodes();
- getBoundingEdges();
+ geoTriangulation = new triangulation(discrete_triangles,gf);
+
+ allNodes = geoTriangulation->vert;
+ _U0 = geoTriangulation->bord.rbegin()->second;
+ _totLength = geoTriangulation->bord.rbegin()->first;
+
+
orderVertices(_totLength, _U0, _coords);
parametrize(false);
buildOct(CAD);
if (!checkOrientationUV()){
+ Msg::Info("discreteDIskFace:: parametrization is not one-to-one; fixing the discrete system.");
parametrize(true);
buildOct(CAD);
}
@@ -300,166 +306,11 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, in
discreteDiskFace::~discreteDiskFace()
{
triangles.clear();
- if (_ddft)delete[] _ddft;
- if (oct)Octree_Delete(oct);
+ if (_ddft) delete[] _ddft;
+ if (oct) Octree_Delete(oct);
+ if(geoTriangulation) delete geoTriangulation;
}
-void discreteDiskFace::getBoundingEdges()
-{
-
- // first of all, all the triangles have to be oriented in the same way
- std::map<MEdge,std::vector<MElement*>,Less_Edge> ed2tri; // edge to 1 or 2 triangle(s)
-
- for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- MElement *e = discrete_triangles[i];
- for(int j = 0; j < e->getNumEdges() ; j++){
- MEdge ed = e->getEdge(j);
- ed2tri[ed].push_back(e);
- }
- }
-
- // element to its neighbors
- std::map<MElement*,std::vector<MElement*> > neighbors;
- for (unsigned int i=0; i<discrete_triangles.size(); ++i){
- MElement* e = discrete_triangles[i];
- for(int j=0; j<e->getNumEdges(); j++){ // #improveme: efficiency could be improved by setting neighbors mutually
- std::vector<MElement*> my_mt = ed2tri[e->getEdge(j)];
- if (my_mt.size() > 1){// my_mt.size() = {1;2}
- MElement* neighTri = my_mt[0] == e ? my_mt[1] : my_mt[0];
- neighbors[e].push_back(neighTri);
- }
- }
- }
-
- // queue: first in, first out
- std::queue<MElement*> checkList; // element for reference orientation
- std::queue< std::vector<MElement*> > checkLists; // corresponding neighbor element to be checked for its orientation
- std::queue<MElement*> my_todo; // todo list
- std::map<MElement*,bool> check_todo; // help to complete todo list
-
- my_todo.push(discrete_triangles[0]);
- check_todo[discrete_triangles[0]]=true;
- while(!my_todo.empty()){
-
- MElement* myMT = my_todo.front();
- my_todo.pop();
-
- std::vector<MElement*> myV = neighbors[myMT];
- std::vector<MElement*> myInsertion;
-
- checkList.push(myMT);
-
- for(unsigned int i=0; i<myV.size(); ++i){
- if (check_todo.find(myV[i]) == check_todo.end()){
- myInsertion.push_back(myV[i]);
- check_todo[myV[i]] = true;
- my_todo.push(myV[i]);
- }
- }
- checkLists.push(myInsertion);
- }// end while
-
-
- while(!checkList.empty() && !checkLists.empty()){
- MElement* current = checkList.front();
- checkList.pop();
- std::vector<MElement*> neigs = checkLists.front();
- checkLists.pop();
- for (unsigned int i=0; i<neigs.size(); i++){
- bool myCond = false;
- for (unsigned int k=0; k<3; k++){
- for (unsigned int j=0; j<3; j++){
- if (current->getVertex(k) == neigs[i]->getVertex(j)){
- myCond = true;
- if (!(current->getVertex(k!=2 ?k+1 : 0 ) == neigs[i]->getVertex(j!=0 ? j-1 : 2) ||
- current->getVertex(k!=0 ?k-1 : 2 ) == neigs[i]->getVertex(j!=2 ? j+1 : 0))){
- neigs[i]->reverse();
- // std::cout << "discreteDiskFace: triangle " << neigs[i]->getNum() << " has been reoriented." << std::endl;
- }
- break;
- }
- }
- if (myCond)
- break;
- }
- } // end for unsigned int i
- } // end while
-
-
- // now, it is possible to identify the border(s) of the triangulation
- // allEdges contains all edges of the boundary
- std::set<MEdge,Less_Edge> allEdges;
-
- for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- MElement *e = discrete_triangles[i];
- for(int j = 0; j < e->getNumEdges() ; j++){
- MEdge ed = e->getEdge(j);
- std::set<MEdge,Less_Edge>::iterator it = allEdges.find(ed);
- if (it == allEdges.end()) allEdges.insert(ed);
- else allEdges.erase(it);
- }
- }
-
- std::map<MVertex*,std::vector<MVertex*> > firstNode2Edge;
- for (std::set<MEdge>::iterator ie = allEdges.begin(); ie != allEdges.end() ; ++ie) {
- MEdge ed = *ie;
-
- std::vector<MElement *> vectri = ed2tri[ed];
- MElement* tri = vectri[0]; // boundary edge: only one triangle :-)
-
- std::vector<MVertex*> vecver;
- tri->getEdgeVertices(edgeLocalNum(tri,ed),vecver);
- MVertex *first = vecver[0];
- MVertex *last = vecver[1];
- vecver.erase(vecver.begin());
- 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", tag());
-
- firstNode2Edge[first] = vecver;
- firstNode2Edge[first].push_back(last);
- }
-
- while (!firstNode2Edge.empty()) { // quid 'firstNode2Edge' is not empty but 'in' within the map ?
- std::vector<MVertex*> loop;
- double length = 0.;
-
- std::map<MVertex*,std::vector<MVertex*> >::iterator in = firstNode2Edge.begin();
- while(in != firstNode2Edge.end()) {
- MVertex *first = in->first;
- std::vector<MVertex*> myV = in->second;
-
- loop.push_back(first);
-
- MVertex* previous = first;
-
- for(unsigned int i=0; i<=myV.size()-1; i++){
-
- MVertex* current = myV[i];
-
- loop.push_back(current);
-
- 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()) );
-
- _loops.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 ?
-
- previous = current;
-
- }
- firstNode2Edge.erase(in);
- in = firstNode2Edge.find(previous);
- }// end while in
- } // end while firstNode2Edge
-
- // dirichlet BC
- _U0 = _loops.rbegin()->second;
- _totLength = _loops.rbegin()->first;
-}
-
-
void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
{
if (oct)Octree_Delete(oct);
@@ -468,22 +319,22 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
double ssize[3] = {2.02,2.02,2.0};
const int maxElePerBucket = 15;
oct = Octree_Create(maxElePerBucket, origin, ssize, discreteDiskFaceBB,
- discreteDiskFaceCentroid, discreteDiskFaceInEle);
-
+ discreteDiskFaceCentroid, discreteDiskFaceInEle);
+
_ddft = new discreteDiskFaceTriangle[discrete_triangles.size()];
// if (CAD) printf("-------------> %ld %ld\n",CAD->size(),discrete_triangles.size());
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- MTriangle *t = discrete_triangles[i];
+ MElement *t = discrete_triangles[i];
discreteDiskFaceTriangle* my_ddft = &_ddft[i];
-
+ 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);
@@ -491,7 +342,7 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
bool discreteDiskFace::parametrize(bool one2one) const
-{ // mark, to improve
+{ // #improveme
if (one2one)
Msg::Info("Parametrizing surface %d with 'one-to-one map'", tag());
@@ -500,34 +351,34 @@ 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){
- MTriangle *t = discrete_triangles[i];
+ MElement *t = discrete_triangles[i];
for(int j=0; j<t->getNumVertices(); j++){
myAssemblerU.numberVertex(t->getVertex(j), 0, 1);
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);
@@ -552,14 +403,12 @@ bool discreteDiskFace::parametrize(bool one2one) const
}
}
-
-
double t2 = Cpu();
Msg::Debug("Assembly done in %8.3f seconds", t2 - t1);
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;
@@ -577,21 +426,21 @@ 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);
@@ -602,29 +451,50 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
bool discreteDiskFace::checkOrientationUV(){
- double initial, current; // initial and current orientation
discreteDiskFaceTriangle *ct;
- ct = &_ddft[0];
- double p1[2] = {ct->p[0].x(), ct->p[0].y()};
- double p2[2] = {ct->p[1].x(), ct->p[1].y()};
- double p3[2] = {ct->p[2].x(), ct->p[2].y()};
- initial = robustPredicates::orient2d(p1, p2, p3);
- unsigned int i=1;
- for (; i<discrete_triangles.size(); i++){
- ct = &_ddft[i];
- p1[0] = ct->p[0].x(); p1[1] = ct->p[0].y();
- p2[0] = ct->p[1].x(); p2[1] = ct->p[1].y();
- p3[0] = ct->p[2].x(); p3[1] = ct->p[2].y();
- current = robustPredicates::orient2d(p1, p2, p3);
- if(initial*current < 0.) break;
+
+ if(_order==1){
+ double initial, current; // initial and current orientation
+ ct = &_ddft[0];
+ double p1[2] = {ct->p[0].x(), ct->p[0].y()};
+ double p2[2] = {ct->p[1].x(), ct->p[1].y()};
+ double p3[2] = {ct->p[2].x(), ct->p[2].y()};
+ initial = robustPredicates::orient2d(p1, p2, p3);
+ unsigned int i=1;
+ for (; i<discrete_triangles.size(); i++){
+ ct = &_ddft[i];
+ p1[0] = ct->p[0].x(); p1[1] = ct->p[0].y();
+ p2[0] = ct->p[1].x(); p2[1] = ct->p[1].y();
+ p3[0] = ct->p[2].x(); p3[1] = ct->p[2].y();
+ current = robustPredicates::orient2d(p1, p2, p3);
+ if(initial*current < 0.) {
+ return false;
+ break;
+ }
+ }
+ return true;
}
- if(i < discrete_triangles.size())
- return false;
- // Msg::Error("discreteDiskFace: The parametrization is not one-to-one :-(");
- 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(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];
+ if (min < 0){
+ return false;
+ break;
+ }
+ }
return true;
+ }
}
-
+
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
@@ -637,11 +507,9 @@ 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;
- }
+ }
- //polynomialBasis myPolynomial = polynomialBasis(dt->tri->getTypeForMSH());
std::vector<double> eval(_N);
- //mynodalbasis->f(U,V,0.,eval);//#mark
functionShapes(_order,Xi,&eval[0]);
double X=0,Y=0,Z=0;
for(int io=0; io<_N; io++){
@@ -651,7 +519,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){
@@ -660,7 +528,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
@@ -681,82 +549,81 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
}
Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
}
- else if (v->onWhat()->dim()==0){
+ else if (v->onWhat()->dim()==0)
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
+ 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);
-
- MTriangle* tri = NULL;
+
+ 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];
- //mynodalbasis->df(U,V,0.,df);
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.;
@@ -765,81 +632,119 @@ 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
-
-
- return;
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{ // cf Sauvage's thesis
+ return;
}
-void discreteDiskFace::buildAllNodes()
+
+void discreteDiskFace::putOnView()
{
- // allNodes contains all mesh-nodes
- for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- MElement *e = discrete_triangles[i];
- for(int j = 0; j < e->getNumVertices() ; j++){
- MVertex* ev = e->getVertex(j);
- std::set<MVertex*>::iterator it = allNodes.find(ev); // several times the same nodes !
- if (it == allNodes.end()) allNodes.insert (ev);
- }
- }
-}
+ char mybuffer [64];
-void discreteDiskFace::putOnView()
-{
+ snprintf(mybuffer,sizeof(mybuffer),"param_u_part%d_order%d.pos",initialTriangulation->idNum,_order);
+ FILE* view_u = Fopen(mybuffer,"w");
- FILE* view_u = Fopen("param_u.pos","w");
- FILE* view_v = Fopen("param_v.pos","w");
- FILE* UVxyz = Fopen("UVxyz.pos","w");
- if(view_u && view_v){
+ 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);
+ FILE* UVx = Fopen(mybuffer,"w");
+
+ 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);
+ FILE* UVz = Fopen(mybuffer,"w");
+
+ if(view_u && view_v && UVx && UVy && UVz){
fprintf(view_u,"View \"paramU\"{\n");
fprintf(view_v,"View \"paramV\"{\n");
- fprintf(UVxyz,"View \"Uxyz\"{\n");
+ fprintf(UVx,"View \"Ux\"{\n");
+ fprintf(UVy,"View \"Uy\"{\n");
+ fprintf(UVz,"View \"Uz\"{\n");
for (unsigned int i=0; i<discrete_triangles.size(); i++){
discreteDiskFaceTriangle* my_ddft = &_ddft[i];
if (_order == 1){
fprintf(view_u,"ST(");
fprintf(view_v,"ST(");
- fprintf(UVxyz,"ST(");
+ fprintf(UVx,"ST(");
+ fprintf(UVy,"ST(");
+ fprintf(UVz,"ST(");
}
else{
fprintf(view_u,"ST%d(",_order);
fprintf(view_v,"ST%d(",_order);
- fprintf(UVxyz,"ST%d(",_order);
+ fprintf(UVx,"ST%d(",_order);
+ fprintf(UVy,"ST%d(",_order);
+ 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(UVxyz,"%g,%g,%g,",my_ddft->p[j].x(),my_ddft->p[j].y(),0.);
+ 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(UVxyz,"%g,%g,%g){",my_ddft->p[_N-1].x(),my_ddft->p[_N-1].y(),0.);
+ 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.);
for (int j=0; j<_N-1; j++){
fprintf(view_u,"%g,",my_ddft->p[j].x());
fprintf(view_v,"%g,",my_ddft->p[j].y());
- fprintf(UVxyz,"%d,",my_ddft->gf->tag());
+ fprintf(UVx,"%g,",my_ddft->tri->getVertex(j)->x());
+ fprintf(UVy,"%g,",my_ddft->tri->getVertex(j)->y());
+ fprintf(UVz,"%g,",my_ddft->tri->getVertex(j)->z());
}
fprintf(view_u,"%g};\n",my_ddft->p[_N-1].x());
fprintf(view_v,"%g};\n",my_ddft->p[_N-1].y());
- fprintf(UVxyz,"%d};\n",my_ddft->gf->tag());
+ fprintf(UVx,"%g};\n",my_ddft->tri->getVertex(_N-1)->x());
+ fprintf(UVy,"%g};\n",my_ddft->tri->getVertex(_N-1)->y());
+ fprintf(UVz,"%g};\n",my_ddft->tri->getVertex(_N-1)->z());
}
fprintf(view_u,"};\n");
fprintf(view_v,"};\n");
- fprintf(UVxyz,"};\n");
+ fprintf(UVx,"};\n");
+ fprintf(UVy,"};\n");
+ fprintf(UVz,"};\n");
fclose(view_u);
fclose(view_v);
- fclose(UVxyz);
+ fclose(UVx);
+ fclose(UVy);
+ fclose(UVz);
}
+ /*
+ #ifdef HAVE_POST
+ std::map<int, std::vector<double> > u;
+ std::map<int, std::vector<double> > v;
+ for(std::set<MVertex*>::iterator iv = allNodes.begin(); iv!=allNodes.end(); ++iv){
+ MVertex *p = *iv;
+ u[p->getNum()].push_back(coordinates[p].x());
+ v[p->getNum()].push_back(coordinates[p].y());
+ }
+ PView* view_u = new PView("u", "NodeData", GModel::current(), u);
+ PView* view_v = new PView("v", "NodeData", GModel::current(), v);
+ view_u->setChanged(true);
+ view_v->setChanged(true);
+ view_u->write("param_u.msh", 5);
+ view_v->write("param_v.msh", 5);
+ delete view_u;
+ delete view_v;
+ #else
+ Msg::Error("discreteDiskFace: cannot export without post module")
+ #endif
+ */
}
-
+
// useful for mesh generators ----------------------------------------
// Intersection of a circle and a plane
GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
@@ -894,7 +799,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;
@@ -912,14 +817,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);
@@ -943,6 +848,6 @@ 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 1fb99a2df7..56fbb1fb36 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,33 +23,159 @@
#include "PView.h"
#include "robustPredicates.h"
+/*inline utilities*/
+inline int nodeLocalNum(MElement* e, MVertex* v) {// const
+ 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
+ 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 {
+
+ public:
+
+ // 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*
+ 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(){
+ return ( -vert.size() + ed2tri.size() - tri.size() + 2 - bord.size() )/2;
+ }
+
+ void assignVert(){
+ for(unsigned int i = 0; i < tri.size(); ++i){
+ MElement* t = tri[i];
+ for(int j = 0; j < t->getNumVertices() ; j++){
+ MVertex* tv = t->getVertex(j);
+ std::set<MVertex*>::iterator it = vert.find(tv);
+ if (it == vert.end()) vert.insert (tv);
+ }
+ }
+ }
+
+ void assignEd2tri(){
+ for(unsigned int i = 0; i < tri.size(); ++i){
+ MElement *t = tri[i];
+ for(int j = 0; j < 3 ; j++){
+ MEdge ed = t->getEdge(j);
+ ed2tri[ed].push_back(i);
+ }
+ }
+ }
+
+ void assignBord(){
+ for(unsigned int i = 0; i < tri.size(); ++i){
+ MElement *t = tri[i];
+ for(int j = 0; j < t->getNumEdges() ; j++){
+ MEdge ed = t->getEdge(j);
+ std::set<MEdge,Less_Edge>::iterator it = borderEdg.find(ed);
+ if (it == borderEdg.end()) borderEdg.insert(ed);
+ else borderEdg.erase(it);
+ }
+ }
+
+ std::map<MVertex*,std::vector<MVertex*> > firstNode2Edge;
+ for (std::set<MEdge>::iterator ie = borderEdg.begin(); ie != borderEdg.end() ; ++ie) {
+ MEdge ed = *ie;
+ std::vector<int> nT = ed2tri[ed];
+ MElement* t = tri[nT[0]];
+
+ std::vector<MVertex*> vecver;
+ t->getEdgeVertices(edgeLocalNum(t,ed),vecver);
+ MVertex *first = vecver[0];
+ MVertex *last = vecver[1];
+ vecver.erase(vecver.begin());
+ 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());
+ firstNode2Edge[first] = vecver;
+ firstNode2Edge[first].push_back(last);
+ }
+ while (!firstNode2Edge.empty()) {
+ std::vector<MVertex*> loop;
+ double length = 0.;
+
+ 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];
+
+ 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 ?
+ } // end while firstNode2Edge
+ }// end method
+
+ void assign(){
+ assignVert();
+ assignEd2tri();
+ assignBord();
+ }
+
+ //builder
+ 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
class discreteDiskFaceTriangle {
public:
- SPoint3 p[6]; // vertices in (u;v) #improveme std::vector instead
+ std::vector<SPoint3> p; // vertices in (u;v) #improveme std::vector instead
//SPoint2 gfp[6]; // CAD model
GFace *gf; // GFace tag
- MTriangle *tri; // mesh triangle in (x;y;z)
+ MElement *tri; // mesh triangle in (x;y;z)
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 buildAllNodes() ;
- void getBoundingEdges();
void putOnView();
bool checkOrientationUV();
public:
- discreteDiskFace(GFace *parent, std::vector<MTriangle*> &mesh, int p=1, std::vector<GFace*> *CAD = NULL);// MTriangle -> MTriangle 6
+ discreteDiskFace(GFace *parent, triangulation* diskTriangulation, int p=1, std::vector<GFace*> *CAD = NULL);// BUILDER
virtual ~discreteDiskFace();
void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v) const;
GPoint point(double par1, double par2) const;
@@ -59,7 +185,7 @@ class discreteDiskFace : public GFace {
double curvatures(const SPoint2&,SVector3*,SVector3*,double*,double*) const;
virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 ¶m,
- SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
GEntity::GeomType geomType() const { return DiscreteDiskSurface; }
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
@@ -67,38 +193,28 @@ class discreteDiskFace : public GFace {
protected:
//------------------------------------------------
// a copy of the mesh that should not be destroyed
- std::vector<MTriangle*> discrete_triangles;
+ triangulation* initialTriangulation;
+ triangulation* geoTriangulation;// parametrized triangulation
+ std::vector<MElement*> discrete_triangles;
std::vector<MVertex*> discrete_vertices;
//------------------------------------------------
- int nodeLocalNum(MElement* e, MVertex* v) const{
- for(int i=0; i<e->getNumVertices(); i++)
- if (v == e->getVertex(i))
- return i;
- return -1;
- }
- int edgeLocalNum(MElement* e, MEdge ed) const{
- for(int i=0; i<e->getNumEdges(); i++)
- if (ed == e->getEdge(i))
- return i;
- return -1;
- }
- //------------------------------------------------
+
+ //-----------------------------------------------
int _order;
- int _N;
+ int _N;// number of dof's for a triangle
double _totLength;
std::map<double,std::vector<MVertex*> > _loops;
std::vector<MVertex*> _U0; // dirichlet's bc's
mutable std::set<MVertex*> allNodes;
mutable std::map<MVertex*, SPoint3> coordinates;
- mutable std::map<SPoint3,SPoint3 > _coordPoints; // ?
- mutable v2t_cont adjv; // ? v2t_cont ? ?????
+ mutable std::map<SPoint3,SPoint3 > _coordPoints;
+ mutable v2t_cont adjv;
mutable std::map<MVertex*, Pair<SVector3,SVector3> > firstDerivatives;
mutable discreteDiskFaceTriangle *_ddft;
- mutable Octree *oct;
+ mutable Octree *oct = NULL;
mutable std::vector<double> _coords;
- // const nodalBasis* mynodalbasis;
};
-
+
#endif
-
+
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 039681c756..e9a69c2605 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -3,23 +3,26 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
-#include <stdlib.h>
+
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "discreteFace.h"
#include "discreteDiskFace.h"
-#include "MTriangle.h"
-#include "MEdge.h"
#include "Geo.h"
#include "GFaceCompound.h"
#include "Context.h"
#include "OS.h"
+#include <stack>
+#include <queue>
+extern "C" {
+#include <metis.h>
+}
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
- Tree_Add(model->getGEOInternals()->Surfaces, &s);
- meshStatistics.status = GFace::DONE;
+ Tree_Add(model->getGEOInternals()->Surfaces, &s);
+ meshStatistics.status = GFace::DONE;
}
void discreteFace::setBoundEdges(GModel *gm, std::vector<int> tagEdges)
@@ -77,19 +80,18 @@ SPoint2 discreteFace::parFromPoint(const SPoint3 &p, bool onSurface) const
SVector3 discreteFace::normal(const SPoint2 ¶m) const
{
- return SVector3();
+ return SVector3();
}
double discreteFace::curvatureMax(const SPoint2 ¶m) const
{
-
- return false;
+ return false;
}
double discreteFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const
{
- return false;
+ return false;
}
Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
@@ -106,13 +108,52 @@ void discreteFace::secondDer(const SPoint2 ¶m,
// FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
void discreteFace::createGeometry()
{
+ checkAndFixOrientation();
+
+ int order = 1;
+ int nPart = 2;
+ std::vector<triangulation*> part;
+ part.resize(nPart);
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
+
if (!_atlas.empty())return;
- // parametrization is done here !!!
- // if (_CAD != empty) --> _triangles[i] is classified on _CAD[i]
- discreteDiskFace *df = new discreteDiskFace (this, triangles,1,(_CAD.empty() ? NULL : &_CAD));
- df->replaceEdges(l_edges);
- _atlas.push_back(df);
+
+ int id=1;
+ std::stack<triangulation*> toSplit;
+ std::vector<triangulation*> toParam;
+ std::vector<MElement*> tem(triangles.begin(),triangles.end());
+ toSplit.push(new triangulation(tem,this));
+ if((toSplit.top())->genus()!=0){
+
+ while( !toSplit.empty()){
+
+ triangulation* tosplit = toSplit.top();
+ toSplit.pop();
+
+ split(tosplit,part,nPart);
+ delete tosplit; // #mark
+
+ for(int i=0; i<nPart; i++){
+ if(part[i]->genus()!=0)
+ toSplit.push(part[i]);
+ else{
+ toParam.push_back(part[i]);
+ part[i]->idNum=id++;
+ }
+ }// end for i
+ }
+
+ }// end if it is not disk-like
+ else{
+ toParam.push_back(toSplit.top());
+ toSplit.top()->idNum=id++;
+ }
+
+ for(unsigned int i=0; i<toParam.size(); i++){
+ discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
+ df->replaceEdges(l_edges);// #FIXME
+ _atlas.push_back(df);
+ }
#endif
}
@@ -169,6 +210,145 @@ void discreteFace::mesh(bool verbose)
#endif
}
+
+void discreteFace::checkAndFixOrientation(){
+
+ // first of all, all the triangles have to be oriented in the same way
+ std::map<MEdge,std::vector<MElement*>,Less_Edge> ed2tri; // edge to 1 or 2 triangle(s)
+
+ for(unsigned int i = 0; i < triangles.size(); ++i){
+ MElement *e = triangles[i];
+ for(int j = 0; j < e->getNumEdges() ; j++){
+ MEdge ed = e->getEdge(j);
+ ed2tri[ed].push_back(e);
+ }
+ }
+
+ // element to its neighbors
+ std::map<MElement*,std::vector<MElement*> > neighbors;
+ for (unsigned int i=0; i<triangles.size(); ++i){
+ MElement* e = triangles[i];
+ for(int j=0; j<e->getNumEdges(); j++){ // #improveme: efficiency could be improved by setting neighbors mutually
+ std::vector<MElement*> my_mt = ed2tri[e->getEdge(j)];
+ if (my_mt.size() > 1){// my_mt.size() = {1;2}
+ MElement* neighTri = my_mt[0] == e ? my_mt[1] : my_mt[0];
+ neighbors[e].push_back(neighTri);
+ }
+ }
+ }
+
+ // queue: first in, first out
+ std::queue<MElement*> checkList; // element for reference orientation
+ std::queue< std::vector<MElement*> > checkLists; // corresponding neighbor element to be checked for its orientation
+ std::queue<MElement*> my_todo; // todo list
+ std::map<MElement*,bool> check_todo; // help to complete todo list
+
+ my_todo.push(triangles[0]);
+
+ check_todo[triangles[0]]=true;
+ while(!my_todo.empty()){
+
+ MElement* myMT = my_todo.front();
+ my_todo.pop();
+
+ std::vector<MElement*> myV = neighbors[myMT];
+ std::vector<MElement*> myInsertion;
+
+ checkList.push(myMT);
+
+ for(unsigned int i=0; i<myV.size(); ++i){
+ if (check_todo.find(myV[i]) == check_todo.end()){
+ myInsertion.push_back(myV[i]);
+ check_todo[myV[i]] = true;
+ my_todo.push(myV[i]);
+ }
+ }
+ checkLists.push(myInsertion);
+ }// end while
+
+
+ while(!checkList.empty() && !checkLists.empty()){
+ MElement* current = checkList.front();
+ checkList.pop();
+ std::vector<MElement*> neigs = checkLists.front();
+ checkLists.pop();
+ for (unsigned int i=0; i<neigs.size(); i++){
+ bool myCond = false;
+ for (unsigned int k=0; k<3; k++){
+ for (unsigned int j=0; j<3; j++){
+ if (current->getVertex(k) == neigs[i]->getVertex(j)){
+ myCond = true;
+ if (!(current->getVertex(k!=2 ?k+1 : 0 ) == neigs[i]->getVertex(j!=0 ? j-1 : 2) ||
+ current->getVertex(k!=0 ?k-1 : 2 ) == neigs[i]->getVertex(j!=2 ? j+1 : 0))){
+ neigs[i]->reverse();
+ Msg::Info("discreteDiskFace: triangle %d has been reoriented.",neigs[i]->getNum());
+ }
+ break;
+ }
+ }
+ if (myCond)
+ break;
+ }
+ } // end for unsigned int i
+ } // end while
+}
+
+void discreteFace::split(triangulation* trian,std::vector<triangulation*> &partition,int nPartitions)
+{
+
+ int nVertex = trian->tri.size(); // number of elements
+ int nEdge = trian->ed2tri.size() - trian->borderEdg.size();// number of edges, (without the boundary ones)
+
+ std::vector<int> idx;
+ idx.resize(nVertex+1);
+
+ std::vector<int> nbh;
+ nbh.resize(2*nEdge);
+
+ idx[0]=0;
+ for(int i=0; i<nVertex; i++){// triangle by triangle
+
+ MElement* current = trian->tri[i];
+
+ int temp = 0;
+ for(int j=0; j<3; j++){ // edge by edge
+
+ MEdge ed = current->getEdge(j);
+ int nEd = trian->ed2tri[ed].size();
+
+ if (nEd > 1){
+ std::vector<int> local = trian->ed2tri[ed];
+ nbh[idx[i]+temp] = local[0] == i ? local[1] : local[0];
+ temp++;
+ }
+
+ }// end for j
+
+ idx[i+1] = idx[i]+temp;
+
+ }// end for i
+
+
+ int edgeCut;
+ std::vector<int> part;
+ part.resize(nVertex);
+ int zero = 0;
+ METIS_PartGraphRecursive(&nVertex,&idx[0],&nbh[0],NULL,NULL,&zero,&zero,&nPartitions,&zero,&edgeCut,&part[0]);
+ // CONNEC
+ std::vector<std::vector<MElement*> > elem;
+ elem.resize(nPartitions);
+
+ for(int i=0; i<nVertex; i++)
+ elem[part[i]].push_back(trian->tri[i]);
+
+ for(int i=0; i<nPartitions; i++)
+ partition[i] = new triangulation(elem[i],this);
+
+
+ return;
+}
+
+
// delete all discrete disk faces
//void discreteFace::deleteAtlas() {
//}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index c489a934cf..8b4abfb2d4 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -10,17 +10,27 @@
#include "GFace.h"
#include "discreteEdge.h"
#include "MEdge.h"
+#include "meshPartitionObjects.h"
+#include "meshPartitionOptions.h"
+#include "meshPartition.h"
+#include "MTriangle.h"
+#include "MEdge.h"
+#include <stdlib.h>
class discreteDiskFace;
+class triangulation;
+
class discreteFace : public GFace {
// FIXME we should at the end use a
- // mesh() functioon that is specific to
+ // mesh() function that is specific to
// discreteFace
// we should also SAVE those data's
public:
discreteFace(GModel *model, int num);
virtual ~discreteFace() {}
+ void checkAndFixOrientation();
+ void split(triangulation*,std::vector<triangulation*>&,int);
GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
SVector3 normal(const SPoint2 ¶m) const;
--
GitLab