metric.cpp

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

  Real8 hmin = Max(hmin1,MinimalHmin());
  Real8 hmax = Min(hmax1,MaximalHmax());
  MeshIstream f_metrix(fmetrix);
  Int4 k,j;
  f_metrix >>  k >> j ;
  if(verbosity>1)
    cout << " metrix: open " << fmetrix 
	 << ", le coef = " << coef
	 << ", hmin = " << hmin 
	 << ", hmax = " << hmax 
	 << (  (j == 1)? " Iso " : " AnIso " )<< endl;
  
  if (k != nbv || !(j == 1 || j == 3)) 
    {
      cerr << " Error Pb metrix " << k << " <> " 
	   <<  nbv << " or  1 or 3 <> " << j << endl;
      MeshError(1002);
    }
  
  cout << " j = " << j << endl;
  //  Int4 nberr = 0;
  for (Int4 iv=0;iv<nbv;iv++)
    {
      Real8 h;
      if (j == 1) 
	{
	f_metrix >>  h ;
	vertices[iv].m=Metric(Max(hmin,Min(hmax, h*coef)));
	}
      else if (j==3) 
	{
	  Real8 a,b,c;	     
	  f_metrix >>  a >> b >> c  ;
	  MetricAnIso M(a,b,c);
	  MatVVP2x2 Vp(M/coef);
	  
	  Vp.Maxh(hmin);
	  Vp.Minh(hmax);
	  vertices[iv].m = Vp;
	  
	}
    }
 
}

void Triangles::WriteMetric(ostream & f,int iso)
{
  if (iso)
    {
      f <<  nbv <<" " << 1 << endl ;
      for (Int4 iv=0;iv<nbv;iv++)
	{
	  MatVVP2x2 V=vertices[iv].m;
	  f <<  V.hmin()  << endl;
	}
    }
else
  {
    f <<  nbv <<" " << 3 << endl ;
    for (Int4 iv=0;iv<nbv;iv++)
      f <<  vertices[iv].m.a11 << " " 
	<<  vertices[iv].m.a21 << " " 
	<<  vertices[iv].m.a22 << endl;
  }
}
void  Triangles::MaxSubDivision(Real8 maxsubdiv)
{
const  Real8 maxsubdiv2 = maxsubdiv*maxsubdiv;
#ifdef DRAWING2
  inquire();
#endif	    
  if(verbosity>1)
    cout << " -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv <<   endl  ;
  // for all the edges 
  // if the len of the edge is to long 
  Int4 it,nbchange=0;    
  Real8 lmax=0;
  for (it=0;it<nbt;it++)
    {
      Triangle &t=triangles[it];
      for (int j=0;j<3;j++)
	{
	  Triangle &tt = *t.TriangleAdj(j);
	  if ( ! &tt ||  it < Number(tt) && ( tt.link || t.link)) 
	    {
		Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
		Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
		R2 AB= (R2) v1-(R2) v0;
		Metric M = v0;
		Real8 l = M(AB,AB);
		lmax = Max(lmax,l);
		if(l> maxsubdiv2)
		  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
		    Real8 lc = M(AC,AC);
		    D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
		    D2xD2 Rt1(Rt.inv());
		    D2xD2 D(maxsubdiv2,0,0,lc);
		    D2xD2 MM = Rt1*D*Rt1.t();
#ifdef DRAWING1
		    v0.m.Draw(v0);
#endif	    
		    v0.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
#ifdef DRAWING1
		    v0.m.Draw(v0);
#endif	    
		    //		    cout << " M(AB,AB) = " << M(AB,AB) << " == " << maxsubdiv 
		    //	 << " M(AC,AC) = " << M(AC,AC) << " == " << lc << endl; 
		    nbchange++;
		  }
		M = v1;
		l = M(AB,AB);
		lmax = Max(lmax,l);
		if(l> maxsubdiv2)
		  { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
		    Real8 lc = M(AC,AC);
		    D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
		    D2xD2 Rt1(Rt.inv());
		    D2xD2 D(maxsubdiv2,0,0,lc);
		    D2xD2  MM = Rt1*D*Rt1.t();
#ifdef DRAWING1
		    v1.m.Draw(v1);
#endif	    
		    v1.m =  M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y);
#ifdef DRAWING1
		    v1.m.Draw(v1);
		    inquire();
#endif	    
		    // cout << " M(AB,AB) = " << M(AB,AB) << " == " << maxsubdiv 
		    //	 << " M(AC,AC) = " << M(AC,AC) << " == " << lc << endl; 
		    nbchange++;
		  }
		
		
	    }
	}
    }
  if(verbosity>3)
  cout << "    Nb of metric change = " << nbchange 
       << " Max  of the subdivision of a edges before change  = " << sqrt(lmax) << endl;
