← Back to branch summary

~ubuntu-branches/ubuntu/wily/freefem++/wily

« back to all changes in this revision

Viewing changes to examples++-load/msh3.cpp

  • Committer: Package Import Robot
  • Author(s): Dimitrios Eftaxiopoulos, Dimitrios Eftaxiopoulos, Christophe Trophime
  • Date: 2013-09-12 00:02:58 UTC
  • mfrom: (1.2.1) (11.1.1 experimental)
  • Revision ID: package-import@ubuntu.com-20130912000258-aclq2zfa1svt0p3x
Tags: 3.25-1
http://bugs.debian.org/701161
[ Dimitrios Eftaxiopoulos ]
* Imported Upstream version 3.25 (Closes: #701161 #706714)
* Change installation directory of header-like *.idp files
  from /usr/lib/freefem++ to /usr/include/freefem++, in order
  to fix a lintian warning
* Update patch to examples++-load/Makefile.am in order to enable
  functioning of load *.so and include *.idp commands in *.edp
  scripts
* Delete patches to src/Graphics/sansgraph.cpp and
  src/Graphics/xglrgraph.cpp because they are not needed any more
* Fix lintian warning about missing LDFLAGS
* Override dh_auto_test in debian/rules, such that in case it is 
  used, it completes executing all *.edp example files, regardless
  of aborting on some of them
* Add libmetis-dev to build-deps in d/control
* Remove libparmetis-dev from build deps
* Add --parallel option to dh $@ in debian/rules
* Add hardening compilation flags to mpic++
* Allow testing of compiling and running the example files after build

[ Christophe Trophime ]
* update C. Trophime email
* add support for nlopt, ipopt - simplify debian/rules
* upload CT changes to 3.20
* add patch for configure
* add patch for examples++-mpi
* fix bamg install
* add corrected scripts to build plugins
* add patch for properly build examples++-load
* add lintian overrides for libfreefem++
* add some missing files
* update patches
* update rules
* reorder BuildDepends - comment out unsupported libs

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples++-load/msh3.cpp
610
610
    return sqrt(pow(x,2)+pow(y,2));
611
611
    break;
612
612
  default :
613
 
    cout << "zmin_func pas d�finis" << endl;
 
613
    cout << "zmin_func no defined" << endl;
614
614
    return 0.;
615
615
  }
616
616
}
628
628
    return 3.+sqrt(pow(x,2)+pow(y,2));
629
629
    break;
630
630
  default :
631
 
    cout << "zmaxfunc pas d�finis" << endl;
 
631
    cout << "zmaxfunc no defined" << endl;
632
632
    return 0.;
633
633
  }
634
634
}
660
660
    return int(multi*(3+sqrt(pow(x,2)+pow(y,2))));
661
661
    break;
662
662
  default :
663
 
    cout << "Ni_func pas d�finis" << endl;
 
663
    cout << "Ni_func no defined" << endl;
664
664
    return 0;
665
665
  }
666
666
}
1677
1677
  R3 Pn(1e100,1e100,1e100),Px(-1e100,-1e100,-1e100);
1678
1678
  const list<Mesh3 *> lth(*lst.lth);
1679
1679
  Mesh3 * th0=0;
1680
 
  
 
1680
  int kk=0; 
1681
1681
  for(list<Mesh3 *>::const_iterator i=lth.begin();i != lth.end();i++)
1682
1682
    {
1683
 
 
 
1683
      if( ! *i) continue ;
 
1684
      kk++;
1684
1685
      Mesh3 &Th3(**i);  // definis ???
1685
1686
      th0=&Th3;
1686
1687
      if(verbosity>1)  cout << " determination of hmin : GluMesh3D + "<< Th3.nv << " " << Th3.nt << " "<< Th3.nbe << endl;
1707
1708
        Px=Maxc(P,Px);     
1708
1709
      }
1709
1710
    } 
1710
 
  
 
1711
  if(kk==0) return 0; // no mesh ....
 
1712
 
1711
1713
  if(verbosity > 1) cout << "      - hmin =" <<  hmin << " ,  Bounding Box: " << Pn << " "<< Px << endl;
1712
1714
  
1713
1715
  // probleme memoire
1731
1733
  //int nbv0=0;
1732
1734
  for(list<Mesh3 *>::const_iterator i=lth.begin(); i!=lth.end();i++)
