← Back to branch summary

~ubuntu-branches/ubuntu/karmic/gmsh/karmic

« back to all changes in this revision

Viewing changes to Mesh/meshGFaceDelaunayInsertion.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Daniel Leidert
  • Date: 2008-05-18 12:46:05 UTC
  • mfrom: (1.2.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20080518124605-716xqbqeo07o497k
Tags: 2.2.0-2
[Christophe Prud'homme]
* Bug fix: "gmsh ships no .desktop", thanks to Vassilis Pandis (Closes:
  #375770). Applied the Ubuntu patch.

[Daniel Leidert]
* debian/control (Vcs-Svn): Fixed.
  (Build-Depends): Use texlive instead of tetex-bin.
* debian/gmsh.doc-base (Section): Fixed accordingly to doc-base (>= 0.8.10).
* debian/rules: Removed some variable declarations, that lead to double
  configuration and seem to be useless.
  (build/gmsh): Try to avoid multiple runs by using a stamp.
  (orig-tarball): Renamed to get-orig-source and changed to use uscan.
* debian/watch: Added.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Mesh/meshGFaceDelaunayInsertion.cpp
1
 
// $Id: meshGFaceDelaunayInsertion.cpp,v 1.2 2007-05-24 14:44:06 remacle Exp $
 
1
// $Id: meshGFaceDelaunayInsertion.cpp,v 1.25 2008-04-17 09:07:01 remacle Exp $
2
2
//
3
 
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 
3
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
4
4
//
5
5
// This program is free software; you can redistribute it and/or modify
6
6
// it under the terms of the GNU General Public License as published by
22
22
#include "BDS.h"
23
23
#include "BackgroundMesh.h"
24
24
#include "meshGFaceDelaunayInsertion.h"
 
25
#include "meshGFaceOptimize.h"
 
26
#include "meshGFace.h"
25
27
#include "GFace.h"
26
28
#include "Numeric.h"
27
29
#include "Message.h"
29
31
#include <map>
30
32
#include <algorithm>
31
33
 
32
 
// computes the center of the circum circle in the tangent plane
33
 
// the metric given by a b d is supposed to be constant
34
 
void MTri3::Center_Circum_Aniso(double a, double b, double d, double &x, double &y, double &r) const
35
 
{
36
 
  double sys[2][2], X[2];
 
34
#include "Context.h"
 
35
extern Context_T CTX;
 
36
 
 
37
const double LIMIT_ = 0.5 * sqrt(2.0);
 
38
 
 
39
/*  This function controls the frontal algorithm  */
 
40
static bool isActive ( MTri3 *t , double limit_, int &active){
 
41
  if (t->isDeleted()) return false;
 
42
  for (active=0;active<3;active++){
 
43
    MTri3 *neigh = t->getNeigh(active);
 
44
    if (!neigh || neigh->getRadius() < limit_)return true;
 
45
        //if (!neigh || neigh->done)return true;
 
46
  }
 
47
  return false;
 
48
}
 
49
 
 
50
 
 
51
void circumCenterXY(double *p1, double *p2, double *p3, double *res)
 
52
{
 
53
  double d, a1, a2, a3;
 
54
 
 
55
  const double x1 = p1[0];
 
56
  const double x2 = p2[0];
 
57
  const double x3 = p3[0];
 
58
  const double y1 = p1[1];
 
59
  const double y2 = p2[1];
 
60
  const double y3 = p3[1];
 
61
 
 
62
  d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
 
63
  if(d == 0.0) {
 
64
    Msg(WARNING, "Colinear points in circum circle computation");
 
65
    res[0] = res[1] = -99999.;
 
66
    return ;
 
67
  }
 
68
 
 
69
  a1 = x1 * x1 + y1 * y1;
 
70
  a2 = x2 * x2 + y2 * y2;
 
71
  a3 = x3 * x3 + y3 * y3;
 
72
  res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
 
73
  res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
 
74
}
 
75
 
 
76
void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv)
 
77
{
 
78
  double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
 
79
  double v2[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
 
80
  double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
 
81
  double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
 
82
  double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
 
83
  norme(vx); norme(vy); norme(vz);
 
84
  double p1P[2] = {0.0, 0.0};
 
85
  double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
 
86
  double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
 
87
  double resP[2];
 
88
 
 
89
  circumCenterXY(p1P, p2P, p3P,resP);
 
90
 
 
91
  if(uv){
 
92
    double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
 
93
                        {p2P[1] - p1P[1], p3P[1] - p1P[1]}};
 
94
    double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
 
95
    sys2x2(mat, rhs, uv);
 
96
  }
 
97
  
 
98
  res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
 
99
  res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
 
100
  res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
 
101
}
 
102
 
 
103
bool circumCenterMetricInTriangle(MTriangle *base, 
 
104
                                  const double *metric,
 
105
                                  const std::vector<double> &Us,
 
106
                                  const std::vector<double> &Vs)
 
107
{
 
108
  double R, x[2], uv[2];
 
109
  circumCenterMetric(base, metric, Us, Vs, x, R);
 
110
  return invMapUV(base, x, Us, Vs, uv, 1.e-8);
 
111
}
 
112
 
 
113
void circumCenterMetric(double *pa, double *pb, double *pc,
 
114
                        const double *metric,
 
115
                        double *x, double &Radius2) 
 
116
{
 
117
  // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1 
 
118
  double sys[2][2];
37
119
  double rhs[2];
38
 
  double x1, y1, x2, y2, x3, y3;
39
 
 
40
 
  x1 = base->getVertex(0)->x();
41
 
  y1 = base->getVertex(0)->y();
42
 
  x2 = base->getVertex(1)->x();
43
 
  y2 = base->getVertex(1)->y();
44
 
  x3 = base->getVertex(2)->x();
45
 
  y3 = base->getVertex(2)->y();
46
 
 
47
 
  sys[0][0] = 2. * a * (x1 - x2) + 2. * b * (y1 - y2);
48
 
  sys[0][1] = 2. * d * (y1 - y2) + 2. * b * (x1 - x2);
49
 
  sys[1][0] = 2. * a * (x1 - x3) + 2. * b * (y1 - y3);
50
 
  sys[1][1] = 2. * d * (y1 - y3) + 2. * b * (x1 - x3);
 
120
 
 
121
  const double a = metric[0];
 
122
  const double b = metric[1];
 
123
  const double d = metric[2];
 
124
 
 
125
  sys[0][0] = 2. * a * (pa[0] - pb[0]) + 2. * b * (pa[1] - pb[1]);
 
126
  sys[0][1] = 2. * d * (pa[1] - pb[1]) + 2. * b * (pa[0] - pb[0]);
 
127
  sys[1][0] = 2. * a * (pa[0] - pc[0]) + 2. * b * (pa[1] - pc[1]);
 
128
  sys[1][1] = 2. * d * (pa[1] - pc[1]) + 2. * b * (pa[0] - pc[0]);
51
129
 
52
130
  rhs[0] =
53
 
    a * (x1 * x1 - x2 * x2) + d * (y1 * y1 - y2 * y2) + 2. * b * (x1 * y1 -
54
 
                                                                  x2 * y2);
 
131
    a * (pa[0] * pa[0] - pb[0] * pb[0]) + 
 
132
    d * (pa[1] * pa[1] - pb[1] * pb[1]) + 
 
133
    2. * b * (pa[0] * pa[1] - pb[0] * pb[1]);
55
134
  rhs[1] =
56
 
    a * (x1 * x1 - x3 * x3) + d * (y1 * y1 - y3 * y3) + 2. * b * (x1 * y1 -
57
 
                                                                  x3 * y3);
58
 
  sys2x2(sys, rhs, X);
59
 
  
60
 
  x = X[0];
61
 
  y = X[1];
62
 
  
63
 
  r = sqrt((X[0] - x1) * (X[0] - x1) * a
64
 
           + (X[1] - y1) * (X[1] - y1) * d
65
 
           + 2. * (X[0] - x1) * (X[1] - y1) * b);
66
 
}
67
 
 
68
 
 
69
 
// find a common basis for 2 metrics in 2D 
70
 
void simultaneousMetricReduction( const gmsh2dMetric &M1,  const gmsh2dMetric &M2,
71
 
                                  double & l1,double & l2, double V[2][2]) 
72
 
{
73
 
  double a11=M1.a11,a21=M1.a21,a22=M1.a22;
74
 
  double b11=M2.a11,b21=M2.a21,b22=M2.a22;
75
 
  const double /*c11 = a11*a11,*/ c21= a21*a21;
76
 
  const double /*d11 = b11*b11,*/ d21= b21*b21;
77
 
  const double a=b11*b22 - d21;
78
 
  const double b=-a11*b22-a22*b11+2*a21*b21;
79
 
  const double c=-c21+a11*a22;
80
 
  const double bb = b*b,ac= a*c;
81
 
  const double delta = bb - 4 * ac;
82
 
  if (bb + fabs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) )
83
 
    {
84
 
      if (fabs(a) < 1.e-30 )
85
 
        l1 = l2 = 0;
86
 
      else 
87
 
        l1=l2=-b/(2*a); 
88
 
      V[0][0] = 1.;
89
 
      V[1][0] = 0.;
90
 
      V[0][1] = 0.;
91
 
      V[1][1] = 1.;
92
 
    }
93
 
  else {
94
 
    const double delta2 = sqrt(delta);
95
 
    l1= (-b - delta2)/(2*a);
96
 
    l2= (-b + delta2)/(2*a);
97
 
    double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
98
 
    double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
99
 
    double vp1x,vp1y,vp2x,vp2y;   
100
 
    if(s1 < s0)
101
 
      s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
102
 
    else
103
 
      s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
104
 
    
105
 
    v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
106
 
    s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
107
 
    if(s1 < s0)
108
 
      s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
109
 
    else
110
 
      s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
111
 
      V[0][0] = vp1x;
112
 
      V[0][1] = vp2x;
113
 
      V[1][0] = vp1y;
114
 
      V[1][1] = vp2y;
115
 
  }
116
 
}
117
 
 
118
 
gmsh2dMetric metricIntersection(const gmsh2dMetric &Ma,const gmsh2dMetric &Mb) 
119
 
{
120
 
  double M[2][2],M1[2][2];
121
 
  double l1,l2;
122
 
  simultaneousMetricReduction(Ma,Mb,l1,l2,M);
123
 
  inv2x2 ( M, M1 );
124
 
  const double M1ev1 = Ma.eval(M[0][0],M[0][1]);
125
 
  const double M1ev2 = Ma.eval(M[0][0],M[0][1]);
126
 
  const double M2ev1 = Mb.eval(M[1][0],M[1][1]);
127
 
  const double M2ev2 = Mb.eval(M[1][0],M[1][1]);
128
 
 
129
 
  const double D0 = std::max( M1ev1,M1ev2);
130
 
  const double D1 = std::max( M2ev1,M2ev2);
131
 
 
132
 
  return gmsh2dMetric(M1[0][0] * D0 * M1[0][0] + M1[0][1] * D1 * M1[0][1] ,
133
 
                      0.5 * ( M1[0][0] * D0 * M1[1][0] + M1[0][1] * D1 * M1[1][1] +
134
 
                              M1[1][0] * D0 * M1[0][0] + M1[1][1] * D1 * M1[0][1] ), 
135
 
                      M1[1][0] * D0 * M1[1][0] + M1[1][1] * D1 * M1[1][1] );                        
136
 
}
137
 
 
138
 
gmsh2dMetric metricInterpolationTriangle(const gmsh2dMetric &M1,
139
 
                                         const gmsh2dMetric &M2,
140
 
                                         const gmsh2dMetric &M3, 
141
 
                                         const double u, 
142
 
                                         const double v) 
143
 
{
144
 
  double m1[2][2] = {{M1.a11,0.5*M1.a21},{0.5*M1.a21,M1.a22}};
145
 
  double m2[2][2] = {{M2.a11,0.5*M2.a21},{0.5*M2.a21,M2.a22}};
146
 
  double m3[2][2] = {{M3.a11,0.5*M3.a21},{0.5*M3.a21,M3.a22}};
147
 
  double invm1 [2][2]; inv2x2 ( m1, invm1 );
148
 
  double invm2 [2][2]; inv2x2 ( m2, invm2 );
149
 
  double invm3 [2][2]; inv2x2 ( m3, invm3 ); 
150
 
  double invm [2][2] ;
151
 
  for (int i=0;i<2;i++)
152
 
    for (int j=0;j<2;j++)
153
 
      invm[i][j] = (1.-u-v) * invm1[i][j] + u * invm2[i][j] + v * invm3[i][j];
154
 
  double m [2][2]; inv2x2 ( invm, m );
155
 
  return gmsh2dMetric ( m[0][0] , m[0][1] , m[1][1] );
156
 
}
157
 
 
158
 
gmsh2dMetric :: gmsh2dMetric ( double lc )
159
 
  : a11 ( 1./(lc*lc) ) ,a21 ( 0.0 ) ,a22 ( 1./(lc*lc) )
160
 
{  
161
 
}
162
 
 
163
 
 
164
 
MTri3::MTri3 ( MTriangle * t, std::vector<gmsh2dMetric> & sizes) : deleted(false), base (t)
 
135
    a * (pa[0] * pa[0] - pc[0] * pc[0]) + 
 
136
    d * (pa[1] * pa[1] - pc[1] * pc[1]) + 
 
137
    2. * b * (pa[0] * pa[1] - pc[0] * pc[1]);
 
138
 
 
139
  sys2x2(sys, rhs, x);
 
140
 
 
141
  Radius2 = 
 
142
    (x[0] - pa[0]) * (x[0] - pa[0]) * a +
 
143
    (x[1] - pa[1]) * (x[1] - pa[1]) * d +
 
144
    2. * (x[0] - pa[0]) * (x[1] - pa[1]) * b;
 
145
}
 
146
 
 
147
void circumCenterMetric(MTriangle *base, 
 
148
                        const double *metric,
 
149
                        const std::vector<double> &Us,
 
150
                        const std::vector<double> &Vs,
 
151
                        double *x, double &Radius2) 
 
152
{
 
153
  // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1 
 
154
  double pa[2] = {Us[base->getVertex(0)->getNum()],
 
155
                  Vs[base->getVertex(0)->getNum()]};
 
156
  double pb[2] = {Us[base->getVertex(1)->getNum()],
 
157
                  Vs[base->getVertex(1)->getNum()]};
 
158
  double pc[2] = {Us[base->getVertex(2)->getNum()],
 
159
                  Vs[base->getVertex(2)->getNum()]};
 
160
  circumCenterMetric(pa, pb, pc, metric, x, Radius2);
 
161
}
 
162
 
 
163
void buildMetric(GFace *gf, double *uv, double *metric)
 
164
{
 
165
  Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
 
166
  metric[0] = dot(der.first(), der.first());
 
167
  metric[1] = dot(der.second(), der.first());
 
168
  metric[2] = dot(der.second(), der.second());
 
169
} 
 
170
 
 
171
int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3, 
 
172
                        double *uv, double *metric) 
 