#ifdef DRAWING2
  inquire();
#endif	    

}

void Triangles::SmoothMetric(Real8 raisonmax) 
{ 
  if(raisonmax<1.1) return;
  if(verbosity > 1)
     cout << " -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl;
  ReMakeTriangleContainingTheVertex();
  Int4 i,j,kch,kk,ip;
  Int4 *first_np_or_next_t0 = new Int4[nbv];
  Int4 *first_np_or_next_t1 = new Int4[nbv];
  Int4 Head0 =0,Head1=-1;
  Real8 logseuil= log(raisonmax);

  for(i=0;i<nbv-1;i++)
    first_np_or_next_t0[i]=i+1; 
  first_np_or_next_t0[nbv-1]=-1;// end;
  for(i=0;i<nbv;i++)
    first_np_or_next_t1[i]=-1;
  kk=0;
  while (Head0>=0&& kk++<100) {
    kch=0;
    for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1)
      {  //  pour tous les triangles autour du sommet s
	// 	cout << kk << " i = " << i << " " << ip << endl;
	register Triangle * t= vertices[i].t;
	assert(t);
	Vertex & vi = vertices[i];
	TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);
	Vertex *pvj0 = ta.EdgeVertex(0);
	while (1) {
	  //	  cout << i << " " <<  Number(ta.EdgeVertex(0)) << " "
	  //      << Number(ta.EdgeVertex(1)) << " ---> " ;
	  ta=Previous(Adj(ta));
	  // cout <<  Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl;
	  assert(vertices+i == ta.EdgeVertex(1));
	  Vertex & vj = *(ta.EdgeVertex(0));
	  if ( &vj ) {
	    j= &vj-vertices;
	    assert(j>=0 && j < nbv);
	    R2 Aij = (R2) vj - (R2) vi;
	    Real8 ll =  Norme2(Aij);
	    if (0) {  
	      Real8 hi = ll/vi.m(Aij);
	      Real8 hj = ll/vj.m(Aij);
	      if(hi < hj)
		{
		  Real8 dh=(hj-hi)/ll;
		  //cout << " dh = " << dh << endl;
		  if (dh>logseuil) {
		    vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
		    if(first_np_or_next_t1[j]<0)
		      kch++,first_np_or_next_t1[j]=Head1,Head1=j;
		  }
		}
	    } 
	    else
	      {
		Real8 li = vi.m(Aij);
		//Real8 lj = vj.m(Aij);
		//		if ( i == 2 || j == 2)
		//  cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) <<  endl;
	      	if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
		  //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) )
		  if(first_np_or_next_t1[j]<0) // if the metrix change 
		    kch++,first_np_or_next_t1[j]=Head1,Head1=j;
	      }
	  }
	  if  ( &vj ==  pvj0 ) break;
	}
      }
    Head0 = Head1;
    Head1 = -1;
    Exchange(first_np_or_next_t0,first_np_or_next_t1);
    if(verbosity>5)
    cout << "     Iteration " << kk << " Nb de  vertices with change  " << kch << endl;
  }
  if(verbosity>2 && verbosity < 5) 
    cout << "    Nb of Loop " << kch << endl;
  delete [] first_np_or_next_t0;
  delete [] first_np_or_next_t1;
}