1733
1735
    {
 
1736
      if( ! *i) continue ;
1734
1737
      const Mesh3 &Th3(**i);
1735
1738
      if(verbosity>1)  cout << " loop over mesh for create new mesh "<< endl;
1736
1739
      if(verbosity>1)  cout << " GluMesh3D + "<< Th3.nv << " " << Th3.nt <<" " << Th3.nbe << endl;
1777
1780
  
1778
1781
  double hseuil_border = hseuil/3.;
1779
1782
  //nbv0=0;
1780
 
  for(list<Mesh3 *>::const_iterator i=lth.begin();i != lth.end();i++){
1781
 
    const Mesh3 &Th3(**i);
 
1783
  for(list<Mesh3 *>::const_iterator i=lth.begin();i != lth.end();i++)
 
1784
    {
 
1785
      if( ! *i) continue ;
 
1786
      const Mesh3 &Th3(**i);
1782
1787
      
1783
1788
    for (int k=0;k<Th3.nbe;k++)
1784
1789
      {
1860
1865
template<class RR,class AA=RR,class BB=AA> 
1861
1866
struct Op3_addmesh: public binary_function<AA,BB,RR> { 
1862
1867
  static RR f(Stack s,const AA & a,const BB & b)  
1863
 
  { cout << "Op3_addmesh" << endl; return RR(s, a, b );} 
 
1868
  { return RR(s, a, b );} 
1864
1869
};
1865
1870
 
1866
1871
template<bool INIT,class RR,class AA=RR,class BB=AA> 
1870
1875
    ffassert(a );
1871
1876
    pmesh3  p=GluMesh3(b);
1872
1877
    
1873
 
    cout << "INIT=" << INIT << endl;
1874
1878
    if(!INIT && *a){
1875
1879
      //Add2StackOfPtr2FreeRC(stack,*a);
1876
1880
       (**a).destroy();
2193
2197
public:
2194
2198
  Expression a; 
2195
2199
  
2196
 
  static const int n_name_param =2+2+2; // 
 
2200
  static const int n_name_param =2+2+2+2; //
2197
2201
  static basicAC_F0::name_and_type name_param[] ;
2198
2202
  Expression nargs[n_name_param];
2199
2203
  KN_<long>  arg(int i,Stack stack,KN_<long> a ) const
2200
2204
    { ffassert( !(nargs[i] && nargs[i+2]));
2201
2205
      i = nargs[i] ? i : i+2;
2202
2206
      return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
 
2207
    long  arg(int i,Stack stack, long  a ) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
 
2208
    bool   arg(int i,Stack stack, bool   a ) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
2203
2209
 
2204
2210
  
2205
2211
public:
2220
2226
  {  "refface", &typeid(KN_<long> )},
2221
2227
  {  "region", &typeid(KN_<long> )},
2222
2228
  {  "label", &typeid(KN_<long> )},
2223
 
    {  "fregion", &typeid(long)},
2224
 
    {  "flabel", &typeid(long )}
2225
 
   
 
2229
  {  "fregion", &typeid(long)},
 
2230
  {  "flabel", &typeid(long )},
 
2231
  {  "rmlfaces", &typeid(long)},
 
2232
  {  "rmInternalFaces", &typeid(bool)}
 
2233
 
2226
2234
};
2227
2235
//  besoin en cas de fichier 2D / fichier 3D 
2228
2236
 
2239
2247
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
2240
2248
  Mesh3 * pTh= GetAny<Mesh3 *>((*a)(stack));
2241
2249
  Mesh3 & Th=*pTh;
 
2250
  if(!pTh) return pTh;
2242
2251
  Mesh3 *m= pTh;
2243
2252
  int nbv=Th.nv; // nombre de sommet 
2244
2253
  int nbt=Th.nt; // nombre de triangles
2249
2258
  KN<long> nrface (arg(1,stack,zz));  
2250
2259
  Expression freg = nargs[4];
2251
2260
  Expression flab = nargs[5];
 
2261
    bool  rm_faces = nargs[6];
 
2262
    long  rmlabfaces (arg(6,stack,0L));
 
2263
    bool  rm_i_faces (arg(7,stack,false));
 
2264
 
2252
2265
 // cout << " Chnage " << freg << " " << flab << endl;   
2253
 
  if(nrface.N() <=0 && nrtet.N() <=0 && (!freg) && (!flab)  ) return m; // modf J.M. oct 2010 
 
2266
  if(nrface.N() <=0 && nrtet.N() <=0 && (!freg) && (!flab) && !rmlabfaces && !rm_i_faces ) return m; // modf J.M. oct 2010
2254
2267
  ffassert( nrtet.N() %2 ==0);
2255
2268
  ffassert( nrface.N() %2 ==0);
2256
2269
  
2331
2344
  
2332
2345
  Triangle3 * bb=b;
2333
2346
     R2 PtHat2(1./3.,1./3.);
 
2347
    int nrmf=0; 
2334
2348
  for (int i=0;i<nbe;i++)
2335
2349
    { 
2336
2350
      const Triangle3 &K( Th.be(i) );
2337
2351
      int fk,ke = Th.BoundaryElement(i,fk); // element co
 
2352
      int fkk,kke = Th.ElementAdj(ke,fkk=fk); // element co
 
2353
      bool onborder = (kke==ke) || (kke <0) ;
2338
2354
      const Tet & KE(Th[ke]);
2339
2355
      R3 B= KE.PBord(fk,PtHat2);
2340
2356
      int iv[3];       
2341
 
    
 
2357
      bool  rmf = rm_i_faces && ! onborder; 
2342
2358
      iv[0] = Th.operator()(K[0]);
2343
2359
      iv[1] = Th.operator()(K[1]);
2344
2360
      iv[2] = Th.operator()(K[2]);
2345
2361
      
2346
2362
      int l0,l1=ChangeLab3D(mapface,l0=K.lab) ;
2347
 
      (*bb).set( v, iv, l1);   
2348
 
        if(flab)
 
2363
      if(flab)
2349
2364
          {//      R3 B(1./4.,1./4.,1./4.);  // 27/09/10 : J.Morice error in msh3.cpp
2350
2365
              mp->set(Th,KE(B),B,KE,K.lab);
2351
 
              bb->lab =GetAny<long>( (* freg)(stack)) ;  
 
2366
              l1 =GetAny<long>( (* flab)(stack)) ;
2352
2367
              lmn= min (lmn,bb->lab);
2353
2368
              lmx= max (lmx,bb->lab);
2354
2369
          }
 
2370
        if( !rmf && rm_faces)
 
2371
          rmf =  !onborder &&  ( l1 == rmlabfaces  );
 
2372
        if(rmf)
 
2373
            nrmf++;
 
2374
        else
 
2375
         (*bb++).set( v, iv, l1);
 
2376
 
2355
2377
        
2356
 
        bb++;
2357
2378
    }
 
2379
    if(nrmf && verbosity > 2) cout << "   change  mesh3 : number of removed  internal faces " << nrmf << " == " << nbe - (bb-b) << endl;
2358
2380
    
 
2381
  nben -= nrmf;
 
2382
  nbe -= nrmf;
2359
2383
  assert(nben==bb-b);
2360
2384
  *mp=mps; 
2361
2385
  if(nbt != 0)
2381
2405
      
2382
2406
      return mpq;
2383
2407
    }
 
2408
 
 
2409
  Mesh3 *mpq = NULL;
 
2410
  return mpq;
2384
2411
}
2385
2412
 
2386
2413
 
4007
4034
 
4008
4035
  //hseuil = hseuil/10.;
4009
4036
        
4010
 
  //Vertex3  *v= new Vertex3[tab_nv];
4011
 
  Vertex3  v[tab_nv];
 
4037
  Vertex3  *v= new Vertex3[tab_nv];
 
4038
  //Vertex3  v[tab_nv];
4012
4039
  
4013
4040
  EF23::GTree<Vertex3> *gtree = new EF23::GTree<Vertex3>(v,bmin,bmax,0);
4014
4041
        
4049
4076
  }
4050
4077
 
4051
4078
  delete gtree;
4052
 
 
 
4079
  delete [] v;
4053
4080
        
4054
4081
  if(verbosity >1) cout << "hseuil=" << hseuil <<endl;
4055
4082
  if(verbosity >1) cout << "nv_t = " << nv_t << " / " << "nv_t(anc)" << tab_nv <<endl;   
4081
4108
                              const R3 & bmin, const R3 & bmax, const double &hmin, int * ind_np, int * ind_label, int & np){
4082
4109
        
4083
4110
  double hseuil =hmin/10.; 
4084
 
  //Vertex3  *v= new Vertex3[NbPoints];
4085
 
  Vertex3 v[NbPoints];
 
4111
  Vertex3  *v= new Vertex3[NbPoints];
 
4112
  //Vertex3 v[NbPoints];
4086
4113
  EF23::GTree<Vertex3> *gtree = new EF23::GTree<Vertex3>(v,bmin,bmax,0);
4087
4114
        
4088
4115
  if(verbosity >1) cout<< "verif hmin vertex3 GTree switch" << point_confondus_ok << endl;
4167
4194
  }
4168
4195
        
4169
4196
  delete gtree;
 
4197
  delete [] v; 
4170
4198
 
4171
4199
  /*
4172
4200
    int z_verifnumberofpoints;
4891
4919
 };
4892
4920
 
4893
4921
 
4894
 
Mesh3 * truncmesh(const Mesh3 &Th,const long &kksplit,int *split, bool kk, const int newbelabel){
4895
 
  
4896
 
  static const int FaceTriangle[4]={3,0,1,2};  //={{3,2,1}, {0,2,3},{ 3,1,0},{ 0,1,2}}
4897
 
 
4898
 
 
4899
 
  // computation of number of border elements and vertex without split
4900
 
  int nbe = 0;
4901
 
  int nt  = 0;
4902
 
  int nv  = 0;
4903
 
  int nvtrunc =0;
4904
 
  int nbedge=0;
4905
 
  double hmin=1e100;
4906
 
  R3 bmin,bmax;
4907
 
    int nbeee=0;
4908
 
  const int kksplit2 = kksplit*kksplit; 
4909
 
  const int kksplit3 = kksplit2*kksplit; 
4910
 
  int tagb[4]={1,2,4,8} ; 
4911
 
  KN<int> tagTonB(Th.nt);
4912
 
  tagTonB=0; 
4913
 
    
4914
 
  for( int ibe=0; ibe < Th.nbe; ibe++)
4915
 
    {  
 
4922
Mesh3 * truncmesh(const Mesh3 &Th,const long &kksplit,int *split, bool kk, const int newbelabel)
 
4923
{
 
4924
    
 
4925
    static const int FaceTriangle[4]={3,0,1,2};  //={{3,2,1}, {0,2,3},{ 3,1,0},{ 0,1,2}}
 
4926
    
 
4927
    
 
4928
    // computation of number of border elements and vertex without split
 
4929
    int nbe = 0;
 
4930
    int nt  = 0;
 
4931
    int nv  = 0;
 
4932
    int nvtrunc =0;
 
4933
    int nedge=0;
 
4934
    int nface=0;
 
4935
    double hmin=1e100;
 
4936
    R3 bmin,bmax;
 
4937
    int nbeee=0,nbfi=0;
 
4938
    const int kksplit2 = kksplit*kksplit;
 
4939
    const int kksplit3 = kksplit2*kksplit;
 
4940
    int ntsplit =0;
 
4941
    int tagb[4]={1,2,4,8} ;
 
4942
    KN<int> tagTonB(Th.nt);
 
4943
    tagTonB=0;
 
4944
    
 
4945
    for( int ibe=0; ibe < Th.nbe; ibe++)
 
4946
    {
4916
4947
        int iff;
4917
4948
        int it=Th.BoundaryElement(ibe,iff);
4918
4949
        tagTonB[it]|= tagb[iff];
4919
 
         int ifff=iff,itt=Th.ElementAdj(it,ifff);
 
4950
        int ifff=iff,itt=Th.ElementAdj(it,ifff);
4920
4951
        if(itt >=0 &&  itt != it)
4921
 
             tagTonB[itt]|= tagb[ifff];        
4922
 
    }
4923
 
    
4924
 
  for (int i=0;i<Th.nt;i++)
4925
 
    if(split[i]) 
4926
 
      { 
4927
 
        // computation of number of tetrahedrons 
4928
 
        nt=nt+kksplit3;  
4929
 
        // computation of number of border elements
4930
 
        for (int j=0;j<4;j++)
4931
 
          {
4932
 
            int jt=j,it=Th.ElementAdj(i,jt);
4933
 
            if ( (it==i || it <0)  || ! split[it]) nbeee++;   
4934
 
            if(it==i || it <0) nbe += kksplit2;  //on est sur la frontiere
4935
 
            else if (!split[it]) nbe += kksplit2; //le voisin ne doit pas etre decoupe
4936
 
            else if ( (tagTonB[i]&tagb[j] ) != 0 && i<it) nbe += kksplit2; // internal boundary ..
4937
 
            //else{
4938
 
            // rien a faire
4939
 
            //}
4940
 
              
4941
 
          }
4942
 
 
4943
 
        for (int e=0;e<6;e++){
4944
 
          hmin=min(hmin,Th[i].lenEdge(e));   // calcul de .lenEdge pour un Mesh3
4945
 
        }
4946
 
      }
4947
 
  if(verbosity>5) 
4948
 
  cout << "  number of  not intern boundary faces = " << nbeee << ",  all faces  =  " << nbe << endl;
4949
 
  double hseuil = (hmin/kksplit)/10.; 
4950
 
 
4951
 
  /* determination de bmin, bmax et hmin */ 
4952
 
 
4953
 
  KN<int> takevertex(Th.nv);
4954
 
  for(int i=0; i<Th.nv; i++){
4955
 
    takevertex[i]=-1;
4956
 
  } 
4957
 
  for(int i=0; i<Th.nt; i++){
4958
 
    if(split[i]){
4959
 
      const Tet &K(Th.elements[i]);
4960
 
          
4961
 
      for(int ii=0; ii<4; ii++){
4962
 
        int iv= Th.operator()( K[ii] );
4963
 
        if( takevertex[iv] == -1 ){
4964
 
          R3 P( Th.vertices[iv].x, Th.vertices[iv].y, Th.vertices[iv].z);
4965
 
          bmin=Minc(P,bmin);
4966
 
          bmax=Maxc(P,bmax); 
4967
 
          takevertex[iv]=nvtrunc;
4968
 
          nvtrunc++;
4969
 
        }
4970
 
      }
4971
 
      
4972
 
    }
4973
 
  }
4974
 
  
4975
 
  if( kksplit > 1 ){
4976
 
    KN<int> NbFois(nvtrunc,0);
4977
 
    for(int i=0; i<nvtrunc; i++)
4978
 
      NbFois[i] = 0;
4979
 
    int majedge=0;
4980
 
    for(int i=0; i<Th.nt; i++){
4981
 
      if(split[i]){
4982
 
        const Tet &K(Th.elements[i]);
4983
 
        for(int e=0; e<6; e++){
4984
 
          int e1 = Th.operator()( K[ Th[i].nvedge[e][0] ] );
4985
 
          int e2 = Th.operator()( K[ Th[i].nvedge[e][1] ] );
4986
 
          
4987
 
          if( takevertex[e1] < takevertex[e2] ){
4988
 
            NbFois[ takevertex[e1] ]++;
4989
 
            majedge++;
4990
 
          }
4991
 
          else{
4992
 
            NbFois[ takevertex[e2] ]++;
4993
 
            majedge++;
4994
 
          }
4995
 
        }
4996
 
 
4997
 
      }
4998
 
    }
4999
 
 
5000
 
    KN<int> first(nvtrunc+1,0);
5001
 
    KN<int> current(nvtrunc,0);
5002
 
    for(int i=1; i<nvtrunc; i++)
5003
 
      first[i] = first[i-1] + NbFois[i-1];
5004
 
    
5005
 
    first[nvtrunc] = majedge;
5006
 
    for(int i=0; i<nvtrunc; i++)
5007
 
      current[i] = first[i];
5008
 
    if(verbosity>9)
5009
 
      cout << "    trunc (3d) majedge =" << majedge << endl;
5010
 
 
5011
 
    KN<int> tableau(majedge);
5012
 
    for(int i=0; i<Th.nt; i++){
5013
 
      if(split[i]){
5014
 
        const Tet &K(Th.elements[i]);
5015
 
        for(int e=0; e<6; e++){
5016
 
          int e1 = Th.operator()( K[ Th[i].nvedge[e][0] ] );
5017
 
          int e2 = Th.operator()( K[ Th[i].nvedge[e][1] ] );
5018
 
          //int e1 = Th[i].nvedge[e][0];
5019
 
          //int e2 = Th[i].nvedge[e][1];
5020
 
          
5021
 
          if( takevertex[e1] < takevertex[e2] ){
5022
 
            tableau[ current[ takevertex[e1] ] ] = takevertex[e2]; 
5023
 
            current[ takevertex[e1] ]++;
5024
 
          }
5025
 
          else{
5026
 
            tableau[ current[ takevertex[e2] ] ] = takevertex[e1];
5027
 
            current[ takevertex[e2] ]++;
5028
 
          }
5029
 
          
5030
 
        }
5031
 
      }
5032
 
    }
5033
 
    for(int i=0; i<nvtrunc; i++)
5034
 
      assert(current[i] == first[i+1]);
5035
 
 
5036
 
    // determination du nombre d'edge
5037
 
    
5038
 
    for(int i=0; i<nvtrunc; i++){
5039
 
      list<int> list1;
5040
 
      list<int>::const_iterator ilist;
5041
 
      //cout << "i = "<< i  << " nvtrunc = " <<  nvtrunc << endl;
5042
 
      for(int jj=first[i]; jj < first[i+1]; jj++){
5043
 
        int labOk = 0;
5044
 
        //cout << "tableauu[ jj ] = " <<  tableau[ jj ] << endl;
5045
 
        for( ilist=list1.begin(); ilist!=list1.end(); ilist++){
5046
 
          if( *ilist == tableau[ jj ] ){ labOk = 1;   break; }    
5047
 
        }
5048
 
        if(labOk == 0){
5049
 
          list1.push_back(tableau[jj]);
5050
 
          //cout << "push back :: tableauu[ jj ] = " <<  tableau[ jj ] << endl;
5051
 
          nbedge++;
5052
 
        }
5053
 
 
5054
 
      }
5055
 
      //cout << " nbedge="<< nbedge  << endl;
5056
 
    }
5057
 
  }
5058
 
 
5059
 
 
5060
 
 
5061
 
  /* determination des vertex, triangles et tetrah�dre obtenue apr�s splitting dans le Simplex */ 
5062
 
 
5063
 
  int nfacesub = kksplit2;
5064
 
  int ntetsub = kksplit3;
5065
 
  int nvsub = (kksplit+1)*(kksplit+2)*(kksplit+3)/6;
5066
 
  int ntrisub = 4*kksplit2;
5067
 
   
5068
 
  R3 *vertexsub; //[nvsub];
5069
 
  int *tetsub;   //[4*ntetsub];
5070
 
  int *trisub;   //[4*kksplit*kksplit];
5071
 
 
5072
 
  SplitSimplex<R3>( kksplit, nvsub, vertexsub, ntetsub, tetsub);
5073
 
  SplitSurfaceSimplex( kksplit, ntrisub, trisub);
5074
 
  
5075
 
  /*
5076
 
    for( int iii=0; iii< nvsub; iii++){
5077
 
    cout << "vertexsub["<< iii <<"]=" << " "<< vertexsub[iii].x << " "<< vertexsub[iii].y << " "<< vertexsub[iii].z << endl;
5078
 
    }
5079
 
    for( int iii=0; iii< ntetsub; iii++){
5080
 
    cout << "tetsub=" << tetsub[4*iii] << " "<< tetsub[4*iii+1] << " "<< tetsub[4*iii+2] << " "<< tetsub[4*iii+3] <<endl;
5081
 
    }
5082
 
 
5083
 
    for( int iii=0; iii< 4*kksplit2; iii++){
5084
 
    cout << iii << " tetsub=" << trisub[3*iii] << " "<< trisub[3*iii+1] << " " << trisub[3*iii+2] << endl;
5085
 
    }
5086
 
  */
5087
 
  if(verbosity>3) 
5088
 
    cout << "  -- trunc (3d) : Th.nv= " << Th.nv << "kksplit="<< kksplit << endl;
5089
 
  // determination de nv 
5090
 
  /*if(kksplit == 1)
5091
 
    nv=nvtrunc;
5092
 
  else{
5093
 
  */
 
4952
            tagTonB[itt]|= tagb[ifff];
 
4953
    }
 