173
{
 
174
  double x[2], Radius2;
 
175
  circumCenterMetric(p1, p2, p3, metric, x, Radius2);
 
176
  const double a = metric[0];
 
177
  const double b = metric[1];
 
178
  const double d = metric[2];
 
179
  double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
 
180
    + (x[1] - uv[1]) * (x[1] - uv[1]) * d
 
181
    + 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
 
182
  return d2 < Radius2;  
 
183
}
 
184
 
 
185
int inCircumCircleAniso(GFace *gf, MTriangle *base, 
 
186
                        const double *uv, const double *metricb,
 
187
                        const std::vector<double> &Us,
 
188
                        const std::vector<double> &Vs) 
 
189
{
 
190
  double x[2], Radius2, metric[3];
 
191
  double pa[2] = {(Us[base->getVertex(0)->getNum()] + Us[base->getVertex(1)->getNum()] +
 
192
                   Us[base->getVertex(2)->getNum()]) / 3.,
 
193
                  (Vs[base->getVertex(0)->getNum()] + Vs[base->getVertex(1)->getNum()] + 
 
194
                   Vs[base->getVertex(2)->getNum()]) / 3.};
 
195
  
 
196
  buildMetric(gf, pa, metric);
 
197
  circumCenterMetric(base, metric, Us, Vs, x, Radius2);
 
198
 
 
199
  const double a = metric[0];
 
200
  const double b = metric[1];
 
201
  const double d = metric[2];
 
202
 
 
203
  double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
 
204
    + (x[1] - uv[1]) * (x[1] - uv[1]) * d
 
205
    + 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
 
206
  
 
207
  return d2 < Radius2;  
 
208
}
 
