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

4aefbf052ee710f7c5bb18cd0c41db4a38f34511
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   
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 XXX(struct primalMesh * myPrimalMesh)
723  * 
724  * 
725  * Input : myPrimalMesh struct
726  * 
727  * Use :  
728  *        
729  * 
730  * Allocate : 
731  * Produce : 
732  * 
733  */
735   
736   clock_t t=startFunction(__FUNCTION__);
737   affiche("\n");
738   
739   connectivity_int selectedFace = 0;
740   
741   
742   myPrimalMesh->faceToCells  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->faceNumber);
743   for (connectivity_int ii=0;ii<myPrimalMesh->faceNumber;ii++) {
744       myPrimalMesh->faceToCells[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
745     }
746   
747   myPrimalMesh->faceToCellsNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->faceNumber);
748   memset(myPrimalMesh->faceToCellsNumber,0,myPrimalMesh->faceNumber*sizeof (connectivity_short) );
749   
750   
751   
752   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
753     {
754       
755       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
756         {
757           //          myPrimalMesh->cellToFacesOwner[celli][faceAi];
758           //          debug_print("O: %ld %ld\n",celli, myPrimalMesh->cellToFacesOwner[celli][faceAi]);
759           
760           selectedFace = myPrimalMesh->cellToFacesOwner[celli][faceAi];
761           
762           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
763           
764           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
765           myPrimalMesh->faceToCellsNumber[selectedFace]++;
766           
767         }
768       
769       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesNeighbourNumber[celli];faceAi++)
770         {
771           //          myPrimalMesh->cellToFacesOwner[celli][faceAi];
772           //          debug_print("O: %ld %ld\n",celli, myPrimalMesh->cellToFacesNeighbour[celli][faceAi]);
773           
774           selectedFace = myPrimalMesh->cellToFacesNeighbour[celli][faceAi];
775           
776           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
777           
778           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
779           myPrimalMesh->faceToCellsNumber[selectedFace]++;
780           
781           
782         }
783       
784       
785     }
786   
787   
788 #ifdef DEBUG
789   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
790     {
791       debug_print( "%-40s : %ld ( ","myPrimalMesh->faceToCells", faceAi);
792       
793       for(connectivity_int celli=0;celli<myPrimalMesh->faceToCellsNumber[faceAi];celli++)
794         {
795           debug_print( "%ld " , myPrimalMesh->faceToCells[faceAi][celli] );
796         }
797       debug_print(")\n");
798     }
799   
800 #endif
801   
802   //  free(copyFace);
803   //  free(vertexList);
804   //  free(vertexListInverted);
805   
806   
807   endFunction(__FUNCTION__, t);
808   
813 void setHexahedreSegments(struct primalMesh * myPrimalMesh)
814 /*
815  * void XXX(struct primalMesh * myPrimalMesh)
816  * 
817  * 
818  * Input : myPrimalMesh struct
819  * 
820  * Use :  
821  *        
822  * 
823  * Allocate : 
824  * Produce : 
825  * 
826  */
828   
829   clock_t t=startFunction(__FUNCTION__);
830   affiche("\n");
831   //  affiche("%s : REFAIRE CETTE FONCTION EN UTILISANT CELLTOFACE PUIS REDUIRE LA RECHERCHE DANS SEGMENTS EN INTRODUISANT PAR EXEMPLE CELLTOSEGMENTS\n",__FUNCTION__);
832   
833   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
834   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
835   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
836   
837   connectivity_int testSegment = 0;
838   
839   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(0));
840   
841   /// MODIF
842   
843   char ** cellSegmentTest  = (char ** ) malloc(sizeof (char*)*myPrimalMesh->cellNumber);
844   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
845     {
846       cellSegmentTest[celli]  = (char * ) malloc(sizeof (char)*HEXAHEDRON_SEGMENTS);
847       memset(cellSegmentTest[celli],0,sizeof (char)*HEXAHEDRON_SEGMENTS);
848     }
849   
850   myPrimalMesh->cellToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
851   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
852     {
853       myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
854     }
855   
856   myPrimalMesh->cellToSegmentOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
857   memset(myPrimalMesh->cellToSegmentOwnerNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_short));
858   
859   myPrimalMesh->cellToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
860   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
861     {
862       myPrimalMesh->cellToSegmentNeighbour[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
863     }
864   
865   myPrimalMesh->cellToSegmentNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->cellNumber);
866   memset(myPrimalMesh->cellToSegmentNeighbourNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_short));
867   
868   /// MODIF
869   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
870     {
871       for (connectivity_int segmentAi=0;segmentAi<HEXAHEDRON_SEGMENTS;segmentAi++)
872         { 
873           if(cellSegmentTest[celli][segmentAi] ==0)
874             {
875               connectivity_int localPointA1 = hexaedron_localSegmentList[segmentAi][PT1];
876               connectivity_int localPointA2 = hexaedron_localSegmentList[segmentAi][PT2];
877               
878               cellSegmentTest[celli][segmentAi] = 1;
879               
880               myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->segmentNumber+1));
881               myPrimalMesh->segments[myPrimalMesh->segmentNumber]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
882               
883               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT1] = myPrimalMesh->cellToVertex[celli][localPointA1];
884               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT2] = myPrimalMesh->cellToVertex[celli][localPointA2];
885               
886               myPrimalMesh->segmentNumber++;
887               
888               myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentOwner[celli], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentOwnerNumber[celli]+1));
889               
890               myPrimalMesh->cellToSegmentOwner[celli][myPrimalMesh->cellToSegmentOwnerNumber[celli]] = myPrimalMesh->segmentNumber -1;
891               
892               myPrimalMesh->cellToSegmentOwnerNumber[celli]++;
893               
894               if(myPrimalMesh->segmentNumber % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && myPrimalMesh->segmentNumber != 0)
895                 release_print( "%-40s :  (%ld)\n","SEGMENTS", myPrimalMesh->segmentNumber);
896               
897               for (connectivity_int cellAi=0;cellAi<myPrimalMesh->cellToCellsNumbers[celli];cellAi++)
898                 {      
899                   connectivity_int selectedCell = myPrimalMesh->cellToCells[celli][cellAi];
900                   
901                   
902                   if(celli!=selectedCell) {
903                       
904                       //                      affiche("selectedCell %ld:%ld\n",celli,selectedCell);
905                       
906                       for (connectivity_int segmentBi=0;segmentBi<HEXAHEDRON_SEGMENTS;segmentBi++)
907                         { 
908                           if(cellSegmentTest[selectedCell][segmentBi] ==0)
909                             {
910                               
911                               connectivity_int localPointB1 = hexaedron_localSegmentList[segmentBi][PT1];
912                               connectivity_int localPointB2 = hexaedron_localSegmentList[segmentBi][PT2];
913                               
914                               if(
915                                  myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
916                                  &&
917                                  myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
918                                  )
919                                 {
920                                   //                                  affiche("A: %ld %ld\n",segmentAi,segmentBi);
921                                   cellSegmentTest[selectedCell][segmentBi] = 1;
922                                   
923                                 }
924                               else 
925                                 if(
926                                    myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
927                                    &&
928                                    myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
929                                    )
930                                   {
931                                     //                                    affiche("B: %ld %ld\n",segmentAi,segmentBi);
932                                     cellSegmentTest[selectedCell][segmentBi] = -1;
933                                   }
934                               
935                               
936                               if(cellSegmentTest[selectedCell][segmentBi] !=0)
937                                 {
938                                   
939                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentNeighbour[selectedCell], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]+1));
940                                   
941                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell][myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]] = myPrimalMesh->segmentNumber -1;
942                                   
943                                   myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]++;
944                                   
945                                 }
946                               
947                             }
948                           
949                         }
950                     }  
951                   else 
952                     {
953                       //                      affiche("Noeud voisin \n");
954                     }
955                   
956                 }
957             }
958         }
959       
960     }
961   
962 #ifdef DEBUG  
963   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
964     {
965       debug_print( "%-40s : %ld ( ","cellToSegmentOwner", celli);
966       
967       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli];segmenti++)
968         {
969           debug_print( "%ld " , myPrimalMesh->cellToSegmentOwner[celli][segmenti] );
970         }
971       debug_print(")\n");
972     }  
973   
974   
975   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
976     {
977       debug_print( "%-40s : %ld ( ","cellToSegmentNeighbour", celli);
978       
979       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli];segmenti++)
980         {
981           debug_print( "%ld " , myPrimalMesh->cellToSegmentNeighbour[celli][segmenti] );
982         }
983       debug_print(")\n");
984     }  
985   
986   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
987     {
988       debug_print( "%-40s : %ld ( ","myPrimalMesh->segments", segmentAi);
989       
990       for(connectivity_int vertexAi=0;vertexAi<SEGMENTVERTEX;vertexAi++)
991         {
992           debug_print( "%ld " , myPrimalMesh->segments[segmentAi][vertexAi] );
993         }
994       debug_print(")\n");
995     }
996   
997 #endif
998   
999   free(copyFace);
1000   free(vertexList);
1001   free(vertexListInverted);
1002   
1003   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1004     free(cellSegmentTest[celli]);
1005   free(cellSegmentTest);  
1006   
1007   
1008   
1009   
1010   endFunction(__FUNCTION__, t);
1011   
1016 void setHexahedreSegmentsCentres(struct primalMesh * myPrimalMesh)
1017 /*
1018  * void XXX(struct primalMesh * myPrimalMesh)
1019  * 
1020  * 
1021  * Input : myPrimalMesh struct
1022  * 
1023  * Use :  
1024  *        
1025  * 
1026  * Allocate : 
1027  * Produce : 
1028  * 
1029  */
1031   
1032   clock_t t=startFunction(__FUNCTION__);
1033   affiche("\n");
1034   
1035   myPrimalMesh->segmentCentres = (dataType ** ) malloc(sizeof (dataType*)*myPrimalMesh->segmentNumber);
1036   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1037     {
1038       myPrimalMesh->segmentCentres[segmenti]  = (dataType * ) malloc(sizeof (dataType)*DIM3D);
1039     }
1040   
1041   
1042   for (connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1043     {
1044       
1045       if(segmenti % 100000 == 0 && segmenti != 0)
1046         release_print( "%-40s :  (%ld)\n","SEGMENT", segmenti);
1047       
1048       myPrimalMesh->segmentCentres[segmenti][0] = (
1049             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][0]
1050           +
1051           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][0]
1052           )/2.0;
1053       myPrimalMesh->segmentCentres[segmenti][1] = (
1054             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][1]
1055           +
1056           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][1]
1057           )/2.0;
1058       myPrimalMesh->segmentCentres[segmenti][2] = (
1059             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][2]
1060           +
1061           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][2]
1062           )/2.0;
1063     }
1064   
1065 #ifdef DEBUG
1066   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1067     {
1068       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentCentres", segmentAi);
1069       
1070       for(connectivity_int vertexAi=0;vertexAi<DIM3D;vertexAi++)
1071         {
1072           debug_print( "%lf " , myPrimalMesh->segmentCentres[segmentAi][vertexAi] );
1073         }
1074       debug_print(")\n");
1075     }
1076   
1077 #endif  
1078   
1079   endFunction(__FUNCTION__, t);
1080   
1085 void setHexahedreSegments2(struct primalMesh * myPrimalMesh)
1086 /*
1087  * void XXX(struct primalMesh * myPrimalMesh)
1088  * 
1089  * 
1090  * Input : myPrimalMesh struct
1091  * 
1092  * Use :  
1093  *        
1094  * 
1095  * Allocate : 
1096  * Produce : 
1097  * 
1098  */
1100   
1101   clock_t t=startFunction(__FUNCTION__);
1102   affiche("\n");
1103   affiche("%s : Nouvelle fonction setHexahedreSegments pour réduire le temps de calcul\n",__FUNCTION__);
1104   
1105   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1106   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1107   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1108   
1109   connectivity_int testSegment = 0;
1110   
1111   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->cellNumber*6));
1112   for (connectivity_int segmentAi=0;segmentAi<(myPrimalMesh->cellNumber*6);segmentAi++)
1113     myPrimalMesh->segments[segmentAi]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1114   
1115   for (connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1116     //  for (long faceAi=(long)(myPrimalMesh->faceNumber-1);faceAi>=0;faceAi--)
1117     {
1118       
1119       if(faceAi % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && faceAi != 0)
1120         release_print( "%-40s :  (%ld)\n","FACE", faceAi);
1121       
1122       
1123       memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1124       copyFace[QUAD] = copyFace[0];
1125       
1126       for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1127         {
1128           
1129           vertexList[0] = copyFace[vertexAi];
1130           vertexList[1] = copyFace[(vertexAi+1) % QUAD];
1131           
1132           vertexListInverted[0] = vertexList[1];
1133           vertexListInverted[1] = vertexList[0];
1134           
1135           testSegment = 0;
1136           
1137           /*
1138            * AMELIORATION du temps de calcul :
1139            * Augmenter la vitesse de recherche
1140            * en incluant une fonction faceToCells puis cellToSegments pour réduire l'impact de la recherche
1141            */
1142           for (connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber  && testSegment == 0;segmentAi++)
1143             {
1144               
1145               if( vertexList[0] == myPrimalMesh->segments[segmentAi][0])
1146                 {
1147                   if( vertexList[1] == myPrimalMesh->segments[segmentAi][1])
1148                     {
1149                       testSegment = 1;
1150                     }
1151                 }
1152               else
1153                 if( vertexListInverted[0] == myPrimalMesh->segments[segmentAi][0])
1154                   {
1155                     if( vertexListInverted[1] == myPrimalMesh->segments[segmentAi][1])
1156                       
1157                       {
1158                         testSegment = 2;
1159                       }
1160                   }
1161             }
1162           
1163           if(testSegment == 0)
1164             {
1165               if(myPrimalMesh->cellNumber*HEXAHEDRON_FACES<=myPrimalMesh->segmentNumber)
1166                 {
1167                   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->segmentNumber+1));
1168                   myPrimalMesh->segments[myPrimalMesh->segmentNumber]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1169                 }
1170               memcpy(myPrimalMesh->segments[myPrimalMesh->segmentNumber], vertexList, sizeof(connectivity_int)*SEGMENTVERTEX);
1171               
1172               myPrimalMesh->segmentNumber++;
1173             }
1174           
1175         }
1176     }
1177   
1178 #ifdef DEBUG
1179   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1180     {
1181       debug_print( "%-40s : %ld ( ","myPrimalMesh->segments", segmentAi);
1182       
1183       for(connectivity_int vertexAi=0;vertexAi<SEGMENTVERTEX;vertexAi++)
1184         {
1185           debug_print( "%ld " , myPrimalMesh->segments[segmentAi][vertexAi] );
1186         }
1187       debug_print(")\n");
1188     }
1189   
1190 #endif
1191   
1192   free(copyFace);
1193   free(vertexList);
1194   free(vertexListInverted);
1195   
1196   
1197   endFunction(__FUNCTION__, t);
1198   
1204 void setHexahedreVertexToSegments(struct primalMesh * myPrimalMesh)
1205 /*
1206  * void XXX(struct primalMesh * myPrimalMesh)
1207  * 
1208  * 
1209  * Input : myPrimalMesh struct
1210  * 
1211  * Use :  
1212  *        
1213  * 
1214  * Allocate : 
1215  * Produce : 
1216  * 
1217  */
1219   
1220   clock_t t=startFunction(__FUNCTION__);
1221   affiche("\n");
1222   
1223   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1224   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1225   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1226   
1227   //connectivity_int testSegment = 0;
1228   
1229   myPrimalMesh->vertexToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
1230   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
1231       myPrimalMesh->vertexToSegmentOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1232     }
1233   
1234   myPrimalMesh->vertexToSegmentOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->vertexNumber);
1235   memset(myPrimalMesh->vertexToSegmentOwnerNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
1236   
1237   myPrimalMesh->vertexToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
1238   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
1239       myPrimalMesh->vertexToSegmentNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1240     }
1241   
1242   myPrimalMesh->vertexToSegmentNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->vertexNumber);
1243   memset(myPrimalMesh->vertexToSegmentNeighbourNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
1244   
1245   
1246   connectivity_int pointOwner = 0;
1247   connectivity_int pointNeighbour = 0;
1248   
1249   for (connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1250     {
1251       if(segmentAi % (connectivity_int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && segmentAi!=0)
1252         release_print( "%-40s :  (%ld)\n","SEGMENTS", segmentAi);
1253       
1254       pointOwner = myPrimalMesh->segments[segmentAi][1];
1255       
1256       myPrimalMesh->vertexToSegmentOwner[pointOwner] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentOwner[pointOwner],(myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]+1)*sizeof (connectivity_int));
1257       myPrimalMesh->vertexToSegmentOwner[pointOwner][myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]] = segmentAi;
1258       myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]++;
1259       
1260       pointNeighbour = myPrimalMesh->segments[segmentAi][0];
1261       
1262       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour],(myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]+1)*sizeof (connectivity_int));
1263       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour][myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]] = segmentAi;
1264       myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]++;
1265       
1266     }
1267   
1268 #ifdef DEBUG
1269   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
1270     {
1271       
1272       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentOwner", vertexAi);
1273       
1274       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentOwnerNumber[vertexAi];segmentAi++)
1275         {
1276           debug_print( "%ld " , myPrimalMesh->vertexToSegmentOwner[vertexAi][segmentAi] );
1277         }
1278       
1279       debug_print(")\n");
1280       
1281     }
1282   
1283   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
1284     {
1285       
1286       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentNeighbour", vertexAi);
1287       
1288       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentNeighbourNumber[vertexAi];segmentAi++)
1289         {
1290           debug_print( "%ld " , myPrimalMesh->vertexToSegmentNeighbour[vertexAi][segmentAi] );
1291         }
1292       
1293       debug_print(")\n");
1294       
1295     }
1296 #endif
1297   
1298   free(copyFace);
1299   free(vertexList);
1300   free(vertexListInverted);
1301   
1302   
1303   endFunction(__FUNCTION__, t);
1304   
1308 void setHexahedreSegmentToFaces(struct primalMesh * myPrimalMesh)
1309 /*
1310  * void XXX(struct primalMesh * myPrimalMesh)
1311  * 
1312  * 
1313  * Input : myPrimalMesh struct
1314  * 
1315  * Use :  
1316  *        
1317  * 
1318  * Allocate : 
1319  * Produce : 
1320  * 
1321  */
1323   
1324   clock_t t=startFunction(__FUNCTION__);
1325   affiche("\n");
1326   
1327   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1328   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1329   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1330   connectivity_int segmentAi = 0;
1331   connectivity_int testSegment=0;
1332   connectivity_int faceAi =0;
1333   
1334   //connectivity_int testSegment = 0;
1335   
1336   //  affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
1337   
1338   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1339   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1340       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1341     }
1342   
1343   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1344   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1345   
1346   
1347   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1348   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1349       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1350     }
1351   
1352   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1353   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1354   
1355   
1356   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1357     {
1358       if(celli % 100000 == 0 && celli != 0)
1359         release_print( "%-40s :  (%ld)\n","CELL", celli);
1360       
1361       
1362       
1363       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesOwnerNumber[celli];facei++)
1364         {
1365           faceAi = myPrimalMesh->cellToFacesOwner[celli][facei];
1366           
1367           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1368           copyFace[QUAD] = copyFace[0];  
1369           
1370           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1371             {
1372               testSegment=0;
1373               
1374               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
1375                 {
1376                   
1377                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
1378                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1379                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1380                   
1381                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1382                     {
1383                       
1384                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1385                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1386                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1387                       testSegment=1;
1388                     }
1389                   else
1390                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1391                       {
1392                         
1393                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1394                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1395                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1396                         testSegment=2;
1397                       }
1398                   
1399                 }
1400               
1401               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
1402                 {
1403                   
1404                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
1405                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1406                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1407                   
1408                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1409                     {
1410                       
1411                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1412                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1413                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1414                       testSegment=1;
1415                     }
1416                   else
1417                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1418                       {
1419                         
1420                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1421                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1422                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1423                         testSegment=2;
1424                       }
1425                   
1426                   
1427                   
1428                 }
1429             }
1430         }
1431       
1432       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesNeighbourNumber[celli];facei++)
1433         {
1434           faceAi = myPrimalMesh->cellToFacesNeighbour[celli][facei];
1435           
1436           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1437           copyFace[QUAD] = copyFace[0];  
1438           
1439           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1440             {
1441               testSegment=0;
1442               
1443               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
1444                 {
1445                   
1446                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
1447                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1448                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1449                   
1450                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1451                     {
1452                       
1453                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1454                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1455                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1456                       testSegment=1;
1457                     }
1458                   else
1459                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1460                       {
1461                         
1462                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1463                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1464                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1465                         testSegment=2;
1466                       }
1467                 }
1468               
1469               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
1470                 {
1471                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
1472                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1473                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1474                   
1475                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1476                     {
1477                       
1478                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1479                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1480                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1481                       testSegment=1;
1482                     }
1483                   else
1484                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1485                       {
1486                         
1487                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1488                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1489                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1490                         testSegment=2;
1491                       }
1492                   
1493                   
1494                 }
1495             }
1496         }
1497       
1498       
1499       
1500       
1501       
1502     }
1503   
1504 #ifdef DEBUG
1505   
1506   
1507   
1508   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1509     {
1510       
1511       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
1512       
1513       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
1514         {
1515           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
1516         }
1517       
1518       debug_print(")\n");
1519       
1520     }
1521   
1522   
1523   
1524   
1525   
1526   
1527   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1528     {
1529       
1530       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
1531       
1532       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
1533         {
1534           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
1535         }
1536       
1537       debug_print(")\n");
1538       
1539     }
1540   
1541 #endif
1542   
1543   free(copyFace);
1544   free(vertexList);
1545   free(vertexListInverted);
1546   
1547   
1548   endFunction(__FUNCTION__, t);
1549   
1554 void setHexahedreSegmentToFaces2(struct primalMesh * myPrimalMesh)
1555 /*
1556  * void XXX(struct primalMesh * myPrimalMesh)
1557  * 
1558  * 
1559  * Input : myPrimalMesh struct
1560  * 
1561  * Use :  
1562  *        
1563  * 
1564  * Allocate : 
1565  * Produce : 
1566  * 
1567  */
1569   
1570   clock_t t=startFunction(__FUNCTION__);
1571   affiche("\n");
1572   
1573   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1574   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1575   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1576   
1577   //connectivity_int testSegment = 0;
1578   
1579   affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
1580   
1581   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1582   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1583       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1584     }
1585   
1586   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1587   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1588   
1589   
1590   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
1591   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
1592       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1593     }
1594   
1595   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_short * ) malloc(sizeof (connectivity_short)*myPrimalMesh->segmentNumber);
1596   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_short));
1597   
1598   
1599   for (connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1600     {
1601       if(faceAi % ((int)(floor(myPrimalMesh->cellNumber*10/NBAFFICHE))) == 0 && faceAi!=0)
1602         release_print( "%-40s :  (%ld)\n","FACES", faceAi);
1603       
1604       memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1605       copyFace[QUAD] = copyFace[0];
1606       
1607       for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1608         {
1609           connectivity_int testSegment=0;
1610           /*
1611            * AMELIORATION du temps de calcul :
1612            * Augmenter la vitesse de recherche
1613            * en incluant une fonction faceToCells puis cellToSegments pour réduire l'impact de la recherche
1614            * de segmentToFaceOwner et segmentToFaceNeighbour
1615            */
1616           for (connectivity_int segmentAi=0 ; segmentAi<myPrimalMesh->segmentNumber && testSegment==0 ;segmentAi++)
1617             {
1618               
1619               vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
1620               vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
1621               
1622               if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1623                 {
1624                   
1625                   myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
1626                   myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
1627                   myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
1628                   testSegment=1;
1629                 }
1630               else
1631                 if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
1632                   {
1633                     
1634                     myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
1635                     myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
1636                     myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
1637                     testSegment=2;
1638                   }
1639             }
1640         }
1641     }
1642   
1643   
1644   
1645 #ifdef DEBUG
1646   
1647   
1648   
1649   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1650     {
1651       
1652       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
1653       
1654       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
1655         {
1656           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
1657         }
1658       
1659       debug_print(")\n");
1660       
1661     }
1662   
1663   
1664   
1665   
1666   
1667   
1668   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1669     {
1670       
1671       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
1672       
1673       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
1674         {
1675           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
1676         }
1677       
1678       debug_print(")\n");
1679       
1680     }
1681   
1682 #endif
1683   
1684   free(copyFace);
1685   free(vertexList);
1686   free(vertexListInverted);
1687   
1688   
1689   endFunction(__FUNCTION__, t);
1690   
1696 void setHexahedreSortFaceToVertex(struct primalMesh * myPrimalMesh)
1697 /*
1698  * void XXX(struct primalMesh * myPrimalMesh)
1699  * 
1700  * 
1701  * Input : myPrimalMesh struct
1702  * 
1703  * Use :  
1704  *        
1705  * 
1706  * Allocate : 
1707  * Produce : 
1708  * 
1709  */
1711   //////////////////// TRIE LES NOEUDS
1712   
1713   clock_t t=startFunction(__FUNCTION__);
1714   
1715   //    connectivity_int actualPoint;
1716   //    connectivity_int nextPoint;
1717   
1718   // on trie les noeund dans les faces
1719   for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
1720     {
1721       connectivity_int              position=0;
1722       
1723       
1724       for(connectivity_int ii=1;ii<myPrimalMesh->faceToVertexNumber[k];ii++)
1725         {
1726           if(myPrimalMesh->faceToVertex[k][position]>myPrimalMesh->faceToVertex[k][ii])
1727             {
1728               position=ii;
1729             }
1730         }
1731       
1732       for(connectivity_int ii=0;ii<position;ii++)
1733         rotate(myPrimalMesh->faceToVertex[k], myPrimalMesh->faceToVertexNumber[k]);
1734       
1735       
1736     }
1737   
1738  
1739   endFunction(__FUNCTION__, t);
1740   
1744 void setHexahedreAllocVertexToSegment(struct primalMesh * myPrimalMesh)
1746   
1747   
1748   clock_t t=startFunction(__FUNCTION__);
1749   
1750   
1751   
1752   myPrimalMesh->vertexToSegments = (connectivity_int **) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int*));
1753   for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
1754     {
1755       myPrimalMesh->vertexToSegments[i] = (connectivity_int *) malloc( myPrimalMesh->vertexToSegmentNumber[i] * sizeof(connectivity_int));
1756       memset(myPrimalMesh->vertexToSegments[i], 0.0, QUAD*sizeof(connectivity_int));
1757     }
1758   
1759   
1760 #ifdef DEBUG  
1761   for(connectivity_int i=0;i<myPrimalMesh->vertexNumber;i++)
1762     {
1763       
1764       debug_print("%-40s : %ld (%d)\n",
1765                   "myPrimalMesh->vertexToSegmentNumber",
1766                   i,
1767                   myPrimalMesh->vertexToSegmentNumber[i]
1768                   );
1769     }
1770 #endif
1771   endFunction(__FUNCTION__, t);
1774 void setHexahedreFaceCentersAreas(struct primalMesh * myPrimalMesh)
1775 /*
1776  * void XXX(struct primalMesh * myPrimalMesh)
1777  * 
1778  * 
1779  * Input : myPrimalMesh struct
1780  * 
1781  * Use :  
1782  *        
1783  * 
1784  * Allocate : 
1785  * Produce : 
1786  * 
1787  */
1789   
1790   clock_t t=startFunction(__FUNCTION__);
1791   
1792   connectivity_int actualPoint;
1793   connectivity_int nextPoint;
1794   
1795   dataType *fC;
1796   
1797   
1798   fC =(dataType *) calloc(DIM3D, sizeof(dataType));
1799   
1800   myPrimalMesh->faceCentres = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType *));
1801   for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
1802     {
1803       myPrimalMesh->faceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1804       memset(myPrimalMesh->faceCentres[i], 0.0, DIM3D*sizeof(dataType));
1805     }
1806   
1807   myPrimalMesh->faceAreas = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType*));
1808   for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
1809     {
1810       myPrimalMesh->faceAreas[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1811       memset(myPrimalMesh->faceAreas[i], 0.0, DIM3D*sizeof(dataType));
1812     }
1813   
1814   
1815   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1816     {
1817       
1818       zeroVector(&fC);
1819       
1820       for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1821         {
1822           actualPoint = myPrimalMesh->faces[faceAi][vertexAi];
1823           
1824           sumTwoVectors(&fC,fC,myPrimalMesh->vertex[actualPoint]);
1825           
1826         }
1827       
1828       
1829       scalarDivVector(&(myPrimalMesh->faceCentres[faceAi]), QUAD,fC);
1830       
1831 #ifdef DEBUG      
1832       debug_print( "%-40s : %ld %d ( ","myPrimalMesh->faceCentres Estimated", faceAi, QUAD);
1833       
1834       for(connectivity_int ii=0;ii<DIM3D;ii++)
1835         {
1836           debug_print( "%lf " , myPrimalMesh->faceCentres[faceAi][ii]);
1837         }
1838 #endif
1839       debug_print(")\n");
1840       if(QUAD==3) {
1841           
1842           // Triangle
1843           affiche("NON IMPLEMENTED ...");
1844           exit(1);
1845           
1846         }
1847       else
1848         { //https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
1849           dataType *sumN;
1850           dataType sumA       = ZEROSCALAR;
1851           dataType *sumAc;
1852           
1853           dataType *c;
1854           dataType a          = ZEROSCALAR;
1855           dataType * n       ;
1856           
1857           dataType *result1;
1858           dataType *result2;
1859           
1860           
1861           result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
1862           sumN =(dataType *) calloc(DIM3D, sizeof(dataType));
1863           sumAc =(dataType *) calloc(DIM3D, sizeof(dataType));
1864           c =(dataType *) calloc(DIM3D, sizeof(dataType));
1865           n =(dataType *) calloc(DIM3D, sizeof(dataType));
1866           result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
1867           
1868           
1869           
1870           for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++) {
1871               
1872               actualPoint = myPrimalMesh->faces[faceAi][vertexAi];
1873               nextPoint = myPrimalMesh->faces[faceAi][(vertexAi + 1) % QUAD];
1874               
1875               sumThreeVectors(&c, (myPrimalMesh->vertex[actualPoint]), (myPrimalMesh->vertex[nextPoint]), (myPrimalMesh->faceCentres[faceAi]) );
1876               
1877               subTwoVectors(&result1, (myPrimalMesh->vertex[nextPoint]), (myPrimalMesh->vertex[actualPoint]));
1878               subTwoVectors(&result2, (myPrimalMesh->faceCentres[faceAi]), (myPrimalMesh->vertex[actualPoint]));
1879               
1880               crossProduct(&n, result1, result2);
1881               mag(&(a), n);
1882               
1883               sumTwoVectors(&sumN, sumN, n);
1884               
1885               sumA+=a;
1886               scalarDotVector(&result1, a, c);
1887               sumTwoVectors(&sumAc, sumAc, result1);
1888             }
1889           
1890           if (sumA < SMALL)
1891             {
1892               zeroVector( &(myPrimalMesh->faceAreas[faceAi]) );
1893             }
1894           else
1895             {
1896               scalarDotVector(&(myPrimalMesh->faceCentres[faceAi]), 1.0/(3.0*sumA), sumAc); //correct faceCentres after estimating them
1897               scalarDotVector(&(myPrimalMesh->faceAreas[faceAi]), 0.5, sumN);
1898             }
1899           
1900           free(sumN);
1901           free(sumAc);
1902           free(c);
1903           free(n)       ;
1904           free(result1);
1905           free(result2);
1906           
1907         }
1908       
1909     }
1910   
1911 #ifdef DEBUG
1912   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1913     {
1914       debug_print( "%-40s : %ld %d ( ","myPrimalMesh->faceCentres", faceAi, QUAD);
1915       
1916       for(connectivity_int ii=0;ii<DIM3D;ii++)
1917         {
1918           debug_print( "%lf " , myPrimalMesh->faceCentres[faceAi][ii]);
1919         }
1920       
1921       debug_print(")\n");
1922       
1923       debug_print( "%-40s : %ld ( ","myPrimalMesh->faceAreas", faceAi);
1924       
1925       for(connectivity_int ii=0;ii<DIM3D;ii++)
1926         {
1927           debug_print( "%lf " , myPrimalMesh->faceAreas[faceAi][ii]);
1928         }
1929       
1930       debug_print(")\n");
1931     }
1932 #endif
1933   
1934   free(fC);
1935   endFunction(__FUNCTION__, t);
1936   
1940 void setHexahedreEstimateVolumeCentroid(struct primalMesh * myPrimalMesh)
1941 /*
1942  * void XXX(struct primalMesh * myPrimalMesh)
1943  * 
1944  * 
1945  * Input : myPrimalMesh struct
1946  * 
1947  * Use :  
1948  *        
1949  * 
1950  * Allocate : 
1951  * Produce : 
1952  * 
1953  */
1955   
1956   clock_t t=startFunction(__FUNCTION__);
1957   
1958   
1959   dataType *result;
1960   
1961   // vertex number on faces
1962   
1963   connectivity_int it_faceOwner,it_faceNeighbour;
1964   
1965   result =(dataType *) calloc(DIM3D, sizeof(dataType));
1966   
1967   myPrimalMesh->volumeCentroid = (dataType **) malloc(myPrimalMesh->cellNumber * sizeof(dataType*));
1968   for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
1969     {
1970       myPrimalMesh->volumeCentroid[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
1971       memset(myPrimalMesh->volumeCentroid[i], 0.0, DIM3D*sizeof(dataType));
1972     }
1973   
1974   //    myPrimalMesh->cellToFacesNumber = (connectivity_short*) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
1975   
1976   
1977   myPrimalMesh->volume  = (dataType *) malloc(myPrimalMesh->cellNumber * sizeof(dataType));
1978   memset(myPrimalMesh->volume, 0.0, myPrimalMesh->cellNumber*sizeof(dataType));
1979   
1980   
1981   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1982     {
1983       
1984       zeroVector(&result);
1985       
1986       for (connectivity_int ii=0;ii<myPrimalMesh->cellToFacesOwnerNumber[celli];ii++)
1987         {
1988           it_faceOwner = myPrimalMesh->cellToFacesOwner[celli][ii];
1989           
1990           
1991           sumTwoVectors(&result, result, (myPrimalMesh->faceCentres[it_faceOwner]) );
1992           
1993         }
1994       
1995       for (connectivity_int ii=0;ii<myPrimalMesh->cellToFacesNeighbourNumber[celli];ii++)
1996         {
1997           it_faceNeighbour = myPrimalMesh->cellToFacesNeighbour[celli][ii];
1998           
1999           sumTwoVectors(&result, result, (myPrimalMesh->faceCentres[it_faceNeighbour]) );
2000           
2001         }
2002       scalarDivVector(&(myPrimalMesh->volumeCentroid[celli]) , myPrimalMesh->cellToFacesNeighbourNumber[celli]+myPrimalMesh->cellToFacesOwnerNumber[celli], result);
2003       
2004     }
2005   
2006 #ifdef DEBUG  
2007   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
2008     {
2009       
2010       
2011       debug_print("%-40s : %ld (%lf %lf %lf)\n",
2012                   "myPrimalMesh->volumeCentroid Estimated",
2013                   celli,
2014                   myPrimalMesh->volumeCentroid[celli][0],
2015           myPrimalMesh->volumeCentroid[celli][1],
2016           myPrimalMesh->volumeCentroid[celli][2]
2017           );
2018     }
2019 #endif  
2020   
2021   free(result);
2022   endFunction(__FUNCTION__, t);
2023   
2027 void setHexahedreVolumeCentroid(struct primalMesh * myPrimalMesh)
2028 /*
2029  * void XXX(struct primalMesh * myPrimalMesh)
2030  * 
2031  * 
2032  * Input : myPrimalMesh struct
2033  * 
2034  * Use :  
2035  *        
2036  * 
2037  * Allocate : 
2038  * Produce : 
2039  * 
2040  */
2042   clock_t t=startFunction(__FUNCTION__);
2043   
2044   //  connectivity_int  it_cell=0;
2045   dataType        *cellCentroid;
2046   dataType        pyr3Vol=0.0;
2047   dataType        *result1;
2048   dataType        *result2;
2049   dataType        *result3;
2050   
2051   result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
2052   result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
2053   result3 =(dataType *) calloc(DIM3D, sizeof(dataType));
2054   cellCentroid  =(dataType *) calloc(DIM3D, sizeof(dataType));
2055   
2056   // vertex number on faces
2057   connectivity_int it_face;
2058   
2059   // Estimate centroids
2060   setHexahedreEstimateVolumeCentroid(myPrimalMesh);
2061   
2062   
2063   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
2064     {
2065       
2066       equalVector(&cellCentroid, myPrimalMesh->volumeCentroid[celli]);
2067       
2068       zeroVector(&(myPrimalMesh->volumeCentroid[celli]));
2069       
2070       for (connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
2071         {
2072           it_face = myPrimalMesh->cellToFacesOwner[celli][faceAi];
2073           
2074           subTwoVectors( &result1, myPrimalMesh->faceCentres[it_face], cellCentroid );
2075           
2076           dotProduct(&pyr3Vol,myPrimalMesh->faceAreas[it_face], result1);
2077           
2078           scalarDotVector(&result1, 3.0/4.0, myPrimalMesh->faceCentres[it_face]); //
2079           scalarDotVector(&result2, 1.0/4.0, cellCentroid); //
2080           
2081           sumTwoVectors(&result3,result1,result2);
2082           
2083           scalarDotVector(&result1,pyr3Vol,result3);
2084           
2085           sumTwoVectors(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volumeCentroid[celli],result1);
2086           
2087           (myPrimalMesh->volume[celli])+=pyr3Vol;
2088           
2089         }
2090       
2091       for (connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesNeighbourNumber[celli];faceAi++)
2092         {
2093           it_face = myPrimalMesh->cellToFacesNeighbour[celli][faceAi];
2094           
2095           subTwoVectors( &result1, cellCentroid , myPrimalMesh->faceCentres[it_face]);
2096           
2097           dotProduct(&pyr3Vol,myPrimalMesh->faceAreas[it_face], result1);
2098           
2099           scalarDotVector(&result1, 3.0/4.0, myPrimalMesh->faceCentres[it_face]); //
2100           scalarDotVector(&result2, 1.0/4.0, cellCentroid); //
2101           
2102           sumTwoVectors(&result3,result1,result2);
2103           
2104           scalarDotVector(&result1,pyr3Vol,result3);
2105           
2106           sumTwoVectors(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volumeCentroid[celli],result1);
2107           
2108           (myPrimalMesh->volume[celli])+=pyr3Vol;
2109           
2110         }
2111       
2112       
2113       
2114       dataType fVol = fabs(myPrimalMesh->volume[celli]);
2115       if (fVol > VSMALL)
2116         {
2117           scalarDivVector(&(myPrimalMesh->volumeCentroid[celli]),myPrimalMesh->volume[celli],myPrimalMesh->volumeCentroid[celli]);
2118         }
2119       else
2120         {
2121           equalVector(&(myPrimalMesh->volumeCentroid[celli]), cellCentroid);
2122         }
2123       
2124       
2125       myPrimalMesh->volume[celli] *= (1.0/3.0);
2126       
2127 #ifdef DEBUG      
2128       debug_print("%-40s : %ld (%lf) (%lf %lf %lf)\n",
2129                   "myPrimalMesh->volumeCentroid",
2130                   celli,
2131                   myPrimalMesh->volume[celli],
2132                   myPrimalMesh->volumeCentroid[celli][0],
2133           myPrimalMesh->volumeCentroid[celli][1],
2134           myPrimalMesh->volumeCentroid[celli][2]
2135           );
2136 #endif      
2137     }
2138   
2139   free(result1);
2140   free(result2);
2141   free(result3);
2142   
2143   free(cellCentroid);
2144   
2145   endFunction(__FUNCTION__, t);
2146   
2150 void setHexahedreBoundaryFaces(struct primalMesh * myPrimalMesh)
2151 /*
2152  * void XXX(struct primalMesh * myPrimalMesh)
2153  * 
2154  * 
2155  * Input : myPrimalMesh struct
2156  * 
2157  * Use :  
2158  *        
2159  * 
2160  * Allocate : 
2161  * Produce : 
2162  * 
2163  */
2165   clock_t t=startFunction(__FUNCTION__);
2166   
2167   
2168   for(connectivity_int facei=0;facei<myPrimalMesh->faceNumber;facei++)
2169     {
2170       if(myPrimalMesh->faceToCellsNumber[facei] == 1)
2171         {
2172           
2173           myPrimalMesh->faceBoundary = (connectivity_int *) realloc( myPrimalMesh->faceBoundary, (myPrimalMesh->faceBoundaryNumber+1) * sizeof(connectivity_int));
2174           
2175           myPrimalMesh->faceBoundary[myPrimalMesh->faceBoundaryNumber] = facei;
2177           myPrimalMesh->faceBoundaryNumber++;
2178         }
2180     }  
2182   for(connectivity_int facei=0;facei<myPrimalMesh->faceBoundaryNumber;facei++)
2183     {
2184 //      myPrimalMesh->faces[myPrimalMesh->faceBoundary[facei]];
2185     }
2186   
2187 //https://codeforwin.org/2015/07/c-program-to-print-all-unique-element-in-array.html  
2190 #ifdef DEBUG      
2191   for(connectivity_int facei=0;facei<myPrimalMesh->faceBoundaryNumber;facei++)
2192     {
2193       debug_print("%-40s : %ld (%ld)\n",
2194                   "myPrimalMesh->faceBoundary",
2195                   facei,
2196                   myPrimalMesh->faceBoundary[facei]
2197                   
2198                   );
2199     }
2200 #endif   
2201   
2202   endFunction(__FUNCTION__, t);
2203