4954
    
 
4955
    for (int i=0;i<Th.nt;i++)
 
4956
        if(split[i])
 
4957
        {
 
4958
            ++ntsplit;
 
4959
            // computation of number of tetrahedrons
 
4960
            nt=nt+kksplit3;
 
4961
            // computation of number of border elements
 
4962
            for (int j=0;j<4;j++)
 
4963
            {
 
4964
                int jt=j,it=Th.ElementAdj(i,jt);
 
4965
                if ( (it==i || it <0)  || ! split[it]) nbeee++;// boundary face ...
 
4966
                else nbfi++; // internal face count 2 times ...
 
4967
                if(it==i || it <0) nbe += kksplit2;  //on est sur la frontiere
 
4968
                else if (!split[it]) nbe += kksplit2; //le voisin ne doit pas etre decoupe
 
4969
                else if ( (tagTonB[i]&tagb[j] ) != 0 && i<it) nbe += kksplit2; // internal boundary ..
 
4970
            }
 
4971
            
 
4972
            for (int e=0;e<6;e++){
 
4973
                hmin=min(hmin,Th[i].lenEdge(e));   // calcul de .lenEdge pour un Mesh3
 
4974
            }
 
4975
        }
 
4976
    ffassert( nbfi %2 ==0) ;
 