209
 
 
210
MTri3::MTri3(MTriangle *t, double lc) : deleted(false), base(t)//,done(0)
165
211
{
166
212
  neigh[0] = neigh[1] = neigh[2] = 0;
167
 
  // compute the metric at the centroid
168
 
  metric = metricInterpolationTriangle(sizes [base->getVertex(0)->getNum()] , 
169
 
                                       sizes [base->getVertex(1)->getNum()] , 
170
 
                                       sizes [base->getVertex(2)->getNum()] ,
171
 
                                       1./3., 1./3.); 
172
 
  Center_Circum_Aniso(metric.a11,metric.a21,metric.a22, xc, yc, circum_radius);
173
 
  //  printf("metric %g %g %g circum_radius = %g\n",metric.a11,metric.a21,metric.a22,circum_radius);
174
 
}
175
 
 
176
 
int inCircumCircle ( MTri3 *t , const double *p , const gmsh2dMetric &metric ) 
177
 
{
178
 
  double xc, yc;
179
 
  double eps, d1, d2, x[2];
180
 
  t->Center_Circum_Aniso(metric.a11,metric.a21,metric.a22,xc, yc, d1);
181
 
 
182
 
  x[0] = xc - p[0];
183
 
  x[1] = yc - p[1];
184
 
 
185
 
  d2 = sqrt(fabs(metric.eval(x[0],x[1])));
186
 
 
187
 
  eps = fabs(d1 - d2) / (d1 + d2);
188
 
  if(eps < 1.e-12) {
189
 
    return 0; 
190
 
  }
191
 
  if(d2 < d1)
192
 
    return 1;
193
 
  return 0;
194
 
 
195
 
}
196
 
 
197
 
int MTri3::inCircumCircle ( const double *p ) const
198
 
{
199
 
 
200
 
    double pa[2] = {base->getVertex(0)->x(),base->getVertex(0)->y()};
201
 
    double pb[2] = {base->getVertex(1)->x(),base->getVertex(1)->y()};
202
 
    double pc[2] = {base->getVertex(2)->x(),base->getVertex(2)->y()};
203
 
    double result = gmsh::incircle(pa, pb, pc, (double*)p) * gmsh::orient2d(pa, pb, pc);  
204
 
    return (result > 0) ? 1 : 0;
205
 
 
206
 
 
207
 
  double eps, d1, d2, x[2];
208
 
 
209
 
  x[0] = xc - p[0];
210
 
  x[1] = yc - p[1];
211
 
 
212
 
  d1 = circum_radius;
213
 
  d2 = sqrt(fabs(metric.eval(x[0],x[1])));
214
 
 
215
 
  eps = fabs(d1 - d2) / (d1 + d2);
216
 
  if(eps < 1.e-12) {
217
 
    return 0; 
218
 
  }
219
 
  if(d2 < d1)
220
 
    return 1;
221
 
  return 0;
222
 
}
223
 
 
224
 
struct edgeXface
225
 
{
226
 
  MVertex *v[2];
227
 
  MTri3 * t1;
228
 
  int i1;
229
 
  edgeXface ( MTri3 *_t, int iFac)
230
 
    : t1(_t),i1(iFac)
231
 
  {
232
 
    v[0] = t1->tri()->getVertex ( iFac == 0 ? 2 : iFac-1 );
233
 
    v[1] = t1->tri()->getVertex ( iFac );
234
 
    std::sort ( v, v+2 );
235
 
  }
236
 
  inline bool operator < ( const edgeXface & other) const
237
 
  {
238
 
    if (v[0] < other.v[0])return true;
239
 
    if (v[0] > other.v[0])return false;
240
 
    if (v[1] < other.v[1])return true;
241
 
    return false;
242
 
  }
243
 
};
 
213
 
 
214
  // compute the circumradius of the triangle
 
215
  double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(), base->getVertex(0)->z()};
 
216
  double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
 
217
  double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
 
218
  double center[3];
 
219
  circumCenterXYZ(pa, pb, pc, center);
 
220
  const double dx = base->getVertex(0)->x() - center[0];
 
221
  const double dy = base->getVertex(0)->y() - center[1];
 
222
  const double dz = base->getVertex(0)->z() - center[2];
 
223
  
 
224
  circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
 
225
  circum_radius /= lc;
 
226
}
 
227
 
 
228
int MTri3::inCircumCircle(const double *p) const
 
229
{
 
230
  double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(), base->getVertex(0)->z()};
 
231
  double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
 
232
  double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
 
233
  double fourth[3];
 
234
  fourthPoint(pa, pb, pc, fourth);
 
235
 
 
236
  double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) * 
 
237
    gmsh::orient3d(pa, pb, pc,fourth);  
 
238
  return (result > 0) ? 1 : 0;  
 
239
}
 
240
 
 
241
int inCircumCircle(MTriangle *base, 
 
242
                   const double *p, const double *param ,
 
243
                   std::vector<double> &Us, std::vector<double> &Vs) 
 
244
{
 
245
  double pa[2] = {Us[base->getVertex(0)->getNum()],
 
246
                  Vs[base->getVertex(0)->getNum()]};
 
247
  double pb[2] = {Us[base->getVertex(1)->getNum()],
 
248
                  Vs[base->getVertex(1)->getNum()]};
 
249
  double pc[2] = {Us[base->getVertex(2)->getNum()],
 
250
                  Vs[base->getVertex(2)->getNum()]};
 
251
 
 
252
  double result = gmsh::incircle(pa, pb, pc, (double*)param) * gmsh::orient2d(pa, pb, pc);  
 
253
  return (result > 0) ? 1 : 0;  
 
254
}
244
255
 
245
256
template <class ITER>
246
 
void connectTris ( ITER beg, ITER end)
 
257
void connectTris(ITER beg, ITER end)
247
258
{
248
259
  std::set<edgeXface> conn;
249
 
  {
250
 
    while (beg != end)
251
 
      {
252
 
        if (!(*beg)->isDeleted())
253
 
          {
254
 
            for (int i=0;i<3;i++)
255
 
              {
256
 
                edgeXface fxt ( *beg ,  i );
257
 
                std::set<edgeXface>::iterator found = conn.find (fxt);
258
 
                if (found == conn.end())
259
 
                  conn.insert ( fxt );
260
 
                else if (found->t1 != *beg)
261
 
                  {
262
 
                    found->t1->setNeigh(found->i1,*beg);
263
 
                    (*beg)->setNeigh ( i, found->t1);
264
 
                  }
265
 
              }
266
 
          }
267
 
        ++beg;
 
260
  while (beg != end){
 
261
    if (!(*beg)->isDeleted()){
 
262
      for (int i = 0; i < 3; i++){
 
263
        edgeXface fxt(*beg, i);
 
264
        std::set<edgeXface>::iterator found = conn.find(fxt);
 
265
        if (found == conn.end())
 
266
          conn.insert(fxt);
 
267
        else if (found->t1 != *beg){
 
268
          found->t1->setNeigh(found->i1, *beg);
 
269
          (*beg)->setNeigh(i, found->t1);
 
270
        }
268
271
      }
269
 
  }
270
 
}
271
 
 
272
 
 
273
 
void recurFindCavity (const gmsh2dMetric &metric, 
274
 
                      std::list<edgeXface> & shell, 
275
 
                      std::list<MTri3*> & cavity, 
276
 
                      double *v, 
277
 
                      MTri3 *t)
278
 
{
279
 
  t->setDeleted(true);
280
 
  // the cavity that has to be removed
281
 
  // because it violates delaunay criterion
282
 
  cavity.push_back(t);
283
 
 
284
 
  for (int i=0;i<3;i++)
285
 
    {
286
 
      MTri3 *neigh =  t->getNeigh(i) ;
287
 
      if (!neigh)
288
 
          shell.push_back ( edgeXface ( t, i ) );
289
 
      else  if (!neigh->isDeleted())
290
 
        {
291
 
          int circ =  inCircumCircle ( neigh, v, metric );
292
 
          if (circ)
293
 
            recurFindCavity ( metric,shell, cavity,v , neigh);
294
 
          else
295
 
            shell.push_back ( edgeXface ( t, i ) );
296
 
        }
297
 
    }
298
 
}
299
 
 
300
 
bool insertVertex (MVertex *v , 
301
 
                   MTri3   *t ,
302
 
                   std::set<MTri3*,compareTri3Ptr> &allTets,
303
 
                   std::vector<gmsh2dMetric> & vSizes,
304
 
                   const gmsh2dMetric &metric)
305
 
{
306
 
  std::list<edgeXface>  shell;
307
 
  std::list<MTri3*>  cavity; 
308
 
  std::list<MTri3*>  new_cavity;
309
 
 
310
 
  double p[2] = {v->x(),v->y()};
311
 
 
312
 
  recurFindCavity ( metric,shell,  cavity, p , t);  
313
 
 
 
272
    }
 
273
    ++beg;
 
274
  }
 
275
}
 
276
 
 
277
void connectTriangles(std::list<MTri3*> &l)
 
278
{
 
279
  connectTris(l.begin(), l.end());
 
280
}
 
281
 
 
282
void connectTriangles(std::vector<MTri3*> &l)
 
283
{
 
284
  connectTris(l.begin(), l.end());
 
285
}
 
286
 
 
287
void connectTriangles(std::set<MTri3*, compareTri3Ptr> &l)
 
288
{
 
289
  connectTris(l.begin(), l.end());
 
290
}
 
291
 
 
292
void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity, 
 
293
                     double *v, double *param, MTri3 *t,
 
294
                     std::vector<double> &Us, std::vector<double> &Vs)
 
