← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Numeric/Numeric.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-02 18:12:15 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090902181215-yla8zvcas2ucvkm9
[Christophe Prud'homme]
* New upstream release
  + fixed surface mesh orientation bug introduced in 2.4.0;
  + mesh and graphics code refactoring;
  + small usability enhancements and bug fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Numeric/Numeric.cpp
267
267
  return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
268
268
}
269
269
 
 
270
void circumCenterXY(double *p1, double *p2, double *p3, double *res)
 
271
{
 
272
  double d, a1, a2, a3;
 
273
 
 
274
  const double x1 = p1[0];
 
275
  const double x2 = p2[0];
 
276
  const double x3 = p3[0];
 
277
  const double y1 = p1[1];
 
278
  const double y2 = p2[1];
 
279
  const double y3 = p3[1];
 
280
 
 
281
  d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
 
282
  if(d == 0.0) {
 
283
    Msg::Warning("Colinear points in circum circle computation");
 
284
    res[0] = res[1] = -99999.;
 
285
    return ;
 
286
  }
 
287
 
 
288
  a1 = x1 * x1 + y1 * y1;
 
289
  a2 = x2 * x2 + y2 * y2;
 
290
  a3 = x3 * x3 + y3 * y3;
 
291
  res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
 
292
  res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
 
293
}
 
294
 
 
295
void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv)
 
296
{
 
297
  double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
 
298
  double v2[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
 
299
  double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
 
300
  double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
 
301
  double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
 
302
  norme(vx); norme(vy); norme(vz);
 
303
  double p1P[2] = {0.0, 0.0};
 
304
  double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
 
305
  double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
 
306
  double resP[2];
 
307
 
 
308
  circumCenterXY(p1P, p2P, p3P,resP);
 
309
 
 
310
  if(uv){
 
311
    double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
 
312
                        {p2P[1] - p1P[1], p3P[1] - p1P[1]}};
 
313
    double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
 
314
    sys2x2(mat, rhs, uv);
 
315
  }
 
316
  
 
317
  res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
 
318
  res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
 
319
  res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
 
320
}
 
321
 