4977
    nface = nbeee + nbfi/2;
 
4978
    double hseuil = (hmin/kksplit)/1000.;
 
4979
    if(verbosity>5)
 
4980
        cout << "  number of  not intern boundary faces = " << nbeee << ",  all faces  =  " << nbe << ", hseuil=" << hseuil <<endl;
 
4981
    
 
4982
    /* determination de bmin, bmax et hmin */
 
4983
    
 
4984
    KN<int> takevertex(Th.nv,-1);
 
4985
    
 
4986
    for(int i=0; i<Th.nt; i++){
 
4987
        if(split[i])
 
4988
        {
 
4989
            const Tet &K(Th.elements[i]);
 
4990
            for(int ii=0; ii<4; ii++)
 
4991
            {
 
4992
                int iv= Th.operator()( K[ii] );
 
4993
                if( takevertex[iv] == -1 )
 
4994
                {
 
4995
                    bmin=Minc(Th.vertices[iv],bmin);
 
4996
                    bmax=Maxc(Th.vertices[iv],bmax);
 
4997
                    takevertex[iv]=nvtrunc++;
 
4998
                }
 
4999
            }
 
5000
            
 
5001
        }
 
5002
    }
 
5003
    
 
5004
    if( kksplit > 1 )
 
5005
    { // compute the number of slip edge ...
 
5006
        nedge=0;
 
5007
        HashTable<SortArray<int,2>,int> edges(3*nface,nface);
 
5008
        for(int i=0; i<Th.nt; i++){
 
5009
            if(split[i])
 
5010
            {
 
5011
                const Tet &K(Th.elements[i]);
 
5012
                for(int e=0;e<6;++e)
 
5013
                {
 
5014
                    
 
5015
                    int e1 = Th( K[ Th[i].nvedge[e][0] ] );
 
5016
                    int e2 = Th( K[ Th[i].nvedge[e][1] ] );
 
5017
                    SortArray<int,2> key(e1,e2);
 
5018
                    if(!edges.find(key) )
 
5019
                        edges.add(key,nedge++);
 
5020
                }
 
5021
            }
 
5022
        }
 
5023
    }
 