295
{
 
296
  t->setDeleted(true);
 
297
  // the cavity that has to be removed
 
298
  // because it violates delaunay criterion
 
299
  cavity.push_back(t);
 
300
 
 
301
  for (int i = 0; i < 3; i++){
 
302
    MTri3 *neigh =  t->getNeigh(i) ;
 
303
    if (!neigh)
 
304
      shell.push_back(edgeXface(t, i));
 
305
    else if (!neigh->isDeleted()){
 
306
      int circ =  inCircumCircle(neigh->tri(), v , param, Us, Vs);
 
307
      if (circ)
 
308
        recurFindCavity(shell, cavity, v, param, neigh, Us, Vs);
 
309
      else
 
310
        shell.push_back(edgeXface(t, i));
 
311
    }
 
312
  }
 
313
}
 
314
 
 
315
void recurFindCavityAniso (GFace *gf,
 
316
                           std::list<edgeXface> &shell, std::list<MTri3*> &cavity, 
 
317
                           double *metric, double *param,  MTri3 *t,
 
318
                           std::vector<double> &Us, std::vector<double> &Vs)
 
319
{
 
320
  t->setDeleted(true);
 
321
  // the cavity that has to be removed
 
322
  // because it violates delaunay criterion
 
323
  cavity.push_back(t);
 
324
 
 
325
  for (int i = 0; i < 3; i++){
 
326
    MTri3 *neigh = t->getNeigh(i) ;
 
327
    if (!neigh)
 
328
      shell.push_back(edgeXface(t, i));
 
329
    else  if (!neigh->isDeleted()){
 
330
      int circ =  inCircumCircleAniso(gf, neigh->tri(), param, metric, Us, Vs);
 
331
      if (circ)
 
332
        recurFindCavityAniso(gf, shell, cavity,metric, param, neigh, Us, Vs);
 
333
      else
 
334
        shell.push_back(edgeXface(t, i));
 
335
    }
 
336
  }
 
337
}
 
338
 
 
339
bool circUV(MTriangle *t, std::vector<double> & Us, std::vector<double> &Vs,
 
340
            double *res, GFace *gf)
 
