meshquad.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// E-MAIL   : [e-mail]
//
// ********** DO NOT REMOVE THIS BANNER **********
/*
 
 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
 */

//
// SUMMARY: Bamg: Bidimensional Anisotrope Mesh Generator
// RELEASE: 0 
// AUTHOR:   F. Hecht,    
// ORG    :  UMPC
// E-MAIL :   [e-mail]   
//
// ORIG-DATE:     frev 98
//---------------------------------------------------------
//  to make quad 
// -------------------

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "Meshio.h"
#include "Mesh2.h"
#include "QuadTree.h"
#include "SetOfE4.h"

namespace bamg {

static const  Direction NoDirOfSearch=Direction();

#ifdef __MWERKS__
#pragma global_optimizer on
#pragma optimization_level 1
//#pragma inline_depth 0
#endif

class DoubleAndInt4 {
public: double q; Int4 i3j;
  int operator<(DoubleAndInt4 a){return q > a.q;}
};// to sort by decroissant



template<class T>  inline void  HeapSort(T *c,long n)
{
  long l,j,r,i;
  T crit;
  c--; // on decale de 1 pour que le tableau commence a 1
  if( n <= 1) return;
  l = n/2 + 1;
  r = n;
  while (1) { // label 2
    if(l <= 1 ) { // label 20
      crit = c[r];
      c[r--] = c[1];
    if ( r == 1 ) { c[1]=crit; return;}
    } else  crit = c[--l]; 
    j=l;
    while (1) {// label 4
      i=j;
      j=2*j;
      if  (j>r) {c[i]=crit;break;} // L8 -> G2
      if ((j<r) && (c[j] < c[j+1])) j++; // L5
      if (crit < c[j]) c[i]=c[j]; // L6+1 G4
      else {c[i]=crit;break;} //L8 -> G2
    }
  }
}

Triangle * swapTest(Triangle *t1,Int2 a);
// 


double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d)
{
  // calcul de 4 angles --
  R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);
  R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);
  //  Move(A),Line(B),Line(C),Line(D),Line(A);
  const Metric & Ma  = a;
  const Metric & Mb  = b;
  const Metric & Mc  = c;
  const Metric & Md  = d;
    
  double lAB=Norme2(AB);
  double lBC=Norme2(BC);
  double lCD=Norme2(CD);
  double lDA=Norme2(DA);
  AB /= lAB;  BC /= lBC;  CD /= lCD;  DA /= lDA;
  // version aniso 
  double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);
  double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);
  double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);
  double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);
  double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));
  // cout << A << B << C << D ;
  // cout << " sinmin " << sinmin << " " << cosDAB << " " << cosABC << " " << cosBCD << " " << cosCDA << endl;
  // rattente(1);
  if (sinmin<=0) return sinmin;
  return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));
}

 GeometricalEdge *   Triangles::ProjectOnCurve( Edge & BhAB, Vertex &  vA, Vertex & vB,
						Real8 theta,
						Vertex & R,VertexOnEdge &  BR,VertexOnGeom & GR) 
                      
{
   void *pA=0,*pB=0;
   Real8 tA=0,tB=0;
   R2 A=vA,B=vB;
    Vertex * pvA=&vA, * pvB=&vB;
  if (vA.vint == IsVertexOnVertex)
    {
      //  cout << " Debut vertex = " << BTh.Number(vA.onbv) ;
      pA=vA.onbv;
    }
  else if (vA.vint == IsVertexOnEdge)
    {
      pA=vA.onbe->be;
      tA=vA.onbe->abcisse;
      // cout << " Debut edge = " << BTh.Number(vA.onbv) << " " << tA ;

    }
  else
    {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vA) << " " << endl;
     cerr << " forget call to SetVertexFieldOnBTh" << endl;
     MeshError(-1);
    } 
    
  if (vB.vint == IsVertexOnVertex)
    {
      // cout << " Fin vertex = " << BTh.Number(vB.onbv) << endl;
      pB=vB.onbv;
    }
  else if(vB.vint == IsVertexOnEdge)
    {
      pB=vB.onbe->be;
      tB=vB.onbe->abcisse;
      // cout << " Fin edge = " << BTh.Number(vB.onbe->be) << " " <<  tB ;

    }
  else
    {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vB) << " " << endl;
     cerr << " forget call to SetVertexFieldOnBTh" << endl;
     MeshError(-1);
    } 
  Edge * e = &BhAB;
 assert( pA && pB && e);
 // be carefull the back ground edge e is on same geom edge 
 // of the initiale edge def by the 2 vertex A B;
 assert(e>=BTh.edges && e<BTh.edges+BTh.nbe);// Is a background Mesh;   
// walk on BTh edge 
 //assert(0 /* not finish ProjectOnCurve with BackGround Mesh*/);
// 1 first find a back ground edge contening the vertex A
// 2 walk n back gound boundary to find the final vertex B