5024
    if(verbosity>10) cout << "    -- nvertex  " << nvtrunc << ", nedges = "<< nedge
 
5025
        << ", nfaces = " << nface << " ntet =" << ntsplit
 
5026
        << endl
 
5027
        << "    -- Euler/Poincare constante = " << nvtrunc-nedge+nface-ntsplit
 
5028
        << endl;
 
5029
    
 
5030
    /* determination des vertex, triangles et tetrahedre obtenue apres splitting dans le Simplex */
 
5031
    
 
5032
    int nfacesub = kksplit2;
 
5033
    int ntetsub = kksplit3;
 
5034
    int nvsub = (kksplit+1)*(kksplit+2)*(kksplit+3)/6;
 
5035
    int ntrisub = 4*kksplit2;
 
5036
    
 
5037
    R3 *vertexsub; //[nvsub];
 
5038
    int *tetsub;   //[4*ntetsub];
 
5039
    int *trisub;   //[4*kksplit*kksplit];
 
5040
    
 
5041
    SplitSimplex<R3>( kksplit, nvsub, vertexsub, ntetsub, tetsub);
 
5042
    SplitSurfaceSimplex( kksplit, ntrisub, trisub);
 
5043
    
 
5044
    if(verbosity>3)
 
5045
        cout << "  -- trunc (3d) : Th.nv= " << Th.nv << "kksplit="<< kksplit << endl;
 
5046
    
5094
5047
    int ntnosplit  = nt/kksplit3;
5095
5048
    int nbenosplit = nbe/kksplit2;
5096
5049
    int nfacenosplit = (4*ntnosplit+nbenosplit)/2;
5098
5051
    if(verbosity>100) cout << "       1) nv= " << nv << endl;
5099
5052
    nv = nv + nfacenosplit*( (kksplit+1)*(kksplit+2)/2 - 3*(kksplit-1) -3 );
5100
5053
    if(verbosity>100) cout << "       2) nv= " << nv << endl;
5101
 
    nv = nv + nbedge*( kksplit-1 );
 
5054
    nv = nv + nedge*( kksplit-1 );
5102
5055
    if(verbosity>100) cout << "       3) nv= " << nv << endl;
5103
5056
    nv = nv + nvtrunc;
5104
5057
    if(verbosity>100) cout << "       4) nv= " << nv << endl;
5105
 
    //}
5106
 
  
5107
 
 
5108
 
  int itt=0; 
5109
 
  int ie=0;
5110
 
 
5111
 
 
5112
 
  Vertex3 *v=new Vertex3[nv];
5113
 
  Tet *t = new Tet[nt];
5114
 
  Tet *tt = t;
5115
 
  
5116
 
  Triangle3 *b  = new Triangle3[nbe];
5117
 
  Triangle3 *bb = b;
5118
 
 
5119
 
  EF23::GTree<Vertex3> *gtree = new EF23::GTree<Vertex3>(v,bmin,bmax,0);
5120
 
 
5121
 
  int np=0;
5122
 
  for(int i=0; i<Th.nt; i++){
5123
 
    if(split[i]){      
5124
 
      const Tet &K(Th.elements[i]);
5125
 
      int ivv[4]; 
5126
 
      for(int ii=0; ii< 4; ii++){
5127
 
        ivv[ii] =Th.operator()(K[ii]);
5128
 
      }
5129
 
 
5130
 
      R3 vertextetsub[nvsub];
5131
 
      for( int iv=0; iv<nvsub; iv++){
5132
 
        double alpha=vertexsub[iv].x;
5133
 
        double beta=vertexsub[iv].y;
5134
 
        double gamma=vertexsub[iv].z;
5135
 
 
5136
 
        vertextetsub[iv].x = (1-alpha-beta-gamma)*Th.vertices[ivv[0]].x + alpha*Th.vertices[ivv[1]].x + beta*Th.vertices[ivv[2]].x + gamma*Th.vertices[ivv[3]].x;  
5137
 
        vertextetsub[iv].y = (1-alpha-beta-gamma)*Th.vertices[ivv[0]].y + alpha*Th.vertices[ivv[1]].y + beta*Th.vertices[ivv[2]].y + gamma*Th.vertices[ivv[3]].y;
5138
 
        vertextetsub[iv].z = (1-alpha-beta-gamma)*Th.vertices[ivv[0]].z + alpha*Th.vertices[ivv[1]].z + beta*Th.vertices[ivv[2]].z + gamma*Th.vertices[ivv[3]].z;  
5139
 
      }
5140
 
    
5141
 
      int newindex[nvsub];
5142
 
      for( int iv=0; iv<nvsub; iv++){
5143
 
        const R3 viR3(vertextetsub[iv].x,vertextetsub[iv].y,vertextetsub[iv].z);
5144
 
        const Vertex3 &vi( viR3 );
5145
 
        Vertex3 * pvi=gtree->ToClose(vi,hseuil);
5146
 
        
5147
 
        if(!pvi){
5148
 
          v[np].x   = vi.x;
5149
 
          v[np].y   = vi.y;
5150
 
          v[np].z   = vi.z;
5151
 
          v[np].lab = K.lab;
5152
 
          newindex[iv] = np;
5153
 
          gtree->Add( v[np] );
5154
 
          np++;
5155
 
          }
5156
 
        else{
5157
 
          newindex[iv] = pvi-v;
5158
 
        }
5159
 
        
5160
 
          //assert(pvi);
5161
 
          //newindex[iv] = pvi-v;
5162
 
        if(np>nv) cout << "np=" << np << " nv=" << nv << endl; 
5163
 
        ffassert( np <= nv );
5164
 
      }
5165
 
 
5166
 
      for( int ii=0; ii<ntetsub; ii++){
5167
 
        int ivt[4];
5168
 
        for( int jj=0; jj< 4; jj++){
5169
 
          ivt[jj] = newindex[tetsub[4*ii+jj]];
5170
 
          assert( tetsub[4*ii+jj] < nvsub );
5171
 
          assert( ivt[jj] < np );
5172
 
        }
5173
 
        (tt++)->set( v, ivt, K.lab);
5174
 
        itt++;
5175
 
        assert( itt <= nt );
5176
 
      }
5177
 
      
5178
 
      for (int j=0;j<4;j++)
5179
 
      {
5180
 
          int jt=j,it=Th.ElementAdj(i,jt);          
5181
 
          
5182
 
          if ( ( (tagTonB[i]&tagb[j]) ==0 ) &&  !(it==i || it <0)  && !split[it]) 
5183
 
          {   
5184
 
              // new border not on boundary 
5185
 
              int ivb[3];
5186
 
              
5187
 
              for( int ii=0; ii<nfacesub; ii++)
5188
 
              {
5189
 
                  int iface = 3*FaceTriangle[j]*nfacesub+3*ii;
5190
 
                  
5191
 
                  for( int jjj=0; jjj<3; jjj++)
5192
 
                  {
5193
 
                      ivb[jjj] = newindex[ trisub[iface+jjj] ];
5194
 
                      assert( trisub[ iface+jjj ] < nvsub );
5195
 
                      assert( ivb[jjj] < np );
5196
 
                  }
5197
 
                  
5198
 
                  (bb++)->set( v, ivb, newbelabel);
5199
 
                  ie++;
5200
 
              } 
5201
 
          }
5202
 
          assert( ie <= nbe);
5203
 
      }
5204
 
      
5205
 
    }
5206
 
  }