341
{
 
342
  double u1 [3] = {Us[t->getVertex(0)->getNum()], Vs[t->getVertex(0)->getNum()], 0};
 
343
  double u2 [3] = {Us[t->getVertex(1)->getNum()], Vs[t->getVertex(1)->getNum()], 0};
 
344
  double u3 [3] = {Us[t->getVertex(2)->getNum()], Vs[t->getVertex(2)->getNum()], 0};
 
345
  circumCenterXY(u1, u2, u3, res);
 
346
  return true;
 
347
  double p1 [3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
 
348
  double p2 [3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
 
349
  double p3 [3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
 
350
  double resxy[3], uv[2];
 
351
  circumCenterXYZ(p1, p2, p3, resxy,uv);
 
352
  return true;
 
353
}
 
354
 
 
355
bool invMapUV(MTriangle *t, double *p,
 
356
              const std::vector<double> &Us, const std::vector<double> &Vs,
 
357
              double *uv, double tol)
 
358
{
 
359
  double mat[2][2];
 
360
  double b[2];
 
361
 
 
362
  double u0 = Us[t->getVertex(0)->getNum()];
 
363
  double v0 = Vs[t->getVertex(0)->getNum()];
 
364
  double u1 = Us[t->getVertex(1)->getNum()];
 
365
  double v1 = Vs[t->getVertex(1)->getNum()];
 
366
  double u2 = Us[t->getVertex(2)->getNum()];
 
367
  double v2 = Vs[t->getVertex(2)->getNum()];
 
368
 
 
369
  mat[0][0] = u1 - u0;
 
370
  mat[0][1] = u2 - u0;
 
371
  mat[1][0] = v1 - v0;
 
372
  mat[1][1] = v2 - v0;
 
373
 
 
374
  b[0] = p[0] - u0;
 
375
  b[1] = p[1] - v0;
 
376
  sys2x2(mat, b, uv);
 
377
 
 
378
  if(uv[0] >= -tol && 
 
379
     uv[1] >= -tol && 
 
380
     uv[0] <= 1. + tol && 
 
381
     uv[1] <= 1. + tol && 
 
382
     1. - uv[0] - uv[1] > -tol) {
 
383
    return true;
 
384
  }
 
385
  return false; 
 
386
}
 
387
 
 
388
double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
 
389
{
 
390
  double u1 = Us[t->getVertex(0)->getNum()];
 
391
  double v1 = Vs[t->getVertex(0)->getNum()];
 
392
  double u2 = Us[t->getVertex(1)->getNum()];
 
393
  double v2 = Vs[t->getVertex(1)->getNum()];
 
394
  double u3 = Us[t->getVertex(2)->getNum()];
 
395
  double v3 = Vs[t->getVertex(2)->getNum()];
 
396
  const double vv1 [2] = {u2 - u1, v2 - v1};
 
397
  const double vv2 [2] = {u3 - u1, v3 - v1};
 
398
  double s = vv1[0] * vv2[1] - vv1[1] * vv2[0]; 
 
399
  return s * 0.5;
 
400
}
 
401
 
 
402
bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
 
403
                  std::set<MTri3*, compareTri3Ptr> &allTets,
 
404
                  std::set<MTri3*, compareTri3Ptr> *activeTets,
 
405
                  std::vector<double> &vSizes, std::vector<double> &vSizesBGM,
 
406
                  std::vector<double> &Us, std::vector<double> &Vs,
 
407
                  double *metric = 0)
 
408
{
 
409
  std::list<edgeXface> shell;
 
410
  std::list<MTri3*> cavity; 
 
411
  std::list<MTri3*> new_cavity;
 
412
 
 
413
  if (!metric){
 
414
    double p[3] = {v->x(), v->y(), v->z()};
 
415
    recurFindCavity(shell, cavity, p, param, t, Us, Vs);
 
416
  }
 
417
  else{
 
418
    recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);  
 
419
  }
 
420
  
314
421
  // check that volume is conserved
315
422
  double newVolume = 0;
316
423
  double oldVolume = 0;
317
424
 
318
425
  std::list<MTri3*>::iterator ittet = cavity.begin();
319
426
  std::list<MTri3*>::iterator ittete = cavity.end();
320
 
  while ( ittet != ittete )
321
 
    {
322
 
      oldVolume += fabs((*ittet)->getSurfaceXY());
323
 
      ++ittet;
324
 
    }
325
 
 
326
 
//    char name[245];
327
 
//    int III = 1;
328
 
//    sprintf(name,"test%d.pos",III);
329
 
//    FILE *ff = fopen (name,"w");
330
 
//    fprintf(ff,"View\"test\"{\n");
331
 
 
332
 
 
333
 
  MTri3** newTris = new MTri3*[shell.size()];;
 
427
  while(ittet != ittete){
 
428
    oldVolume += fabs(getSurfUV((*ittet)->tri(),Us,Vs));
 
429
    ++ittet;
 
430
  }
 
431
  
 
432
  MTri3 **newTris = new MTri3*[shell.size()];
334
433
  int k = 0;
335
 
  std::list<edgeXface>::iterator it =   shell.begin();
336
 
 
337
 
  while (it != shell.end())
338
 
    {
339
 
      MTriangle *t = new MTriangle ( it->v[0],
340
 
                                     it->v[1],
341
 
                                     v);
342
 
      MTri3 *t4 = new MTri3 ( t , vSizes); 
343
 
//       fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
344
 
//            it->v[0]->x(),
345
 
//            it->v[0]->y(),
346
 
//            it->v[0]->z(),
347
 
//            it->v[1]->x(),
348
 
//            it->v[1]->y(),
349
 
//            it->v[1]->z(),
350
 
//            v->x(),
351
 
//            v->y(),
352
 
//            v->z());
353
 
      newTris[k++]=t4;
354
 
      // all new triangles are pushed front in order to
355
 
      // ba able to destroy them if the cavity is not
356
 
      // star shaped around the new vertex.
357
 
      new_cavity.push_back (t4);
358
 
      MTri3 *otherSide = it->t1->getNeigh(it->i1);
359
 
 
360
 
      if (otherSide)
361
 
        new_cavity.push_back (otherSide);
362
 
      //      if (!it->t1->isDeleted())throw;
363
 
      double ss = fabs(t4->getSurfaceXY());
364
 
      if (ss < 1.e-12)ss = 1256172121;
365
 
      newVolume += ss;
366
 
      ++it;
367
 
    }
368
 
//    fprintf(ff,"};\n");
369
 
//    fclose (ff);
370
 
//   printf("%d %d newVolume %g oldVolume %g\n",cavity.size(),new_cavity.size(),newVolume,oldVolume);
371
 
//    getchar();
372
 
 
373
 
 
374
 
  if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume )
375
 
    {      
376
 
      connectTris ( new_cavity.begin(),new_cavity.end() );      
377
 
      allTets.insert(newTris,newTris+shell.size());
378
 
      delete [] newTris;
379
 
      return true;
380
 
    }
 
434
  std::list<edgeXface>::iterator it = shell.begin();
 
435
 
 
436
 
 
437
  bool onePointIsTooClose = false;
 
438
  while (it != shell.end()){
 
439
    MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
 
440
    double lc = 0.3333333333 * (vSizes[t->getVertex(0)->getNum()] +
 
441
                                vSizes[t->getVertex(1)->getNum()] +
 
442
                                vSizes[t->getVertex(2)->getNum()]);
 
443
    double lcBGM = 0.3333333333 * (vSizesBGM[t->getVertex(0)->getNum()] +
 
444
                                   vSizesBGM[t->getVertex(1)->getNum()] +
 
445
                                   vSizesBGM[t->getVertex(2)->getNum()]);
 
446
    double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
 
447
    MTri3 *t4 = new MTri3(t, LL); 
 
448
 
 
449
     double d1 = sqrt((it->v[0]->x()-v->x())*(it->v[0]->x()-v->x())+
 
450
                     (it->v[0]->y()-v->y())*(it->v[0]->y()-v->y())+
 
451
                     (it->v[0]->z()-v->z())*(it->v[0]->z()-v->z()));
 
452
     double d2 = sqrt((it->v[1]->x()-v->x())*(it->v[1]->x()-v->x())+
 
453
                     (it->v[1]->y()-v->y())*(it->v[1]->y()-v->y())+
 
454
                     (it->v[1]->z()-v->z())*(it->v[1]->z()-v->z()));
 
455
     if (d1 < LL*.25 ||d2 < LL*.25)onePointIsTooClose = true;
 
456
                     
 
457
    //    if (t4->getRadius () < LIMIT_ / 2)onePointIsTooClose = true;
 
458
 
 
459
    newTris[k++] = t4;
 
460
    // all new triangles are pushed front in order to
 
461
    // ba able to destroy them if the cavity is not
 
462
    // star shaped around the new vertex.
 
463
    new_cavity.push_back (t4);
 
464
    MTri3 *otherSide = it->t1->getNeigh(it->i1);
 
465
    if (otherSide)
 
466
      new_cavity.push_back (otherSide);
 
467
    double ss = fabs(getSurfUV(t4->tri(),Us,Vs));
 
468
    if (ss < 1.e-25)ss = 1256172121;
 
469
    newVolume += ss;
 
470
    ++it;
 
471
  }
 
472
  if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 && !onePointIsTooClose){      
 
473
    connectTris(new_cavity.begin(),new_cavity.end());      
 
474
    allTets.insert(newTris,newTris+shell.size());
 
475
    if (activeTets){
 
476
      for (std::list<MTri3*>::iterator i = new_cavity.begin();i!=new_cavity.end();++i){
 
477
        int active_edge;
 
478
        if(isActive(*i,LIMIT_,active_edge) && (*i)->getRadius() > LIMIT_){
 
479
          if ((*activeTets).find(*i) == (*activeTets).end())
 
480
            (*activeTets).insert(*i);
 
481
        }
 
482
      }
 
483
    }
 
484
    delete [] newTris;
 
485
    return true;
 
486
  }
381
487
  // The cavity is NOT star shaped
382
 
  else
383
 
    {      
384
 
      for (unsigned int i=0;i<shell.size();i++)delete newTris[i];
385
 
      delete [] newTris;      
386
 
      ittet = cavity.begin();
387
 
      ittete = cavity.end();  
388
 
      while ( ittet != ittete )
389
 
        {
390
 
          (*ittet)->setDeleted(false);
391
 
          ++ittet;
392
 
        }
393
 
      return false;
394
 
    }
395
 
}
396
 
 
397
 
static gmsh2dMetric evalMetricTensor ( GFace *gf , MVertex *mv , double lc )
398
 
{
399
 
  Pair<SVector3,SVector3> der = gf->firstDer(SPoint2(mv->x(),mv->y())) ;
400
 
  const double oneoverlc2 = 1./(lc*lc);
401
 
  const double a11 = dot(der.first() ,der.first() ) * oneoverlc2; 
402
 
  const double a21 = dot(der.second(),der.first() ) * oneoverlc2; 
403
 
  const double a22 = dot(der.second(),der.second()) * oneoverlc2; 
404
 
 
405
 
  //  printf("metric : %g %g %g - %g %g\n",a11,a21,a22,oneoverlc2,lc);
406
 
 
407
 
  return gmsh2dMetric ( a11, a21 , a22);
408
 
}
409
 
 
410
 
static void setLcs ( MTriangle *t, std::map<MVertex*,double> &vSizes)
411
 
{
412
 
  for (int i=0;i<3;i++)
413
 
    {
414
 
      for (int j=i+1;j<3;j++)
415
 
        {
416
 
          MVertex *vi = t->getVertex(i);
417
 
          MVertex *vj = t->getVertex(j);
418
 
          double dx = vi->x()-vj->x();
419
 
          double dy = vi->y()-vj->y();
420
 
          double dz = vi->z()-vj->z();
421
 
          double l = sqrt(dx*dx+dy*dy+dz*dz);
422
 
          std::map<MVertex*,double>::iterator iti = vSizes.find(vi);      
423
 
          std::map<MVertex*,double>::iterator itj = vSizes.find(vj);      
424
 
          if (iti==vSizes.end() || iti->second > l)vSizes[vi] = l;
425
 
          if (itj==vSizes.end() || itj->second > l)vSizes[vj] = l;
426
 
        }
427
 
    }
428
 
}
429
 
 
430
 
void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
431
 
{
432
 
 
 
488
  else{      
 
489
    for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
 
490
    delete [] newTris;      
 
491
    ittet = cavity.begin();
 
492
    ittete = cavity.end();  
 
493
    while(ittet != ittete){
 
494
      (*ittet)->setDeleted(false);
 
495
      ++ittet;
 
496
    }
 
497
    return false;
 
498
  }
 
499
}
 
500
 
 
501
void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
 
502
                std::vector<double> &Us, std::vector<double> &Vs, bool param=true)
 
