Changeset 2782


Ignore:
Timestamp:
01/07/10 14:31:05 (15 years ago)
Author:
Mathieu Morlighem
Message:

back to initial version of FreeFem++

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/Bamgx/Mesh2.cpp

    r2740 r2782  
    41284128                }
    41294129
    4130                 void Triangles::FillHoleInMesh(){
    4131 
     4130                void Triangles::FillHoleInMesh() {
    41324131                        Triangles * OldCurrentTh =CurrentTh;
    4133                         long int verbosity=200;
    4134                        
    41354132                        CurrentTh=this;
    41364133                        //  Int4 NbTold = nbt;
    41374134                        // generation of the integer coor
    4138 
    4139                         //  double coef = coefIcoor;
    4140                         // recherche des extrema des vertices pmin,pmax
    4141                         Int4 i;
    4142                         if(verbosity>2)
    4143                                 cout << "  -- FillHoleInMesh: Nb of vertices =" << nbv
     4135                          {
     4136
     4137                                //  double coef = coefIcoor;
     4138                                // recherche des extrema des vertices pmin,pmax
     4139                                Int4 i;
     4140                                if(verbosity>2)
     4141                                 cout << "  -- FillHoleInMesh: Nb of vertices =" << nbv
    41444142                                        << " Pmin = "<< pmin << " Pmax = "<< pmax << endl;
    41454143
    4146                         assert(ordre);
    4147                         for (i=0;i<nbv;i++)
    4148                                 ordre[i]= 0 ;
    4149 
    4150 
    4151                         NbSubDomains =0;
    4152 
    4153                         // generation of the adjacence of the triangles
    4154                         SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
    4155                         Int4 * st = new Int4[nbt*3];
    4156                         for (i=0;i<nbt*3;i++)
    4157                                 st[i]=-1;
    4158                         Int4 kk =0;
    4159                         for (i=0;i<nbe;i++)
    4160                                 kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
    4161                         if (kk != nbe)
    4162                         {
    4163                                 cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
    4164                                 MeshError(1002,this);
    4165                         }
    4166                         for (i=0;i<nbt;i++)
    4167                                 for (int j=0;j<3;j++)
    4168                                 {
    4169                                         // Int4 i0,i1;
    4170                                         Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
    4171                                                         Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
    4172                                         Int4 invisible = triangles[i].Hidden(j);
    4173                                         if(st[k]==-1)
    4174                                                 st[k]=3*i+j;
    4175                                         else if(st[k]>=0) {
    4176                                                 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));
    4177 
    4178                                                 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
    4179                                                 if (invisible)  triangles[i].SetHidden(j);
    4180                                                 if (k<nbe) {
    4181                                                         triangles[i].SetLocked(j);
    4182                                                 }
    4183                                                 st[k]=-2-st[k]; }
    4184                                         else {
    4185                                                 cerr << " The edge ("
     4144                                assert(ordre);
     4145                                for (i=0;i<nbv;i++)
     4146                                 ordre[i]= 0 ;
     4147
     4148
     4149                                NbSubDomains =0;
     4150
     4151                                // generation of the adjacence of the triangles
     4152                                SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
     4153                                Int4 * st = new Int4[nbt*3];
     4154                                for (i=0;i<nbt*3;i++)
     4155                                 st[i]=-1;
     4156                                Int4 kk =0;
     4157                                for (i=0;i<nbe;i++)
     4158                                 kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
     4159                                if (kk != nbe)
     4160                                  {
     4161                                        cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
     4162                                        MeshError(1002,this);
     4163                                  }
     4164                                for (i=0;i<nbt;i++)
     4165                                 for (int j=0;j<3;j++)
     4166                                        {
     4167                                         // Int4 i0,i1;
     4168                                         Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
     4169                                                                 Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
     4170                                         Int4 invisible = triangles[i].Hidden(j);
     4171                                         if(st[k]==-1)
     4172                                          st[k]=3*i+j;
     4173                                         else if(st[k]>=0) {
     4174                                                 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));
     4175
     4176                                                 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
     4177                                                 if (invisible)  triangles[i].SetHidden(j);
     4178                                                 if (k<nbe) {
     4179                                                         triangles[i].SetLocked(j);
     4180                                                 }
     4181                                                 st[k]=-2-st[k]; }
     4182                                         else {
     4183                                                 cerr << " The edge ("
    41864184                                                        << Number(triangles[i][VerticesOfTriangularEdge[j][0]])
    41874185                                                        << " , "
    41884186                                                        << Number(triangles[i][VerticesOfTriangularEdge[j][1]])
    41894187                                                        << " ) is in more than 2 triangles " <<k <<endl;
    4190                                                 cerr << " Edge " << j << " Of Triangle " << i << endl;
    4191                                                 cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3  << endl;
    4192                                                 cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
     4188                                                 cerr << " Edge " << j << " Of Triangle " << i << endl;
     4189                                                 cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3  << endl;
     4190                                                 cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
    41934191                                                        << " Of Triangle " <<  Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << endl;
    4194                                                 MeshError(9999,this);} 
    4195 
    4196 
    4197                                 }
    4198                        
    4199                         if(verbosity>5) {
    4200                                 cout << "    On Mesh " << name << endl;
    4201                                 cout << "    - The number of Vertices  = " << nbv << endl;
    4202                                 cout << "    - The number of Triangles = " << nbt << endl;
    4203                                 cout << "    - The number of given edge = " << nbe << endl;
    4204                                 cout << "    - The number of all edges = " << edge4->nb() << endl;
    4205                                 cout << "    - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl;
    4206                         }
    4207 
    4208 
    4209                         // check the consistant of edge[].adj and the geometrical required  vertex
    4210                         Int4 k=0;
    4211                         for (i=0;i<edge4->nb();i++)
    4212                                 if (st[i] >=0) // edge alone
    4213                                         if (i < nbe)
    4214                                         {
    4215                                                 Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
    4216                                                 Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
     4192                                                 MeshError(9999,this);}
     4193
     4194
    42174195                                        }
    4218                                         else {
    4219                                                 k++;
    4220                                                 if (verbosity>20 && k <20)
    4221                                                 {
    4222                                                         Int4 i0=edge4->i(i);
    4223                                                         Int4 i1=edge4->j(i);
    4224                                                         cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
    4225                                                 }
     4196                                if(verbosity>5) {
     4197                                        cout << "    On Mesh " << name << endl;
     4198                                        cout << "    - The number of Vertices  = " << nbv << endl;
     4199                                        cout << "    - The number of Triangles = " << nbt << endl;
     4200                                        cout << "    - The number of given edge = " << nbe << endl;
     4201                                        cout << "    - The number of all edges = " << edge4->nb() << endl;
     4202                                        cout << "    - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; }
     4203
     4204
     4205                                        // check the consistant of edge[].adj and the geometrical required  vertex
     4206                                        Int4 k=0;
     4207                                        for (i=0;i<edge4->nb();i++)
     4208                                         if (st[i] >=0) // edge alone
     4209                                          if (i < nbe)
     4210                                                 {
     4211                                                  Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
     4212                                                  Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
     4213                                                 }
     4214                                          else {
     4215                                                  k++;
     4216                                                  if (verbosity>20 && k <20)
     4217                                                         {
     4218                                                          Int4 i0=edge4->i(i);
     4219                                                          Int4 i1=edge4->j(i);
     4220                                                          cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
     4221                                                         }
     4222                                          }
     4223
     4224                                        if(k != 0) {
     4225                                                if (verbosity>20)
     4226                                                  {
     4227                                                        cout << " The given edge are " << endl;
     4228                                                        for (int i=0;i< nbe;i++)
     4229                                                         cout <<  " Edge " << i << " : " <<  Number(edges[i][0]) << " " <<  Number(edges[i][1])
     4230                                                                << " " << edges[i].ref << endl;
     4231                                                  }
     4232                                                cerr << k << " boundary edges  are not defined as edges " << endl;
     4233                                                MeshError(9998,this);
    42264234                                        }
    4227 
    4228                         if(k != 0) {
    4229                                 if (verbosity>20)
    4230                                 {
    4231                                         cout << " The given edge are " << endl;
    4232                                         for (int i=0;i< nbe;i++)
    4233                                                 cout <<  " Edge " << i << " : " <<  Number(edges[i][0]) << " " <<  Number(edges[i][1])
    4234                                                         << " " << edges[i].ref << endl;
    4235                                 }
    4236                                 cerr << k << " boundary edges  are not defined as edges " << endl;
    4237                                 MeshError(9998,this);
    4238                         }
    4239                         printf("salope1\n");
    4240 
    4241                         // generation of the mesh with boundary points   
    4242                         Int4 nbvb = 0;
    4243                         for (i=0;i<nbv;i++)
    4244                         {
    4245                                 vertices[i].t=0;
    4246                                 vertices[i].vint=0;
    4247                                 if (ordre[i])
    4248                                         ordre[nbvb++] = ordre[i];
    4249                         }
    4250 
    4251                         Triangle *savetriangles= triangles;
    4252                         Int4 savenbt=nbt;
    4253                         Int4 savenbtx=nbtx;
    4254                         SubDomain * savesubdomains = subdomains;
    4255                         subdomains = 0;
    4256 
    4257                         Int4  Nbtriafillhole = 2*nbvb;
    4258                         Triangle * triafillhole =new Triangle[Nbtriafillhole];
    4259                         if (verbosity>9)
    4260                                 cout << " Nbtriafillhole triafillhole*" << triafillhole << endl;
    4261                         triangles =  triafillhole;
    4262 
    4263                         nbt=2;
    4264                         nbtx= Nbtriafillhole;
    4265 
    4266                         for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){
    4267                                 if  ( ++i >= nbvb) {
    4268                                         cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
    4269                                         MeshError(998,this);
    4270                                 }
    4271                         }
    4272                         Exchange( ordre[2], ordre[i]);
    4273 
    4274                         Vertex *  v0=ordre[0], *v1=ordre[1];
    4275 
    4276 
    4277                         triangles[0](0) = 0; // sommet pour infini
    4278                         triangles[0](1) = v0;
    4279                         triangles[0](2) = v1;
    4280 
    4281                         triangles[1](0) = 0;// sommet pour infini
    4282                         triangles[1](2) = v0;
    4283                         triangles[1](1) = v1;
    4284                         const int e0 = OppositeEdge[0];
    4285                         const int e1 = NextEdge[e0];
    4286                         const int e2 = PreviousEdge[e0];
    4287                         triangles[0].SetAdj2(e0, &triangles[1] ,e0);
    4288                         triangles[0].SetAdj2(e1, &triangles[1] ,e2);
    4289                         triangles[0].SetAdj2(e2, &triangles[1] ,e1);
    4290 
    4291                         triangles[0].det = -1;  // faux triangles
    4292                         triangles[1].det = -1;  // faux triangles
    4293 
    4294                         triangles[0].SetTriangleContainingTheVertex();
    4295                         triangles[1].SetTriangleContainingTheVertex();
    4296 
    4297                         triangles[0].link=&triangles[1];
    4298                         triangles[1].link=&triangles[0];
    4299 
    4300                         printf("salope2\n");
    4301                         //  nbtf = 2;
    4302                         if (  !quadtree ) delete  quadtree; // ->ReInitialise();
    4303 
    4304                         quadtree = new QuadTree(this,0);
    4305                         quadtree->Add(*v0);
    4306                         quadtree->Add(*v1);
    4307 
    4308                         // on ajoute les sommets un a un
    4309                         Int4 NbSwap=0;
    4310                         for (Int4 icount=2; icount<nbvb; icount++) {
    4311 
    4312                                 Vertex *vi  = ordre[icount];
    4313                                 //        cout << " Add vertex " <<  Number(vi) << endl;
    4314                                 Icoor2 dete[3];
    4315                                 Triangle *tcvi = FindTriangleContening(vi->i,dete);
    4316                                 quadtree->Add(*vi);
    4317                                 Add(*vi,tcvi,dete);
    4318                                 NbSwap += vi->Optim(1,1);
    4319 
    4320                         }// end loop on  icount
    4321                         printf("salope3\n");
    4322 
    4323                         //Int4 nbtfillhole = nbt;
    4324                         // inforce the boundary
    4325                         TriangleAdjacent ta(0,0);
    4326                         Int4 nbloss = 0,knbe=0;
    4327                         for ( i = 0; i < nbe; i++)
    4328                                 if (st[i] >=0)  // edge alone => on border ...  FH oct 2009
    4329                                 {
    4330                                         Vertex & a=edges[i][0], & b =    edges[i][1];
    4331                                         if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
    4332                                         {
    4333                                                 knbe++;
    4334                                                 if (ForceEdge(a,b,ta)<0)
    4335                                                         nbloss++;
    4336                                         }
    4337                                 }
    4338                         if(nbloss)
    4339                         {
    4340                                 cerr << " we loss some  " << nbloss << " "  << " edges other " << knbe << endl;
    4341                                 MeshError(1100,this);
    4342                         }
    4343                         FindSubDomain(1);
    4344                         // remove all the hole
    4345                         // remove all the good sub domain
    4346                         Int4 krm =0;
    4347                         for (i=0;i<nbt;i++)
    4348                                 if (triangles[i].link) // remove triangles
    4349                                 {
    4350                                         krm++;
    4351                                         for (int j=0;j<3;j++)
    4352                                         {
    4353                                                 TriangleAdjacent ta =  triangles[i].Adj(j);
    4354                                                 Triangle & tta = * (Triangle *) ta;
    4355                                                 if(! tta.link) // edge between remove and not remove
    4356                                                 { // change the link of ta;
    4357                                                         int ja = ta;
    4358                                                         Vertex *v0= ta.EdgeVertex(0);
    4359                                                         Vertex *v1= ta.EdgeVertex(1);
    4360                                                         Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
    4361                                                         assert(st[k] >=0);
    4362                                                         tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
    4363                                                         ta.SetLock();
    4364                                                         st[k]=-2-st[k];
    4365                                                 }
    4366                                         }
    4367                                 }
    4368                         Int4 NbTfillHoll =0;
    4369                         for (i=0;i<nbt;i++)
    4370                                 if (triangles[i].link) {
    4371                                         triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
    4372                                         triangles[i].color=-1;
    4373                                 }
    4374                                 else
    4375                                 {
    4376                                         triangles[i].color= savenbt+ NbTfillHoll++;
    4377                                 }
    4378                         // cout <<      savenbt+NbTfillHoll << " " <<  savenbtx  << endl;
    4379                         assert(savenbt+NbTfillHoll <= savenbtx );
    4380                         // copy of the outside triangles in saveTriangles
    4381                         for (i=0;i<nbt;i++)
    4382                                 if(triangles[i].color>=0)
    4383                                 {
    4384                                         savetriangles[savenbt]=triangles[i];
    4385                                         savetriangles[savenbt].link=0;
    4386                                         savenbt++;
    4387                                 }
    4388                         // gestion of the adj
    4389                         k =0;
    4390                         Triangle * tmax = triangles + nbt;
    4391                         for (i=0;i<savenbt;i++) 
    4392                         {
    4393                                 Triangle & ti = savetriangles[i];
    4394                                 for (int j=0;j<3;j++)
    4395                                 {
    4396                                         Triangle * ta = ti.TriangleAdj(j);
    4397                                         int aa = ti.NuEdgeTriangleAdj(j);
    4398                                         int lck = ti.Locked(j);
    4399                                         if (!ta) k++; // bug
    4400                                         else if ( ta >= triangles && ta < tmax)
    4401                                         {
    4402                                                 ta= savetriangles + ta->color;
    4403                                                 ti.SetAdj2(j,ta,aa);
    4404                                                 if(lck) ti.SetLocked(j);
    4405                                         }
    4406                                 }
    4407                         }
    4408                         //       OutSidesTriangles = triangles;
    4409                         //      Int4 NbOutSidesTriangles = nbt;
    4410 
    4411                         // restore triangles;
    4412                         nbt=savenbt;
    4413                         nbtx=savenbtx;
    4414                         delete [] triangles;
    4415                         delete [] subdomains;
    4416                         triangles = savetriangles;
    4417                         subdomains = savesubdomains;
    4418                         //       cout <<  triangles << " <> " << OutSidesTriangles << endl;
    4419                         /*       k=0;
    4420                                  for (i=0;i<nbt;i++)
    4421                                  for (int j=0;j<3;j++)
    4422                                  if (!triangles[i].TriangleAdj(j))
    4423                                  k++;
    4424 
    4425 */
    4426 
    4427                         printf("salope4\n");
    4428                         if (k) {
    4429                                 cerr << "Error Nb of triangles edge alone = " << k << endl;
    4430                                 MeshError(9997,this);
    4431                         }
    4432                         printf("salope4a\n");
    4433                         FindSubDomain();
    4434                         // cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
    4435 
    4436                         //
    4437 
    4438                         printf("salope4b\n");
    4439                         delete edge4;
    4440                         printf("salope4ba\n");
    4441                         delete [] st;
    4442                         printf("salope4bb %i\n",nbv);
    4443                         for (i=0;i<nbv;i++){
    4444                                 printf("%i\n",i);
    4445                                 quadtree->Add(vertices[i]);
    4446                         }
    4447                         printf("salope4c\n");
    4448 
    4449                         SetVertexFieldOn();
    4450                         printf("salope5a\n");
    4451 
    4452                         for (i=0;i<nbe;i++){
    4453                                 if(edges[i].on){
    4454                                         for(int j=0;j<2;j++){
    4455                                                 if (!edges[i].adj[j]){
    4456                                                         if(!edges[i][j].on->IsRequiredVertex()) {
    4457                                                                 cerr << " Erreur adj et sommet requis edges [" << i <<  "][ " << j << "]= "
    4458                                                                         <<  Number(edges[i][j]) << " : "  << " on = " << Gh.Number(edges[i].on) ;
    4459                                                                 if (edges[i][j].on->OnGeomVertex())
     4235                                        // generation of the mesh with boundary points   
     4236                                        Int4 nbvb = 0;
     4237                                        for (i=0;i<nbv;i++)
     4238                                          {
     4239                                                vertices[i].t=0;
     4240                                                vertices[i].vint=0;
     4241                                                if (ordre[i])
     4242                                                 ordre[nbvb++] = ordre[i];
     4243                                          }
     4244
     4245                                        Triangle *savetriangles= triangles;
     4246                                        Int4 savenbt=nbt;
     4247                                        Int4 savenbtx=nbtx;
     4248                                        SubDomain * savesubdomains = subdomains;
     4249                                        subdomains = 0;
     4250
     4251                                        Int4  Nbtriafillhole = 2*nbvb;
     4252                                        Triangle * triafillhole =new Triangle[Nbtriafillhole];
     4253                                        if (verbosity>9)
     4254                                         cout << " Nbtriafillhole triafillhole*" << triafillhole << endl;
     4255                                        triangles =  triafillhole;
     4256
     4257                                        nbt=2;
     4258                                        nbtx= Nbtriafillhole;
     4259
     4260                                        for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
     4261                                         if  ( ++i >= nbvb) {
     4262                                                 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
     4263                                                 MeshError(998,this); }
     4264                                                 Exchange( ordre[2], ordre[i]);
     4265
     4266                                                 Vertex *  v0=ordre[0], *v1=ordre[1];
     4267
     4268
     4269                                                 triangles[0](0) = 0; // sommet pour infini
     4270                                                 triangles[0](1) = v0;
     4271                                                 triangles[0](2) = v1;
     4272
     4273                                                 triangles[1](0) = 0;// sommet pour infini
     4274                                                 triangles[1](2) = v0;
     4275                                                 triangles[1](1) = v1;
     4276                                                 const int e0 = OppositeEdge[0];
     4277                                                 const int e1 = NextEdge[e0];
     4278                                                 const int e2 = PreviousEdge[e0];
     4279                                                 triangles[0].SetAdj2(e0, &triangles[1] ,e0);
     4280                                                 triangles[0].SetAdj2(e1, &triangles[1] ,e2);
     4281                                                 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
     4282
     4283                                                 triangles[0].det = -1;  // faux triangles
     4284                                                 triangles[1].det = -1;  // faux triangles
     4285
     4286                                                 triangles[0].SetTriangleContainingTheVertex();
     4287                                                 triangles[1].SetTriangleContainingTheVertex();
     4288
     4289                                                 triangles[0].link=&triangles[1];
     4290                                                 triangles[1].link=&triangles[0];
     4291
     4292#ifdef DEBUG
     4293                                                 triangles[0].check();
     4294                                                 triangles[1].check();
     4295#endif 
     4296                                                 //  nbtf = 2;
     4297                                                 if (  !quadtree )
     4298                                                  delete  quadtree; // ->ReInitialise();
     4299
     4300                                                 quadtree = new QuadTree(this,0);
     4301                                                 quadtree->Add(*v0);
     4302                                                 quadtree->Add(*v1);
     4303
     4304                                                 // on ajoute les sommets un a un
     4305                                                 Int4 NbSwap=0;
     4306                                                 for (Int4 icount=2; icount<nbvb; icount++) {
     4307
     4308                                                         Vertex *vi  = ordre[icount];
     4309                                                         //       cout << " Add vertex " <<  Number(vi) << endl;
     4310                                                         Icoor2 dete[3];
     4311                                                         Triangle *tcvi = FindTriangleContening(vi->i,dete);
     4312                                                         quadtree->Add(*vi);
     4313                                                         Add(*vi,tcvi,dete);
     4314                                                         NbSwap += vi->Optim(1,1);
     4315
     4316#ifdef DRAWING2
     4317                                                         cout << Number(vi) << " " <<  NbSwap <<  endl;
     4318                                                         reffecran();
     4319                                                         Draw();
     4320                                                         vi->Draw();
     4321                                                         inquire();
     4322#endif
     4323                                                 }// end loop on  icount       
     4324#ifdef DRAWING1
     4325                                                 inquire();
     4326#endif
     4327
     4328                                                 //Int4 nbtfillhole = nbt;
     4329                                                 // inforce the boundary
     4330                                                 TriangleAdjacent ta(0,0);
     4331                                                 Int4 nbloss = 0,knbe=0;
     4332                                                 for ( i = 0; i < nbe; i++)
     4333                                                  if (st[i] >=0)  // edge alone => on border ...  FH oct 2009
     4334                                                         {
     4335                                                          Vertex & a=edges[i][0], & b =    edges[i][1];
     4336                                                          if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
     4337                                                                 {
     4338                                                                  knbe++;
     4339                                                                  if (ForceEdge(a,b,ta)<0)
     4340                                                                        nbloss++;
     4341                                                                 }
     4342                                                         }
     4343                                                 if(nbloss)
     4344                                                        {
     4345                                                         cerr << " we loss some  " << nbloss << " "  << " edges other " << knbe << endl;
     4346                                                         MeshError(1100,this);
     4347                                                        }
     4348                                                 FindSubDomain(1);
     4349                                                 // remove all the hole
     4350                                                 // remove all the good sub domain
     4351                                                 Int4 krm =0;
     4352                                                 for (i=0;i<nbt;i++)
     4353                                                  if (triangles[i].link) // remove triangles
     4354                                                         {
     4355                                                          krm++;
     4356                                                          for (int j=0;j<3;j++)
     4357                                                                 {
     4358                                                                  TriangleAdjacent ta =  triangles[i].Adj(j);
     4359                                                                  Triangle & tta = * (Triangle *) ta;
     4360                                                                  if(! tta.link) // edge between remove and not remove
     4361                                                                         { // change the link of ta;
     4362                                                                          int ja = ta;
     4363                                                                          Vertex *v0= ta.EdgeVertex(0);
     4364                                                                          Vertex *v1= ta.EdgeVertex(1);
     4365                                                                          Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
     4366                                                                          assert(st[k] >=0);
     4367                                                                          tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
     4368                                                                          ta.SetLock();
     4369                                                                          st[k]=-2-st[k];
     4370                                                                         }
     4371                                                                 }
     4372                                                         }
     4373                                                 Int4 NbTfillHoll =0;
     4374                                                 for (i=0;i<nbt;i++)
     4375                                                  if (triangles[i].link) {
     4376                                                          triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);
     4377                                                          triangles[i].color=-1;
     4378                                                  }
     4379                                                  else
     4380                                                         {
     4381                                                          triangles[i].color= savenbt+ NbTfillHoll++;
     4382#ifdef DEBUG
     4383                                                          triangles[i].check();
     4384#endif
     4385                                                         }
     4386                                                 // cout <<      savenbt+NbTfillHoll << " " <<  savenbtx  << endl;
     4387                                                 assert(savenbt+NbTfillHoll <= savenbtx );
     4388                                                 // copy of the outside triangles in saveTriangles
     4389                                                 for (i=0;i<nbt;i++)
     4390                                                  if(triangles[i].color>=0)
     4391                                                         {
     4392                                                          savetriangles[savenbt]=triangles[i];
     4393                                                          savetriangles[savenbt].link=0;
     4394                                                          savenbt++;
     4395                                                         }
     4396                                                 // gestion of the adj
     4397                                                 k =0;
     4398                                                 Triangle * tmax = triangles + nbt;
     4399                                                 for (i=0;i<savenbt;i++) 
     4400                                                        {
     4401                                                         Triangle & ti = savetriangles[i];
     4402                                                         for (int j=0;j<3;j++)
     4403                                                                {
     4404                                                                 Triangle * ta = ti.TriangleAdj(j);
     4405                                                                 int aa = ti.NuEdgeTriangleAdj(j);
     4406                                                                 int lck = ti.Locked(j);
     4407                                                                 if (!ta) k++; // bug
     4408                                                                 else if ( ta >= triangles && ta < tmax)
     4409                                                                        {
     4410                                                                         ta= savetriangles + ta->color;
     4411                                                                         ti.SetAdj2(j,ta,aa);
     4412                                                                         if(lck) ti.SetLocked(j);
     4413                                                                        }
     4414                                                                }
     4415                                                        }
     4416                                                 //      OutSidesTriangles = triangles;
     4417                                                 //     Int4 NbOutSidesTriangles = nbt;
     4418
     4419                                                 // restore triangles;
     4420                                                 nbt=savenbt;
     4421                                                 nbtx=savenbtx;
     4422                                                 delete [] triangles;
     4423                                                 delete [] subdomains;
     4424                                                 triangles = savetriangles;
     4425                                                 subdomains = savesubdomains;
     4426                                                 //      cout <<  triangles << " <> " << OutSidesTriangles << endl;
     4427                                                 /*      k=0;
     4428                                                                 for (i=0;i<nbt;i++)
     4429                                                                 for (int j=0;j<3;j++)
     4430                                                                 if (!triangles[i].TriangleAdj(j))
     4431                                                                 k++;
     4432                                                                 */
     4433                                                 if (k) {
     4434                                                         cerr << "Error Nb of triangles edge alone = " << k << endl;
     4435                                                         MeshError(9997,this);
     4436                                                 }
     4437                                                 FindSubDomain();
     4438                                                 // cout << " NbTOld = " << NbTold << " ==  " << nbt - NbOutT << " " << nbt << endl;
     4439
     4440                                                 //
     4441
     4442                                                 delete edge4;
     4443                                                 delete [] st;
     4444                                                 for (i=0;i<nbv;i++)
     4445                                                  quadtree->Add(vertices[i]);
     4446
     4447                                                 SetVertexFieldOn();
     4448
     4449                                                 for (i=0;i<nbe;i++)
     4450                                                  if(edges[i].on)
     4451                                                        for(int j=0;j<2;j++)
     4452                                                         if (!edges[i].adj[j])
     4453                                                          if(!edges[i][j].on->IsRequiredVertex()) {
     4454                                                                  cerr << " Erreur adj et sommet requis edges [" << i <<  "][ " << j << "]= "
     4455                                                                         <<  Number(edges[i][j]) << " : "  << " on = " << Gh.Number(edges[i].on) ;
     4456                                                                  if (edges[i][j].on->OnGeomVertex())
    44604457                                                                        cerr << " vertex " << Gh.Number(edges[i][j].on->gv);
    4461                                                                 else if (edges[i][j].on->OnGeomEdge())
     4458                                                                  else if (edges[i][j].on->OnGeomEdge())
    44624459                                                                        cerr << "Edges " << Gh.Number(edges[i][j].on->ge);
    4463                                                                 else
     4460                                                                  else
    44644461                                                                        cerr << " = " << edges[i][j].on ;
    4465                                                                 cerr << endl;
    4466                                                         }
    4467                                                 }
    4468                                         }
    4469                                 }
    4470                         }
    4471 
    4472                         printf("salope5\n");
     4462                                                                  cerr << endl;
     4463                                                          }
     4464                          }
    44734465                        CurrentTh=OldCurrentTh;
    44744466                }
Note: See TracChangeset for help on using the changeset viewer.