#ifdef DEBUG
// we suppose if the vertex A is on a background edge and
// the abcisse is 0 or 1 when the related point is not
// a end of curve <=>  !IsRequiredVertex
    if (vA.vint == IsVertexOnEdge)
      if (tA<=0)
	assert(! (*vA.onbe->be)[0].on->IsRequiredVertex());
      else if (tA>=1) 
      	assert(!(*vA.onbe->be)[1].on->IsRequiredVertex());
#endif

     if( vA.vint == IsVertexOnEdge) 
       { // find the start edge 
	     e = vA.onbe->be;	 

       } 
     else if (vB.vint == IsVertexOnEdge) 
       {
         theta = 1-theta;
	 Exchange(tA,tB);
	 Exchange(pA,pB);
	 Exchange(pvA,pvB);
	 Exchange(A,B);
	 e =  vB.onbe->be;

	 // cout << " EXCHANGE  A et B) " << endl;
       } 
     else
       { // do the search by walking 
	 assert(0 /* A FAIRE */);
       }

     // find the direction of walking with sens of edge and pA,PB;
   R2 AB=B-A;

   Real8 cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
   int kkk=0;
   int sens = (cosE01AB>0) ? 1 : 0;
 
   //   Real8 l=0; // length of the edge AB
   Real8 abscisse = -1;
 
   for (int cas=0;cas<2;cas++)
     {// 2 times algo:
       //    1 for computing the length l
       //    2 for find the vertex 
       int  iii;
        Vertex  *v0=pvA,*v1; 
       Edge *neee,*eee;
       Real8 lg =0; // length of the curve 
       Real8 te0;
       // we suppose take the curve's abcisse 
       // cout << kkk << " e = " << BTh.Number(e) << "  v0=  " 
       //    << BTh.Number(v0) << " v1 = " << BTh.Number((*e)[sens]) << endl;
       for ( eee=e,iii=sens,te0=tA;
             eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
             neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii )   
	      { 
		//	cout << kkk << " eee = " << BTh.Number(eee) << "  v0=  " 
		//     << BTh.Number(v0) << " v1 = " << BTh.Number(v1) << endl;
		
		assert(kkk++<100);
		assert(eee);
		Real8 lg0 = lg;
		Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
	     	lg += dp;
		if (cas && abscisse <= lg)
		  { // ok we find the geom edge 
		    Real8 sss  =   (abscisse-lg0)/dp;
		    Real8 thetab = te0*(1-sss)+ sss*iii;
		    assert(thetab>=0 && thetab<=1);
		    BR = VertexOnEdge(&R,eee,thetab);

		    // cout << Number(R) << " = " <<  thetab << " on  " <<  BTh.Number(eee)
		    //	 << " = " << R << endl;

		    return  Gh.ProjectOnCurve(*eee,thetab,R,GR);

		  }
	      }
       // we find the end 
       if (v1 != pvB) 
	 {
	   if (( void*) v1 == pB)
	      tB = iii;
	     
	   Real8 lg0 = lg;
	   assert(eee);
	   v1 = pvB;
	   Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
	   lg += dp;	
	   abscisse = lg*theta;
	   if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
	     { // ok we find the geom edge 
	       Real8 sss  =   (abscisse-lg0)/dp;
	       Real8 thetab = te0*(1-sss)+ sss*tB;
	       assert(thetab>=0 && thetab<=1);
	       BR = VertexOnEdge(&R,eee,thetab);
	      
	       //	cout << kkk << " eee = " << BTh.Number(eee) << "  v0=  " 
	       //     << BTh.Number(v0) << " " << te0
	       //      << " v1 = " << BTh.Number(v1) <<  " " << tB  << endl;

	       //out << Number(R) << " Opt  = " <<  thetab << " on  " <<  BTh.Number(eee) 
	       //	    << " = " << R << endl;

	       return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
	     }
	 }
       abscisse = lg*theta;
       
     }
   cerr << " Big Bug" << endl;
   MeshError(678);
   return 0; // just for the compiler 
      
}                  


void Triangles::MakeQuadrangles(double costheta)
{
  if (verbosity>2) 
    cout << " -- MakeQuadrangles costheta = " << costheta << endl;
  if (verbosity>5)  
    cout << "    (in)  Nb of Quadrilaterals = " << NbOfQuad 
	 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 
         << " Nb of outside triangles = " << NbOutT << endl;

  if (costheta >1) {
    if (verbosity>5)
      cout << "     do nothing costheta >1 "<< endl;
    return;}

  
  Int4 nbqq = (nbt*3)/2;
  DoubleAndInt4  *qq = new DoubleAndInt4[nbqq];

  Int4 i,ij;
  int j;
  Int4 k=0;
  for (i=0;i<nbt;i++)
    for (j=0;j<3;j++)
      if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)
 	   qq[k++].i3j=i*3+j;