503
{
 
504
  FILE *ff = fopen (name,"w");
 
505
  fprintf(ff,"View\"test\"{\n");
 
506
  std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
 
507
  while (it != AllTris.end() ){
 
508
    MTri3 *worst = *it;
 
509
    if (!worst->isDeleted()){
 
510
      if (param)
 
511
        fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
 
512
                Us [(worst)->tri()->getVertex(0)->getNum()],
 
513
                Vs [(worst)->tri()->getVertex(0)->getNum()],
 
514
                0.0,
 
515
                Us [(worst)->tri()->getVertex(1)->getNum()],
 
516
                Vs [(worst)->tri()->getVertex(1)->getNum()],
 
517
                0.0,
 
518
                Us [(worst)->tri()->getVertex(2)->getNum()],
 
519
                Vs [(worst)->tri()->getVertex(2)->getNum()],
 
520
                0.0);
 
521
      else
 
522
        fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
 
523
                (worst)->tri()->getVertex(0)->x(),
 
524
                (worst)->tri()->getVertex(0)->y(),
 
525
                (worst)->tri()->getVertex(0)->z(),
 
526
                (worst)->tri()->getVertex(1)->x(),
 
527
                (worst)->tri()->getVertex(1)->y(),
 
528
                (worst)->tri()->getVertex(1)->z(),
 
529
                (worst)->tri()->getVertex(2)->x(),
 
530
                (worst)->tri()->getVertex(2)->y(),
 
531
                (worst)->tri()->getVertex(2)->z(),
 
532
                (worst)->getRadius(),
 
533
                (worst)->getRadius(),
 
534
                (worst)->getRadius()
 
535
                );
 
536
    }
 
537
    ++it;
 
538
  }
 
539
  fprintf(ff,"};\n");
 
540
  fclose (ff);
 
541
}
 
542
 
 
543
 
 
544
static void insertAPoint(GFace *gf,                      
 
545
                         std::set<MTri3*,compareTri3Ptr>::iterator it, 
 
546
                         double center[2],
 
547
                         double metric[3], 
 
548
                         std::vector<double> &Us, 
 
549
                         std::vector<double> &Vs,
 
550
                         std::vector<double> &vSizes, 
 
551
                         std::vector<double> &vSizesBGM,
 
552
                         std::set<MTri3*,compareTri3Ptr> &AllTris,
 
553
                         std::set<MTri3*,compareTri3Ptr> * ActiveTris = 0,
 
554
                         MTri3 *worst = 0){
 
555
  if (worst){
 
556
    it = AllTris.find(worst);
 
557
    if (worst != *it)throw;
 
558
  }
 
559
  else worst = *it;
 
560
 
 
561
  MTri3 *ptin = 0;
 
562
  double uv[2];
 
563
  MTriangle *base = worst->tri();
 
564
  bool inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);    
 
565
  if (inside)ptin = worst;
 
566
  if (!inside && worst->getNeigh(0)){
 
567
    inside |= invMapUV(worst->getNeigh(0)->tri(), center, Us, Vs, uv, 1.e-8);
 
568
    if (inside)ptin = worst->getNeigh(0);
 
569
  }
 
570
  if (!inside && worst->getNeigh(1)){
 
571
    inside |= invMapUV(worst->getNeigh(1)->tri(), center, Us, Vs, uv, 1.e-8);
 
572
    if (inside)ptin = worst->getNeigh(1);
 
573
  }
 
574
  if (!inside && worst->getNeigh(2)){
 
575
    inside |= invMapUV(worst->getNeigh(2)->tri(), center, Us, Vs, uv, 1.e-8);
 
576
    if (inside)ptin = worst->getNeigh(2);
 
577
  }
 
578
  if (inside) {
 
579
    // we use here local coordinates as real coordinates
 
580
    // x,y and z will be computed hereafter
 
581
    //        Msg(INFO,"Point is inside");
 
582
    GPoint p = gf->point(center[0], center[1]);
 
583
    //    printf("the new point is %g %g\n",p.x(),p.y());
 
584
    MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
 
585
    v->setNum(Us.size());
 
586
    double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] + 
 
587
                  uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] + 
 
588
                  uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]); 
 
589
    // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
 
590
    double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
 
591
    vSizesBGM.push_back(lc);
 
592
    vSizes.push_back(lc1);
 
593
    Us.push_back(center[0]);
 
594
    Vs.push_back(center[1]);
 
595
    
 
596
    if (!insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM, 
 
597
                      Us, Vs, metric)) {
 
598
      Msg(DEBUG2,"2D Delaunay : a cavity is not star shaped");
 
599
      AllTris.erase(it);
 
600
      worst->forceRadius(-1);
 
601
      AllTris.insert(worst);                    
 
602
      delete v;
 
603
    }
 
604
    else 
 
605
      gf->mesh_vertices.push_back(v);
 
606
  }
 
607
  else {
 
608
    Msg(DEBUG2,"Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
 
609
        center[0], center[1],
 
610
        Us[base->getVertex(0)->getNum()], 
 
611
        Vs[base->getVertex(0)->getNum()], 
 
612
        Us[base->getVertex(1)->getNum()], 
 
613
        Vs[base->getVertex(1)->getNum()], 
 
614
        Us[base->getVertex(2)->getNum()], 
 
615
        Vs[base->getVertex(2)->getNum()], 
 
616
        metric[0], metric[1], metric[2]);
 
617
    AllTris.erase(it);
 
618
    worst->forceRadius(0);
 
619
    AllTris.insert(worst);
 
620
  }
 
621
}
 
622
 
 
623
void gmshBowyerWatson(GFace *gf)
 
624
{
 
625
  if (CTX.mesh.algo2d == ALGO_2D_FRONTAL){
 
626
    gmshBowyerWatsonFrontal(gf);
 
627
    return;
 
628
  }
433
629
  std::set<MTri3*,compareTri3Ptr> AllTris;
434
 
  std::map<MVertex*,double> vSizesMap;
435
 
  std::map<const BDS_Point*,MVertex*> recoverMap;
436
 
  std::vector<gmsh2dMetric> vMetrics;
437
 
  std::vector<double> vSizes;
438
 
 
439
 
 
440
 
  // compute edge sizes on the contour and propagate thos sizes inside the surface
441
 
  for (unsigned int i=0;i<gf->triangles.size();i++)setLcs ( gf->triangles[i] , vSizesMap);
442
 
 
443
 
  // set vertex positions to local coordinates of the face  
444
 
  // compute metric tensors of the surface
445
 
  int NUM=0;
446
 
  for (std::map<MVertex*,double>::iterator it = vSizesMap.begin();it!=vSizesMap.end();++it)
447
 
    {
448
 
      BDS_Point *p;
449
 
      p = bds->find_point(it->first->getNum());
450
 
      it->first->x() = p->u;
451
 
      it->first->y() = p->v;
452
 
      it->first->z() = 0;
453
 
      it->first->setNum(NUM++);      
454
 
      vSizes.push_back(it->second);
455
 
      recoverMap[p]= (it->first);
456
 
      vMetrics.push_back(evalMetricTensor ( gf , it->first, it->second));
457
 
    }
458
 
  
459
 
  for (unsigned int i=0;i<gf->triangles.size();i++)
460
 
    AllTris.insert ( new MTri3 ( gf->triangles[i] ,vMetrics ) );
461
 
 
462
 
  gf->triangles.clear();
463
 
  connectTris ( AllTris.begin(), AllTris.end() );      
464
 
  
465
 
  Msg(DEBUG,"All %d tris were connected",AllTris.size());
466
 
 
467
 
  // here the classification should be done
 
630
  std::vector<double> vSizes, vSizesBGM, Us, Vs;
 
631
 
 
632
  buidMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
 
633
 
 
634
  // _printTris ("before.pos", AllTris, Us,Vs);
 
635
  int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
 
636
  // _printTris ("after2.pos", AllTris, Us,Vs);
 
637
  Msg(DEBUG2,"Delaunization of the initial mesh done (%d swaps)", nbSwaps);
468
638
 
469
639
  int ITER = 0;
470
 
 
471
 
  while (1)
472
 
    {
473
 
      MTri3 *worst = *AllTris.begin();
474
 
 
475
 
      if (worst->isDeleted())
476
 
        {
477
 
          delete worst->tri();
478
 
          delete worst;
479
 
          AllTris.erase(AllTris.begin());
480
 
          //      Msg(INFO,"Worst tet is deleted");
481
 
        }
482
 
      else
483
 
        {
484
 
          if(ITER++%5000 ==0)
485
 
            Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
486
 
          if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
487
 
          double center[2],uv[2];
488
 
          worst->getCenter(center);
489
 
          //      worst->tri()->circumcenterXY(center);
490
 
          bool inside = worst->tri()->invertmappingXY(center,uv);
491
 
          if (inside || 1)
492
 
            {
493
 
              // we use here local coordinates as real coordinates
494
 
              // x,y and z will be computed hereafter
495
 
              //              Msg(INFO,"Point is inside");
496
 
              MVertex *v = new MFaceVertex (center[0],center[1],0.0,gf,center[0],center[1]);
497
 
              v->setNum(NUM++);
498
 
              double lc1 = (1.-uv[0]-uv[1]) * vSizes [worst->tri()->getVertex(0)->getNum()] + 
499
 
                uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] + 
500
 
                uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()] ; 
501
 
              GPoint p = gf->point (center[0],center[1]);
502
 
              double lc = std::min(lc1,BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z()));
503
 
              vSizes.push_back(lc);
504
 
              gmsh2dMetric metr = evalMetricTensor ( gf , v, lc); 
505
 
              vMetrics.push_back(metr);
506
 
              
507
 
              // compute mesh spacing there
508
 
              if (!insertVertex ( v  , worst, AllTris,vMetrics,metr))
509
 
                {
510
 
                  AllTris.erase(AllTris.begin());
511
 
                  worst->forceRadius(0);
512
 
                  AllTris.insert(worst);                  
513
 
                  delete v;
514
 
                }
515
 
              else 
516
 
                {
517
 
                  gf->mesh_vertices.push_back(v);
518
 
                }
519
 
            }
