metric.cpp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· C++ 代码 · 共 1,214 行 · 第 1/3 页

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//

/*
 
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include <stdio.h>
#include "Meshio.h"
#include "Mesh2.h"


namespace bamg {


inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3])
{ return    A[0] * ( B[1]*C[2]-B[2]*C[1])
          - A[1] * ( B[0]*C[2]-B[2]*C[0])
          + A[2] * ( B[0]*C[1]-B[1]*C[0]);
}

SaveMetricInterpole  LastMetricInterpole;

void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) 
{
  double a11=M1.a11,a21=M1.a21,a22=M1.a22;
  double b11=M2.a11,b21=M2.a21,b22=M2.a22;
  //  M1 v = l M2 v
  // (M1 - l M2) v =0
  // det (M1 - l M2) =0
  // det (M1 - l M2) = a l^2 + b l + c;
  // = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
  //  const double eps = 1.e-5;
  const double /*c11 = a11*a11,*/ c21= a21*a21;
  const double /*d11 = b11*b11,*/ d21= b21*b21;
  const double a=b11*b22 - d21;
  const double b=-a11*b22-a22*b11+2*a21*b21;
  const double c=-c21+a11*a22;
  const double bb = b*b,ac= a*c;
  const double delta = bb - 4 * ac;
  //  const double kk=c11+c22+c21+d11+d21+d22;
  // modif F Hecht feb 1998 
  // cerr.precision(14);
  //cerr  <<  bb << " " << ac << " " <<  bb <<  " " <<a << endl;
  // cerr << a11 << " " << a21 << " " << a22 << endl;
  //cerr << b11 << " " << b21 << " " << b22 << endl;
  if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) )
   {
   // racine double;
     // cerr << "double " << endl ;
    if (Abs(a) < 1.e-30 )
     l1 = l2 = 0;
    else 
     l1=l2=-b/(2*a); 
    V= D2xD2(1,0,0,1);
   }
  else {
    // cerr << " -- " << a << endl ;
    const double delta2 = sqrt(delta);
    l1= (-b - delta2)/(2*a);
    l2= (-b + delta2)/(2*a);
    // M1 v = l M2 v
    //  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
    // so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
      double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
      double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
      double vp1x,vp1y,vp2x,vp2y;

      if(s1 < s0)
	s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
      else
	s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;

      v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
      s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
      if(s1 < s0)
	s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
      else
	s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
#ifdef DEBUG
      assert(Abs(vp1y)+Abs(vp2y)>0);
#endif
      V=D2xD2(vp1x,vp2x,vp1y,vp2y);
  }
  return;

}

MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) 
{
      D2xD2 M;
      double l1,l2;
      ReductionSimultanee(M1,M2,l1,l2,M);
      R2 v0(M.x.x,M.y.x);
      R2 v1(M.x.y,M.y.y);
      D2xD2 M_1(M.inv());
      D2xD2 D(Max(M1(v0,v0),M2(v0,v0)),0,0,Max(M1(v1,v1),M2(v1,v1)));
      D2xD2 Mi(M_1.t()*D*M_1);
      return MetricAnIso(Mi.x.x,0.5*(Mi.x.y+Mi.y.x),Mi.y.y);
}

MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
	   const  MetricAnIso m1,const  MetricAnIso m2 )
{
  MetricAnIso mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11,
		  a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21,
		  a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22);
  
  MatVVP2x2 vab(mab);
 
  R2 v1(vab.v.x,vab.v.y);
  R2 v2(-v1.y,v1.x);
  
  Real8 h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1);
  Real8 h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2);

  vab.lambda1 =  1 / (h1*h1);
  vab.lambda2 =  1 / (h2*h2);
  *this = vab;
}

 MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma,
	                   Real8  b,const  MetricAnIso mb)
{ 
  MetricAnIso mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
  MatVVP2x2 vab(mab);
  
  R2 v1(vab.v.x,vab.v.y);
  R2 v2(-v1.y,v1.x);
  

  Real8 h1 = a / ma(v1) + b / mb(v1);
  Real8 h2 = a / ma(v2) + b / mb(v2);
  vab.lambda1 =  1 / (h1*h1);
  vab.lambda2 =  1 / (h2*h2);
  *this = vab;
}



 MatVVP2x2::MatVVP2x2(const MetricAnIso M) 
{
  double a11=M.a11,a21=M.a21,a22=M.a22;
  const double eps = 1.e-5;
  double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
  double b=-a11-a22,c=-c21+a11*a22;
  double   delta = b*b - 4 * c ;
  double n2=(c11+c22+c21);
  if ( n2 < 1e-30) 
     lambda1=lambda2=0,v.x=1,v.y=0;
  else if (delta < eps*n2)
    { 
      lambda1=lambda2=-b/2, v.x=1,v.y=0;
    }
  else 
    {  //    ---  construction  de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 
      delta = sqrt(delta);
      lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
      double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
      double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
 
      if(s1 < s0)
	s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
      else
	s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
    };
}


 int MetricAnIso::IntersectWith(const MetricAnIso M2) 
{
  //cerr << " - " << *this << M2 <<  endl;
      int r=0;
      MetricAnIso & M1 = *this;
      D2xD2 M;
      double l1,l2;
      
      ReductionSimultanee(*this,M2,l1,l2,M);
      // cerr << M << endl;
      R2 v1(M.x.x,M.y.x);
      R2 v2(M.x.y,M.y.y);
      double l11=M1(v1,v1);
      double l12=M1(v2,v2);
      double l21=M2(v1,v1);
      double l22=M2(v2,v2);
      if ( l11 < l21 )  r=1,l11=l21;
      if ( l12 < l22 )  r=1,l12=l22; 
      // cerr << r << endl;
      if (r) { // change
        D2xD2 M_1(M.inv());
        D2xD2 D(l11,0,0,l12); 
        D2xD2 Mi(M_1.t()*D*M_1);
        a11=Mi.x.x;
        a21=0.5*(Mi.x.y+Mi.y.x);
        a22=Mi.y.y; }
      return r;
}
void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0)