void Geometry::ReadMetric(const char * fmetrix,Real8 hmin=1.0e-30,Real8 hmax=1.0e30,Real8 coef=1)
{
  hmin = Max(hmin,MinimalHmin());
  MeshIstream f_metrix(fmetrix);
  Int4 k,j;
  f_metrix >>  k >> j ;
  if(verbosity>1)
    cout << " -- ReadMetric  " << fmetrix 
	 << ",  coef = " << coef
	 << ", hmin = " << hmin 
	 << ", hmax = " << hmax 
	 << (  (j == 1)? " Iso " : " AnIso " ) << endl;
  
  if (k != nbv ||  !(j == 1 || j == 3)) {
    cerr << " Error Pb metrix " << k << " <> " 
	 <<  nbv << " or  1 or 3  <> " << j << endl;
    MeshError(1003);}
	 
  
  //  Int4 nberr = 0;
  for (Int4 iv=0;iv<nbv;iv++)
    {
    Real8 h;
    if (j == 1) 
      {
      f_metrix >>  h ;
      vertices[iv].m=Metric(Max(hmin,Min(hmax, h*coef)));
      }
    else if (j==3) 
      {
	Real8 a,b,c;	     
	f_metrix >>  a >> b >> c  ;
	MetricAnIso M(a,b,c);
	MatVVP2x2 Vp(M/coef);
	Vp.Maxh(hmin);
      Vp.Minh(hmax);
      vertices[iv].m = Vp;
      }
    }
  
}


Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB)
{
  Real8 k=1./2.;
  int level=0;
  static int kkk=0;
  static  Metric Ms1[32],Ms2[32];
  static Real8 lMs1[32],lMs2[32];
  static double K[32];
  Real8 l=0,sss=0;
  Ms1[level]=Ma;
  Ms2[level]=Mb;
  Real8 sa =  Ma(AB);
  Real8 sb =  Mb(AB);
  lMs1[level]=sa;
  lMs2[level]=sb;
  K[level]=k;
  level++;
  int i=0;
  Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
  Real8  sstop = 0.1; // Max(0.6,(sa+sb)/5000);
  while (level) {
    level--;
    Metric M1=Ms1[level];
    Metric M2=Ms2[level];
    k=K[level];
    Real8 s1=  lMs1[level];
    Real8 s2=  lMs2[level];

    Real8 s= (s1+s2)*k;
//    if (level >20  && i < 2030-level)
//    cout << "                  level " << level << " " << i << " " << s << " " << k <<endl;
    if( s > sstop   && level < 30 && i < 500-level ) {
      Metric Mi(0.5,M1,0.5,M2);
      Real8 si = Mi(AB);
      if( Abs((s1+s2)-(si+si)) > s1*0.001) 
	{
	  k=k/2;
	  // we begin by the end to walk in the correct sens from a to b
	       // due to the stack 
	       Ms1[level]=Mi;
	  Ms2[level]=M2;
	  lMs1[level]=si;
	  lMs2[level]=s2;
	  K[level]=k;
	  level++;
	  Ms1[level]=M1;
	  Ms2[level]=Mi;
	  lMs1[level]=s1;
	  lMs2[level]=si;
	  K[level]=k;
	  level++;
	}
      else
	L[i]= l += s,S[i]=sss+=k,i++;
    }
    else 
      L[i]= l += s,S[i]=sss+=k,i++;//cout << i << " l = " << l << " sss = " << sss << endl;
  }
  // warning for optimisation S is in [0:0.5] not in [0:1]
  assert(i<512);
  LastMetricInterpole.lab=l;
  LastMetricInterpole.opt=i;
  if (i>200 && kkk++<10)
     cout << "Warning LengthInterpole: ( i = " << i << " l = " << l << " sss " << sss << " ) " << sstop <<endl;
  return l;
}

Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim)
{ 
  if(!optim)  LengthInterpole(Ma,Mb,AB);
  Real8 l  = s* LastMetricInterpole.lab,r;
  int j=LastMetricInterpole.opt-1,i=0,k;
  
  Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
  // warning for optimisation S is the abcisse in [0:0.5]
  // and L is le lenght 
  if(l<=L[0])
    r=2*S[0]*l/L[0];
  else if (l>=L[j])
    r=1;
  else 
    {
      while (j-i>1)
	{
	  k= (i+j)/2;
	  if(l<=L[k])
	    j=k;// l<=L[j] 
	  else
	    i=k; //  L[i]<l
	};
      //   cout  << i << " " << j  <<" " << L[i] << " " << L[j] << " " << S[i] << " " <<  S[j]  << " l=" << l << endl;
      if (i==j)
	r = 2*S[i];
      else
	r =  2*(S[i]*(L[j]-l)+ S[j]*(l-L[i]))/(L[j]-L[i]);
    }
  assert(r<=1 && r>=0);
  return r ;
    
}


#ifdef DRAWING 

void MetricAnIso::Draw(R2 c) const 
{ 
  float x= c.x,y= c.y;
  if (InPtScreen(x,y)) {
  R2 X(cos(0.0),sin(0.0));
  X = X / operator()(X);
  rmoveto(x+X.x,y+X.y);
  
  for (int i=1;i<=100;i++)
   { double t= 2*Pi*i/100.0;
     R2 X(cos(t),sin(t));
     X = X / Max(operator()(X),1.0e-5);
     rlineto(x+X.x,y+X.y); } 
  }
}

void MetricIso::Draw(R2 c) const 
{ 
  float x= c.x,y= c.y;
  if (InPtScreen(x,y)) {
  rmoveto(x+h,y);
  for (int i=1;i<=40;i++)
   { double t= Pi*i/20.0;
     rlineto(x+h*cos(t),y+h*sin(t)); } 
  }
}


#endif
}   // end of namespace bamg