5207
 
  if(verbosity>8)
5208
 
    cout << "   -- Number of new  border face not on Border " << ie << endl;
5209
 
  delete [] vertexsub; //[nvsub];
5210
 
  delete [] tetsub;   //[4*ntetsub];
5211
 
  delete [] trisub;   //[4*kksplit*kksplit];
5212
 
 
5213
 
  // split border elements
5214
 
  int nv2Dsub   = (kksplit+1)*(kksplit+2)/4;
5215
 
  int ntri2Dsub = kksplit2;
5216
 
  R2 *vertex2Dsub; //[nvsub];  
5217
 
  int *tri2Dsub;   //[4*kksplit*kksplit];
5218
 
 
5219
 
  SplitSimplex<R2>( kksplit, nv2Dsub, vertex2Dsub, ntri2Dsub, tri2Dsub);
5220
 
 
5221
 
 
5222
 
  for( int ibe=0; ibe < Th.nbe; ibe++)
5223
 
  {  
5224
 
      int iff;
5225
 
      int it=Th.BoundaryElement(ibe,iff);
5226
 
      int ifff=iff,itt=Th.ElementAdj(it,ifff);
5227
 
      if(itt<0) itt=it; 
5228
 
      if( split[it] == 0 && split[itt] == 0) continue; // boundary not on one element 
5229
 
 
5230
 
      const Triangle3 &K(Th.be(ibe));
5231
 
      int ivv[3];
5232
 
      
5233
 
      ivv[0] = Th.operator()(K[0]);
5234
 
      ivv[1] = Th.operator()(K[1]);
5235
 
      ivv[2] = Th.operator()(K[2]);
5236
 
      
5237
 
      R3 vertextrisub[nv2Dsub];
5238
 
      for( int iv=0; iv<nv2Dsub; iv++)
5239
 
      {
5240
 
          double alpha=vertex2Dsub[iv].x;
5241
 
          double beta=vertex2Dsub[iv].y;
5242
 
          
5243
 
          vertextrisub[iv].x = (1-alpha-beta)*Th.vertices[ivv[0]].x + alpha*Th.vertices[ivv[1]].x + beta*Th.vertices[ivv[2]].x;  
5244
 
          vertextrisub[iv].y = (1-alpha-beta)*Th.vertices[ivv[0]].y + alpha*Th.vertices[ivv[1]].y + beta*Th.vertices[ivv[2]].y;
5245
 
          vertextrisub[iv].z = (1-alpha-beta)*Th.vertices[ivv[0]].z + alpha*Th.vertices[ivv[1]].z + beta*Th.vertices[ivv[2]].z;  
5246
 
          
5247
 
      }
5248
 
      int newindex[nv2Dsub];
5249
 
      for( int iv=0; iv<nv2Dsub; iv++)
5250
 
      {
5251
 
          const Vertex3 &vi( vertextrisub[iv] );
5252
 
          Vertex3 * pvi=gtree->ToClose(vi,hseuil);
5253
 
          assert(pvi);
5254
 
          newindex[iv] = pvi-v;
5255
 
      }
5256
 
      
5257
 
      for( int ii=0; ii<nfacesub; ii++)
5258
 
      {
5259
 
          int ivb[3];
5260
 
          for( int jjj=0; jjj<3; jjj++)
5261
 
          {
5262
 
              ivb[jjj] = newindex[ tri2Dsub[3*ii+jjj] ]; 
5263
 
              assert( tri2Dsub[ 3*ii+jjj  ] < nvsub );
5264
 
              if(verbosity > 4 ) cout << "        " << ivb[jjj] << " np:" << np<< endl;
5265
 
              assert( ivb[jjj] < np );
5266
 
          }
5267
 
          
5268
 
          (bb++)->set( v, ivb, K.lab);
5269
 
          ie++;
5270
 
          assert( ie <= nbe);
5271
 
      }
5272
 
      
5273
 
      
5274
 
  }
5275
 
 
5276
 
  delete [] vertex2Dsub;   //[4*ntetsub];
5277
 
  delete [] tri2Dsub;   //[4*kksplit*kksplit];
5278
 
 
5279
 
  
5280
 
  if(verbosity>99)
5281
 
    {
5282
 
      cout << "nbofv initial" << Th.nv << endl; 
5283
 
      cout << "nv=" << nv << " np=" << np << endl; 
5284
 
      cout << "itt=" << itt << " nt=" << nt << endl;
5285
 
      cout << "ie=" << ie << " nbe=" << nbe << endl;
5286
 
    }
5287
 
  ffassert( nv == np );
5288
 
  ffassert( ie ==nbe);
5289
 
  ffassert( itt == nt );
