Skip to content
Snippets Groups Projects
Select Git revision
  • 8ec22431a23af2496b4bc2837b5843e557098e41
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

compute_normal_form.c

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    delaunay_refinement.cpp 13.73 KiB
    #ifdef _OPENMP
    #include <omp.h>
    #endif
    #include <stack>
    #include <set>
    #include <vector>
    #include <algorithm>
    #include <math.h>
    #include "GmshMessage.h"
    #include "OS.h"
    #include "SPoint3.h"
    #include "delaunay3d_private.h"
    #include "delaunay3d.h"
    #include "rtree.h"
    #include "MVertex.h"
    #include "MTetrahedron.h"
    #include "MTriangle.h"
    #include "GRegion.h"
    #include "GFace.h"
    #include "SBoundingBox3d.h"
    
    typedef  std::set< Edge > edgeContainer2;
    
    long int AVGSEARCH;
    
    struct edgeContainer
    {
      std::set< Edge > _hash2;
      std::vector<std::vector<Edge> > _hash;
      edgeContainer (unsigned int N = 1000000) {
        _hash.resize(N);
      }
      bool addNewEdge2 (const Edge &e) {
        std::set< Edge >::iterator it = _hash2.find(e);
        if (it != _hash2.end())return false;
        _hash2.insert(e);
        return true;
      }
      bool addNewEdge (const Edge &e)
      {
        size_t h = ((size_t) e.first >> 3) ;
        std::vector<Edge> &v = _hash[h %_hash.size()];
        AVGSEARCH+=v.size();
        for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
        v.push_back(e);
        return true;
      }
    };
    
    struct IPT {
      double _x1,_x2,_x3,_x4;
      IPT(double x1, double x2, double x3, double x4) :
        _x1(x1),_x2(x2),_x3(x3),_x4(x4){};
    };
    
    double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
    				double lc1 , double lc2 ,
    				double (*f)(const SPoint3 &p, void *),
    				void *data, std::vector< IPT > & _result,
    				double &dl, std::stack<IPT> &_stack, double epsilon = 1.e-5)
    {
      //  _stack.clear();
      _result.clear();
      // local parameters on the edge
      double t1 = 0.0;
      double t2 = 1.0;
      // edge vector
      SPoint3 dp = p2-p1;
    
      // value of f on both sides
      double f1 = lc1; //f(p1 + dp*t1,data);
      double f2 = lc2; //f(p1 + dp*t2,data);
    
      dl = p1.distance(p2);
    
      //  printf ("edge length %g lc %g %g\n",dl,f1,f2);
    
      // add one subsegment on the stack
      IPT top (t1,t2,f1,f2);
      _stack.push(top);
      // store total value of the integral
      double totalValue = 0.0;
      while(!_stack.empty()){
        IPT pp = _stack.top();
        _stack.pop();
        t1 = pp._x1;
        t2 = pp._x2;
        f1 = pp._x3;
        f2 = pp._x4;
        // mid point
        double t12 = 0.5* (t1+t2);
        SPoint3 pmid = p1 + dp*t12;
        double dt = t2-t1;
        // average should be compared to mid value
        double f12 = 0.5* (f1+f2);
        if (fabs (f12 - 0.5*(f1+f2)) > epsilon*dt ) {
          IPT left  (t1,t12,f1,f12);
          IPT right (t12,t2,f12,f2);
          _stack.push(left);
          _stack.push(right);
        }
        else {
          _result.push_back (pp);
          // compute the integral using trapezoidal rule on both sides
          totalValue += 1./((0.5*f12+0.25*(f1+f2)))*dt;
        }
      }
      // take into account the real length of the edge
      totalValue *= dl;
      //  printf("adimensional length %g\n",totalValue);
      return totalValue;
    }
    
    
    void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &p, void *), void *data, std::stack<IPT> &temp) {
      std::vector< IPT > _result;
      double dl;
      SPoint3 p1 = e.first->point();
      SPoint3 p2 = e.second->point();
      const double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl, temp);
      const int N = (int) (dN+0.1);
      const double interval = dN/N;
      double L = 0.0;
    
      //  printf("edge length %g %d intervals of size %g (%d results)\n",dl,N,interval,_result.size());
      const unsigned int Nr = _result.size();
      for (unsigned int i=0; i< Nr ; i++) {
        const IPT & rr = _result[i];
        const double t1 = rr._x1;
        const double t2 = rr._x2;
        const double f1 = rr._x3;
        const double f2 = rr._x4;
        const double dL = 2.*(t2-t1) * dl / (f1+f2);
    
        //    printf("%g --> %g for %g --> %g\n",L,dL,t1,t2);
        double L0 = L;
        while (1) {
          const double t = t1 + (L+interval-L0)*(t2-t1) / dL;
          if (t >= t2*.999) {
    	break;
          }
          else {
    	//	printf("%g ",t);
    	SPoint3 p = p1 * (1.-t) + p2*t;
    	double lc = e.first->lc() * (1.-t) + e.second->lc()*t;
    	const double dx = 0;//1.e-12 * (double) rand() / RAND_MAX;
    	const double dy = 0;//1.e-12 * (double) rand() / RAND_MAX;
    	const double dz = 0;//1.e-12 * (double) rand() / RAND_MAX;
    	S.push_back(new Vertex(p.x()+dx,p.y()+dy,p.z()+dz,lc));
    	L += interval;
          }
        }
      }
      //  printf(" press enter\n");
      //  getchar();
      //  printf("%d points added\n",S.size());
    
      //  exit(1);
    }
    
    void saturateEdges ( edgeContainer &ec,
    		     tetContainer &T,
    		     int nbThreads,
    		     std::vector<Vertex*> &S,
    		     double (*f)(const SPoint3 &p, void *), void *data) {
      std::stack<IPT> temp;
      AVGSEARCH= 0;
      // FIXME
      const int N = T.size(0);
      for (int i=0;i<N;i++){
        Tet *t = T(0,i);
       if (t->V[0] && t->_modified){
          t->_modified = false;
          for (int iEdge=0;iEdge<6;iEdge++){
    	Edge ed = t->getEdge(iEdge);
    	bool isNew = ec.addNewEdge(ed);
    	if (isNew){
    	  saturateEdge (ed, S, f, data, temp);
    	}
          }
        }
      }
    }
    
    /////////////////////////   F I L T E R I N G ////////////////////////////////////////////////////
    
    #define SQR(X) (X)*(X)
    
    class volumePointWithExclusionRegion {
    public :
      Vertex *_v;
      volumePointWithExclusionRegion (Vertex *v) : _v(v){
      }
    
      inline bool inExclusionZone (volumePointWithExclusionRegion *p) const
      {
        const double FACTOR = 0.8;
        const double K = FACTOR*p->_v->lc();
        const double d =
          SQR(p->_v->x() - _v->x())+
          SQR(p->_v->y() - _v->y())+
          SQR(p->_v->z() - _v->z());
        //    printf(" %g %g\n",p-//>_v->lc(),d);
        return d < SQR(K);
      }
      void minmax (double _min[3], double _max[3]) const
      {
        _min[0] = _v->x() - _v->lc();
        _min[1] = _v->y() - _v->lc();
        _min[2] = _v->z() - _v->lc();
        _max[0] = _v->x() + _v->lc();
        _max[1] = _v->y() + _v->lc();
        _max[2] = _v->z() + _v->lc();
      }
    };
    
    struct my_wrapper_3D {
      bool _tooclose;
      volumePointWithExclusionRegion *_p;
      my_wrapper_3D (volumePointWithExclusionRegion *sp) :
        _tooclose (false), _p(sp) {}
    };
    
    
    bool rtree_callback(volumePointWithExclusionRegion *neighbour,void* point)
    {
      my_wrapper_3D *w = static_cast<my_wrapper_3D*>(point);
    
      if (neighbour->inExclusionZone(w->_p)){
        w->_tooclose = true;
        return false;
      }
      return true;
    }
    
    
    
    class vertexFilter {
      RTree<volumePointWithExclusionRegion*,double,3,double> _rtree;
    public:
      void insert (Vertex * v) {
        volumePointWithExclusionRegion *sp = new volumePointWithExclusionRegion (v);
        double _min[3],_max[3];
        sp->minmax(_min,_max);
        _rtree.Insert (_min,_max,sp);
      }
    
      bool inExclusionZone  (volumePointWithExclusionRegion *p)
      {
        my_wrapper_3D w (p);
        double _min[3] = {p->_v->x()-1.e-8, p->_v->y()-1.e-8, p->_v->z()-1.e-8};
        //    double _max[3] = {p->_v->x()+1.e-8, p->_v->y()+1.e-8, p->_v->z()+1.e-8};
        _rtree.Search(_min,_min,rtree_callback,&w);
        return w._tooclose;
      }
    };
    
    void filterVertices (const int numThreads,
    		     vertexFilter &_filter,
    		     std::vector<Vertex*> &add,
    		     double (*f)(const SPoint3 &p, void *),
    		     void *data) {
    
      std::vector<int> indices;
      SortHilbert(add, indices);
      std::vector<Vertex*> _add=add;
    
      // std::vector<Vertex*> _add;
      // Vertex *current = add[0];
      // printf("before %d\n",add.size());
      // for (unsigned int i=1;i<add.size();i++){
      //   const double d = sqrt (SQR(add[i]->x()-current->x())  +
      // 			   SQR(add[i]->y()-current->y())  +
      // 			   SQR(add[i]->z()-current->z())  );
      //   if (0.8*current->lc() > d){
      //     delete add[i];
      //   }
      //   else {
      //     current = add[i];
      //     _add.push_back(add[i]);
      //   }
      // }
      //  printf("after %d\n",_add.size());
    
      add.clear();
      for (unsigned int i=0;i<_add.size();i++){
        SPoint3 p (_add[i]->x(),_add[i]->y(),_add[i]->z());
        volumePointWithExclusionRegion v (_add[i]);
        if (! _filter. inExclusionZone (&v)){
          _filter.insert( _add[i]);
          add.push_back(_add[i]);
        }
        else
          delete _add[i];
      }
    }
    
    
    double _fx (const SPoint3 &p, void *){
      return fabs(0.0125 + .02*p.x());
    }
    
    
    static void _print (const char *name, std::vector<Vertex*> &T){
      FILE *f = fopen(name,"w");
      fprintf(f,"View \"\"{\n");
      for (unsigned int i=0;i<T.size();i++){
        fprintf(f,"SP(%g,%g,%g){%d};\n",
    	    T[i]->x(),T[i]->y(),T[i]->z(),i);
      }
      fprintf(f,"};\n");
      fclose(f);
    }
    
    typedef std::set<conn>   connSet;
    
    void computeAdjacencies (Tet *t, int iFace, connSet &faceToTet){
      conn c (t->getFace(iFace), iFace, t);
      connSet::iterator it = faceToTet.find(c);
      if (it == faceToTet.end()){
        faceToTet.insert(c);
      }
      else{
        t->T[iFace] = it->t;
        it->t->T[it->i] =t;
        faceToTet.erase(it);
      }
    }
    
    bool edgeSwaps(tetContainer &T, int myThread)
    {
      // TODO
      return false;
    }
    
    
    void edgeBasedRefinement (const int numThreads,
    			  const int nptsatonce,
    			  GRegion *gr) {
    
      // fill up old Datastructures
    
      tetContainer allocator (numThreads,1000000);
    
      SBoundingBox3d bb;
      std::vector<Vertex *> _vertices;
      edgeContainer ec;
      std::map<Vertex*,MVertex*> _ma;
    
      {
        std::vector<MTetrahedron*> &T = gr->tetrahedra;
        std::set<MVertex *> all;
        for (unsigned int i=0;i<T.size();i++){
          for (unsigned int j=0;j<4;j++){
    	all.insert(T[i]->getVertex(j));
          }
        }
    
    
        //    FILE *f = fopen ("pts_init.dat","w");
        //    fprintf(f,"%d\n",all.size());
        //    for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
        //      MVertex *mv = *it;
        //      fprintf(f,"%12.5E %12.5E %12.5E\n",mv->x(),mv->y(),mv->z());
        //    }
        //    fclose(f);
    
    
        _vertices.resize(all.size());
        int counter=0;
        for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
          MVertex *mv = *it;
          mv->setIndex(counter);
          Vertex *v = new Vertex (mv->x(),mv->y(),mv->z(),1.e22, counter);
          _vertices[counter] = v;
          bb += SPoint3(v->x(),v->y(),v->z());
          _ma[v] = mv;
          counter++;
        }
        bb *= 1.1;
        {
          connSet faceToTet;
          // FIXME MULTITHREADING
          for (unsigned int i=0;i<T.size();i++){
    	MTetrahedron  *tt = T[i];
    	int i0 = tt->getVertex(0)->getIndex();
    	int i1 = tt->getVertex(1)->getIndex();
    	int i2 = tt->getVertex(2)->getIndex();
    	int i3 = tt->getVertex(3)->getIndex();
    	Tet *t = allocator.newTet(0) ; t->setVertices (_vertices[i0],_vertices[i1],_vertices[i2],_vertices[i3]);
    	computeAdjacencies (t,0,faceToTet);
    	computeAdjacencies (t,1,faceToTet);
    	computeAdjacencies (t,2,faceToTet);
    	computeAdjacencies (t,3,faceToTet);
    	delete tt;
          }
          T.clear();
        }
      }
    
      // do not allow to saturate boundary edges
      {
        for (unsigned int i=0;i< allocator.size(0);i++) {
          Tet  *tt = allocator (0,i);
          for (int j=0;j<4;j++){
    	if (!tt->T[j]){
    	  Face f = tt->getFace(j);
    	  for (int k=0;k<3;k++){
    	    Vertex *vi = f.V[k];
    	    Vertex *vj = f.V[(k+1)%3];
    	    double l = sqrt ((vi->x()-vj->x())*(vi->x()-vj->x())+
    			     (vi->y()-vj->y())*(vi->y()-vj->y())+
    			     (vi->z()-vj->z())*(vi->z()-vj->z()));
    	    ec.addNewEdge(Edge(vi,vj));
    	    vi->lc() = std::min (l,vi->lc() );
    	    vj->lc() = std::min (l,vj->lc() );
    	  }
    	}
          }
        }
        for (unsigned int i=0;i< allocator.size(0);i++) {
          Tet  *tt = allocator (0,i);
          for (int j=0;j<6;j++){
    	Edge e = tt->getEdge(j);
    	if(e.first->lc() == 1.e22){/*printf("coucou\n");*/e.first->lc() = e.second->lc();}
    	else if(e.second->lc() == 1.e22){/*printf("coucou\n");*/e.second->lc() = e.first->lc();}
          }
        }
      }
    
      std::vector<Vertex*> add_all;
      {
        //   vertexFilter _filter (bb, 20);
        vertexFilter _filter;
        for (unsigned int i=0;i<_vertices.size();i++){
          _filter.insert( _vertices[i] );
        }
    
        int iter = 1;
    
        Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME  TETS");
    
        double __t__ = Cpu();
        //    Tet::in_sphere_counter = 0;
        while(1){
          std::vector<Vertex*> add;
          double t1 = Cpu();
          saturateEdges (ec, allocator, numThreads, add, _fx, NULL);
          double t2 = Cpu();
          filterVertices (numThreads, _filter, add, _fx, NULL);
          double t3 = Cpu();
          if (add.empty())break;
          // randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
          std::random_shuffle(add.begin(), add.end());
          // sort them using BRIO
          std::vector<int> indices;
          SortHilbert(add, indices);
          double t4 = Cpu();
          delaunayTrgl (1,1,add.size(), &add,allocator,1.e-28);
          double t5 = Cpu();
          add_all.insert (add_all.end(), add.begin(), add.end());
          Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(),
    		(t2-t1),
    		(t3-t2),
    		(t4-t3),
    		(t5-t4),
    		(t5-__t__),
    		allocator.size(0));
          iter++;
        }
      }
    
    
    
      for (unsigned int i=0; i< allocator.size(0);i++){
        Tet  *tt = allocator (0,i);
        MVertex *mvs[4];
        if (tt->V[0]){
          for (int j=0;j<4;j++){
    	Vertex *v = tt->V[j];
    	std::map<Vertex*,MVertex*>::iterator it = _ma.find(v);
    	if (it == _ma.end()){
    	  MVertex *mv = new MVertex (v->x(),v->y(),v->z(),gr);
    	  gr->mesh_vertices.push_back(mv);
    	  _ma[v] = mv;
    	  mvs[j] = mv;
    	}
    	else mvs[j] = it->second;
          }
          gr->tetrahedra.push_back(new MTetrahedron(mvs[0],mvs[1],mvs[2],mvs[3]));
        }
      }
    
      if (Msg::GetVerbosity() == 99) {
        std::map<Edge,double> _sizes;
        for (unsigned int i=0; i< allocator.size(0);i++){
          Tet  *tt = allocator (0,i);
          if (tt->V[0]){
    	for (int j=0;j<6;j++){
    	  Edge e =  tt->getEdge(j);
    	  std::map<Edge,double>::iterator it = _sizes.find(e);
    	  if (it == _sizes.end()){
    	    double l = sqrt ((e.first->x() -  e.second->x()) * (e.first->x() -  e.second->x()) +
    			     (e.first->y() -  e.second->y()) * (e.first->y() -  e.second->y()) +
    			     (e.first->z() -  e.second->z()) * (e.first->z() -  e.second->z()));
    	    _sizes[e]= 2* l / (e.first->lc() + e.second->lc());
    	  }
    	}
          }
        }
        std::map<Edge,double>::iterator it = _sizes.begin();
        double sum = 0;
        int nbBad = 0;
        for (; it !=_sizes.end();++it){
          double d = it->second;
          double tau = d < 1 ? d - 1 : 1./d - 1;
          if (d > 2.)nbBad++;
          sum += tau;
        }
        Msg::Info("MESH EFFICIENCY : %22.15E %6d edges among %d are out of range",exp (sum / _sizes.size()),nbBad,_sizes.size());
      }
    }