270
322
char float2char(float f)
271
323
{
272
324
  // float normalized in [-1, 1], char in [-127, 127]
495
547
  gmshVector<double> f(N), feps(N), dx(N);
496
548
  
497
549
  for (int iter = 0; iter < MAXIT; iter++){
498
 
    func(x, f, data);
 
550
     func(x, f, data);
 
551
     
 
552
     bool isZero = false;
 
553
     for (int k=0; k<N; k++) {
 
554
         if (f(k) == 0. ) isZero = true;
 
555
         else isZero = false;
 
556
         if (isZero == false) break;
 
557
       }
 
558
     if (isZero) break;
 
559
    //printf("**** fffffff0=%g %g %g %g %g %g %g %g %g\n", f(0), f(1), f(2),  f(3), f(4),  f(5), f(6),  f(7), f(8));
499
560
 
500
561
    for (int j = 0; j < N; j++){
501
562
      double h = EPS * fabs(x(j));
519
580
  }
520
581
  return false;
521
582
}
 
583
 
 
584
/*
 
585
min_a f(x+a*d);
 
586
 
 
587
f(x+a*d) = f(x) + f'(x) (
 
588
 
 
589
*/
 
590
 
 
591
void gmshLineSearch(double (*func)(gmshVector<double> &, void *), void* data, 
 
592
                    gmshVector<double> &x, gmshVector<double> &p,  
 
593
                    gmshVector<double> &g, double &f, 
 
594
                    double stpmax, int &check)
 
595
{
 
596
  int i;
 
597
  double alam, alam2 = 1., alamin, f2 = 0., fold2 = 0., rhs1, rhs2, temp, tmplam = 0.;
 
598
 
 
599
  const double ALF = 1.0e-4;
 
600
  const double TOLX = 1.0e-9;
 
601
 
 
602
  gmshVector<double> xold(x);
 
603
  const double fold = (*func)(xold,data);
 
604
  
 
605
  check=0;
 
606
  int n = x.size();
 
607
  double norm = p.norm();
 
608
  if (norm > stpmax)p.scale(stpmax/norm);
 
609
  double slope=0.0;
 
610
  for (i=0;i<n;i++)slope += g(i)*p(i);
 
611
  double test=0.0;
 
612
  for (i=0;i<n;i++) {
 
613
    temp=fabs(p(i))/std::max(fabs(xold(i)),1.0);
 
614
    if (temp > test) test=temp;
 
615
  }
 
616
  /*
 
617
  for (int j=0;j<100;j++){
 
618
    double sx = (double)j/99;
 
619
    for (i=0;i<n;i++) x(i)=xold(i)+10*sx*p(i);    
 
620
    double jzede = (*func)(x,data); 
 
621
  }
 
622
  */
 
623
 
 
624
  alamin=TOLX/test;
 
625
  alam=1.0;
 
626
  while(1) {
 
627
    for (i=0;i<n;i++) x(i)=xold(i)+alam*p(i);
 
628
    f=(*func)(x,data);
 
629
    //    printf("f = %g x = %g %g alam = %g p = %g %g\n",f,x(0),x(1),alam,p(0),p(1));
 
630
   if (alam < alamin) {
 
631
      for (i=0;i<n;i++) x(i)=xold(i);
 
632
      //      printf("ALERT : alam %g alamin %g\n",alam,alamin);
 
633
      check=1;
 
634
      return;
 
635
    } 
 
636
    else if (f <= fold + ALF*alam*slope) return;
 
637
    else {
 
638
      if (alam == 1.0)
 
639
        tmplam = -slope/(2.0*(f-fold-slope));
 
640
      else {
 
641
        rhs1 = f-fold-alam*slope;
 
642
        rhs2=f2-fold2-alam2*slope;
 
643
        const double a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
 
644
        const double b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
 
645
        if (a == 0.0) tmplam = -slope/(2.0*b);
 
646
        else {
 
647
          const double disc=b*b-3.0*a*slope;
 
648
          if (disc<0.0) Msg::Error("Roundoff problem in gmshLineSearch.");
 
649
          else tmplam=(-b+sqrt(disc))/(3.0*a);
 
650
        }
 
651
        if (tmplam>0.5*alam)
 
652
          tmplam=0.5*alam;
 
653
      }
 
654
    }
 
655
    alam2=alam;
 
656
    f2 = f;
 
657
    fold2=fold;
 
658
    alam=std::max(tmplam,0.1*alam);
 
659
  }
 
660
}
 
661
 
 
662
double minimize_grad_fd (double (*func)(gmshVector<double> &, void *),
 
663
                       gmshVector<double> &x, void *data)
 
664
{
 
665
  const int MAXIT = 3;
 
666
  const double EPS = 1.e-4;
 
667
  const int N = x.size();
 
668
  
 
669
  gmshVector<double> grad(N);
 
670
  gmshVector<double> dir(N);
 
671
  double f,feps,finit;
 
672
 
 
673
  for (int iter = 0; iter < MAXIT; iter++){
 
674
    // compute gradient of func
 
675
    f = func(x,data);
 
676
    if (iter == 0)finit = f;
 
677
    //    printf("Opti iter %d x = (%g %g) f = %g\n",iter,x(0),x(1),f);
 
678
    //    printf("grad = (");
 
679
    for (int j = 0; j < N; j++){
 
680
      double h = EPS * fabs(x(j));
 
681
      if(h == 0.) h = EPS;
 
682
      x(j) += h;
 
683
      feps = func(x, data);
 
684
      grad(j) = (feps - f) / h;
 
685
      //      printf("%g ",grad(j));
 
686
      dir(j) = -grad(j);
 
687
      x(j) -= h;
 
688
    }
 
689
    //    printf(")\n ");
 
690
    // do a 1D line search to fine the minimum
 
691
    // of f(x - \alpha \nabla f)
 
692
    double f, stpmax=100000;
 
693
    int check;
 
694
    gmshLineSearch (func, data, x,dir,grad,f,stpmax,check);
 
695
    //    printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
 
696
    if (check == 1)break;
 
697
  }
 
698
  
 
699
  return f;
 
700
}
 
701