5290
 
 
5291
 
  //delete gtree;
5292
 
 
5293
 
  Mesh3 *Tht = new Mesh3( nv, nt, nbe, v, t, b); 
5294
 
  Tht->BuildGTree(); // Add JM. Oct 2010 
5295
 
  delete gtree;
5296
 
  
5297
 
 
5298
 
  return Tht;
 
5058
    
 
5059
    
 
5060
    
 
5061
    int itt=0;
 
5062
    int ie=0;
 
5063
    
 
5064
    
 
5065
    Vertex3 *v=new Vertex3[nv];
 
5066
    Tet *t = new Tet[nt];
 
5067
    Tet *tt = t;
 
5068
    
 
5069
    Triangle3 *b  = new Triangle3[nbe];
 
5070
    Triangle3 *bb = b;
 
5071
    R3 hh = (bmax-bmax)/10.;
 
5072
    EF23::GTree<Vertex3> *gtree = new EF23::GTree<Vertex3>(v,bmin-hh,bmax+hh,0);
 
5073
    const R3 * pP[4];
 
5074
    int np=0; // nb of new points .. 
 
5075
 
 
5076
    {
 
5077
        KN<R3>  vertextetsub(nvsub);
 
5078
        KN<int> newindex (nvsub);
 
5079
        
 
5080
        for(int i=0; i<Th.nt; i++)
 
5081
            if(split[i])
 
5082
            {
 
5083
                const Tet &K(Th.elements[i]);
 
5084
                
 
5085
                for(int ii=0; ii< 4; ii++)
 
5086
                    pP[ii] = & K[ii];
 
5087
                
 
5088
                for( int iv=0; iv<nvsub; iv++)
 
5089
                    (R3&) vertextetsub[iv]= vertexsub[iv].Bary(pP);
 
5090
                
 
5091
                for( int iv=0; iv<nvsub; iv++)
 
5092
                {
 
5093
                    Vertex3 * pvi=gtree->ToClose(vertextetsub[iv],hseuil);
 
5094
                    
 
5095
                    if(!pvi)
 
5096
                    {
 
5097
                        (R3&) v[np]   = vertextetsub[iv];
 
5098
                        v[np].lab = K.lab;
 
5099
                        newindex[iv] = np;
 
5100
                        gtree->Add( v[np] );
 
5101
                        np++;
 
5102
                    }
 
5103
                    else
 
5104
                        newindex[iv] = pvi-v;
 
5105
                    
 
5106
                    ffassert( np <= nv );
 
5107
                }
 
5108
                
 
5109
                for( int ii=0; ii<ntetsub; ii++)
 
5110
                {
 
5111
                    int ivt[4];
 
5112
                    for( int jj=0; jj< 4; jj++)
 
5113
                    {
 
5114
                        ivt[jj] = newindex[tetsub[4*ii+jj]];
 
5115
                        assert( tetsub[4*ii+jj] < nvsub );
 
5116
                        assert( ivt[jj] < np );
 
5117
                    }
 
5118
                    (tt++)->set( v, ivt, K.lab);
 
5119
                    itt++;
 
5120
                    assert( itt <= nt );
 
5121
                }
 
5122
                
 
5123
                for (int j=0;j<4;j++)
 
5124
                {
 
5125
                    int jt=j,it=Th.ElementAdj(i,jt);
 
5126
                    
 
5127
                    if ( ( (tagTonB[i]&tagb[j]) ==0 ) &&  !(it==i || it <0)  && !split[it])
 
5128
                    {
 
5129
                        // new border not on boundary
 
5130
                        int ivb[3];
 
5131
                        
 
5132
                        for( int ii=0; ii<nfacesub; ii++)
 
5133
                        {
 
5134
                            int iface = 3*FaceTriangle[j]*nfacesub+3*ii;
 
5135
                            
 
5136
                            for( int jjj=0; jjj<3; jjj++)
 
5137
                            {
 
5138
                                ivb[jjj] = newindex[ trisub[iface+jjj] ];
 
5139
                                assert( trisub[ iface+jjj ] < nvsub );
 
5140
                                assert( ivb[jjj] < np );
 
5141
                            }
 
5142
                            
 
5143
                            (bb++)->set( v, ivb, newbelabel);
 
5144
                            ie++;
 
5145
                        }
 
5146
                    }
 
5147
                    assert( ie <= nbe);
 
5148
                    
 
5149
                }
 
5150
            }
 
5151
    }
 
5152
    if(verbosity>10)
 
5153
        cout  << "    ++ np=" << np << "==  nv=" << nv << endl;
 
5154
    ffassert( np == nv); 
 
5155
    if(verbosity>8)
 
5156
        cout << "   -- Number of new  border face not on Border " << ie << endl;
 
5157
    delete [] vertexsub; //[nvsub];
 
5158
    delete [] tetsub;   //[4*ntetsub];
 
5159
    delete [] trisub;   //[4*kksplit*kksplit];
 
5160
    
 
5161
    // split border elements
 
5162
    int nv2Dsub   = (kksplit+1)*(kksplit+2)/4;
 
5163
    int ntri2Dsub = kksplit2;
 
5164
    R2 *vertex2Dsub; //[nvsub];
 
5165
    int *tri2Dsub;   //[4*kksplit*kksplit];
 
5166
    
 
5167
    SplitSimplex<R2>( kksplit, nv2Dsub, vertex2Dsub, ntri2Dsub, tri2Dsub);
 
5168
    
 
5169
    
 
5170
    for( int ibe=0; ibe < Th.nbe; ibe++)
 