{
  if(verbosity>1)
    cout << " -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso "  ) << endl;
  Real8 ss[2]={0.00001,0.99999};
  Real8 errC = 2*sqrt(2*err);
  Real8 hmax = Gh.MaximalHmax();
  Real8 hmin = Gh.MinimalHmin();
  Real8 maxaniso = 1e6;
  assert(hmax>0);
  SetVertexFieldOn();
  if (errC > 1) errC = 1;
  for (Int4  i=0;i<nbe;i++)
   for (int j=0;j<2;j++)
    {
      
      Vertex V;
      VertexOnGeom GV;
      // cerr << Number(edges[i]) << " " << ss[j] << endl;
      Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
	{
	  GeometricalEdge * eg = GV;
	  Real8 s = GV;
	  R2 tg;
	  //	   cerr << i << " " << j << " " << Number(V) << " on = " 
	  //	<< Gh.Number(eg) << " at s = " << s << " " << endl;
	  Real8  R1= eg->R1tg(s,tg);
	  // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x " 
	  //    << V.r << errC/ Max(R1,1e-20) <<  " hold=" <<V.m(tg) << " "  << endl;
	  Real8 ht = hmax;
          if (R1>1.0e-20) 
	    {  // err relative to the length of the edge
	      ht = Min(Max(errC/R1,hmin),hmax);
	    }
	  Real8 hn = iso? ht : Min(hmax,ht*maxaniso);
	  //cerr << ht << " " << hn << "m=" << edges[i][j].m <<  endl;
	  assert(ht>0 && hn>0);
	  MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg);
	  //cerr << " : " ;
	  Metric MVp(Vp);
	  // cerr << " : "  << MVp  << endl;
	  edges[i][j].m.IntersectWith(MVp);
	  //cerr << " . " << endl;
	}

    }
  // the problem is for the vertex on vertex 
     
}
/*
void  Triangles::BoundAnisotropy(Real8 anisomax)
{
  if (verbosity > 1) 
    cout << " -- BoundAnisotropy by  " << anisomax << endl; 
  Real8 h1=1.e30,h2=1e-30,rx=0;
  Real8 coef = 1./(anisomax*anisomax);
  Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
  for (Int4 i=0;i<nbv;i++)
    {

      MatVVP2x2 Vp(vertices[i]);
      
      h1=Min(h1,Vp.lmin());
      h2=Max(h2,Vp.lmax());
      rx = Max(rx,Vp.Aniso2());
      
      Vp.BoundAniso2(coef);
      
      hn1=Min(hn1,Vp.lmin());
      hn2=Max(hn2,Vp.lmax());
      rnx = Max(rnx,Vp.Aniso2());

      
      vertices[i].m = Vp;

    }

  if (verbosity>2)
    {
      cout << "     input :  Hmin = " << sqrt(1/h2)  << " Hmax = " << sqrt(1/h1) 
	   << " factor of anisotropy max  = " << sqrt(rx) << endl;
      cout << "     output:  Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 
	   << " factor of anisotropy max  = " << sqrt(rnx) << endl;
    }
}
*/
void  Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso)
{
  double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100));
  if (verbosity > 1) 
    cout << " -- BoundAnisotropy by  " << anisomax << endl; 
  Real8 h1=1.e30,h2=1e-30,rx=0;
  Real8 coef = 1./(anisomax*anisomax);
  Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
  for (Int4 i=0;i<nbv;i++)
    {

      MatVVP2x2 Vp(vertices[i]);
      double lmax=Vp.lmax();
      h1=Min(h1,Vp.lmin());
      h2=Max(h2,Vp.lmax());
      rx = Max(rx,Vp.Aniso2());

      Vp *= Min(lminaniso,lmax)/lmax;
      
      Vp.BoundAniso2(coef);
      
      hn1=Min(hn1,Vp.lmin());
      hn2=Max(hn2,Vp.lmax());
      rnx = Max(rnx,Vp.Aniso2());

      
      vertices[i].m = Vp;

    }

  if (verbosity>2)
    {
      cout << "     input :  Hmin = " << sqrt(1/h2)  << " Hmax = " << sqrt(1/h1) 
	   << " factor of anisotropy max  = " << sqrt(rx) << endl;
      cout << "     output:  Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 
	   << " factor of anisotropy max  = " << sqrt(rnx) << endl;
    }
}
void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols,
				    const  Real8 hmin1,const Real8 hmax1,const Real8 coef,
				    const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,
				    const int DoNormalisation,const double power,const int choice)
{ //  the array of solution s is store    
  // sol0,sol1,...,soln    on vertex 0
  //  sol0,sol1,...,soln   on vertex 1
  //  etc.
  //  choise = 0 =>  H is computed with green formule
  //   otherwise  => H is computed from P2 on 4T 
  const int dim = 2;
  
  int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 

  // computation of the nb of field 
  Int4 ntmp = 0;
  if (typsols)
    {
      for (Int4 i=0;i<nbsol;i++)
	     ntmp += sizeoftype[typsols[i]];
    }
  else
    ntmp = nbsol;

  // n is the total number of fields

  const Int4 n = ntmp;

  Int4 i,k,iA,iB,iC,iv;
  R2 O(0,0);
  int RelativeMetric = CutOff>1e-30;
  Real8 hmin = Max(hmin1,MinimalHmin());
  Real8 hmax = Min(hmax1,MaximalHmax());
  Real8 coef2 = 1/(coef*coef);

  if(verbosity>1) 
    {
      cout << " -- Construction of Metric: Nb of field. " << n << " nbt = " 
	   << nbt << " nbv= " << nbv 
	   << " coef = " << coef << endl
	   << "     hmin = " << hmin << " hmax=" << hmax 
	   << " anisomax = " << anisomax <<  " Nb Jacobi " << NbJacobi << " Power = " << power ;
      if (RelativeMetric)
	cout << " RelativeErr with CutOff= "  <<  CutOff << endl;
      else
	cout << " Absolute Err" <<endl;
    }
  double *ss=(double*)s;//, *ssiii = ss;

  double sA,sB,sC;

  Real8 *detT = new Real8[nbt];
  Real8 *Mmass= new Real8[nbv];
  Real8 *Mmassxx= new Real8[nbv];
  Real8 *dxdx= new Real8[nbv];
  Real8 *dxdy= new Real8[nbv];
  Real8 *dydy= new Real8[nbv];
  Real8 *workT= new Real8[nbt];
  Real8 *workV= new Real8[nbv];
  int *OnBoundary = new int[nbv];
  for (iv=0;iv<nbv;iv++)