The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings MMD
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

1 /*---------------------------------------------------------------------------*\
2     \o/\o/\o/
3     MMD
4     Version : 0.0.3
5     Web : https://github.com/alainbastide/MMD
6 -------------------------------------------------------------------------------
7 License
9     MMD is free software: you can redistribute it and/or modify it
10     under the terms of the GNU General Public License as published by the
11     Free Software Foundation, either version 3 of the License, or (at your
12     option) any later version.
13     
14     MMD is distributed in the hope that it will be useful, but
15     WITHOUT ANY WARRANTY; without even the implied warranty of
16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17     General Public License for more details.
18     
19     You should have received a copy of the GNU General Public License
20     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
21     
22 Application
23     explain
24     
25 Description
26     explain
27     
28 Author
29     Alain Bastide, Université de La Réunion, FRANCE.  All rights reserved
30     
31 \*---------------------------------------------------------------------------*/
33 #include "mmd.h"
34 #include "primal.h"
35 #include "memory.h"
36 #include "tensor.h"
37 #include "listop.h"
40 DEF_hexaedron_localNodeList
41 DEF_hexaedron_localNodeListNumbers
42 DEF_hexaedron_localSegmentList
43 DEF_hexaedron_localSegmentsListNumbers
49 void setHexahedreVertex(struct primalMesh * myPrimalMesh)
50 /*
51  * void setHexahedreVertex(struct primalMesh * myPrimalMesh)
52  * 
53  * Build vertex for each point without overlap
54  * 
55  * Input : myPrimalMesh struct
56  * 
57  * Use : L= myPrimalMesh->L
58  *          myPrimalMesh->M
59  *          myPrimalMesh->l
60  *          myPrimalMesh->N
61  *          myPrimalMesh->H
62  *          myPrimalMesh->P
63  * 
64  * Allocate : myPrimalMesh->vertex
65  * Produce : myPrimalMesh->vertex
66  * 
67  */
68 {
69   
70   //    clock_t t;
71   //    t = clock();
72   
73   clock_t t=startFunction(__FUNCTION__);
74   affiche("\n");
75   
76   // definition des vertex
77   connectivity_int it_vertex=0;
78   //    connectivity_int it_cell=0;
79   
80   dataType L=myPrimalMesh->L/myPrimalMesh->M;
81   dataType l=myPrimalMesh->l/myPrimalMesh->N;
82   dataType H=myPrimalMesh->H/myPrimalMesh->P;
83   
84   allocMemoryForMatrixDataType(&(myPrimalMesh->vertex) , myPrimalMesh->vertexNumber,DIM3D);
85   
86   for(connectivity_int i=0;i<myPrimalMesh->M+1;i++)
87     for(connectivity_int j=0;j<myPrimalMesh->N+1;j++)
88       for(connectivity_int k=0;k<myPrimalMesh->P+1;k++)
89         {
90           it_vertex = k + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
91           myPrimalMesh->vertex[it_vertex][0] = i * L;
92           myPrimalMesh->vertex[it_vertex][1] = j * l;
93           myPrimalMesh->vertex[it_vertex][2] = k * H;
94 #ifdef DEBUG  
95           debug_print("%-40s : %ld (%f,%f,%f)\n","myPrimalMesh->vertex", it_vertex, myPrimalMesh->vertex[it_vertex][0], myPrimalMesh->vertex[it_vertex][1], myPrimalMesh->vertex[it_vertex][2]);
96 #endif
97         }
98   
99   endFunction(__FUNCTION__, t);
100   
101   
104 void setHexahedreCellToVertexNumber(struct primalMesh * myPrimalMesh)
105 /*
106  * void setHexahedreCellToVertexNumber(struct primalMesh * myPrimalMesh)
107  * 
108  * This function builds a vector with the number of vertices of each cell; 
109  * here, we only use the hexahedron (HEXAHEDRON) but the code could be 
110  * updated to use several cell types
111  * 
112  * Input : myPrimalMesh struct
113  * 
114  * Use : 
115  * 
116  * Allocate : myPrimalMesh->cellToVertexNumber
117  * Produce : myPrimalMesh->cellToVertexNumber
118  * 
119  */
121   clock_t t=startFunction(__FUNCTION__);
122   
123   // definition des vertex
124   //    connectivity_int it_vertex=0;
125   connectivity_int it_cell=0;
126   
127   
128   myPrimalMesh->cellToVertexNumber = (connectivity_short *) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
129   
130   // hexaedron number of vertex on hexaedron
131   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
132     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
133       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
134         {
135           it_cell     = k + j*myPrimalMesh->P + i*myPrimalMesh->P*myPrimalMesh->N;
136           myPrimalMesh->cellToVertexNumber[it_cell]=HEXAHEDRON;
137           debug_print("%-40s : %ld (%d)\n","myPrimalMesh->cellToVertexNumber", it_cell, myPrimalMesh->cellToVertexNumber[it_cell] );
138         }
139   endFunction(__FUNCTION__, t);
140   
145 void setHexahedreCellToVertex(struct primalMesh * myPrimalMesh)
146 /*
147  * void setHexahedreCellToVertex(struct primalMesh * myPrimalMesh)
148  * 
149  * This function builds a vector by assigning for each cell the 
150  * number of vertices. The order of the vertices is ordered with 
151  * a positive sense of membership.
152  * 
153  * Input : myPrimalMesh struct
154  * 
155  * Use : 
156  * 
157  * Allocate : myPrimalMesh->cellToVertex
158  * Produce : myPrimalMesh->cellToVertex
159  * 
160  */
162   
163   clock_t t=startFunction(__FUNCTION__);
164   
165   connectivity_int it_vertex=0;
166   connectivity_int it_vertex2=0;
167   connectivity_int it_vertex3=0;
168   connectivity_int it_vertex4=0;
169   connectivity_int it_vertex5=0;
170   connectivity_int it_vertex6=0;
171   connectivity_int it_vertex7=0;
172   connectivity_int it_vertex8=0;
173   connectivity_int it_cell=0;
174   // cell 0 : 0 1 9 8 / 2 3 11 10
175   // hexaedron vertex number to define hexaedron
176   
177   connectivity_int vN1 = (myPrimalMesh->P+1) ;
178   connectivity_int vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1) ;
179   connectivity_int cN1 = myPrimalMesh->P*myPrimalMesh->N;
180   
181   
182   
183   myPrimalMesh->cellToVertex = (connectivity_int **) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_int *));
184   for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
185     {
186       myPrimalMesh->cellToVertex[i] = (connectivity_int *) malloc(myPrimalMesh->cellToVertexNumber[i] * sizeof(connectivity_int));
187       memset(myPrimalMesh->cellToVertex[i], 0, myPrimalMesh->cellToVertexNumber[i]*sizeof(connectivity_int));
188     }
189   
190   
191   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
192     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
193       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
194         {
195           it_cell     = k + j*myPrimalMesh->P + i*cN1; //cN1 = myPrimalMesh->P*myPrimalMesh->N;
196           
197           /*
198 *                 it_vertex =   k      + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
199             it_vertex2 = (k + 1) + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
200             it_vertex3 = (k + 1) + j*(myPrimalMesh->P+1) +(i+1)*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
201             it_vertex4 =  k      + j*(myPrimalMesh->P+1) +(i+1)*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
202 */
203           it_vertex =   k      + j*vN1 +i*vN2; // vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1);
204           it_vertex2 = it_vertex  + 1;
205           it_vertex3 = it_vertex2 + vN2; // vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1);
206           it_vertex4 = it_vertex3 - 1;
207           
208           myPrimalMesh->cellToVertex[it_cell][0]=it_vertex;
209           myPrimalMesh->cellToVertex[it_cell][1]=it_vertex2;
210           myPrimalMesh->cellToVertex[it_cell][2]=it_vertex3;
211           myPrimalMesh->cellToVertex[it_cell][3]=it_vertex4;
212           
213           
214           
215           it_vertex5 =  k      + (j+1)*vN1 +i*vN2;
216           it_vertex6 = (k + 1) + (j+1)*vN1 +i*vN2;
217           it_vertex7 = (k + 1) + (j+1)*vN1 +(i+1)*vN2;
218           it_vertex8 =  k      + (j+1)*vN1 +(i+1)*vN2;
219           
220           myPrimalMesh->cellToVertex[it_cell][4]=it_vertex5;
221           myPrimalMesh->cellToVertex[it_cell][5]=it_vertex6;
222           myPrimalMesh->cellToVertex[it_cell][6]=it_vertex7;
223           myPrimalMesh->cellToVertex[it_cell][7]=it_vertex8;
224           
225           
226 #ifdef DEBUG  
227           
228           debug_print("%-40s : %ld (%ld %ld %ld %ld %ld %ld %ld %ld)\n",
229                       "myPrimalMesh->cellToVertex",
230                       it_cell,
231                       myPrimalMesh->cellToVertex[it_cell][0],
232               myPrimalMesh->cellToVertex[it_cell][1],
233               myPrimalMesh->cellToVertex[it_cell][2],
234               myPrimalMesh->cellToVertex[it_cell][3],
235               myPrimalMesh->cellToVertex[it_cell][4],
236               myPrimalMesh->cellToVertex[it_cell][5],
237               myPrimalMesh->cellToVertex[it_cell][6],
238               myPrimalMesh->cellToVertex[it_cell][7] );
239 #endif
240           
241         }
242   endFunction(__FUNCTION__, t);
247 void setHexahedreVertexToCellNumbers(struct primalMesh * myPrimalMesh)
248 /*
249  * void setHexahedreVertexToCellNumbers(struct primalMesh * myPrimalMesh)
250  * 
251  * This function constructs a vector by assigning for each vertex 
252  * the number of cells attached to it. 
253  * 
254  * Input : myPrimalMesh struct
255  * 
256  * Use :  myPrimalMesh->cellNumber
257  *        myPrimalMesh->vertexNumber 
258  *        myPrimalMesh->cellToVertex
259  *        myPrimalMesh->cellToVertexNumber
260  * 
261  * Allocate : myPrimalMesh->vertexToCellNumber
262  * Produce : myPrimalMesh->vertexToCellNumber
263  * 
264  */
266   
267   clock_t t=startFunction(__FUNCTION__);
268   affiche("\n");
269   
270   myPrimalMesh->vertexToCellNumber = (connectivity_short *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_short));
271   memset(myPrimalMesh->vertexToCellNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
272   
273   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
274     {
275       for(connectivity_int i=0;i<myPrimalMesh->cellToVertexNumber[k];i++)
276         {
277           myPrimalMesh->vertexToCellNumber[myPrimalMesh->cellToVertex[k][i]]++;
278         }
279     }
280   
281 #ifdef DEBUG  
282   
283   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
284     {
285       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToCellNumbers", k);
286       
287       {
288         debug_print( "%d " , myPrimalMesh->vertexToCellNumber[k]);
289       }
290       debug_print(")\n");
291     }
292 #endif  
293   endFunction(__FUNCTION__, t);
294   
300 void setHexahedreVertexToCells(struct primalMesh * myPrimalMesh)
301 /*
302  * void setHexahedreVertexToCells(struct primalMesh * myPrimalMesh)
303  * 
304  * 
305  * Input : myPrimalMesh struct
306  * 
307  * Use :  myPrimalMesh->vertexNumber
308  *        myPrimalMesh->vertexToCellNumber
309  *        myPrimalMesh->cellToVertex
310  *        
311  * 
312  * Allocate : myPrimalMesh->vertexToCells
313  * Produce : myPrimalMesh->vertexToCells
314  * 
315  */
317   
318   clock_t t=startFunction(__FUNCTION__);
319   affiche("\n");
320   
321   connectivity_int * tmp ;
322   
323   myPrimalMesh->vertexToCells = (connectivity_int **) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int *));
324   
325   for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
326     {
327       myPrimalMesh->vertexToCells[i] = (connectivity_int *) malloc(myPrimalMesh->vertexToCellNumber[i] * sizeof(connectivity_int));
328       memset(myPrimalMesh->vertexToCells[i], 0, myPrimalMesh->vertexToCellNumber[i]*sizeof(connectivity_int));
329     }
330   
331   tmp = (connectivity_int *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int));
332   memset(tmp,0, myPrimalMesh->vertexNumber * sizeof(connectivity_int));
333   
334   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
335     {
336       for(connectivity_int i=0;i<myPrimalMesh->cellToVertexNumber[k];i++)
337         {
338           myPrimalMesh->vertexToCells[myPrimalMesh->cellToVertex[k][i]][tmp[myPrimalMesh->cellToVertex[k][i]]] = k;
339           tmp[myPrimalMesh->cellToVertex[k][i]]++;
340         }
341     }
342   
343   
344   
345   
346 #ifdef DEBUG  
347   for(connectivity_int k=0;k<myPrimalMesh->vertexNumber;k++)
348     {
349       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToCells", k);
350       
351       for(connectivity_int ii=0;ii<myPrimalMesh->vertexToCellNumber[k];ii++)
352         {
353           debug_print( "%ld " , myPrimalMesh->vertexToCells[k][ii]);
354         }
355       debug_print(")\n");
356     }
357 #endif
358   free(tmp);
359   
360   endFunction(__FUNCTION__, t);
361   
362   
367 void setHexahedreCellToCells(struct primalMesh * myPrimalMesh)
368 /*
369  * void setHexahedreCellToCells(struct primalMesh * myPrimalMesh)
370  * 
371  * This function searches the cells around a cell using vertices 
372  * to find them.
373  * 
374  * Input : myPrimalMesh struct
375  * 
376  * Use :  myPrimalMesh->M
377  *        myPrimalMesh->N
378  *        myPrimalMesh->P
379  *        myPrimalMesh->cellNumber
380  *        myPrimalMesh->cellToVertexNumber
381  *        myPrimalMesh->cellToVertex
382  *        
383  * 
384  * Allocate : myPrimalMesh->cellToCells
385  *            myPrimalMesh->cellToCellsNumbers
386  * Produce : myPrimalMesh->cellToCells
387  *           myPrimalMesh->cellToCellsNumbers
388  * 
389  */
391   
392   clock_t t=startFunction(__FUNCTION__);
393   affiche("\n");
394   
395   connectivity_int it_cell=0;
396   
397   // vertex number on faces
398   connectivity_int cN1 =  myPrimalMesh->P*myPrimalMesh->N;
399   
400   connectivity_int * vertexList= (connectivity_int *) malloc(SEGMENTVERTEX * sizeof(connectivity_int));
401   
402   myPrimalMesh->cellToCells = (connectivity_int **) malloc( (myPrimalMesh->cellNumber) * sizeof(connectivity_int *));
403   myPrimalMesh->cellToCellsNumbers = (connectivity_int *) malloc( (myPrimalMesh->cellNumber) * sizeof(connectivity_int));
404   
405   connectivity_int * cellList  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024*10);
406   
407   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
408     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
409       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
410         {
411           
412           it_cell     = k + j*myPrimalMesh->P + i*cN1; //cN1 = myPrimalMesh->P*myPrimalMesh->N;
413           
414           if(myPrimalMesh->cellToVertexNumber[it_cell] == HEXAHEDRON)
415             {
416               
417               connectivity_int cellNumbers    = 0 ;
418               connectivity_int vertexSelect   = 0 ;
419               
420               for (connectivity_short vertexi=0;vertexi<HEXAHEDRON;vertexi++)
421                 {
422                   
423                   vertexSelect = myPrimalMesh->cellToVertex[it_cell][vertexi];
424                   
425                   for (connectivity_short celli=0;celli<myPrimalMesh->vertexToCellNumber[vertexSelect];celli++)
426                     {
427                       
428                       if(
429                          findInList(cellList,cellNumbers, myPrimalMesh->vertexToCells[vertexSelect][celli])<0
430                          )
431                         {
432                           
433                           cellList[cellNumbers]=myPrimalMesh->vertexToCells[vertexSelect][celli];
434                           
435                           cellNumbers++;
436                           
437                         }
438                     }
439                   
440                 }                        //                        cellNumbers+=jj;
441               
442               sortListConnectivity(cellList, cellNumbers);
443               
444               myPrimalMesh->cellToCells[it_cell]=(connectivity_int *) malloc(sizeof (connectivity_int)*cellNumbers);
445               
446               for(connectivity_int kk=0;kk<cellNumbers;kk++)
447                 {
448                   myPrimalMesh->cellToCells[it_cell][kk] = cellList[kk];
449                 }
450               
451               myPrimalMesh->cellToCellsNumbers[it_cell] = cellNumbers;
452               
453             }
454         }
455   
456 #ifdef DEBUG  
457   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
458     {
459       debug_print( "%-40s : %ld (%ld) ( ","myPrimalMesh->cellToCells", k, myPrimalMesh->cellToCellsNumbers[k]);
460       
461       for(connectivity_int ii=0;ii<myPrimalMesh->cellToCellsNumbers[k];ii++)
462         {
463           debug_print( "%ld " , myPrimalMesh->cellToCells[k][ii]);
464         }
465       debug_print(")\n");
466     }
467 #endif
468   
469   
470   free(cellList);
471   free(vertexList);
472   
473   endFunction(__FUNCTION__, t);
474   
480 void setHexahedreCellToFacesOwnerNeighbour(struct primalMesh * myPrimalMesh)
481 /*
482  * void setHexahedreCellToFacesOwnerNeighbour(struct primalMesh * myPrimalMesh)
483  * 
484  * XXX
485  * 
486  * Input : myPrimalMesh struct
487  * 
488  * Use :  myPrimalMesh->cellNumber
489  *        
490  * 
491  * Allocate : myPrimalMesh->cellToFacesOwner
492  *            myPrimalMesh->cellToFacesNeighbour
493  *            myPrimalMesh->cellToFacesOwnerNumber
494  *            myPrimalMesh->cellToFacesNeighbourNumber
495  *            
496  *            
497  *            
498  * Produce : myPrimalMesh->faces 
499  *           myPrimalMesh->faceNumber
500  *           myPrimalMesh->cellToFacesOwner
501  *           myPrimalMesh->cellToFacesOwnerNumber
502  *           myPrimalMesh->cellToFacesNeighbour
503  *           myPrimalMesh->cellToFacesNeighbourNumber
504  *          
505  *           
506  * 
507  */
509   
510   clock_t t=startFunction(__FUNCTION__);
511   affiche("\n");
512   
513   connectivity_int position;
514   connectivity_int selectedCell;
515   connectivity_int numberOfVertex;
516   
517   connectivity_int * vertexList  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024*10);
518   connectivity_int * vertexListBase  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
519   connectivity_int * vertexListInverted  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
520   connectivity_int * vertexListInvertedTwo  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
521   
522   connectivity_int * cellTest  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
523   memset(cellTest,0,sizeof (connectivity_int)*myPrimalMesh->cellNumber);
524   
525   char ** cellFaceTest  = (char ** ) malloc(sizeof (char*)*myPrimalMesh->cellNumber);
526   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
527     {
528       cellFaceTest[celli]  = (char * ) malloc(sizeof (char)*HEXAHEDRON_FACES);
529       memset(cellFaceTest[celli],0,sizeof (char)*HEXAHEDRON_FACES);
530     }
531   
532   myPrimalMesh->cellToFacesOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
533   for (connectivity_int ii=0;ii<myPrimalMesh->cellNumber;ii++) {
534       myPrimalMesh->cellToFacesOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
535     }
536   
537   myPrimalMesh->cellToFacesNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
538   for (connectivity_int ii=0;ii<myPrimalMesh->cellNumber;ii++) {
539       myPrimalMesh->cellToFacesNeighbour[ii]  = NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
540     }
541   
542   myPrimalMesh->cellToFacesOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
543   memset(myPrimalMesh->cellToFacesOwnerNumber,0,myPrimalMesh->cellNumber*sizeof (connectivity_short) );
544   
545   myPrimalMesh->cellToFacesNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
546   memset(myPrimalMesh->cellToFacesNeighbourNumber,0,myPrimalMesh->cellNumber*sizeof (connectivity_short) );
547   
548   
549   //    myPrimalMesh->cellToFacesNeighbour = realloc(myPrimalMesh->cellToFacesNeighbour,sizeof (connectivity_int*)*(myPrimalMesh->cellNumber));
550   //    myPrimalMesh->cellToFacesOwner = realloc(myPrimalMesh->cellToFacesOwner,sizeof (connectivity_int*)*(myPrimalMesh->cellNumber));
551   
552   connectivity_int facesNeighbourCount = 0;
553   connectivity_int facesOwnerCount = 0;
554   
555   
556   for(connectivity_int cellAi=0;cellAi<myPrimalMesh->cellNumber;cellAi++)
557     {
558       if(cellAi % 100000 == 0 && cellAi != 0)
559         release_print( "%-40s :  (%ld)\n","CELLS", cellAi);
560       
561       for (connectivity_int faceAi=0;faceAi<HEXAHEDRON_FACES;faceAi++)
562         {
563           
564           for (connectivity_int vertexAi=0;vertexAi<hexaedron_localNodeListNumbers[faceAi];vertexAi++)
565             {
566               vertexListInverted[vertexAi] = myPrimalMesh->cellToVertex[cellAi][ hexaedron_localNodeList[faceAi][hexaedron_localNodeListNumbers[faceAi] -1 - vertexAi] ];
567               vertexListBase[vertexAi] = myPrimalMesh->cellToVertex[cellAi][ hexaedron_localNodeList[faceAi][vertexAi] ];
568             }
569           
570           position = findMinInList(vertexListInverted,hexaedron_localNodeListNumbers[faceAi]);
571           sortList(&vertexListInverted, hexaedron_localNodeListNumbers[faceAi], position);
572           
573           position = findMinInList(vertexListBase,hexaedron_localNodeListNumbers[faceAi]);
574           sortList(&vertexListBase, hexaedron_localNodeListNumbers[faceAi], position);
575           
576           //            if( cellTest[celli]==0)
577           {
578             
579             for(connectivity_int cellBi=0;cellBi<myPrimalMesh->cellToCellsNumbers[cellAi];cellBi++)
580               {
581                 
582                 selectedCell    = myPrimalMesh->cellToCells[cellAi][cellBi];
583                 numberOfVertex  = myPrimalMesh->cellToVertexNumber[selectedCell];
584                 
585                 
586                 for (connectivity_int faceBi=0;faceBi<HEXAHEDRON_FACES;faceBi++)
587                   {
588                     
589                     for (connectivity_int vertexBi=0;vertexBi<hexaedron_localNodeListNumbers[faceBi];vertexBi++)
590                       {
591                         vertexList[vertexBi] = myPrimalMesh->cellToVertex[selectedCell][ hexaedron_localNodeList[faceBi][vertexBi] ];
592                       }
593                     
594                     position = findMinInList(vertexList,hexaedron_localNodeListNumbers[faceBi]);
595                     sortList(&vertexList, hexaedron_localNodeListNumbers[faceBi], position);
596                     
597                     if(memcmp(vertexListBase,vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int)) == 0)
598                       {
599                         
600                         if(cellFaceTest[selectedCell][faceBi]==0)
601                           {
602                             cellFaceTest[selectedCell][faceBi]=1;
603                             
604                             facesOwnerCount = (connectivity_int) myPrimalMesh->cellToFacesOwnerNumber[selectedCell];
605                             
606                             myPrimalMesh->faces = (connectivity_int**) realloc(myPrimalMesh->faces,sizeof (connectivity_int*)*(myPrimalMesh->faceNumber+1));
607                             myPrimalMesh->faces[myPrimalMesh->faceNumber] =(connectivity_int*)  malloc(hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int) );
608                             
609                             memcpy( myPrimalMesh->faces[myPrimalMesh->faceNumber], vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int) );
610                             
611                             myPrimalMesh->cellToFacesOwner[selectedCell] = (connectivity_int *) realloc(myPrimalMesh->cellToFacesOwner[selectedCell], (facesOwnerCount+1)*sizeof (connectivity_int) );
612                             
613                             myPrimalMesh->cellToFacesOwner[selectedCell][facesOwnerCount] = myPrimalMesh->faceNumber;
614                             myPrimalMesh->cellToFacesOwnerNumber[selectedCell]++;
615                             
616                             myPrimalMesh->faceNumber++;
617                             
618                           }
619                       }
620                     else
621                       {
622                         if(memcmp(vertexListInverted,vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int)) == 0)
623                           {
624                             if(cellFaceTest[selectedCell][faceBi]==0)
625                               {
626                                 cellFaceTest[selectedCell][faceBi]=-1;
627                                 
628                                 facesNeighbourCount = (connectivity_int) myPrimalMesh->cellToFacesNeighbourNumber[selectedCell];
629                                 
630                                 //                                    debug_print("%ld\n",facesNeighbourCount);
631                                 
632                                 if(myPrimalMesh->cellToFacesNeighbour[selectedCell]==NULL)
633                                   myPrimalMesh->cellToFacesNeighbour[selectedCell] = (connectivity_int*) malloc((facesNeighbourCount+1)*sizeof (connectivity_int) );
634                                 else
635                                   myPrimalMesh->cellToFacesNeighbour[selectedCell] = (connectivity_int*) realloc(myPrimalMesh->cellToFacesNeighbour[selectedCell], ((facesNeighbourCount+1)) *sizeof (connectivity_int) );
636                                 
637                                 connectivity_short testFaceCi = 0;
638                                 
639                                 //                                for( connectivity_int faceCi= 0; (faceCi<myPrimalMesh->faceNumber) && (testFaceCi == 0) ; faceCi++ )
640                                 for( connectivity_int faceCi= (connectivity_int) myPrimalMesh->faceNumber-1; (faceCi>0) && (testFaceCi == 0) ; faceCi-- )
641                                   {
642                                     
643                                     if(memcmp(vertexListBase,myPrimalMesh->faces[faceCi],QUAD*sizeof (connectivity_int))==0)
644                                       {
645                                         testFaceCi = 1;
646                                         myPrimalMesh->cellToFacesNeighbour[selectedCell][facesNeighbourCount]=faceCi;
647                                         myPrimalMesh->cellToFacesNeighbourNumber[selectedCell]++;
648                                       }
649                                   }
650                                 
651                               }
652                           }
653                       }
654                   }
655                 
656               }
657           }
658           
659           
660           
661           
662           
663           
664         }
665     }
666   
667   free(vertexList);
668   free(vertexListBase);
669   free(vertexListInverted);
670   free(cellTest);
671   free(vertexListInvertedTwo);
672   
673   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
674     free(cellFaceTest[celli]);
675   free(cellFaceTest);
676   
677 #ifdef DEBUG
678   
679   for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
680     {
681       debug_print( "%-40s : %ld ( ","myPrimalMesh->faces", k);
682       
683       for(connectivity_int ii=0;ii<QUAD;ii++)
684         {
685           debug_print( "%ld " , myPrimalMesh->faces[k][ii]);
686         }
687       debug_print(")\n");
688     }
689   
690   
691   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
692     {
693       debug_print( "%-40s : %ld ( ","myPrimalMesh->cellToFacesOwner", celli);
694       
695       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
696         {
697           debug_print( "%ld " , myPrimalMesh->cellToFacesOwner[celli][faceAi]);
698         }
699       debug_print(")\n");
700     }
701   
702   
703   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
704     {
705       debug_print( "%-40s : %ld ( ","myPrimalMesh->cellToFacesNeighbour", k);
706       
707       for(connectivity_int ii=0;ii<myPrimalMesh->cellToFacesNeighbourNumber[k];ii++)
708         {
709           debug_print( "%ld " , myPrimalMesh->cellToFacesNeighbour[k][ii] );
710         }
711       debug_print(")\n");
712     }
713   
714 #endif
715   
716   endFunction(__FUNCTION__, t);
717   
720 void setHexahedreFaceToCells(struct primalMesh * myPrimalMesh)
721 /*
722  * void setHexahedreFaceToCells(struct primalMesh * myPrimalMesh)
723  * 
724  * 
725  * Input : myPrimalMesh struct
726  * 
727  * Use :  myPrimalMesh->faceNumber
728  *        myPrimalMesh->cellNumber
729  *        myPrimalMesh->cellToFacesOwnerNumber
730  *        myPrimalMesh->cellToFacesNeighbourNumber
731  *        
732  * 
733  * Allocate : 
734  *            myPrimalMesh->faceToCells
735  *            myPrimalMesh->faceToCellsNumber
736  * Produce : 
737  *            myPrimalMesh->faceToCells
738  *            myPrimalMesh->faceToCellsNumber
739  *           
740  * 
741  */
743   
744   clock_t t=startFunction(__FUNCTION__);
745   affiche("\n");
746   
747   connectivity_int selectedFace = 0;
748   
749   
750   myPrimalMesh->faceToCells  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->faceNumber);
751   for (connectivity_int ii=0;ii<myPrimalMesh->faceNumber;ii++) {
752       myPrimalMesh->faceToCells[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
753     }
754   
755   myPrimalMesh->faceToCellsNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->faceNumber);
756   memset(myPrimalMesh->faceToCellsNumber,0,myPrimalMesh->faceNumber*sizeof (connectivity_short) );
757   
758   
759   
760   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
761     {
762       
763       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
764         {
765           
766           selectedFace = myPrimalMesh->cellToFacesOwner[celli][faceAi];
767           
768           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
769           
770           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
771           myPrimalMesh->faceToCellsNumber[selectedFace]++;
772           
773         }
774       
775       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesNeighbourNumber[celli];faceAi++)
776         {
777           
778           selectedFace = myPrimalMesh->cellToFacesNeighbour[celli][faceAi];
779           
780           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
781           
782           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
783           myPrimalMesh->faceToCellsNumber[selectedFace]++;
784           
785           
786         }
787       
788       
789     }
790   
791   
792 #ifdef DEBUG
793   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
794     {
795       debug_print( "%-40s : %ld ( ","myPrimalMesh->faceToCells", faceAi);
796       
797       for(connectivity_int celli=0;celli<myPrimalMesh->faceToCellsNumber[faceAi];celli++)
798         {
799           debug_print( "%ld " , myPrimalMesh->faceToCells[faceAi][celli] );
800         }
801       debug_print(")\n");
802     }
803   
804 #endif
805   
806   
807   endFunction(__FUNCTION__, t);
808   
813 void setHexahedreSegments(struct primalMesh * myPrimalMesh)
814 /*
815  * void setHexahedreSegments(struct primalMesh * myPrimalMesh)
816  * 
817  * 
818  * Input : myPrimalMesh struct
819  * 
820  * Use :  myPrimalMesh->cellNumber
821  *        myPrimalMesh->cellToVertex
822  *      
823  * 
824  *        hexaedron_localSegmentList
825  *        
826  * 
827  * Allocate : myPrimalMesh->segments
828  *            myPrimalMesh->segmentNumber
829  *            myPrimalMesh->cellToSegmentOwner
830  *            myPrimalMesh->cellToSegmentOwnerNumber
831  *            myPrimalMesh->cellToSegmentNeighbour
832  *            myPrimalMesh->cellToSegmentNeighbourNumber
833  * 
834  * Produce :  myPrimalMesh->segments
835  *            myPrimalMesh->segmentNumber
836  *            myPrimalMesh->cellToSegmentOwner
837  *            myPrimalMesh->cellToSegmentOwnerNumber
838  *            myPrimalMesh->cellToSegmentNeighbour
839  *            myPrimalMesh->cellToSegmentNeighbourNumber
840  * 
841  */
843   
844   clock_t t=startFunction(__FUNCTION__);
845   affiche("\n");
846   
847   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
848   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
849   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
850   
851   connectivity_int testSegment = 0;
852   
853   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(0));
854   
855   /// MODIF
856   
857   char ** cellSegmentTest  = (char ** ) malloc(sizeof (char*)*myPrimalMesh->cellNumber);
858   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
859     {
860       cellSegmentTest[celli]  = (char * ) malloc(sizeof (char)*HEXAHEDRON_SEGMENTS);
861       memset(cellSegmentTest[celli],0,sizeof (char)*HEXAHEDRON_SEGMENTS);
862     }
863   
864   myPrimalMesh->cellToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
865   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
866     {
867       myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
868     }
869   
870   myPrimalMesh->cellToSegmentOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
871   memset(myPrimalMesh->cellToSegmentOwnerNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_short));
872   
873   myPrimalMesh->cellToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
874   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
875     {
876       myPrimalMesh->cellToSegmentNeighbour[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
877     }
878   
879   myPrimalMesh->cellToSegmentNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
880   memset(myPrimalMesh->cellToSegmentNeighbourNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_short));
881   
882   /// MODIF
883   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
884     {
885       for (connectivity_int segmentAi=0;segmentAi<HEXAHEDRON_SEGMENTS;segmentAi++)
886         { 
887           if(cellSegmentTest[celli][segmentAi] ==0)
888             {
889               connectivity_int localPointA1 = hexaedron_localSegmentList[segmentAi][PT1];
890               connectivity_int localPointA2 = hexaedron_localSegmentList[segmentAi][PT2];
891               
892               cellSegmentTest[celli][segmentAi] = 1;
893               
894               myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->segmentNumber+1));
895               myPrimalMesh->segments[myPrimalMesh->segmentNumber]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
896               
897               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT1] = myPrimalMesh->cellToVertex[celli][localPointA1];
898               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT2] = myPrimalMesh->cellToVertex[celli][localPointA2];
899               
900               myPrimalMesh->segmentNumber++;
901               
902               myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentOwner[celli], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentOwnerNumber[celli]+1));
903               
904               myPrimalMesh->cellToSegmentOwner[celli][myPrimalMesh->cellToSegmentOwnerNumber[celli]] = myPrimalMesh->segmentNumber -1;
905               
906               myPrimalMesh->cellToSegmentOwnerNumber[celli]++;
907               
908               if(myPrimalMesh->segmentNumber % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && myPrimalMesh->segmentNumber != 0)
909                 release_print( "%-40s :  (%ld)\n","SEGMENTS", myPrimalMesh->segmentNumber);
910               
911               for (connectivity_int cellAi=0;cellAi<myPrimalMesh->cellToCellsNumbers[celli];cellAi++)
912                 {      
913                   connectivity_int selectedCell = myPrimalMesh->cellToCells[celli][cellAi];
914                   
915                   
916                   if(celli!=selectedCell) {
917                       
918                       //                      affiche("selectedCell %ld:%ld\n",celli,selectedCell);
919                       
920                       for (connectivity_int segmentBi=0;segmentBi<HEXAHEDRON_SEGMENTS;segmentBi++)
921                         { 
922                           if(cellSegmentTest[selectedCell][segmentBi] ==0)
923                             {
924                               
925                               connectivity_int localPointB1 = hexaedron_localSegmentList[segmentBi][PT1];
926                               connectivity_int localPointB2 = hexaedron_localSegmentList[segmentBi][PT2];
927                               
928                               if(
929                                  myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
930                                  &&
931                                  myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
932                                  )
933                                 {
934                                   //                                  affiche("A: %ld %ld\n",segmentAi,segmentBi);
935                                   cellSegmentTest[selectedCell][segmentBi] = 1;
936                                   
937                                 }
938                               else 
939                                 if(
940                                    myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
941                                    &&
942                                    myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
943                                    )
944                                   {
945                                     //                                    affiche("B: %ld %ld\n",segmentAi,segmentBi);
946                                     cellSegmentTest[selectedCell][segmentBi] = -1;
947                                   }
948                               
949                               
950                               if(cellSegmentTest[selectedCell][segmentBi] !=0)
951                                 {
952                                   
953                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentNeighbour[selectedCell], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]+1));
954                                   
955                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell][myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]] = myPrimalMesh->segmentNumber -1;
956                                   
957                                   myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]++;
958                                   
959                                 }
960                               
961                             }
962                           
963                         }
964                     }  
965                   else 
966                     {
967                       //                      affiche("Noeud voisin \n");
968                     }
969                   
970                 }
971             }
972         }
973       
974     }
975   
976 #ifdef DEBUG  
977   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
978     {
979       debug_print( "%-40s : %ld ( ","cellToSegmentOwner", celli);
980       
981       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli];segmenti++)
982         {
983           debug_print( "%ld " , myPrimalMesh->cellToSegmentOwner[celli][segmenti] );
984         }
985       debug_print(")\n");
986     }  
987   
988   
989   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
990     {
991       debug_print( "%-40s : %ld ( ","cellToSegmentNeighbour", celli);
992       
993       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli];segmenti++)
994         {
995           debug_print( "%ld " , myPrimalMesh->cellToSegmentNeighbour[celli][segmenti] );
996         }
997       debug_print(")\n");
998     }  
999   
1000   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1001     {
1002       debug_print( "%-40s : %ld ( ","myPrimalMesh->segments", segmentAi);
1003       
1004       for(connectivity_int vertexAi=0;vertexAi<SEGMENTVERTEX;vertexAi++)
1005         {
1006           debug_print( "%ld " , myPrimalMesh->segments[segmentAi][vertexAi] );
1007         }
1008       debug_print(")\n");
1009     }
1010   
1011 #endif
1012   
1013   free(copyFace);
1014   free(vertexList);
1015   free(vertexListInverted);
1016   
1017   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1018     free(cellSegmentTest[celli]);
1019   free(cellSegmentTest);  
1020   
1021   
1022   
1023   
1024   endFunction(__FUNCTION__, t);
1025   
1030 void setHexahedreSegmentsCentres(struct primalMesh * myPrimalMesh)
1031 /*
1032  * void setHexahedreSegmentsCentres(struct primalMesh * myPrimalMesh)
1033  * 
1034  * 
1035  * Input : myPrimalMesh struct
1036  * 
1037  * Use :  myPrimalMesh->segmentNumber
1038  *        
1039  * 
1040  * Allocate : myPrimalMesh->segmentCentres
1041  *            
1042  * Produce :  myPrimalMesh->segmentCentres
1043  * 
1044  */
1046   
1047   clock_t t=startFunction(__FUNCTION__);
1048   affiche("\n");
1049   
1050   myPrimalMesh->segmentCentres = (dataType ** ) malloc(sizeof (dataType*)*myPrimalMesh->segmentNumber);
1051   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1052     {
1053       myPrimalMesh->segmentCentres[segmenti]  = (dataType * ) malloc(sizeof (dataType)*DIM3D);
1054     }
1055   
1056   
1057   for (connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1058     {
1059       
1060       if(segmenti % 100000 == 0 && segmenti != 0)
1061         release_print( "%-40s :  (%ld)\n","SEGMENT", segmenti);
1062       
1063       myPrimalMesh->segmentCentres[segmenti][0] = (
1064             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][0]
1065           +
1066           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][0]
1067           )/2.0;
1068       myPrimalMesh->segmentCentres[segmenti][1] = (
1069             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][1]
1070           +
1071           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][1]
1072           )/2.0;
1073       myPrimalMesh->segmentCentres[segmenti][2] = (
1074             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][2]
1075           +
1076           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][2]
1077           )/2.0;
1078     }
1079   
1080 #ifdef DEBUG
1081   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1082     {
1083       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentCentres", segmentAi);
1084       
1085       for(connectivity_int vertexAi=0;vertexAi<DIM3D;vertexAi++)
1086         {
1087           debug_print( "%lf " , myPrimalMesh->segmentCentres[segmentAi][vertexAi] );
1088         }
1089       debug_print(")\n");
1090     }
1091   
1092 #endif  
1093   
1094   endFunction(__FUNCTION__, t);
1095   
1102 void setHexahedreVertexToSegments(struct primalMesh * myPrimalMesh)
1103 /*
1104  * void setHexahedreVertexToSegments(struct primalMesh * myPrimalMesh)
1105  * 
1106  * 
1107  * Input : myPrimalMesh struct
1108  * 
1109  * Use :  
1110  *        
1111  * 
1112  * Allocate : 
1113  * Produce : 
1114  * 
1115  */
1117   
1118   clock_t t=startFunction(__FUNCTION__);
1119   affiche("\n");
1120   
1121   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1122   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1123   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1124   
1125   //connectivity_int testSegment = 0;
1126   
1127   myPrimalMesh->vertexToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
1128   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
1129       myPrimalMesh->vertexToSegmentOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1130     }
1131   
1132   myPrimalMesh->vertexToSegmentOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->vertexNumber);
1133   memset(myPrimalMesh->vertexToSegmentOwnerNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
1134   
1135   myPrimalMesh->vertexToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
1136   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
1137       myPrimalMesh->vertexToSegmentNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1138     }
1139   
1140   myPrimalMesh->vertexToSegmentNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->vertexNumber);
1141   memset(myPrimalMesh->vertexToSegmentNeighbourNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
1142   
1143   
1144   connectivity_int pointOwner = 0;
1145   connectivity_int pointNeighbour = 0;
1146   
1147   for (connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1148     {
1149       if(segmentAi % (connectivity_int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && segmentAi!=0)
1150         release_print( "%-40s :  (%ld)\n","SEGMENTS", segmentAi);
1151       
1152       pointOwner = myPrimalMesh->segments[segmentAi][1];
1153       
1154       myPrimalMesh->vertexToSegmentOwner[pointOwner] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentOwner[pointOwner],(myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]+1)*sizeof (connectivity_int));
1155       myPrimalMesh->vertexToSegmentOwner[pointOwner][myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]] = segmentAi;
1156       myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]++;
1157       
1158       pointNeighbour = myPrimalMesh->segments[segmentAi][0];
1159       
1160       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour],(myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]+1)*sizeof (connectivity_int));
1161       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour][myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]] = segmentAi;
1162       myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]++;
1163       
1164     }
1165   
1166 #ifdef DEBUG
1167   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
1168     {
1169       
1170       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentOwner", vertexAi);
1171       
1172       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentOwnerNumber[vertexAi];segmentAi++)
1173         {
1174           debug_print( "%ld " , myPrimalMesh->vertexToSegmentOwner[vertexAi][segmentAi] );
1175         }
1176       
1177       debug_print(")\n");
1178       
1179     }
1180   
1181   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
1182     {
1183       
1184       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentNeighbour", vertexAi);
1185       
1186       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentNeighbourNumber[vertexAi];segmentAi++)
1187         {
1188           debug_print( "%ld " , myPrimalMesh->vertexToSegmentNeighbour[vertexAi][segmentAi] );
1189         }
1190       
1191       debug_print(")\n");
1192       
1193     }
1194 #endif
1195   
1196   free(copyFace);
1197   free(vertexList);
1198   free(vertexListInverted);
1199   
1200   
1201   endFunction(__FUNCTION__, t);
1202   
1206 void setHexahedreSegmentToFaces(struct primalMesh * myPrimalMesh)
1207 /*
1208  * void XXX(struct primalMesh * myPrimalMesh)
1209  * 
1210  * 
1211  * Input : myPrimalMesh struct
1212  * 
1213  * Use :  
1214  *        
1215  * 
1216  * Allocate : 
1217  * Produce : 
1218  * 
1219  */
1221   
1222   clock_t t=startFunction(__FUNCTION__);
1223   affiche("\n");
1224   
1225   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1226   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1227   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1228   connectivity_int segmentAi = 0;
1229   connectivity_int testSegment=0;
1230   connectivity_int faceAi =0;
1231   
1232   //connectivity_int testSegment = 0;
1233   
1234   //  affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
1235   
1236   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1237   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1238       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1239     }
1240   
1241   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1242   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1243   
1244   
1245   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1246   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1247       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1248     }
1249   
1250   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1251   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1252   
1253   
1254   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1255     {
1256       if(celli % 100000 == 0 && celli != 0)
1257         release_print( "%-40s :  (%ld)\n","CELL", celli);
1258       
1259       
1260       
1261       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesOwnerNumber[celli];facei++)
1262         {
1263           faceAi = myPrimalMesh->cellToFacesOwner[celli][facei];
1264           
1265           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1266           copyFace[QUAD] = copyFace[0];  
1267           
1268           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1269             {
1270               testSegment=0;
1271               
1272               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
1273                 {
1274                   
1275                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
1276                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1277                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1278                   
1279                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1280                     {
1281                       
1282                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1283                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1284                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1285                       testSegment=1;
1286                     }
1287                   else
1288                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1289                       {
1290                         
1291                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1292                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1293                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1294                         testSegment=2;
1295                       }
1296                   
1297                 }
1298               
1299               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
1300                 {
1301                   
1302                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
1303                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1304                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1305                   
1306                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1307                     {
1308                       
1309                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1310                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1311                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1312                       testSegment=1;
1313                     }
1314                   else
1315                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1316                       {
1317                         
1318                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1319                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1320                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1321                         testSegment=2;
1322                       }
1323                   
1324                   
1325                   
1326                 }
1327             }
1328         }
1329       
1330       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesNeighbourNumber[celli];facei++)
1331         {
1332           faceAi = myPrimalMesh->cellToFacesNeighbour[celli][facei];
1333           
1334           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1335           copyFace[QUAD] = copyFace[0];  
1336           
1337           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1338             {
1339               testSegment=0;
1340               
1341               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
1342                 {
1343                   
1344                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
1345                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1346                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1347                   
1348                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1349                     {
1350                       
1351                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1352                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1353                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1354                       testSegment=1;
1355                     }
1356                   else
1357                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1358                       {
1359                         
1360                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1361                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1362                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1363                         testSegment=2;
1364                       }
1365                 }
1366               
1367               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
1368                 {
1369                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
1370                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1371                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1372                   
1373                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1374                     {
1375                       
1376                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1377                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1378                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1379                       testSegment=1;
1380                     }
1381                   else
1382                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1383                       {
1384                         
1385                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1386                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1387                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1388                         testSegment=2;
1389                       }
1390                   
1391                   
1392                 }
1393             }
1394         }
1395       
1396       
1397       
1398       
1399       
1400     }
1401   
1402 #ifdef DEBUG
1403   
1404   
1405   
1406   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1407     {
1408       
1409       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
1410       
1411       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
1412         {
1413           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
1414         }
1415       
1416       debug_print(")\n");
1417       
1418     }
1419   
1420   
1421   
1422   
1423   
1424   
1425   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1426     {
1427       
1428       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
1429       
1430       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
1431         {
1432           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
1433         }
1434       
1435       debug_print(")\n");
1436       
1437     }
1438   
1439 #endif
1440   
1441   free(copyFace);
1442   free(vertexList);
1443   free(vertexListInverted);
1444   
1445   
1446   endFunction(__FUNCTION__, t);
1447   
1452 void setHexahedreSegmentToFaces2(struct primalMesh * myPrimalMesh)
1453 /*
1454  * void XXX(struct primalMesh * myPrimalMesh)
1455  * 
1456  * 
1457  * Input : myPrimalMesh struct
1458  * 
1459  * Use :  
1460  *        
1461  * 
1462  * Allocate : 
1463  * Produce : 
1464  * 
1465  */
1467   
1468   clock_t t=startFunction(__FUNCTION__);
1469   affiche("\n");
1470   
1471   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1472   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1473   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1474   
1475   //connectivity_int testSegment = 0;
1476   
1477   affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
1478   
1479   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1480   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1481       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1482     }
1483   
1484   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1485   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1486   
1487   
1488   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1489   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1490       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1491     }
1492   
1493   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1494   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1495   
1496   
1497   for (connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1498     {
1499       if(faceAi % ((int)(floor(myPrimalMesh->cellNumber*10/NBAFFICHE))) == 0 && faceAi!=0)
1500         release_print( "%-40s :  (%ld)\n","FACES", faceAi);
1501       
1502       memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1503       copyFace[QUAD] = copyFace[0];
1504       
1505       for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1506         {
1507           connectivity_int testSegment=0;
1508           /*
1509            * AMELIORATION du temps de calcul :
1510            * Augmenter la vitesse de recherche
1511            * en incluant une fonction faceToCells puis cellToSegments pour réduire l'impact de la recherche
1512            * de segmentToFaceOwner et segmentToFaceNeighbour
1513            */
1514           for (connectivity_int segmentAi=0 ; segmentAi<myPrimalMesh->segmentNumber && testSegment==0 ;segmentAi++)
1515             {
1516               
1517               vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1518               vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1519               
1520               if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1521                 {
1522                   
1523                   myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1524                   myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1525                   myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1526                   testSegment=1;
1527                 }
1528               else
1529                 if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1530                   {
1531                     
1532                     myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1533                     myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1534                     myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1535                     testSegment=2;
1536                   }
1537             }
1538         }
1539     }
1540   
1541   
1542   
1543 #ifdef DEBUG
1544   
1545   
1546   
1547   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1548     {
1549       
1550       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
1551       
1552       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
1553         {
1554           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
1555         }
1556       
1557       debug_print(")\n");
1558       
1559     }
1560   
1561   
1562   
1563   
1564   
1565   
1566   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1567     {
1568       
1569       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
1570       
1571       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
1572         {
1573           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
1574         }
1575       
1576       debug_print(")\n");
1577       
1578     }
1579   
1580 #endif
1581   
1582   free(copyFace);
1583   free(vertexList);
1584   free(vertexListInverted);
1585   
1586   
1587   endFunction(__FUNCTION__, t);
1588   
1594 void setHexahedreSortFaceToVertex(struct primalMesh * myPrimalMesh)
1595 /*
1596  * void XXX(struct primalMesh * myPrimalMesh)
1597  * 
1598  * 
1599  * Input : myPrimalMesh struct
1600  * 
1601  * Use :  
1602  *        
1603  * 
1604  * Allocate : 
1605  * Produce : 
1606  * 
1607  */
1609   //////////////////// TRIE LES NOEUDS
1610   
1611   clock_t t=startFunction(__FUNCTION__);
1612   
1613   //    connectivity_int actualPoint;
1614   //    connectivity_int nextPoint;
1615   
1616   // on trie les noeund dans les faces
1617   for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
1618     {
1619       connectivity_int              position=0;
1620       
1621       
1622       for(connectivity_int ii=1;ii<myPrimalMesh->faceToVertexNumber[k];ii++)
1623         {
1624           if(myPrimalMesh->faceToVertex[k][position]>myPrimalMesh->faceToVertex[k][ii])
1625             {
1626               position=ii;
1627             }
1628         }
1629       
1630       for(connectivity_int ii=0;ii<position;ii++)
1631         rotate(myPrimalMesh->faceToVertex[k], myPrimalMesh->faceToVertexNumber[k]);
1632       
1633       
1634     }
1635   
1636   
1637   endFunction(__FUNCTION__, t);
1638   
1642 void setHexahedreAllocVertexToSegment(struct primalMesh * myPrimalMesh)
1644   
1645   
1646   clock_t t=startFunction(__FUNCTION__);
1647   
1648   
1649   
1650   myPrimalMesh->vertexToSegments = (connectivity_int **) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int*));
1651   for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
1652     {
1653       myPrimalMesh->vertexToSegments[i] = (connectivity_int *) malloc( myPrimalMesh->vertexToSegmentNumber[i] * sizeof(connectivity_int));
1654       memset(myPrimalMesh->vertexToSegments[i], 0.0, QUAD*sizeof(connectivity_int));
1655     }
1656   
1657   
1658 #ifdef DEBUG  
1659   for(connectivity_int i=0;i<myPrimalMesh->vertexNumber;i++)
1660     {
1661       
1662       debug_print("%-40s : %ld (%d)\n",
1663                   "myPrimalMesh->vertexToSegmentNumber",
1664                   i,
1665                   myPrimalMesh->vertexToSegmentNumber[i]
1666                   );
1667     }
1668 #endif
1669   endFunction(__FUNCTION__, t);
1672 void setHexahedreFaceCentersAreas(struct primalMesh * myPrimalMesh)
1673 /*
1674  * void XXX(struct primalMesh * myPrimalMesh)
1675  * 
1676  * 
1677  * Input : myPrimalMesh struct
1678  * 
1679  * Use :  
1680  *        
1681  * 
1682  * Allocate : 
1683  * Produce : 
1684  * 
1685  */
1687   
1688   clock_t t=startFunction(__FUNCTION__);
1689   
1690   connectivity_int actualPoint;
1691   connectivity_int nextPoint;
1692   
1693   dataType *fC;
1694   
1695   
1696   fC =(dataType *) calloc(DIM3D, sizeof(dataType));
1697   
1698   myPrimalMesh->faceCentres = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType *));
1699   for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
1700     {
1701       myPrimalMesh->faceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1702       memset(myPrimalMesh->faceCentres[i], 0.0, DIM3D*sizeof(dataType));
1703     }
1704   
1705   myPrimalMesh->faceAreas = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType*));
1706   for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
1707     {
1708       myPrimalMesh->faceAreas[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1709       memset(myPrimalMesh->faceAreas[i], 0.0, DIM3D*sizeof(dataType));
1710     }
1711   
1712   
1713   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1714     {
1715       
1716       zeroVector(&fC);
1717       
1718       for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1719         {
1720           actualPoint = myPrimalMesh->faces[faceAi][vertexAi];
1721           
1722           sumTwoVectors(&fC,fC,myPrimalMesh->vertex[actualPoint]);
1723           
1724         }
1725       
1726       
1727       scalarDivVector(&(myPrimalMesh->faceCentres[faceAi]), QUAD,fC);
1728       
1729 #ifdef DEBUG      
1730       debug_print( "%-40s : %ld %d ( ","myPrimalMesh->faceCentres Estimated", faceAi, QUAD);
1731       
1732       for(connectivity_int ii=0;ii<DIM3D;ii++)
1733         {
1734           debug_print( "%lf " , myPrimalMesh->faceCentres[faceAi][ii]);
1735         }
1736 #endif
1737       debug_print(")\n");
1738       if(QUAD==3) {
1739           
1740           // Triangle
1741           affiche("NON IMPLEMENTED ...");
1742           exit(1);
1743           
1744         }
1745       else
1746         { //https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
1747           dataType *sumN;
1748           dataType sumA       = ZEROSCALAR;
1749           dataType *sumAc;
1750           
1751           dataType *c;
1752           dataType a          = ZEROSCALAR;
1753           dataType * n       ;
1754           
1755           dataType *result1;
1756           dataType *result2;
1757           
1758           
1759           result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
1760           sumN =(dataType *) calloc(DIM3D, sizeof(dataType));
1761           sumAc =(dataType *) calloc(DIM3D, sizeof(dataType));
1762           c =(dataType *) calloc(DIM3D, sizeof(dataType));
1763           n =(dataType *) calloc(DIM3D, sizeof(dataType));
1764           result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
1765           
1766           
1767           
1768           for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++) {
1769               
1770               actualPoint = myPrimalMesh->faces[faceAi][vertexAi];
1771               nextPoint = myPrimalMesh->faces[faceAi][(vertexAi + 1) % QUAD];
1772               
1773               sumThreeVectors(&c, (myPrimalMesh->vertex[actualPoint]), (myPrimalMesh->vertex[nextPoint]), (myPrimalMesh->faceCentres[faceAi]) );
1774               
1775               subTwoVectors(&result1, (myPrimalMesh->vertex[nextPoint]), (myPrimalMesh->vertex[actualPoint]));
1776               subTwoVectors(&result2, (myPrimalMesh->faceCentres[faceAi]), (myPrimalMesh->vertex[actualPoint]));
1777               
1778               crossProduct(&n, result1, result2);
1779               mag(&(a), n);
1780               
1781               sumTwoVectors(&sumN, sumN, n);
1782               
1783               sumA+=a;
1784               scalarDotVector(&result1, a, c);
1785               sumTwoVectors(&sumAc, sumAc, result1);
1786             }
1787           
1788           if (sumA < SMALL)
1789             {
1790               zeroVector( &(myPrimalMesh->faceAreas[faceAi]) );
1791             }
1792           else
1793             {
1794               scalarDotVector(&(myPrimalMesh->faceCentres[faceAi]), 1.0/(3.0*sumA), sumAc); //correct faceCentres after estimating them
1795               scalarDotVector(&(myPrimalMesh->faceAreas[faceAi]), 0.5, sumN);
1796             }
1797           
1798           free(sumN);
1799           free(sumAc);
1800           free(c);
1801           free(n)       ;
1802           free(result1);
1803           free(result2);
1804           
1805         }
1806       
1807     }
1808   
1809 #ifdef DEBUG
1810   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1811     {
1812       debug_print( "%-40s : %ld %d ( ","myPrimalMesh->faceCentres", faceAi, QUAD);
1813       
1814       for(connectivity_int ii=0;ii<DIM3D;ii++)
1815         {
1816           debug_print( "%lf " , myPrimalMesh->faceCentres[faceAi][ii]);
1817         }
1818       
1819       debug_print(")\n");
1820       
1821       debug_print( "%-40s : %ld ( ","myPrimalMesh->faceAreas", faceAi);
1822       
1823       for(connectivity_int ii=0;ii<DIM3D;ii++)
1824         {
1825           debug_print( "%lf " , myPrimalMesh->faceAreas[faceAi][ii]);
1826         }
1827       
1828       debug_print(")\n");
1829     }
1830 #endif
1831   
1832   free(fC);
1833   endFunction(__FUNCTION__, t);
1834   
1838 void setHexahedreEstimateVolumeCentroid(struct primalMesh * myPrimalMesh)
1839 /*
1840  * void XXX(struct primalMesh * myPrimalMesh)
1841  * 
1842  * 
1843  * Input : myPrimalMesh struct
1844  * 
1845  * Use :  
1846  *        
1847  * 
1848  * Allocate : 
1849  * Produce : 
1850  * 
1851  */
1853   
1854   clock_t t=startFunction(__FUNCTION__);
1855   
1856   
1857   dataType *result;
1858   
1859   // vertex number on faces
1860   
1861   connectivity_int it_faceOwner,it_faceNeighbour;
1862   
1863   result =(dataType *) calloc(DIM3D, sizeof(dataType));
1864   
1865   myPrimalMesh->volumeCentroid = (dataType **) malloc(myPrimalMesh->cellNumber * sizeof(dataType*));
1866   for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
1867     {
1868       myPrimalMesh->volumeCentroid[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1869       memset(myPrimalMesh->volumeCentroid[i], 0.0, DIM3D*sizeof(dataType));
1870     }
1871   
1872   //    myPrimalMesh->cellToFacesNumber = (connectivity_short*) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
1873   
1874   
1875   myPrimalMesh->volume  = (dataType *) malloc(myPrimalMesh->cellNumber * sizeof(dataType));
1876   memset(myPrimalMesh->volume, 0.0, myPrimalMesh->cellNumber*sizeof(dataType));
1877   
1878   
1879   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1880     {
1881       
1882       zeroVector(&result);
1883       
1884       for (connectivity_int ii=0;ii<myPrimalMesh->cellToFacesOwnerNumber[celli];ii++)
1885         {
1886           it_faceOwner = myPrimalMesh->cellToFacesOwner[celli][ii];
1887           
1888           
1889           sumTwoVectors(&result, result, (myPrimalMesh->faceCentres[it_faceOwner]) );
1890           
1891         }
1892       
1893       for (connectivity_int ii=0;ii<myPrimalMesh->cellToFacesNeighbourNumber[celli];ii++)
1894         {
1895           it_faceNeighbour = myPrimalMesh->cellToFacesNeighbour[celli][ii];
1896           
1897           sumTwoVectors(&result, result, (myPrimalMesh->faceCentres[it_faceNeighbour]) );
1898           
1899         }
1900       scalarDivVector(&(myPrimalMesh->volumeCentroid[celli]) , myPrimalMesh->cellToFacesNeighbourNumber[celli]+myPrimalMesh->cellToFacesOwnerNumber[celli], result);
1901       
1902     }
1903   
1904 #ifdef DEBUG  
1905   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1906     {
1907       
1908       
1909       debug_print("%-40s : %ld (%lf %lf %lf)\n",
1910                   "myPrimalMesh->volumeCentroid Estimated",
1911                   celli,
1912                   myPrimalMesh->volumeCentroid[celli][0],
1913           myPrimalMesh->volumeCentroid[celli][1],
1914           myPrimalMesh->volumeCentroid[celli][2]
1915           );
1916     }
1917 #endif  
1918   
1919   free(result);
1920   endFunction(__FUNCTION__, t);
1921   
1925 void setHexahedreVolumeCentroid(struct primalMesh * myPrimalMesh)
1926 /*
1927  * void XXX(struct primalMesh * myPrimalMesh)
1928  * 
1929  * 
1930  * Input : myPrimalMesh struct
1931  * 
1932  * Use :  
1933  *        
1934  * 
1935  * Allocate : 
1936  * Produce : 
1937  * 
1938  */
1940   clock_t t=startFunction(__FUNCTION__);
1941   
1942   //  connectivity_int  it_cell=0;
1943   dataType        *cellCentroid;
1944   dataType        pyr3Vol=0.0;
1945   dataType        *result1;
1946   dataType        *result2;
1947   dataType        *result3;
1948   
1949   result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
1950   result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
1951   result3 =(dataType *) calloc(DIM3D, sizeof(dataType));
1952   cellCentroid  =(dataType *) calloc(DIM3D, sizeof(dataType));
1953   
1954   // vertex number on faces
1955   connectivity_int it_face;
1956   
1957   // Estimate centroids
1958   setHexahedreEstimateVolumeCentroid(myPrimalMesh);
1959   
1960   
1961   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1962     {
1963       
1964       equalVector(&cellCentroid, myPrimalMesh->volumeCentroid[celli]);
1965       
1966       zeroVector(&(myPrimalMesh->volumeCentroid[celli]));
1967       
1968       for (connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
1969         {
1970           it_face = myPrimalMesh->cellToFacesOwner[celli][faceAi];
1971           
1972           subTwoVectors( &result1, myPrimalMesh->faceCentres[it_face], cellCentroid );
1973           
1974           dotProduct(&pyr3Vol,myPrimalMesh->faceAreas[it_face], result1);
1975           
1976           scalarDotVector(&result1, 3.0/4.0, myPrimalMesh->faceCentres[it_face]); //
1977           scalarDotVector(&result2, 1.0/4.0, cellCentroid); //
1978           
1979           sumTwoVectors(&result3,result1,result2);
1980           
1981           scalarDotVector(&result1,pyr3Vol,result3);
1982           
1983           sumTwoVectors(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volumeCentroid[celli],result1);
1984           
1985           (myPrimalMesh->volume[celli])+=pyr3Vol;
1986           
1987         }
1988       
1989       for (connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesNeighbourNumber[celli];faceAi++)
1990         {
1991           it_face = myPrimalMesh->cellToFacesNeighbour[celli][faceAi];
1992           
1993           subTwoVectors( &result1, cellCentroid , myPrimalMesh->faceCentres[it_face]);
1994           
1995           dotProduct(&pyr3Vol,myPrimalMesh->faceAreas[it_face], result1);
1996           
1997           scalarDotVector(&result1, 3.0/4.0, myPrimalMesh->faceCentres[it_face]); //
1998           scalarDotVector(&result2, 1.0/4.0, cellCentroid); //
1999           
2000           sumTwoVectors(&result3,result1,result2);
2001           
2002           scalarDotVector(&result1,pyr3Vol,result3);
2003           
2004           sumTwoVectors(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volumeCentroid[celli],result1);
2005           
2006           (myPrimalMesh->volume[celli])+=pyr3Vol;
2007           
2008         }
2009       
2010       
2011       
2012       dataType fVol = fabs(myPrimalMesh->volume[celli]);
2013       if (fVol > VSMALL)
2014         {
2015           scalarDivVector(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volume[celli],myPrimalMesh->volumeCentroid[celli]);
2016         }
2017       else
2018         {
2019           equalVector(&(myPrimalMesh->volumeCentroid[celli]), cellCentroid);
2020         }
2021       
2022       
2023       myPrimalMesh->volume[celli] *= (1.0/3.0);
2024       
2025 #ifdef DEBUG      
2026       debug_print("%-40s : %ld (%lf) (%lf %lf %lf)\n",
2027                   "myPrimalMesh->volumeCentroid",
2028                   celli,
2029                   myPrimalMesh->volume[celli],
2030                   myPrimalMesh->volumeCentroid[celli][0],
2031           myPrimalMesh->volumeCentroid[celli][1],
2032           myPrimalMesh->volumeCentroid[celli][2]
2033           );
2034 #endif      
2035     }
2036   
2037   free(result1);
2038   free(result2);
2039   free(result3);
2040   
2041   free(cellCentroid);
2042   
2043   endFunction(__FUNCTION__, t);
2044   
2048 void setHexahedreBoundaryFaces(struct primalMesh * myPrimalMesh)
2049 /*
2050  * void XXX(struct primalMesh * myPrimalMesh)
2051  * 
2052  * 
2053  * Input : myPrimalMesh struct
2054  * 
2055  * Use :  
2056  *        
2057  * 
2058  * Allocate : 
2059  * Produce : 
2060  * 
2061  */
2063   clock_t t=startFunction(__FUNCTION__);
2064   
2065   
2066   // identifys boundary faces an put face numbers in faceBoundary
2067   for(connectivity_int facei=0;facei<myPrimalMesh->faceNumber;facei++)
2068     {
2069       if(myPrimalMesh->faceToCellsNumber[facei] == 1)
2070         {
2071           
2072           myPrimalMesh->faceBoundary = (connectivity_int *) realloc( myPrimalMesh->faceBoundary, (myPrimalMesh->faceBoundaryNumber+1) * sizeof(connectivity_int));
2073           
2074           myPrimalMesh->faceBoundary[myPrimalMesh->faceBoundaryNumber] = facei;
2075           
2076           myPrimalMesh->faceBoundaryNumber++;
2077         }
2078       
2079     }  
2080   
2081   connectivity_short * tmpVertex = (connectivity_short *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_short));
2082   memset(tmpVertex, (0), myPrimalMesh->vertexNumber * sizeof(connectivity_short));
2083   connectivity_int boundaryVertexNumber = 0;
2084   
2085   for(connectivity_int facei=0;facei<myPrimalMesh->faceBoundaryNumber;facei++)
2086     {
2087       connectivity_int faceNumber = myPrimalMesh->faceBoundary[facei]; // numero de la face frontière
2088       connectivity_int * facej = myPrimalMesh->faces[ faceNumber ]; // récupere le numero des points de cette face
2089       
2090       for(connectivity_short vertexi=0;vertexi<QUAD;vertexi++)
2091         {      
2092           tmpVertex[facej[vertexi]] ++; // compte pour chaque point le nombre de faces frontières
2093         }
2094     }
2095   
2096   connectivity_int * boundaryVertexGlobal = (connectivity_int *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int));
2097   memset(boundaryVertexGlobal, (connectivity_int)(0), myPrimalMesh->vertexNumber * sizeof(connectivity_int));
2098   
2099   for(connectivity_int vertexi=0;vertexi<myPrimalMesh->vertexNumber;vertexi++)
2100     {
2101       if(tmpVertex[vertexi]!=0) boundaryVertexNumber++; // compte le nombre de points frontière
2102     }  
2103   
2104   connectivity_int * boundaryVertex = (connectivity_int *) malloc(boundaryVertexNumber * sizeof(connectivity_int));
2105   connectivity_int count = 0;
2106   for(connectivity_int vertexi=0;vertexi<myPrimalMesh->vertexNumber;vertexi++)
2107     {
2108       if(tmpVertex[vertexi]!=0) { // si le point est frontière
2109           boundaryVertex[count] = vertexi; // vecteur de taille boundaryVertexNumber qui stocke les numéros de points
2110           boundaryVertexGlobal[vertexi] = count; // stocke dans un vecteur global la position locale du point dans boundaryVertex
2111           count++;
2112         }
2113     }
2114   
2115   
2116   // boundaryVertexToBoundaryFaces
2117   count=0;
2118   connectivity_int * boundaryVertexToBoundaryFacesNumber = (connectivity_int *) malloc(boundaryVertexNumber * sizeof(connectivity_int));
2119   for(connectivity_int vertexi=0;vertexi<boundaryVertexNumber;vertexi++)
2120     {
2121       memset(boundaryVertexToBoundaryFacesNumber, (connectivity_int)(0), sizeof(connectivity_int));
2122     }
2123   
2124   connectivity_int ** boundaryVertexToBoundaryFaces = (connectivity_int **) malloc(boundaryVertexNumber * sizeof(connectivity_int*));
2125   for(connectivity_int vertexi=0;vertexi<myPrimalMesh->vertexNumber;vertexi++)
2126     {
2127       if(tmpVertex[vertexi]!=0) {
2128           boundaryVertexToBoundaryFacesNumber[count] = tmpVertex[vertexi];
2129           boundaryVertexToBoundaryFaces[count] = (connectivity_int*) malloc( boundaryVertexToBoundaryFacesNumber[count] * sizeof(connectivity_int*));
2130           count++;
2131         }  
2132     }
2133   
2135   //  for(connectivity_int facei=0;facei<myPrimalMesh->faceBoundaryNumber;facei++)
2136   //    {
2137   //      connectivity_int faceNumber = myPrimalMesh->faceBoundary[facei];
2138   //      connectivity_int * facej = myPrimalMesh->faces[ faceNumber ];
2139   
2140   //      for(connectivity_short vertexi=0;vertexi<QUAD;vertexi++)
2141   //        {      //      myPrimalMesh->faces[myPrimalMesh->faceBoundary[facei]];
2142   //          connectivity_int vertexGlobal = ;
2143   
2144   //          boundaryVertexToBoundaryFaces[facej[vertexi]][boundaryVertexToBoundaryFacesNumber[facej[vertexi]] ] = faceNumber;
2145   //          boundaryVertexToBoundaryFacesNumber[facej[vertexi]]++;
2146   //        }
2147   //    }
2148   
2149   
2150 #ifdef DEBUG      
2151   for(connectivity_int facei=0;facei<myPrimalMesh->faceBoundaryNumber;facei++)
2152     {
2153       debug_print("%-40s : %ld (%ld)\n",
2154                   "myPrimalMesh->faceBoundary",
2155                   facei,
2156                   myPrimalMesh->faceBoundary[facei]
2157                   
2158                   );
2159     }
2160   
2161   debug_print("%-40s : %ld \n","boundaryVertexNumber",boundaryVertexNumber);
2162   
2163   for(connectivity_int vertexi=0;vertexi<boundaryVertexNumber;vertexi++)
2164     {
2165       debug_print("%-40s : %ld (%d)\n",
2166                   "tmpVertex",
2167                   vertexi,
2168                   tmpVertex[vertexi]
2169                   
2170                   );
2171     }
2172   
2173   
2174   for(connectivity_int vertexi=0;vertexi<myPrimalMesh->vertexNumber;vertexi++)
2175     {
2176       debug_print("%-40s : %ld (%ld)\n",
2177                   "boundaryVertexGlobal",
2178                   vertexi,
2179                   boundaryVertexGlobal[vertexi]
2180                   
2181                   );
2182     }
2183   
2184   for(connectivity_int vertexi=0;vertexi<boundaryVertexNumber;vertexi++)
2185     {
2186       debug_print("%-40s : %ld (%ld)\n",
2187                   "boundaryVertex",
2188                   vertexi,
2189                   boundaryVertex[vertexi]
2190                   
2191                   );
2192     }
2193   
2194   
2195   //  for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
2196   //    {
2197   
2198   //      debug_print( "%-40s : %ld","myDualMesh->internalDualFaces", faceAi);
2199   
2200   
2201   //          for(connectivity_int vertexi=0;vertexi<DIM3D;vertexi++)
2202   //            {
2203   //              debug_print( "%lf " , boundaryVertexToBoundaryFaces[faceAi][pointi] );
2204   //            }
2205   
2206   //      debug_print(" \n");
2207   
2208   //    }
2209 #endif   
2210   
2211   free(tmpVertex);
2212   free(boundaryVertexToBoundaryFaces);
2213   //  boundaryVertexGlobal
2214   //  boundaryVertexToBoundaryFacesNumber
2215   
2216   //      boundaryVertexToBoundaryFaces
2217   
2218   endFunction(__FUNCTION__, t);
2219