//  sort  qq
  HeapSort(qq,k);
  
  Int4 kk=0;
  for (ij=0;ij<k;ij++)
    { 
      //      cout << qq[ij].q << " " << endl;
      i=qq[ij].i3j/3;
      j=(int) (qq[ij].i3j%3);
      // optisamition no float computation  
      if (triangles[i].QualityQuad(j,0) >=costheta) 
        triangles[i].SetHidden(j),kk++;
    }
  NbOfQuad = kk;
  if (verbosity>2) 
    {
    cout << "    (out)  Nb of Quadrilaterals = " << NbOfQuad 
	 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 
         << " Nb of outside triangles = " << NbOutT << endl;
    }
  delete [] qq;
#ifdef DRAWING2 
  Draw();
  inquire();
#endif

}
/*
Triangles::BThBoundary(Edge e,Real 8) const 
{
  // pointeur of the background must be on 
  // 
  Edge be = e.on;
}
*/
int  Triangles::SplitElement(int choice)
{
  Direction NoDirOfSearch;
  const  int withBackground = &BTh != this && &BTh;
 if (verbosity>2) 
   cout << " -- SplitElement " << (choice? " Q->4Q and T->4T " : " Q->4Q or T->3Q " ) << endl;;
 if (verbosity>5)
   cout << endl << "    (in)  Nb of Quadrilaterals = " << NbOfQuad 
	<< " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 
	<< " Nb of outside triangles = " << NbOutT << endl;
 
 ReNumberingTheTriangleBySubDomain();
#ifdef DRAWING2 
  Draw();
 inquire();
#endif
 //int nswap =0;
  const Int4 nfortria( choice ? 4 : 6);
    if(withBackground) 
    {
      BTh.SetVertexFieldOn();
      SetVertexFieldOnBTh();
    }
   else
     BTh.SetVertexFieldOn();
   
  Int4 newnbt=0,newnbv=0;
  Int4 * kedge = 0;
 Int4 newNbOfQuad=NbOfQuad;
  Int4 * ksplit= 0, * ksplitarray=0;
  Int4 kkk=0;
  int ret =0;
  if (nbvx<nbv+nbe) return 1;//   
  Triangles *  OCurrentTh= CurrentTh;
  CurrentTh = this;
  // 1) create  the new points by spliting the internal edges 
  // set the 
  Int4 nbvold = nbv;
  Int4 nbtold = nbt;
  Int4 NbOutTold  = NbOutT;
  Int4  NbEdgeOnGeom=0;
  Int4 i;
  
  nbt = nbt - NbOutT; // remove all the  the ouside triangles 
  Int4 nbtsave = nbt;
  Triangle * lastT = triangles + nbt;
  for (i=0;i<nbe;i++)
    if(edges[i].on) NbEdgeOnGeom++;
  Int4 newnbe=nbe+nbe;
  //  Int4 newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
  Int4 newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
  // Int4 newNbVertexOnBThVertex=NbVertexOnBThVertex;
  Int4 newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;

  // do allocation for pointeur to the geometry and background
  VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
  VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
  if (NbVerticesOnGeomEdge)
    memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
  if (NbVertexOnBThEdge)
    memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
  Edge *newedges = new Edge [newnbe];
  //  memcpy(newedges,edges,sizeof(Edge)*nbe);
  SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
#ifdef DEBUG
  for (i=0;i<nbt;i++)
    triangles[i].check();
#endif
#ifdef DRAWING2  
  reffecran();
#endif  
  Int4 k=nbv;
  Int4 kk=0;
  Int4 kvb = NbVertexOnBThEdge;
  Int4 kvg = NbVerticesOnGeomEdge;
  Int4 ie =0;
  Edge ** edgesGtoB=0;
  if (withBackground)
   edgesGtoB= BTh.MakeGeometricalEdgeToEdge();
  Int4 ferr=0;
  for (i=0;i<nbe;i++)
     newedges[ie].on=0;

  for (i=0;i<nbe;i++)
    {
      GeometricalEdge *ong =  edges[i].on;

      newedges[ie]=edges[i];
      newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
      newedges[ie].adj[1]=newedges + ie +1;
      R2 A = edges[i][0],B = edges[i][1];
      // cout << " ie = " << ie <<"  v0 = " <<  Number(newedges[ie][0]) << endl;


      kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
      if (ong) // a geometrical edges 
	{ 
	  if (withBackground)
	    {
	      // walk on back ground mesh 
	      //  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 
	      // a faire -- difficile 
	      // the first PB is to now a background edge between the 2 vertices
	      assert(edgesGtoB); 
	      // cout << " ie = " << ie <<"  v0 = " <<  Number(newedges[ie][0]) << endl;
             ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].on)],
			     edges[i][0],edges[i][1],0.5,vertices[k],
                             newVertexOnBThEdge[kvb],
                             newVerticesOnGeomEdge[kvg++]);
	     vertices[k].ReferenceNumber= edges[i].ref;
	     vertices[k].DirOfSearch =   NoDirOfSearch;        
;