520
 
          else
521
 
            {
522
 
              //              Msg(INFO,"Point is outside");
523
 
              AllTris.erase(AllTris.begin());
524
 
              worst->forceRadius(0);
525
 
              AllTris.insert(worst);              
526
 
            }
527
 
        }
528
 
    }
529
 
 
530
 
  // restore real coordinates of vertices
531
 
  if (1){
532
 
    std::map<const BDS_Point*,MVertex*>::const_iterator itx =  recoverMap.begin();
533
 
    while (itx != recoverMap.end())
534
 
      {
535
 
        MVertex *v = (itx->second);
536
 
        const BDS_Point *p = (itx->first);
537
 
        v->x() = p->X;
538
 
        v->y() = p->Y;
539
 
        v->z() = p->Z;
540
 
        ++itx;
541
 
      }
542
 
    for (unsigned int i=0;i<gf->mesh_vertices.size();i++)
543
 
      {
544
 
        double u = gf->mesh_vertices[i]->x();
545
 
        double v = gf->mesh_vertices[i]->y();
546
 
        GPoint p = gf->point(u,v);
547
 
        gf->mesh_vertices[i]->x() = p.x();
548
 
        gf->mesh_vertices[i]->y() = p.y();
549
 
        gf->mesh_vertices[i]->z() = p.z();
550
 
      }
551
 
  }
552
 
 
553
 
  // fill new gmsh structures with triangles
554
 
  while (1)
555
 
    {
556
 
      if (AllTris.begin() == AllTris.end() ) break;
557
 
      MTri3 *worst = *AllTris.begin();
558
 
      if (worst->isDeleted())
559
 
        {
560
 
          delete worst->tri();
561
 
        }
562
 
      else
563
 
        {
564
 
          gf->triangles.push_back(worst->tri());
565
 
        }
 
640
  while (1){
 
641
    MTri3 *worst = *AllTris.begin();
 
642
    if (worst->isDeleted()){
 
643
      delete worst->tri();
566
644
      delete worst;
567
 
      AllTris.erase(AllTris.begin());      
568
 
    }
569
 
}
 
645
      AllTris.erase(AllTris.begin());
 
646
    }
 
647
    else{
 
648
      if(ITER++ % 5000 == 0)
 
649
        Msg(DEBUG1,"%7d points created -- Worst tri radius is %8.3f",
 
650
            vSizes.size(), worst->getRadius());
 
651
      double center[2],metric[3],r2;
 
652
      if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
 
653
      circUV(worst->tri(), Us, Vs, center, gf);
 
654
      MTriangle *base = worst->tri();
 
655
      double pa[2] = {(Us[base->getVertex(0)->getNum()] + 
 
656
                       Us[base->getVertex(1)->getNum()] + 
 
657
                       Us[base->getVertex(2)->getNum()]) / 3.,
 
658
                      (Vs[base->getVertex(0)->getNum()] + 
 
659
                       Vs[base->getVertex(1)->getNum()] + 
 
660
                       Vs[base->getVertex(2)->getNum()]) / 3.};
 
661
      buildMetric(gf, pa, metric);
 
662
      circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);       
 
663
      insertAPoint(gf,AllTris.begin(),center,metric,Us,Vs,vSizes,vSizesBGM,AllTris);
 
664
    }
 
665
  }    
 
666
  transferDataStructure(gf, AllTris); 
 
667
}
 
668
 
 
669
/*
 
670
  Let's try a frontal delaunay approach now that the delaunay mesher is stable
 
671
 
 
672
  We use the approach of Rebay (JCP 1993) that we extend to anisotropic metrics
 
673
 
 
674
  The point isertion is done on the Voronoi Edges
 
675
  
 
676
*/
 
677
 
 
678
static double length_metric ( const double p[2], const double q[2], const double metric[3])
 
679
{
 
680
  return sqrt (   (p[0]-q[0]) * metric [0] * (p[0]-q[0]) +
 
681
                2*(p[0]-q[0]) * metric [1] * (p[1]-q[1]) +
 
682
                  (p[1]-q[1]) * metric [2] * (p[1]-q[1]) );
 
683
}
 
684
 
 
685
/*
 
686
          /|
 
687
         / |
 
688
        /  |
 
689
       /   |
 
690
   lc /    |  r
 
691
     /     |
 
692
    /      |  
 
693
   /       x
 
694
  /        |    
 
695
 /         |  r/2
 
696
/          |  
 
697
-----------+
 
698
     lc/2
 
699
 
 
700
     (3 r/2)^2 = lc^2 - lc^2/4 
 
701
     -> lc^2 3/4 = 9r^2/4
 
702
     -> lc^2 = 3 r^2
 
703
 
 
704
     r^2 /4 + lc^2/4 = r^2
 
705
     -> lc^2 = 3 r^2
 
706
     
 
707
*/
 