5171
    {
 
5172
        int iff;
 
5173
        int it=Th.BoundaryElement(ibe,iff);
 
5174
        int ifff=iff,itt=Th.ElementAdj(it,ifff);
 
5175
        if(itt<0) itt=it;
 
5176
        if( split[it] == 0 && split[itt] == 0) continue; // boundary not on one element
 
5177
        
 
5178
        const Triangle3 &K(Th.be(ibe));
 
5179
        int ivv[3];
 
5180
        
 
5181
        ivv[0] = Th.operator()(K[0]);
 
5182
        ivv[1] = Th.operator()(K[1]);
 
5183
        ivv[2] = Th.operator()(K[2]);
 
5184
        
 
5185
        R3 *vertextrisub = new R3  [nv2Dsub];
 
5186
        int *newindex = new int[nv2Dsub];
 
5187
        for( int iv=0; iv<nv2Dsub; iv++)
 
5188
        {
 
5189
            double alpha=vertex2Dsub[iv].x;
 
5190
            double beta=vertex2Dsub[iv].y;
 
5191
            
 
5192
            vertextrisub[iv].x = (1-alpha-beta)*Th.vertices[ivv[0]].x + alpha*Th.vertices[ivv[1]].x + beta*Th.vertices[ivv[2]].x;
 
5193
            vertextrisub[iv].y = (1-alpha-beta)*Th.vertices[ivv[0]].y + alpha*Th.vertices[ivv[1]].y + beta*Th.vertices[ivv[2]].y;
 
5194
            vertextrisub[iv].z = (1-alpha-beta)*Th.vertices[ivv[0]].z + alpha*Th.vertices[ivv[1]].z + beta*Th.vertices[ivv[2]].z;
 
5195
            
 
5196
        }
 
5197
        
 
5198
        for( int iv=0; iv<nv2Dsub; iv++)
 
5199
        {
 
5200
            const Vertex3 &vi( vertextrisub[iv] );
 
5201
            Vertex3 * pvi=gtree->ToClose(vi,hseuil);
 
5202
            assert(pvi);
 
5203
            newindex[iv] = pvi-v;
 
5204
        }
 
5205
        
 
5206
        for( int ii=0; ii<nfacesub; ii++)
 
5207
        {
 
5208
            int ivb[3];
 
5209
            for( int jjj=0; jjj<3; jjj++)
 
5210
            {
 
5211
                ivb[jjj] = newindex[ tri2Dsub[3*ii+jjj] ];
 
5212
                assert( tri2Dsub[ 3*ii+jjj  ] < nvsub );
 
5213
                if(verbosity > 199 ) cout << "        " << ivb[jjj] << " np:" << np<< endl;
 
5214
                assert( ivb[jjj] < np );
 
5215
            }
 
5216
            
 
5217
            (bb++)->set( v, ivb, K.lab);
 
5218
            ie++;
 
5219
            assert( ie <= nbe);
 
5220
        }
 
5221
        delete [] vertextrisub;
 
5222
        delete [] newindex;
 
5223
        
 
5224
        
 
5225
    }
 
5226
    
 
5227
    delete [] vertex2Dsub;   //[4*ntetsub];
 
5228
    delete [] tri2Dsub;   //[4*kksplit*kksplit];
 
5229
    
 
5230
    
 
5231
    if(verbosity>99)
 
5232
    {
 
5233
        cout << "nbofv initial" << Th.nv << endl;
 
5234
        cout << "nv=" << nv << " np=" << np << endl;
 
5235
        cout << "itt=" << itt << " nt=" << nt << endl;
 
5236
        cout << "ie=" << ie << " nbe=" << nbe << endl;
 
5237
    }
 
5238
    ffassert( nv == np );
 
5239
    ffassert( ie ==nbe);
 
5240
    ffassert( itt == nt );
 
5241
    
 
5242
    //delete gtree;
 
5243
    
 
5244
    Mesh3 *Tht = new Mesh3( nv, nt, nbe, v, t, b);
 
5245
    Tht->BuildGTree(); // Add JM. Oct 2010
 
5246
    delete gtree;
 
5247
    
 
5248
    
 
5249
    return Tht;
5299
5250
}
5300
5251
 
5301
5252
 
5302
 
AnyType Op_trunc_mesh3::Op::operator()(Stack stack)  const { 
 
5253
AnyType Op_trunc_mesh3::Op::operator()(Stack stack)  const {
5303
5254
    
5304
5255
  Mesh3 *pTh = GetAny<Mesh3 *>((*getmesh)(stack));
5305
5256
  Mesh3 &Th = *pTh;
5323
5274
    cout << "  -- Trunc mesh: Nb of Tetrahedrons = " << kk << " label=" <<label <<endl;
5324
5275
  Mesh3 * Tht = truncmesh(Th,kkksplit,split,false,label);
5325
5276
  
5326
 
  
5327
 
  //  cout << Tht->nv << " " << Tht->nt << " " << Tht->nbe << " " << endl;
5328
 
  //   cout << "==================================" <<  Tht << endl;
5329
 
  //   exit(1);
5330
 
  //   for( int jjj=0; jjj<Tht->nv; jjj++){
5331
 
  //     cout << "apres trunc mesh :: vertex" << jjj+1 <<" " <<  Tht->vertices[jjj].x << " " << Tht->vertices[jjj].y << " " << Tht->vertices[jjj].z << " " << Tht->vertices[jjj].lab << endl;
5332
 
  //if( abs(Tht->vertices[jjj].lab) > 1 ) exit(1); 
5333
 
  //  }
5334
 
  //string filename("Thtpp_res.mesh");
5335
 
  //Tht->Save(filename); 
5336
 
 
5337
 
  //  cout << "==================================" <<  Tht << endl;
5338
 
  
5339
5277
 
5340
5278
  Add2StackOfPtr2FreeRC(stack,Tht);//  07/2008 FH 
5341
5279
  *mp=mps;
5511
5449
                
5512
5450
                
5513
5451
                Mesh *pThnew = new Mesh(nv,nt,ns,v,t,b);  // attention aux composantes connexes.
 
5452
                //Lorenzo
 
5453
                R2 Pn,Px;
 
5454
                pThnew->BoundingBox(Pn,Px);
 
5455
                if(!pThnew->quadtree)
 
5456
                        pThnew->quadtree=new Fem2D::FQuadTree(pTh,Pn,Px,pTh->nv);
 
5457
                //Lorenzo
 
5458
                                
5514
5459
                return pThnew;
5515
5460
                
5516
5461
        }
5761
5706
    
5762
5707
    Mesh3 *pThnew = new Mesh3(nv,nt,ns,v,t,b);  // peut etre a définir ???
5763
5708
    // attention aux composantes connexes.
5764
 
                
 
5709
        pThnew->BuildGTree();  //Lorenzo
5765
5710
    
5766
5711
    return pThnew;
5767
5712
    
5822
5767
    
5823
5768
    cout <<" nv" << nv << " ns " << endl;  
5824
5769
    Mesh3 *pThnew = new Mesh3(nv,ns,v,b);  
5825
 
                
 
5770
        pThnew->BuildGTree();  //Lorenzo        
5826
5771
    return pThnew;    
5827
5772
  }
5828
 
 
 
5773
   
 
5774
  Mesh3 *pThnew = NULL;
 
5775
  return pThnew;    
5829
5776
}
5830
5777
 
5831
5778
 
5846
5793
  typedef Mesh *pmesh;
5847
5794
  typedef Mesh3 *pmesh3;
5848
5795
  
5849
 
  if (verbosity)
 
5796
  if (verbosity && mpirank == 0)
5850
5797
    cout << " load: msh3  " << endl;
5851
 
  //cout << " je suis dans Init " << endl; 
5852
5798
  
5853
5799
  TheOperators->Add("+",new OneBinaryOperator_st< Op3_addmesh<listMesh3,pmesh3,pmesh3>  >      );
5854
5800
  TheOperators->Add("+",new OneBinaryOperator_st< Op3_addmesh<listMesh3,listMesh3,pmesh3>  >      );