708
 
 
709
 
 
710
void gmshBowyerWatsonFrontal(GFace *gf){
 
711
  std::set<MTri3*,compareTri3Ptr> AllTris;
 
712
  std::set<MTri3*,compareTri3Ptr> ActiveTris;
 
713
  std::vector<double> vSizes, vSizesBGM, Us, Vs;
 
714
  buidMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
 
715
 
 
716
  // delaunise the initial mesh
 
717
  int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
 
718
  Msg(DEBUG2,"Delaunization of the initial mesh done (%d swaps)", nbSwaps);
 
719
  
 
720
  int ITER = 0, active_edge;
 
721
  // compute active triangle
 
722
  std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
 
723
  for ( ; it!=AllTris.end();++it){
 
724
    if(isActive(*it,LIMIT_,active_edge))
 
725
      ActiveTris.insert(*it);
 
726
    else if ((*it)->getRadius() < LIMIT_)break;
 
727
  }
 
728
  
 
729
  // insert points
 
730
  while (1){
 
731
    if (!ActiveTris.size())break;
 
732
    MTri3 *worst = (*ActiveTris.begin());
 
733
    ActiveTris.erase(ActiveTris.begin());
 
734
    //    printf("active_tris.size = %d\n",ActiveTris.size());
 
735
 
 
736
    if (!worst->isDeleted() && isActive(worst,LIMIT_,active_edge) && worst->getRadius() > LIMIT_){      
 
737
      if(ITER++ % 5000 == 0)
 
738
        Msg(DEBUG1,"%7d points created -- Worst tri radius is %8.3f",
 
739
            vSizes.size(), worst->getRadius());
 
740
      // compute circum center of that guy
 
741
      double center[2],metric[3],r2;
 
742
      MTriangle *base = worst->tri();
 
743
      circUV(base, Us, Vs, center, gf);
 
744
      double pa[2] = {(Us[base->getVertex(0)->getNum()] + 
 
745
                       Us[base->getVertex(1)->getNum()] + 
 
746
                       Us[base->getVertex(2)->getNum()]) / 3.,
 
747
                      (Vs[base->getVertex(0)->getNum()] + 
 
748
                       Vs[base->getVertex(1)->getNum()] + 
 
749
                       Vs[base->getVertex(2)->getNum()]) / 3.};
 
750
      buildMetric(gf, pa, metric);
 
751
      circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2); 
 
752
      // compute the middle point of the edge
 
753
      int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
 
754
      int ip2 = active_edge;
 
755
      //     printf("the active edge is %d : %g %g -> %g %g\n",active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
 
756
      //           base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
 
757
      double P[2] =  {Us[base->getVertex(ip1)->getNum()],  
 
758
                      Vs[base->getVertex(ip1)->getNum()]};
 
759
      double Q[2] =  {Us[base->getVertex(ip2)->getNum()], 
 
760
                      Vs[base->getVertex(ip2)->getNum()]};
 
761
      double midpoint[2] =  {0.5*(P[0]+Q[0]),0.5*(P[1]+Q[1])};
 
762
      
 
763
      // now we have the edge center and the center of the circumcircle, 
 
764
      // we try to find a point that would produce a perfect triangle while
 
765
      // connecting the 2 points of the active edge
 
766
      
 
767
      double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
 
768
      double norm = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
 
769
      dir[0]/=norm;
 
770
      dir[1]/=norm;
 
771
      const double RATIO = sqrt(  dir[0]*dir[0]*metric[0]+
 
772
                                  2*dir[1]*dir[0]*metric[1]+
 
773
                                  dir[1]*dir[1]*metric[2]);    
 
774
      
 
775
      
 
776
      const double p    = 0.5*length_metric (P,Q,metric);// / RATIO;
 
777
      const double q    = length_metric (center,midpoint,metric);
 
778
      const double rhoM1 = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
779
      const double rhoM2 = 0.5 * (vSizesBGM[base->getVertex(ip1)->getNum()] + vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
780
      const double rhoM  = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
 
781
      //      const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
782
      
 
783
      const double rhoM_hat = std::min(std::max(rhoM,p),(p*p+q*q)/(2*q));
 
784
      const double d = (rhoM_hat + sqrt (rhoM_hat*rhoM_hat - p*p))/RATIO;
 
785
      
 
786
      double newPoint[2] = 
 
787
        {
 
788
          midpoint[0] + d * dir[0],
 
789
          midpoint[1] + d * dir[1]
 
790
        };
 
791
      insertAPoint(gf,AllTris.end(),newPoint,metric,Us,Vs,vSizes,vSizesBGM,AllTris,&ActiveTris,worst);
 
792
    } 
 
793
  }
 
794
 
 
795
  char name[245];
 
796
  sprintf(name,"frontal%d.pos",gf->tag());
 
797
  //  _printTris (name, AllTris, Us,Vs);
 
798
  transferDataStructure(gf, AllTris); 
 
799
} 
 
800
 
 
801
 
 
802
/*
 
803
void gmshBowyerWatsonFrontal(GFace *gf){
 
804
 
 
805
  std::set<MTri3*,compareTri3Ptr> AllTris;
 
806
  std::vector<double> vSizes, vSizesBGM, Us, Vs;
 
807
  buidMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
 
808
 
 
809
  // delaunise the initial mesh
 
810
  int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
 
811
  Msg(DEBUG2,"Delaunization of the initial mesh done (%d swaps)", nbSwaps);
 
812
  
 
813
 
 
814
  int ITER = 0, active_edge;
 
815
  while (1) {
 
816
    // compute active triangle
 
817
    std::set<MTri3*,compareTri3Ptr> ActiveTris;
 
818
    std::set<MTri3*,compareTri3Ptr>::reverse_iterator it = AllTris.rbegin();
 
819
    for ( ; it!=AllTris.rend();++it){
 
820
      if ((*it)->getRadius() < LIMIT_)
 
821
        (*it)->done = true;
 
822
      else if(isActive(*it,LIMIT_,active_edge))
 
823
        ActiveTris.insert(*it);
 
824
    }
 
825
    if (ActiveTris.size() == 0)break;
 
826
    
 
827
    // insert points
 
828
    while (1){
 
829
      if (!ActiveTris.size())break;
 
830
      MTri3 *worst = (*ActiveTris.begin());
 
831
      ActiveTris.erase(ActiveTris.begin());
 
832
      //    printf("active_tris.size = %d\n",ActiveTris.size());
 
833
      
 
834
      if (!worst->isDeleted() && isActive(worst,LIMIT_,active_edge) && worst->getRadius() > LIMIT_){      
 
835
        if(ITER++ % 5000 == 0)
 
836
          Msg(DEBUG1,"%7d points created -- Worst tri radius is %8.3f",
 
837
              vSizes.size(), worst->getRadius());
 
838
        // compute circum center of that guy
 
839
        double center[2],uv[2],metric[3],r2;
 
840
        MTriangle *base = worst->tri();
 
841
        circUV(base, Us, Vs, center, gf);
 
842
        double pa[2] = {(Us[base->getVertex(0)->getNum()] + 
 
843
                         Us[base->getVertex(1)->getNum()] + 
 
844
                         Us[base->getVertex(2)->getNum()]) / 3.,
 
845
                        (Vs[base->getVertex(0)->getNum()] + 
 
846
                         Vs[base->getVertex(1)->getNum()] + 
 
847
                         Vs[base->getVertex(2)->getNum()]) / 3.};
 
848
        buildMetric(gf, pa, metric);
 
849
        circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2); 
 
850
        // compute the middle point of the edge
 
851
        int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
 
852
        int ip2 = active_edge;
 
853
        //     printf("the active edge is %d : %g %g -> %g %g\n",active_edge,base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),
 
854
        //         base->getVertex(ip2)->x(),base->getVertex(ip2)->y());
 
855
        double P[2] =  {Us[base->getVertex(ip1)->getNum()],  
 
856
                        Vs[base->getVertex(ip1)->getNum()]};
 
857
        double Q[2] =  {Us[base->getVertex(ip2)->getNum()], 
 
858
                        Vs[base->getVertex(ip2)->getNum()]};
 
859
        double midpoint[2] =  {0.5*(P[0]+Q[0]),0.5*(P[1]+Q[1])};
 
860
        
 
861
        // now we have the edge center and the center of the circumcircle, 
 
862
        // we try to find a point that would produce a perfect triangle while
 
863
        // connecting the 2 points of the active edge
 
864
        
 
865
        double dir[2] = {center[0]-midpoint[0],center[1]-midpoint[1]};
 
866
        double norm = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
 
867
        dir[0]/=norm;
 
868
        dir[1]/=norm;
 
869
        const double RATIO = sqrt(  dir[0]*dir[0]*metric[0]+
 
870
                                    2*dir[1]*dir[0]*metric[1]+
 
871
                                    dir[1]*dir[1]*metric[2]);    
 
872
        
 
873
        
 
874
        const double p    = 0.5*length_metric (P,Q,metric);// / RATIO;
 
875
        const double q    = length_metric (center,midpoint,metric);
 
876
        const double rhoM1 = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
877
        const double rhoM2 = 0.5 * (vSizesBGM[base->getVertex(ip1)->getNum()] + vSizesBGM[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
878
        const double rhoM  = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
 
879
        //      const double rhoM = 0.5 * (vSizes[base->getVertex(ip1)->getNum()] + vSizes[base->getVertex(ip2)->getNum()] ) / sqrt(3.);// * RATIO;
 
880
        
 
881
        const double rhoM_hat = std::min(std::max(rhoM,p),(p*p+q*q)/(2*q));
 
882
        const double d = (rhoM_hat + sqrt (rhoM_hat*rhoM_hat - p*p))/RATIO;
 
883
        
 
884
        double newPoint[2] = 
 
885
          {
 
886
            midpoint[0] + d * dir[0],
 
887
            midpoint[1] + d * dir[1]
 
888
          };
 
889
        insertAPoint(gf,AllTris.end(),newPoint,metric,Us,Vs,vSizes,vSizesBGM,AllTris,&ActiveTris,worst);
 
890
      } 
 
891
    }
 
892
  }
 
893
  _printTris ("frontal.pos", AllTris, Us,Vs);
 
894
  transferDataStructure(gf, AllTris); 
 
895
} 
 
896
*/
 
897
 
 
898