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

9cdd9c7fded7baee9760fd44c3a881a26c192979
1 #include "mmd.h"
2 #include "tensor.h"
3 #include <signal.h>
5 connectivity_int  hexaedron_localNodeList[HEXAHEDRON_FACES][QUAD]={ \
6   {0,3,2,1}, \
7   {0,1,5,4}, \
8   {0,4,7,3}, \
9   {1,2,6,5}, \
10   {4,5,6,7}, \
11   {3,7,6,2} \
12 };
15 connectivity_int  hexaedron_localNodeListNumbers[HEXAHEDRON_FACES]={ \
16   QUAD, \
17   QUAD, \
18   QUAD, \
19   QUAD, \
20   QUAD, \
21   QUAD \
22 };
24 #define PT1 0
25 #define PT2 1
26 #define PT3 2
27 #define PT4 3
28 #define CELL1 0
29 #define CELL2 1
31 connectivity_int  hexaedron_localSegmentList[HEXAHEDRON_SEGMENTS][SEGMENTVERTEX]={ \
32   {0,1}, \
33   {1,5}, \
34   {5,4}, \
35   {4,0}, \
36   {4,7}, \
37   {1,2}, \
38   {2,3}, \
39   {3,7}, \
40   {7,6}, \
41   {6,2}, \
42   {5,6}, \
43   {0,3} \
44 };
47 connectivity_int  hexaedron_localSegmentsListNumbers[HEXAHEDRON_SEGMENTS]={ \
48   SEGMENTVERTEX, \
49   SEGMENTVERTEX, \
50   SEGMENTVERTEX, \
51   SEGMENTVERTEX, \
52   SEGMENTVERTEX, \
53   SEGMENTVERTEX, \
54   SEGMENTVERTEX, \
55   SEGMENTVERTEX, \
56   SEGMENTVERTEX, \
57   SEGMENTVERTEX, \
58   SEGMENTVERTEX, \
59   SEGMENTVERTEX \
60 };
63 /*!
64  * \fn void rotate(connectivity_int arr[], connectivity_int n)
65  * \brief Circular displacement 
66  * \param arr { vector  }
67  * \param n { Offset number }
68  * \return arr {vector with inverted values}
69  * 
70  * Circular displacement of the elements of a vector arr of n elements
71  */
72 void rotate(connectivity_int arr[], connectivity_int n)
73 {
74   connectivity_int x = arr[0], i;
75   for (i = 0; i <n-1; i++)
76     arr[i] = arr[i+1];
77   arr[n-1] = x;
78 }
80 clock_t startFunction(const char * functionName)
81 {
82   affiche("%-40s : Start\n", functionName);
83   
84   return clock();
85 }
87 void endFunction(const char * functionName, clock_t t)
88 {
89   t = clock() - t;
90   double time_taken = ((double)t) / CLOCKS_PER_SEC; // in seconds
91   
92   affiche("%-40s : %f seconds to execute\n", functionName, time_taken);
93   
94 }
96 void removeDuplicates(connectivity_int arr[], connectivity_int * size)
97 {
98   for(connectivity_int i=0; i<(*size); i++)
99     {
100       for(connectivity_int j=i+1; j<(*size); j++)
101         {
102           /* If any duplicate found */
103           if(arr[i] == arr[j])
104             {
105               /* Delete the current duplicate element */
106               for(connectivity_int k=j; k<(*size)-1; k++)
107                 {
108                   arr[k] = arr[k + 1];
109                 }
110               
111               /* Decrement size after removing duplicate element */
112               (*size)--;
113               
114               /* If shifting of elements occur then don't increment j */
115               j--;
116             }
117         }
118     }
121 void sortListConnectivity(connectivity_int a[], connectivity_int  n)
123   connectivity_int tmp = 0;
124   
125   for (connectivity_int i = 0; i < n; i++)                     //Loop for ascending ordering
126     {
127       for (connectivity_int j = 0; j < n; j++)             //Loop for comparing other values
128         {
129           if (a[j] > a[i])                //Comparing other array elements
130             {
131               tmp = a[i];         //Using temporary variable for storing last value
132               a[i] = a[j];            //replacing value
133               a[j] = tmp;             //storing last value
134             }
135         }
136     }
137   
141 void sortList(connectivity_int ** vector,connectivity_int vectorSize, connectivity_int position)
143   
144   for(connectivity_int ii=0;ii<position;ii++)
145     rotate(*vector, vectorSize);
146   
149 long findInList(connectivity_int * vector,connectivity_int vectorSize, connectivity_int value)
151   
152   for(connectivity_int ii=0;ii<vectorSize;ii++)
153     {
154       if(vector[ii]==value)
155         {
156           return (long)(ii);
157         }
158     }
159   
160   return -1;
161   
163 connectivity_int findMinInList(connectivity_int * vector,connectivity_int vectorSize)
165   
166   connectivity_int              position=0;
167   
168   for(connectivity_int ii=1;ii<vectorSize;ii++)
169     {
170       if(vector[position]>vector[ii])
171         {
172           position=ii;
173         }
174     }
175   
176   return(position);
180 void invertListSort(connectivity_int ** result,connectivity_int * vector,connectivity_int vectorSize)
182   connectivity_int position = 0;
183   
184   for(connectivity_int ii=0;ii<vectorSize;ii++)
185     (*result)[ii] = vector[vectorSize -1 - ii];
186   position = findMinInList((*result),vectorSize);
187   sortList(result, vectorSize, position);
188   
191 void intersectTwoIntegerSet(connectivity_int ** output , connectivity_int * lenOutput , connectivity_int * inputOne , connectivity_int  lenInputOne , connectivity_int * inputTwo, connectivity_int  lenInputTwo)
193   
194   (*lenOutput)=0;
195   //  affiche("%ld\n",lenInputOne);
196   //  affiche("%ld\n",lenInputTwo);
197   
198   for(connectivity_int posAi=0;posAi<lenInputOne;posAi++)
199     for(connectivity_int posBi=0;posBi<lenInputTwo;posBi++)
200       {
201         if(inputOne[posAi]==inputTwo[posBi])
202           {
203             (*output)[(*lenOutput)] = inputOne[posAi];
204             (*lenOutput)++;
205           }
206       }
207   
208   removeDuplicates((*output), lenOutput);
213 void unionTwoIntegerSet(connectivity_int ** output , connectivity_int * lenOutput , connectivity_int * inputOne , connectivity_int  lenInputOne , connectivity_int * inputTwo, connectivity_int  lenInputTwo)
215   
216   (*lenOutput)=0;
217   //  affiche("%ld\n",lenInputOne);
218   //  affiche("%ld\n",lenInputTwo);
219   
220   (*output) = (connectivity_int * ) malloc(sizeof(connectivity_int)*(lenInputOne+lenInputTwo));
221   
222   for(connectivity_int posAi=0;posAi<lenInputOne;posAi++)
223     {
224       (*output)[(*lenOutput)] = inputOne[posAi];
225       (*lenOutput)++;
226     }
227   
228   for(connectivity_int posBi=0;posBi<lenInputTwo;posBi++)
229     {
230       (*output)[(*lenOutput)] = inputTwo[posBi];
231       (*lenOutput)++;
232     }
233   
234   removeDuplicates((*output), lenOutput);
238 void findListFromList(connectivity_int ** matrix, connectivity_int * matrixRowNumbers, connectivity_int size, connectivity_int * whatToFind, connectivity_int sizewhatToFind )
240   
241   //connectivity_int tmp = 0;
242   
243   for (connectivity_int i = 0; i < size; i++)
244     {
245       if(findInList(matrix[i],matrixRowNumbers[i], whatToFind[0])!=-1)
246         {
247           for (connectivity_int j = 0; j < matrixRowNumbers[i]; j++)
248             {
249               
250               for (connectivity_int k=1;k<sizewhatToFind;k++)
251                 {
252                   
253                 }
254               
255             }
256         }
257     }
258   
264 //#define allocMemoryForMatrix(x,y,z) _Generic((x), dataType: allocMemoryForMatrixDataType, char: allocMemoryForMatrixDataType)(x,y,z);
265 void allocMemoryForMatrixDataType(dataType *** result, connectivity_int N, connectivity_int M )
267   
268   //    affiche("\n");
269   connectivity_int sum=0;
270   
271   (*result) = (dataType **) malloc(N * sizeof(dataType *));
272   if((*result)==NULL)
273     {
274       debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
275       exit(1);
276     }
277   for(connectivity_int i=0; i<N; i++)
278     {
279       (*result)[i] = (dataType *) malloc(M * sizeof(dataType));
280       memset((*result)[i], 0.0, M*sizeof(connectivity_int));
281       if((*result)[i]==NULL)
282         {
283           debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
284           debug_print("i=%ld M=%ld N=%ld  \n",i, M, N);
285           exit(1);
286         }
287     }
288   
289   sum = sizeof(dataType *)*N + M * sizeof(dataType) ;
290   
291   affiche("%-40s : (%10ld bytes) : (%8.2f Mo)\n","Allocated Memory",sum, sum/1024000.0);
292   
293   
294   
297 //#define allocMemoryForMatrix(x,y,z) _Generic((x), dataType: allocMemoryForMatrixDataType, char: allocMemoryForMatrixDataType)(x,y,z);
298 void allocMemoryForVectorDataType(dataType ** result, connectivity_int N )
300   
301   //    affiche("\n");
302   connectivity_int sum=0;
303   
304   (*result) = (dataType *) malloc(N * sizeof(dataType ));
305   if((*result)==NULL)
306     {
307       debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
308       exit(1);
309     }
310   
311   sum = sizeof(dataType )*N ;
312   
313   affiche("%-40s : (%10ld bytes) : (%8.2f Mo)\n","Allocated Memory",sum, sum/1024000.0);
314   
315   
316   
320 //#define allocMemoryForMatrix(x,y,z) _Generic((x), dataType: allocMemoryForMatrixDataType, char: allocMemoryForMatrixDataType)(x,y,z);
321 void allocMemoryForVectorConnectivity(connectivity_int ** result, connectivity_int N )
323   
324   //    affiche("\n");
325   connectivity_int sum=0;
326   
327   (*result) = (connectivity_int *) malloc(N * sizeof(connectivity_int ));
328   if((*result)==NULL)
329     {
330       debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
331       exit(1);
332     }
333   
334   sum = sizeof(connectivity_int )*N ;
335   
336   affiche("%-40s : (%10ld bytes) : (%8.2f Mo)\n","Allocated Memory",sum, sum/1024000.0);
337   
338   
339   
343 //#define allocMemoryForMatrix(x,y,z) _Generic((x), dataType: allocMemoryForMatrixDataType, char: allocMemoryForMatrixDataType)(x,y,z);
344 void allocMemoryForMatrixConnectivity(connectivity_int *** result, connectivity_int N, connectivity_int M )
346   
347   //    affiche("\n");
348   connectivity_int sum=0;
349   
350   (*result) = (connectivity_int **) malloc(N * sizeof(connectivity_int *));
351   if((*result)==NULL)
352     {
353       debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
354       exit(1);
355     }
356   for(connectivity_int i=0; i<N; i++)
357     {
358       (*result)[i] = (connectivity_int *) malloc(M * sizeof(connectivity_int));
359       memset((*result)[i], 0.0, M*sizeof(connectivity_int));
360       if((*result)[i]==NULL)
361         {
362           debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
363           debug_print("i=%ld M=%ld N=%ld  \n",i, M, N);
364           exit(1);
365         }
366     }
367   
368   sum = sizeof(connectivity_int *)*N + M * sizeof(connectivity_int) ;
369   
370   affiche("%-40s : (%10ld bytes) : (%8.2f Mo)\n","Allocated Memory",sum, sum/1024000.0);
371   
375 void allocMemoryForSpecialMatrixConnectivity(connectivity_int *** result, connectivity_int N, connectivity_short *M )
377   
378   connectivity_int sum=0;
379   
380   //    affiche("\n");
381   
382   (*result) = (connectivity_int **) malloc(N * sizeof(connectivity_int *));
383   if((*result)==NULL)
384     {
385       debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
386       exit(1);
387     }
388   for(connectivity_int i=0; i<N; i++)
389     {
390       (*result)[i] = (connectivity_int *) malloc(M[i] * sizeof(connectivity_int));
391       memset((*result)[i], 0.0, M[i]*sizeof(connectivity_int));
392       sum+=M[i] * sizeof(connectivity_int);
393       
394       if((*result)[i]==NULL) {
395           debug_print("ERROR : %s %d %s\n",__FUNCTION__, __LINE__, __FILE__);
396           debug_print("i=%ld M=%d N=%ld  \n",i, M[i], N);
397           exit(1);
398         }
399     }
400   sum = sizeof(connectivity_int *)*N + sum ;
401   
402   affiche("%-40s : (%10ld bytes) : (%8.2f Mo)\n","Allocated Memory",sum, sum/1024000.0);
403   
407 //void allocMemory(struct primalMesh * myPrimalMesh)
408 //{
409 //  clock_t t=startFunction(__FUNCTION__);
413 //    myPrimalMesh->vertex = (dataType **) malloc(myPrimalMesh->vertexNumber * sizeof(dataType *));
414 //    //ERROR_TEST( (myPrimalMesh->vertex) );
415 //    for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
416 //    {
417 //        myPrimalMesh->vertex[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
418 //        memset(myPrimalMesh->vertex[i], 0.0, DIM3D*sizeof(dataType));
419 //    }
421 //    myPrimalMesh->cellToVertexNumber = (connectivity_short *) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
425 //    myPrimalMesh->faceToVertexNumber = (connectivity_short *) malloc(myPrimalMesh->faceNumber * sizeof(connectivity_short));
428 //    myPrimalMesh->faceCentres = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType *));
429 //    for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
430 //    {
431 //        myPrimalMesh->faceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
432 //        memset(myPrimalMesh->faceCentres[i], 0.0, DIM3D*sizeof(dataType));
433 //    }
435 //    myPrimalMesh->faceAreas = (dataType **) malloc(myPrimalMesh->faceNumber * sizeof(dataType*));
436 //    for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
437 //    {
438 //        myPrimalMesh->faceAreas[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
439 //        memset(myPrimalMesh->faceAreas[i], 0.0, DIM3D*sizeof(dataType));
440 //    }
442 //    myPrimalMesh->volumeCentroid = (dataType **) malloc(myPrimalMesh->cellNumber * sizeof(dataType*));
443 //    for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
444 //    {
445 //        myPrimalMesh->volumeCentroid[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
446 //        memset(myPrimalMesh->volumeCentroid[i], 0.0, DIM3D*sizeof(dataType));
447 //    }
449 //    //    myPrimalMesh->cellToFacesNumber = (connectivity_short*) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
452 //    myPrimalMesh->volume  = (dataType *) malloc(myPrimalMesh->cellNumber * sizeof(dataType));
453 //    memset(myPrimalMesh->volume, 0.0, myPrimalMesh->cellNumber*sizeof(dataType));
457 //    myPrimalMesh->segmentToVertex = (connectivity_int **) malloc(myPrimalMesh->segmentNumber * sizeof(connectivity_int*));
458 //    for(connectivity_int i=0; i<myPrimalMesh->segmentNumber; i++)
459 //    {
460 //        myPrimalMesh->segmentToVertex[i] = (connectivity_int *) malloc( SEGMENTVERTEX * sizeof(connectivity_int));
461 //        memset(myPrimalMesh->segmentToVertex[i], 0.0, SEGMENTVERTEX*sizeof(connectivity_int));
462 //    }
466 //    myPrimalMesh->faceToSegments = (connectivity_int **) malloc(myPrimalMesh->faceNumber * sizeof(connectivity_int*));
467 //    for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
468 //    {
469 //        myPrimalMesh->faceToSegments[i] = (connectivity_int *) malloc( QUAD * sizeof(connectivity_int));
470 //        memset(myPrimalMesh->faceToSegments[i], 0.0, QUAD*sizeof(connectivity_int));
471 //    }
473 //    myPrimalMesh->vertexToSegmentNumber  = (connectivity_short *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_short));
474 //    memset(myPrimalMesh->vertexToSegmentNumber, 0.0, myPrimalMesh->vertexNumber*sizeof(connectivity_short));
477 //  endFunction(__FUNCTION__, t);
479 //}
483 //void allocMemory2(struct primalMesh * myPrimalMesh)
484 //{
485 //  clock_t t=startFunction(__FUNCTION__);
488 //    myPrimalMesh->cellToVertex = (connectivity_int **) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_int *));
489 //    for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
490 //    {
491 //        myPrimalMesh->cellToVertex[i] = (connectivity_int *) malloc(myPrimalMesh->cellToVertexNumber[i] * sizeof(connectivity_int));
492 //        memset(myPrimalMesh->cellToVertex[i], 0, myPrimalMesh->cellToVertexNumber[i]*sizeof(connectivity_int));
493 //    }
495 //    myPrimalMesh->faceToVertex = (connectivity_int **) malloc(myPrimalMesh->faceNumber * sizeof(connectivity_int *));
496 //    for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
497 //    {
498 //        myPrimalMesh->faceToVertex[i] = (connectivity_int *) malloc(myPrimalMesh->faceToVertexNumber[i] * sizeof(connectivity_int));
499 //        memset(myPrimalMesh->faceToVertex[i], 0, myPrimalMesh->faceToVertexNumber[i]*sizeof(connectivity_int));
500 //    }
502 //    myPrimalMesh->cellToFaces = (connectivity_int **) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_int *));
503 //    for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
504 //    {
505 //        myPrimalMesh->cellToFaces[i] = (connectivity_int *) malloc(myPrimalMesh->cellToFacesNumber[i] * sizeof(connectivity_int));
506 //        memset(myPrimalMesh->cellToFaces[i], 0, myPrimalMesh->cellToFacesNumber[i]*sizeof(connectivity_int));
507 //    }
511 //  endFunction(__FUNCTION__, t);
513 //}
515 void freeMemory(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
517   clock_t t=startFunction(__FUNCTION__);
518   
519   
520   if(myPrimalMesh->vertex!=NULL){
521       for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
522         {
523           if(myPrimalMesh->vertex[i]!=NULL)
524             free(myPrimalMesh->vertex[i]);
525         }
526       free(myPrimalMesh->vertex);
527     }
528   
529   if(myPrimalMesh->cellToVertexNumber!=NULL){
530       free(myPrimalMesh->cellToVertexNumber);
531     }
532   
533   if(myPrimalMesh->cellToVertex!=NULL){
534       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
535         {
536           if(myPrimalMesh->cellToVertex[i]!=NULL)
537             free(myPrimalMesh->cellToVertex[i]);
538         }
539       free(myPrimalMesh->cellToVertex);
540     }
541   
542   if(myPrimalMesh->faceToVertex!=NULL){
543       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
544         {
545           if(myPrimalMesh->faceToVertex[i]!=NULL)
546             free(myPrimalMesh->faceToVertex[i]);
547         }
548       free(myPrimalMesh->faceToVertex);
549     }
550   
551   if(myPrimalMesh->faceToVertexNumber!=NULL){
552       free(myPrimalMesh->faceToVertexNumber);
553     }
554   
555   if(myPrimalMesh->faceCentres!=NULL){
556       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
557         {
558           if(myPrimalMesh->faceCentres[i]!=NULL)
559             free(myPrimalMesh->faceCentres[i]);
560         }
561       free(myPrimalMesh->faceCentres);
562     }
563   
564   
565   if(myPrimalMesh->faceAreas!=NULL){
566       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
567         {
568           if(myPrimalMesh->faceAreas[i]!=NULL)
569             free(myPrimalMesh->faceAreas[i]);
570         }
571       free(myPrimalMesh->faceAreas);
572     }
573   
574   if(myPrimalMesh->cellToFaces!=NULL){
575       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
576         {
577           if(myPrimalMesh->cellToFaces[i]!=NULL)
578             free(myPrimalMesh->cellToFaces[i]);
579         }
580       free(myPrimalMesh->cellToFaces);
581     }
582   
583   if(myPrimalMesh->volumeCentroid!=NULL){
584       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
585         {
586           if(myPrimalMesh->volumeCentroid[i]!=NULL)
587             free(myPrimalMesh->volumeCentroid[i]);
588         }
589       free(myPrimalMesh->volumeCentroid);
590     }
591   
592   if(myPrimalMesh->cellToFacesNumber!=NULL){
593       free(myPrimalMesh->cellToFacesNumber);
594     }
595   
596   if(myPrimalMesh->volume!=NULL){
597       free(myPrimalMesh->volume);
598     }
599   
600   
601   //  if(myPrimalMesh->segmentToVertex!=NULL){
602   //      for(connectivity_int i=0; i<myPrimalMesh->segmentNumber; i++)
603   //        {
604   //          if(myPrimalMesh->segmentToVertex[i]!=NULL)
605   //            free(myPrimalMesh->segmentToVertex[i]);
606   //        }
607   //      free(myPrimalMesh->segmentToVertex);
608   //    }
609   
610   if(myPrimalMesh->faceToSegments!=NULL){
611       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
612         {
613           if(myPrimalMesh->faceToSegments[i]!=NULL)
614             free(myPrimalMesh->faceToSegments[i]);
615         }
616       free(myPrimalMesh->faceToSegments);
617     }
618   
619   
620   if(myPrimalMesh->vertexToSegmentNumber!=NULL){
621       free(myPrimalMesh->vertexToSegmentNumber);
622     }
623   
624   
625   if(myPrimalMesh->vertexToSegments!=NULL){
626       for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
627         {
628           if(myPrimalMesh->vertexToSegments[i]!=NULL)
629             free(myPrimalMesh->vertexToSegments[i]);
630         }
631       free(myPrimalMesh->vertexToSegments);
632     }
633   
634   
635   
636   if(myPrimalMesh->vertexToCells!=NULL){
637       for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
638         {
639           if(myPrimalMesh->vertexToCells[i]!=NULL)
640             free(myPrimalMesh->vertexToCells[i]);
641         }
642       free(myPrimalMesh->vertexToCells);
643     }
644   
645   
646   
647   if(myPrimalMesh->vertexToCells!=NULL){
648       
649       free(myPrimalMesh->vertexToCellNumbers);
650     }
651   
652   
653   if(myPrimalMesh->cellToCellsNumbers!=NULL){
654       
655       free(myPrimalMesh->cellToCellsNumbers);
656     }
657   
658   
659   if(myPrimalMesh->cellToCells!=NULL){
660       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
661         {
662           if(myPrimalMesh->cellToCells[i]!=NULL)
663             free(myPrimalMesh->cellToCells[i]);
664         }
665       free(myPrimalMesh->cellToCells);
666     }
667   
668   
669   
670   if(myPrimalMesh->faces!=NULL){
671       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
672         {
673           if(myPrimalMesh->faces[i]!=NULL)
674             free(myPrimalMesh->faces[i]);
675         }
676       free(myPrimalMesh->faces);
677     }
678   
679   
680   if(myPrimalMesh->cellToFacesOwnerNumber!=NULL)
681     free(myPrimalMesh->cellToFacesOwnerNumber);
682   if(myPrimalMesh->cellToFacesNeighbourNumber!=NULL)
683     free(myPrimalMesh->cellToFacesNeighbourNumber);
684   
685   
686   if(myPrimalMesh->cellToFacesNeighbour!=NULL){
687       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
688         {
689           if(myPrimalMesh->cellToFacesNeighbour[i]!=NULL)
690             free(myPrimalMesh->cellToFacesNeighbour[i]);
691         }
692       free(myPrimalMesh->cellToFacesNeighbour);
693     }
694   
695   
696   if(myPrimalMesh->cellToFacesOwner!=NULL){
697       for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
698         {
699           if(myPrimalMesh->cellToFacesOwner[i]!=NULL)
700             free(myPrimalMesh->cellToFacesOwner[i]);
701         }
702       free(myPrimalMesh->cellToFacesOwner);
703     }
704   
705   
706   if(myPrimalMesh->segments!=NULL){
707       for(connectivity_int i=0; i<myPrimalMesh->segmentNumber; i++)
708         {
709           if(myPrimalMesh->segments[i]!=NULL)
710             free(myPrimalMesh->segments[i]);
711         }
712       free(myPrimalMesh->segments);
713     }
714   
715   
716   
717   if(myPrimalMesh->vertexToSegmentOwner!=NULL){
718       for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
719         {
720           if(myPrimalMesh->vertexToSegmentOwner[i]!=NULL)
721             free(myPrimalMesh->vertexToSegmentOwner[i]);
722         }
723       free(myPrimalMesh->vertexToSegmentOwner);
724     }
725   
726   
727   
728   if(myPrimalMesh->vertexToSegmentNeighbour!=NULL){
729       for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
730         {
731           if(myPrimalMesh->vertexToSegmentNeighbour[i]!=NULL)
732             free(myPrimalMesh->vertexToSegmentNeighbour[i]);
733         }
734       free(myPrimalMesh->vertexToSegmentNeighbour);
735     }
736   
737   
738   if(myPrimalMesh->vertexToSegmentOwnerNumber!=NULL)
739     free(myPrimalMesh->vertexToSegmentOwnerNumber);
740   if(myPrimalMesh->vertexToSegmentNeighbourNumber!=NULL)
741     free(myPrimalMesh->vertexToSegmentNeighbourNumber);
742   
743   
744   
745   
746   if(myPrimalMesh->segmentToFaceOwnerNumber!=NULL)
747     free(myPrimalMesh->segmentToFaceOwnerNumber);
748   
749   
750   if(myPrimalMesh->segmentToFaceOwner!=NULL){
751       for(connectivity_int i=0; i<myPrimalMesh->segmentNumber; i++)
752         {
753           if(myPrimalMesh->segmentToFaceOwner[i]!=NULL)
754             free(myPrimalMesh->segmentToFaceOwner[i]);
755         }
756       free(myPrimalMesh->segmentToFaceOwner);
757     }
758   
759   if(myPrimalMesh->segmentToFaceNeighbourNumber!=NULL)
760     free(myPrimalMesh->segmentToFaceNeighbourNumber);
761   
762   if(myPrimalMesh->segmentToFaceNeighbour!=NULL){
763       for(connectivity_int i=0; i<myPrimalMesh->segmentNumber; i++)
764         {
765           if(myPrimalMesh->segmentToFaceNeighbour[i]!=NULL)
766             free(myPrimalMesh->segmentToFaceNeighbour[i]);
767         }
768       free(myPrimalMesh->segmentToFaceNeighbour);
769     }
770   
771   
772   if(myPrimalMesh->faceToCellsNumber!=NULL)
773     free(myPrimalMesh->faceToCellsNumber);
774   
775   if(myPrimalMesh->faceToCells!=NULL){
776       for(connectivity_int i=0; i<myPrimalMesh->faceNumber; i++)
777         {
778           if(myPrimalMesh->faceToCells[i]!=NULL)
779             free(myPrimalMesh->faceToCells[i]);
780         }
781       free(myPrimalMesh->faceToCells);
782     }
783   
784   if(myPrimalMesh->cellToSegmentOwner!=NULL)
785     {
786       for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
787         free(myPrimalMesh->cellToSegmentOwner[celli]);
788       free(myPrimalMesh->cellToSegmentOwner);  
789     }
790   
791   if(myPrimalMesh->cellToSegmentNeighbour!=NULL)
792     {
793       for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
794         free(myPrimalMesh->cellToSegmentNeighbour[celli]);
795       free(myPrimalMesh->cellToSegmentNeighbour);  
796     }
797   
798   if(myPrimalMesh->cellToSegmentOwnerNumber!=NULL)
799     {
800       free(myPrimalMesh->cellToSegmentOwnerNumber);  
801     }
802   
803   if(myPrimalMesh->cellToSegmentNeighbourNumber!=NULL)
804     {
805       free(myPrimalMesh->cellToSegmentNeighbourNumber); 
806     }
807   
808   
809   
810   /// DUAL MESH
811   
812   
813   if(myDualMesh->internalDualFaces!=NULL)
814     {
815       for(connectivity_int facei=0;facei<myDualMesh->internalDualFacesNumber;facei++)
816         free(myDualMesh->internalDualFaces[facei]);
817       free(myDualMesh->internalDualFaces);  
818     }
819   
820   
821   
822   if(myDualMesh->segmentToInternalDualFace!=NULL)
823     {
824       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
825         free(myDualMesh->segmentToInternalDualFace[segmenti]);
826       free(myDualMesh->segmentToInternalDualFace);  
827     }
828   
829   
830   if(myDualMesh->segmentToInternalDualFaceNumber!=NULL)
831     {
832       free(myDualMesh->segmentToInternalDualFaceNumber); 
833     }
834   
835   
836   
837   
838   
839   endFunction(__FUNCTION__, t);
840   
845 void setHexahedreVertex(struct primalMesh * myPrimalMesh)
847   
848   //    clock_t t;
849   //    t = clock();
850   
851   clock_t t=startFunction(__FUNCTION__);
852   affiche("\n");
853   
854   // definition des vertex
855   connectivity_int it_vertex=0;
856   //    connectivity_int it_cell=0;
857   
858   dataType L=myPrimalMesh->L/myPrimalMesh->M;
859   dataType l=myPrimalMesh->l/myPrimalMesh->N;
860   dataType H=myPrimalMesh->H/myPrimalMesh->P;
861   
862   allocMemoryForMatrixDataType(&(myPrimalMesh->vertex) , myPrimalMesh->vertexNumber,DIM3D);
863   
864   for(connectivity_int i=0;i<myPrimalMesh->M+1;i++)
865     for(connectivity_int j=0;j<myPrimalMesh->N+1;j++)
866       for(connectivity_int k=0;k<myPrimalMesh->P+1;k++)
867         {
868           it_vertex = k + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
869           myPrimalMesh->vertex[it_vertex][0] = i * L;
870           myPrimalMesh->vertex[it_vertex][1] = j * l;
871           myPrimalMesh->vertex[it_vertex][2] = k * H;
872 #ifdef DEBUG  
873           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]);
874 #endif
875         }
876   
877   endFunction(__FUNCTION__, t);
878   
879   
882 void setHexahedreCellToVertexNumber(struct primalMesh * myPrimalMesh)
884   clock_t t=startFunction(__FUNCTION__);
885   
886   // definition des vertex
887   //    connectivity_int it_vertex=0;
888   connectivity_int it_cell=0;
889   
890   
891   myPrimalMesh->cellToVertexNumber = (connectivity_int *) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_int));
892   
893   // hexaedron number of vertex on hexaedron
894   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
895     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
896       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
897         {
898           it_cell     = k + j*myPrimalMesh->P + i*myPrimalMesh->P*myPrimalMesh->N;
899           myPrimalMesh->cellToVertexNumber[it_cell]=HEXAHEDRON;
900           debug_print("%-40s : %ld (%ld)\n","myPrimalMesh->cellToVertexNumber", it_cell, myPrimalMesh->cellToVertexNumber[it_cell] );
901         }
902   endFunction(__FUNCTION__, t);
903   
907 //void setHexahedreFaceToVertexNumber(struct primalMesh * myPrimalMesh)
908 //{
909 //  // vertex number on faces
910 //  clock_t t=startFunction(__FUNCTION__);
912 //  myPrimalMesh->faceToVertexNumber = (connectivity_short *) malloc(myPrimalMesh->faceNumber * sizeof(connectivity_short));
914 //  for(connectivity_int i=0;i<myPrimalMesh->faceNumber;i++)
915 //    {
918 //      myPrimalMesh->faceToVertexNumber[i] = QUAD;
920 //      debug_print( "%-40s : %ld (%d)\n",
921 //                   "myPrimalMesh.faceToVertexNumber",
922 //                   i,
923 //                   myPrimalMesh->faceToVertexNumber[i]
924 //                   );
927 //    }
928 //  endFunction(__FUNCTION__, t);
929 //}
931 //void setHexahedreCellToFacesNumber(struct primalMesh * myPrimalMesh)
932 //{
934 //  clock_t t=startFunction(__FUNCTION__);
936 //  myPrimalMesh->cellToFacesNumber = (connectivity_short*) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_short));
939 //  for (connectivity_int i=0;i<myPrimalMesh->cellNumber;i++)
940 //    {
941 //      if(myPrimalMesh->cellToVertexNumber[i] == HEXAHEDRON)
942 //        myPrimalMesh->cellToFacesNumber[i] = HEXAHEDRON_FACES;
944 //    }
946 //  for(connectivity_int ii=0;ii<myPrimalMesh->cellNumber;ii++)
947 //    {
948 //      debug_print( "%-40s : %ld (%d) \n" ,"myPrimalMesh->cellToFacesNumber", ii, myPrimalMesh->cellToFacesNumber[ii]);
949 //    }
950 //  endFunction(__FUNCTION__, t);
951 //}
954 void setHexahedreCellToVertex(struct primalMesh * myPrimalMesh)
956   
957   clock_t t=startFunction(__FUNCTION__);
958   
959   connectivity_int it_vertex=0;
960   connectivity_int it_vertex2=0;
961   connectivity_int it_vertex3=0;
962   connectivity_int it_vertex4=0;
963   connectivity_int it_vertex5=0;
964   connectivity_int it_vertex6=0;
965   connectivity_int it_vertex7=0;
966   connectivity_int it_vertex8=0;
967   connectivity_int it_cell=0;
968   // cell 0 : 0 1 9 8 / 2 3 11 10
969   // hexaedron vertex number to define hexaedron
970   
971   connectivity_int vN1 = (myPrimalMesh->P+1) ;
972   connectivity_int vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1) ;
973   connectivity_int cN1 = myPrimalMesh->P*myPrimalMesh->N;
974   
975   
976   
977   myPrimalMesh->cellToVertex = (connectivity_int **) malloc(myPrimalMesh->cellNumber * sizeof(connectivity_int *));
978   for(connectivity_int i=0; i<myPrimalMesh->cellNumber; i++)
979     {
980       myPrimalMesh->cellToVertex[i] = (connectivity_int *) malloc(myPrimalMesh->cellToVertexNumber[i] * sizeof(connectivity_int));
981       memset(myPrimalMesh->cellToVertex[i], 0, myPrimalMesh->cellToVertexNumber[i]*sizeof(connectivity_int));
982     }
983   
984   
985   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
986     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
987       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
988         {
989           it_cell     = k + j*myPrimalMesh->P + i*cN1; //cN1 = myPrimalMesh->P*myPrimalMesh->N;
990           
991           /*
992 *                 it_vertex =   k      + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
993             it_vertex2 = (k + 1) + j*(myPrimalMesh->P+1) +i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
994             it_vertex3 = (k + 1) + j*(myPrimalMesh->P+1) +(i+1)*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
995             it_vertex4 =  k      + j*(myPrimalMesh->P+1) +(i+1)*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
996 */
997           it_vertex =   k      + j*vN1 +i*vN2; // vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1);
998           it_vertex2 = it_vertex  + 1;
999           it_vertex3 = it_vertex2 + vN2; // vN2 = (myPrimalMesh->P+1)*(myPrimalMesh->N+1);
1000           it_vertex4 = it_vertex3 - 1;
1001           
1002           myPrimalMesh->cellToVertex[it_cell][0]=it_vertex;
1003           myPrimalMesh->cellToVertex[it_cell][1]=it_vertex2;
1004           myPrimalMesh->cellToVertex[it_cell][2]=it_vertex3;
1005           myPrimalMesh->cellToVertex[it_cell][3]=it_vertex4;
1006           
1007           
1008           
1009           it_vertex5 =  k      + (j+1)*vN1 +i*vN2;
1010           it_vertex6 = (k + 1) + (j+1)*vN1 +i*vN2;
1011           it_vertex7 = (k + 1) + (j+1)*vN1 +(i+1)*vN2;
1012           it_vertex8 =  k      + (j+1)*vN1 +(i+1)*vN2;
1013           
1014           myPrimalMesh->cellToVertex[it_cell][4]=it_vertex5;
1015           myPrimalMesh->cellToVertex[it_cell][5]=it_vertex6;
1016           myPrimalMesh->cellToVertex[it_cell][6]=it_vertex7;
1017           myPrimalMesh->cellToVertex[it_cell][7]=it_vertex8;
1018           
1019           
1020 #ifdef DEBUG  
1021           
1022           debug_print("%-40s : %ld (%ld %ld %ld %ld %ld %ld %ld %ld)\n",
1023                       "myPrimalMesh->cellToVertex",
1024                       it_cell,
1025                       myPrimalMesh->cellToVertex[it_cell][0],
1026               myPrimalMesh->cellToVertex[it_cell][1],
1027               myPrimalMesh->cellToVertex[it_cell][2],
1028               myPrimalMesh->cellToVertex[it_cell][3],
1029               myPrimalMesh->cellToVertex[it_cell][4],
1030               myPrimalMesh->cellToVertex[it_cell][5],
1031               myPrimalMesh->cellToVertex[it_cell][6],
1032               myPrimalMesh->cellToVertex[it_cell][7] );
1033 #endif
1034           
1035         }
1036   endFunction(__FUNCTION__, t);
1040 //void setHexahedreFaceToVertex(struct primalMesh * myPrimalMesh)
1041 //{
1043 //  clock_t t=startFunction(__FUNCTION__);
1044 //  affiche("\n");
1046 //  connectivity_int it_cell=0;
1048 //  // vertex number on faces
1049 //  connectivity_int cN1 =  myPrimalMesh->P*myPrimalMesh->N;
1051 //  connectivity_int it_surface=0;
1052 //  connectivity_int it_face=0;
1055 //  allocMemoryForSpecialMatrixConnectivity(&(myPrimalMesh->faceToVertex),myPrimalMesh->faceNumber, myPrimalMesh->faceToVertexNumber);
1059 //  for(connectivity_int i=0;i<myPrimalMesh->M;i++)
1060 //    for(connectivity_int j=0;j<myPrimalMesh->N;j++)
1061 //      for(connectivity_int k=0;k<myPrimalMesh->P;k++)
1062 //        {
1063 //          //it_surface = k + j*(myPrimalMesh->P+1) + i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
1064 //          it_cell     = k + j*myPrimalMesh->P + i*cN1; //cN1 = myPrimalMesh->P*myPrimalMesh->N;
1065 //          it_face     = 0;
1067 //          if(myPrimalMesh->cellToVertexNumber[it_cell] == HEXAHEDRON)
1068 //            {
1069 //              connectivity_int vertexList[HEXAHEDRON];
1070 //              // FACE 1
1071 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][0];
1072 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][3];
1073 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][2];
1074 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][1];
1076 //              //                    vertexList[0] = myPrimalMesh->cellToVertex[it_cell][0];
1077 //              //                    vertexList[1] = myPrimalMesh->cellToVertex[it_cell][3];
1078 //              //                    vertexList[2] = myPrimalMesh->cellToVertex[it_cell][2];
1079 //              //                    vertexList[3] = myPrimalMesh->cellToVertex[it_cell][1];
1081 //              //                    // FAIRE TRIER LISTE
1082 //              //                    // INVERT LISTE et TRIER LISTE
1084 //              //                    myPrimalMesh->faceToVertex[it_surface][0] = vertexList[0];
1085 //              //                    myPrimalMesh->faceToVertex[it_surface][1] = vertexList[3];
1086 //              //                    myPrimalMesh->faceToVertex[it_surface][2] = vertexList[2];
1087 //              //                    myPrimalMesh->faceToVertex[it_surface][3] = vertexList[1];
1089 //              it_surface+=1;
1091 //              // FACE 2
1092 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][0];
1093 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][1];
1094 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][5];
1095 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][4];
1097 //              it_surface+=1;
1100 //              // FACE 3
1101 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][0];
1102 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][4];
1103 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][7];
1104 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][3];
1106 //              it_surface+=1;
1109 //              // FACE 4
1110 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][1];
1111 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][2];
1112 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][6];
1113 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][5];
1115 //              it_surface+=1;
1118 //              // FACE 5
1119 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][4];
1120 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][5];
1121 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][6];
1122 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][7];
1124 //              it_surface+=1;
1127 //              // FACE 6
1128 //              myPrimalMesh->faceToVertex[it_surface][0] = myPrimalMesh->cellToVertex[it_cell][3];
1129 //              myPrimalMesh->faceToVertex[it_surface][1] = myPrimalMesh->cellToVertex[it_cell][7];
1130 //              myPrimalMesh->faceToVertex[it_surface][2] = myPrimalMesh->cellToVertex[it_cell][6];
1131 //              myPrimalMesh->faceToVertex[it_surface][3] = myPrimalMesh->cellToVertex[it_cell][2];
1133 //              it_surface+=1;
1135 //            }
1137 //        }
1139 //  for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
1140 //    {
1141 //      debug_print( "%-40s : %ld ( ","myPrimalMesh->faceToVertex", k);
1143 //      for(connectivity_int ii=0;ii<myPrimalMesh->faceToVertexNumber[k];ii++)
1144 //        {
1145 //          debug_print( "%ld " , myPrimalMesh->faceToVertex[k][ii]);
1146 //        }
1147 //      debug_print(")\n");
1148 //    }
1151 //  endFunction(__FUNCTION__, t);
1153 //}
1155 void setHexahedreVertexToCellNumbers(struct primalMesh * myPrimalMesh)
1157   
1158   clock_t t=startFunction(__FUNCTION__);
1159   affiche("\n");
1160   
1161   myPrimalMesh->vertexToCellNumbers = (connectivity_int *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int));
1162   memset(myPrimalMesh->vertexToCellNumbers, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_int));
1163   
1164   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
1165     {
1166       for(connectivity_int i=0;i<myPrimalMesh->cellToVertexNumber[k];i++)
1167         {
1168           myPrimalMesh->vertexToCellNumbers[myPrimalMesh->cellToVertex[k][i]]++;
1169         }
1170     }
1171   
1172 #ifdef DEBUG  
1173   
1174   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
1175     {
1176       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToCellNumbers", k);
1177       
1178       {
1179         debug_print( "%ld " , myPrimalMesh->vertexToCellNumbers[k]);
1180       }
1181       debug_print(")\n");
1182     }
1183 #endif  
1184   endFunction(__FUNCTION__, t);
1185   
1191 void setHexahedreVertexToCells(struct primalMesh * myPrimalMesh)
1193   
1194   clock_t t=startFunction(__FUNCTION__);
1195   affiche("\n");
1196   
1197   connectivity_int * tmp ;
1198   
1199   myPrimalMesh->vertexToCells = (connectivity_int **) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int *));
1200   
1201   for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
1202     {
1203       myPrimalMesh->vertexToCells[i] = (connectivity_int *) malloc(myPrimalMesh->vertexToCellNumbers[i] * sizeof(connectivity_int));
1204       memset(myPrimalMesh->vertexToCells[i], 0, myPrimalMesh->vertexToCellNumbers[i]*sizeof(connectivity_int));
1205     }
1206   
1207   tmp = (connectivity_int *) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int));
1208   memset(tmp,0, myPrimalMesh->vertexNumber * sizeof(connectivity_int));
1209   
1210   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
1211     {
1212       for(connectivity_int i=0;i<myPrimalMesh->cellToVertexNumber[k];i++)
1213         {
1214           myPrimalMesh->vertexToCells[myPrimalMesh->cellToVertex[k][i]][tmp[myPrimalMesh->cellToVertex[k][i]]] = k;
1215           tmp[myPrimalMesh->cellToVertex[k][i]]++;
1216         }
1217     }
1218   
1219   
1220   
1221   
1222 #ifdef DEBUG  
1223   for(connectivity_int k=0;k<myPrimalMesh->vertexNumber;k++)
1224     {
1225       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToCells", k);
1226       
1227       for(connectivity_int ii=0;ii<myPrimalMesh->vertexToCellNumbers[k];ii++)
1228         {
1229           debug_print( "%ld " , myPrimalMesh->vertexToCells[k][ii]);
1230         }
1231       debug_print(")\n");
1232     }
1233 #endif
1234   free(tmp);
1235   
1236   endFunction(__FUNCTION__, t);
1237   
1238   
1243 void setHexahedreCellToCells(struct primalMesh * myPrimalMesh)
1245   
1246   clock_t t=startFunction(__FUNCTION__);
1247   affiche("\n");
1248   
1249   connectivity_int it_cell=0;
1250   
1251   // vertex number on faces
1252   connectivity_int cN1 =  myPrimalMesh->P*myPrimalMesh->N;
1253   
1254   //  connectivity_int it_face=0;
1255   //    connectivity_int it_segment=0;
1256   
1257   //  allocMemoryForMatrixConnectivity(&(myPrimalMesh->segmentToVertex),myPrimalMesh->segmentNumber,SEGMENTVERTEX);
1258   connectivity_int * vertexList= (connectivity_int *) malloc(SEGMENTVERTEX * sizeof(connectivity_int));
1259   
1260   myPrimalMesh->cellToCells = (connectivity_int **) malloc( (myPrimalMesh->cellNumber) * sizeof(connectivity_int *));
1261   myPrimalMesh->cellToCellsNumbers = (connectivity_int *) malloc( (myPrimalMesh->cellNumber) * sizeof(connectivity_int));
1262   
1263   connectivity_int * cellList  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024*10);
1264   
1265   for(connectivity_int i=0;i<myPrimalMesh->M;i++)
1266     for(connectivity_int j=0;j<myPrimalMesh->N;j++)
1267       for(connectivity_int k=0;k<myPrimalMesh->P;k++)
1268         {
1269           //it_surface = k + j*(myPrimalMesh->P+1) + i*(myPrimalMesh->P+1)*(myPrimalMesh->N+1);
1270           it_cell     = k + j*myPrimalMesh->P + i*cN1; //cN1 = myPrimalMesh->P*myPrimalMesh->N;
1271           //          it_face     = 0;
1272           
1273           if(myPrimalMesh->cellToVertexNumber[it_cell] == HEXAHEDRON)
1274             {
1275               
1276               connectivity_int cellNumbers    = 0 ;
1277               connectivity_int vertexSelect   = 0 ;
1278               
1279               for (connectivity_short vertexi=0;vertexi<HEXAHEDRON;vertexi++)
1280                 {
1281                   
1282                   vertexSelect = myPrimalMesh->cellToVertex[it_cell][vertexi];
1283                   
1284                   for (connectivity_short celli=0;celli<myPrimalMesh->vertexToCellNumbers[vertexSelect];celli++)
1285                     {
1286                       
1287                       if(
1288                          findInList(cellList,cellNumbers, myPrimalMesh->vertexToCells[vertexSelect][celli])<0
1289                          )
1290                         {
1291                           
1292                           cellList[cellNumbers]=myPrimalMesh->vertexToCells[vertexSelect][celli];
1293                           
1294                           cellNumbers++;
1295                           
1296                         }
1297                     }
1298                   
1299                   
1300                   /*
1301                         debug_print( "%-40s : Celli (%ld)   -> VERTEX (%7ld) -> CELL LIST ( ","cellList", it_cell, vertexSelect);
1302                         
1303                         for(connectivity_int ii=0;ii<cellNumbers;ii++)
1304                         {
1305                             debug_print("%ld ",cellList[ii]);
1306                         }
1307                         
1308                         debug_print(")\n");
1309 */
1310                   
1311                 }                        //                        cellNumbers+=jj;
1312               
1313               sortListConnectivity(cellList, cellNumbers);
1314               
1315               myPrimalMesh->cellToCells[it_cell]=(connectivity_int *) malloc(sizeof (connectivity_int)*cellNumbers);
1316               
1317               for(connectivity_int kk=0;kk<cellNumbers;kk++)
1318                 {
1319                   myPrimalMesh->cellToCells[it_cell][kk] = cellList[kk];
1320                 }
1321               
1322               myPrimalMesh->cellToCellsNumbers[it_cell] = cellNumbers;
1323               
1324             }
1325         }
1326   
1327 #ifdef DEBUG  
1328   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
1329     {
1330       debug_print( "%-40s : %ld (%ld) ( ","myPrimalMesh->cellToCells", k, myPrimalMesh->cellToCellsNumbers[k]);
1331       
1332       for(connectivity_int ii=0;ii<myPrimalMesh->cellToCellsNumbers[k];ii++)
1333         {
1334           debug_print( "%ld " , myPrimalMesh->cellToCells[k][ii]);
1335         }
1336       debug_print(")\n");
1337     }
1338 #endif
1339   
1340   //  for(connectivity_int k=0;k<myPrimalMesh->segmentNumber;k++)
1341   //    {
1342   //      debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToVertex", k);
1343   
1344   //      for(connectivity_int ii=0;ii<SEGMENTVERTEX;ii++)
1345   //        {
1346   //          debug_print( "%ld " , myPrimalMesh->segmentToVertex[k][ii]);
1347   //        }
1348   
1349   //      debug_print(")\n");
1350   //    }
1351   
1352   free(cellList);
1353   free(vertexList);
1354   
1355   endFunction(__FUNCTION__, t);
1356   
1362 void setHexahedreCellToFacesOwnerNeighbour(struct primalMesh * myPrimalMesh)
1364   
1365   clock_t t=startFunction(__FUNCTION__);
1366   affiche("\n");
1367   
1368   connectivity_int position;
1369   connectivity_int selectedCell;
1370   connectivity_int numberOfVertex;
1371   
1372   connectivity_int * vertexList  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024*10);
1373   connectivity_int * vertexListBase  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
1374   connectivity_int * vertexListInverted  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
1375   connectivity_int * vertexListInvertedTwo  = (connectivity_int * ) malloc(sizeof (connectivity_int)*1024);
1376   
1377   connectivity_int * cellTest  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1378   memset(cellTest,0,sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1379   
1380   char ** cellFaceTest  = (char ** ) malloc(sizeof (char*)*myPrimalMesh->cellNumber);
1381   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1382     {
1383       cellFaceTest[celli]  = (char * ) malloc(sizeof (char)*HEXAHEDRON_FACES);
1384       memset(cellFaceTest[celli],0,sizeof (char)*HEXAHEDRON_FACES);
1385     }
1386   
1387   myPrimalMesh->cellToFacesOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
1388   for (connectivity_int ii=0;ii<myPrimalMesh->cellNumber;ii++) {
1389       myPrimalMesh->cellToFacesOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1390     }
1391   
1392   myPrimalMesh->cellToFacesNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
1393   for (connectivity_int ii=0;ii<myPrimalMesh->cellNumber;ii++) {
1394       myPrimalMesh->cellToFacesNeighbour[ii]  = NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1395     }
1396   
1397   myPrimalMesh->cellToFacesOwnerNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1398   memset(myPrimalMesh->cellToFacesOwnerNumber,0,myPrimalMesh->cellNumber*sizeof (connectivity_int) );
1399   
1400   myPrimalMesh->cellToFacesNeighbourNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1401   memset(myPrimalMesh->cellToFacesNeighbourNumber,0,myPrimalMesh->cellNumber*sizeof (connectivity_int) );
1402   
1403   
1404   //    myPrimalMesh->cellToFacesNeighbour = realloc(myPrimalMesh->cellToFacesNeighbour,sizeof (connectivity_int*)*(myPrimalMesh->cellNumber));
1405   //    myPrimalMesh->cellToFacesOwner = realloc(myPrimalMesh->cellToFacesOwner,sizeof (connectivity_int*)*(myPrimalMesh->cellNumber));
1406   
1407   connectivity_int facesNeighbourCount = 0;
1408   connectivity_int facesOwnerCount = 0;
1409   
1410   
1411   for(connectivity_int cellAi=0;cellAi<myPrimalMesh->cellNumber;cellAi++)
1412     {
1413       if(cellAi % 100000 == 0 && cellAi != 0)
1414         release_print( "%-40s :  (%ld)\n","CELLS", cellAi);
1415       
1416       for (connectivity_int faceAi=0;faceAi<HEXAHEDRON_FACES;faceAi++)
1417         {
1418           
1419           for (connectivity_int vertexAi=0;vertexAi<hexaedron_localNodeListNumbers[faceAi];vertexAi++)
1420             {
1421               vertexListInverted[vertexAi] = myPrimalMesh->cellToVertex[cellAi][ hexaedron_localNodeList[faceAi][hexaedron_localNodeListNumbers[faceAi] -1 - vertexAi] ];
1422               vertexListBase[vertexAi] = myPrimalMesh->cellToVertex[cellAi][ hexaedron_localNodeList[faceAi][vertexAi] ];
1423             }
1424           
1425           position = findMinInList(vertexListInverted,hexaedron_localNodeListNumbers[faceAi]);
1426           sortList(&vertexListInverted, hexaedron_localNodeListNumbers[faceAi], position);
1427           
1428           position = findMinInList(vertexListBase,hexaedron_localNodeListNumbers[faceAi]);
1429           sortList(&vertexListBase, hexaedron_localNodeListNumbers[faceAi], position);
1430           
1431           //            if( cellTest[celli]==0)
1432           {
1433             
1434             for(connectivity_int cellBi=0;cellBi<myPrimalMesh->cellToCellsNumbers[cellAi];cellBi++)
1435               {
1436                 
1437                 selectedCell    = myPrimalMesh->cellToCells[cellAi][cellBi];
1438                 numberOfVertex  = myPrimalMesh->cellToVertexNumber[selectedCell];
1439                 
1440                 
1441                 for (connectivity_int faceBi=0;faceBi<HEXAHEDRON_FACES;faceBi++)
1442                   {
1443                     
1444                     for (connectivity_int vertexBi=0;vertexBi<hexaedron_localNodeListNumbers[faceBi];vertexBi++)
1445                       {
1446                         vertexList[vertexBi] = myPrimalMesh->cellToVertex[selectedCell][ hexaedron_localNodeList[faceBi][vertexBi] ];
1447                       }
1448                     
1449                     position = findMinInList(vertexList,hexaedron_localNodeListNumbers[faceBi]);
1450                     sortList(&vertexList, hexaedron_localNodeListNumbers[faceBi], position);
1451                     
1452                     if(memcmp(vertexListBase,vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int)) == 0)
1453                       {
1454                         
1455                         if(cellFaceTest[selectedCell][faceBi]==0)
1456                           {
1457                             cellFaceTest[selectedCell][faceBi]=1;
1458                             
1459                             facesOwnerCount = (connectivity_int) myPrimalMesh->cellToFacesOwnerNumber[selectedCell];
1460                             
1461                             myPrimalMesh->faces = (connectivity_int**) realloc(myPrimalMesh->faces,sizeof (connectivity_int*)*(myPrimalMesh->faceNumber+1));
1462                             myPrimalMesh->faces[myPrimalMesh->faceNumber] =(connectivity_int*)  malloc(hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int) );
1463                             
1464                             memcpy( myPrimalMesh->faces[myPrimalMesh->faceNumber], vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int) );
1465                             
1466                             myPrimalMesh->cellToFacesOwner[selectedCell] = (connectivity_int *) realloc(myPrimalMesh->cellToFacesOwner[selectedCell], (facesOwnerCount+1)*sizeof (connectivity_int) );
1467                             
1468                             myPrimalMesh->cellToFacesOwner[selectedCell][facesOwnerCount] = myPrimalMesh->faceNumber;
1469                             myPrimalMesh->cellToFacesOwnerNumber[selectedCell]++;
1470                             
1471                             myPrimalMesh->faceNumber++;
1472                             
1473                           }
1474                       }
1475                     else
1476                       {
1477                         if(memcmp(vertexListInverted,vertexList,hexaedron_localNodeListNumbers[faceBi]*sizeof (connectivity_int)) == 0)
1478                           {
1479                             if(cellFaceTest[selectedCell][faceBi]==0)
1480                               {
1481                                 cellFaceTest[selectedCell][faceBi]=-1;
1482                                 
1483                                 facesNeighbourCount = (connectivity_int) myPrimalMesh->cellToFacesNeighbourNumber[selectedCell];
1484                                 
1485                                 //                                    debug_print("%ld\n",facesNeighbourCount);
1486                                 
1487                                 if(myPrimalMesh->cellToFacesNeighbour[selectedCell]==NULL)
1488                                   myPrimalMesh->cellToFacesNeighbour[selectedCell] = (connectivity_int*) malloc((facesNeighbourCount+1)*sizeof (connectivity_int) );
1489                                 else
1490                                   myPrimalMesh->cellToFacesNeighbour[selectedCell] = (connectivity_int*) realloc(myPrimalMesh->cellToFacesNeighbour[selectedCell], ((facesNeighbourCount+1)) *sizeof (connectivity_int) );
1491                                 
1492                                 connectivity_short testFaceCi = 0;
1493                                 
1494                                 //                                for( connectivity_int faceCi= 0; (faceCi<myPrimalMesh->faceNumber) && (testFaceCi == 0) ; faceCi++ )
1495                                 for( connectivity_int faceCi= (connectivity_int) myPrimalMesh->faceNumber-1; (faceCi>0) && (testFaceCi == 0) ; faceCi-- )
1496                                   {
1497                                     
1498                                     if(memcmp(vertexListBase,myPrimalMesh->faces[faceCi],QUAD*sizeof (connectivity_int))==0)
1499                                       {
1500                                         testFaceCi = 1;
1501                                         myPrimalMesh->cellToFacesNeighbour[selectedCell][facesNeighbourCount]=faceCi;
1502                                         myPrimalMesh->cellToFacesNeighbourNumber[selectedCell]++;
1503                                       }
1504                                   }
1505                                 
1506                               }
1507                           }
1508                       }
1509                   }
1510                 
1511               }
1512           }
1513           
1514           
1515           
1516           
1517           
1518           
1519         }
1520     }
1521   
1522   free(vertexList);
1523   free(vertexListBase);
1524   free(vertexListInverted);
1525   free(cellTest);
1526   free(vertexListInvertedTwo);
1527   
1528   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1529     free(cellFaceTest[celli]);
1530   free(cellFaceTest);
1531   
1532 #ifdef DEBUG
1533   
1534   for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
1535     {
1536       debug_print( "%-40s : %ld ( ","myPrimalMesh->faces", k);
1537       
1538       for(connectivity_int ii=0;ii<QUAD;ii++)
1539         {
1540           debug_print( "%ld " , myPrimalMesh->faces[k][ii]);
1541         }
1542       debug_print(")\n");
1543     }
1544   
1545   
1546   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1547     {
1548       debug_print( "%-40s : %ld ( ","myPrimalMesh->cellToFacesOwner", celli);
1549       
1550       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
1551         {
1552           debug_print( "%ld " , myPrimalMesh->cellToFacesOwner[celli][faceAi]);
1553         }
1554       debug_print(")\n");
1555     }
1556   
1557   
1558   for(connectivity_int k=0;k<myPrimalMesh->cellNumber;k++)
1559     {
1560       debug_print( "%-40s : %ld ( ","myPrimalMesh->cellToFacesNeighbour", k);
1561       
1562       for(connectivity_int ii=0;ii<myPrimalMesh->cellToFacesNeighbourNumber[k];ii++)
1563         {
1564           debug_print( "%ld " , myPrimalMesh->cellToFacesNeighbour[k][ii] );
1565         }
1566       debug_print(")\n");
1567     }
1568   
1569 #endif
1570   
1571   endFunction(__FUNCTION__, t);
1572   
1575 void setHexahedreFaceToCells(struct primalMesh * myPrimalMesh)
1577   
1578   clock_t t=startFunction(__FUNCTION__);
1579   affiche("\n");
1580   
1581   connectivity_int selectedFace = 0;
1582   
1583   
1584   myPrimalMesh->faceToCells  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->faceNumber);
1585   for (connectivity_int ii=0;ii<myPrimalMesh->faceNumber;ii++) {
1586       myPrimalMesh->faceToCells[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
1587     }
1588   
1589   myPrimalMesh->faceToCellsNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->faceNumber);
1590   memset(myPrimalMesh->faceToCellsNumber,0,myPrimalMesh->faceNumber*sizeof (connectivity_int) );
1591   
1592   
1593   
1594   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1595     {
1596       
1597       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesOwnerNumber[celli];faceAi++)
1598         {
1599           //          myPrimalMesh->cellToFacesOwner[celli][faceAi];
1600           //          debug_print("O: %ld %ld\n",celli, myPrimalMesh->cellToFacesOwner[celli][faceAi]);
1601           
1602           selectedFace = myPrimalMesh->cellToFacesOwner[celli][faceAi];
1603           
1604           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
1605           
1606           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
1607           myPrimalMesh->faceToCellsNumber[selectedFace]++;
1608           
1609         }
1610       
1611       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->cellToFacesNeighbourNumber[celli];faceAi++)
1612         {
1613           //          myPrimalMesh->cellToFacesOwner[celli][faceAi];
1614           //          debug_print("O: %ld %ld\n",celli, myPrimalMesh->cellToFacesNeighbour[celli][faceAi]);
1615           
1616           selectedFace = myPrimalMesh->cellToFacesNeighbour[celli][faceAi];
1617           
1618           myPrimalMesh->faceToCells[selectedFace] = (connectivity_int*) realloc(myPrimalMesh->faceToCells[selectedFace],sizeof (connectivity_int)*(myPrimalMesh->faceToCellsNumber[selectedFace]+1));
1619           
1620           myPrimalMesh->faceToCells[selectedFace][myPrimalMesh->faceToCellsNumber[selectedFace]] = celli;
1621           myPrimalMesh->faceToCellsNumber[selectedFace]++;
1622           
1623           
1624         }
1625       
1626       
1627     }
1628   
1629   
1630 #ifdef DEBUG
1631   for(connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1632     {
1633       debug_print( "%-40s : %ld ( ","myPrimalMesh->faceToCells", faceAi);
1634       
1635       for(connectivity_int celli=0;celli<myPrimalMesh->faceToCellsNumber[faceAi];celli++)
1636         {
1637           debug_print( "%ld " , myPrimalMesh->faceToCells[faceAi][celli] );
1638         }
1639       debug_print(")\n");
1640     }
1641   
1642 #endif
1643   
1644   //  free(copyFace);
1645   //  free(vertexList);
1646   //  free(vertexListInverted);
1647   
1648   
1649   endFunction(__FUNCTION__, t);
1650   
1655 void setHexahedreSegments(struct primalMesh * myPrimalMesh)
1657   
1658   clock_t t=startFunction(__FUNCTION__);
1659   affiche("\n");
1660   //  affiche("%s : REFAIRE CETTE FONCTION EN UTILISANT CELLTOFACE PUIS REDUIRE LA RECHERCHE DANS SEGMENTS EN INTRODUISANT PAR EXEMPLE CELLTOSEGMENTS\n",__FUNCTION__);
1661   
1662   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1663   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1664   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1665   
1666   connectivity_int testSegment = 0;
1667   
1668   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(0));
1669   
1670   /// MODIF
1671   
1672   char ** cellSegmentTest  = (char ** ) malloc(sizeof (char*)*myPrimalMesh->cellNumber);
1673   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1674     {
1675       cellSegmentTest[celli]  = (char * ) malloc(sizeof (char)*HEXAHEDRON_SEGMENTS);
1676       memset(cellSegmentTest[celli],0,sizeof (char)*HEXAHEDRON_SEGMENTS);
1677     }
1678   
1679   myPrimalMesh->cellToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
1680   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1681     {
1682       myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
1683     }
1684   
1685   myPrimalMesh->cellToSegmentOwnerNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1686   memset(myPrimalMesh->cellToSegmentOwnerNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_int));
1687   
1688   myPrimalMesh->cellToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->cellNumber);
1689   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1690     {
1691       myPrimalMesh->cellToSegmentNeighbour[celli]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*0);
1692     }
1693   
1694   myPrimalMesh->cellToSegmentNeighbourNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->cellNumber);
1695   memset(myPrimalMesh->cellToSegmentNeighbourNumber, 0, myPrimalMesh->cellNumber*sizeof(connectivity_int));
1696   
1697   /// MODIF
1698   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1699     {
1700       for (connectivity_int segmentAi=0;segmentAi<HEXAHEDRON_SEGMENTS;segmentAi++)
1701         { 
1702           if(cellSegmentTest[celli][segmentAi] ==0)
1703             {
1704               connectivity_int localPointA1 = hexaedron_localSegmentList[segmentAi][PT1];
1705               connectivity_int localPointA2 = hexaedron_localSegmentList[segmentAi][PT2];
1706               
1707               cellSegmentTest[celli][segmentAi] = 1;
1708               
1709               myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->segmentNumber+1));
1710               myPrimalMesh->segments[myPrimalMesh->segmentNumber]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1711               
1712               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT1] = myPrimalMesh->cellToVertex[celli][localPointA1];
1713               myPrimalMesh->segments[myPrimalMesh->segmentNumber][PT2] = myPrimalMesh->cellToVertex[celli][localPointA2];
1714               
1715               myPrimalMesh->segmentNumber++;
1716               
1717               myPrimalMesh->cellToSegmentOwner[celli]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentOwner[celli], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentOwnerNumber[celli]+1));
1718               
1719               myPrimalMesh->cellToSegmentOwner[celli][myPrimalMesh->cellToSegmentOwnerNumber[celli]] = myPrimalMesh->segmentNumber -1;
1720               
1721               myPrimalMesh->cellToSegmentOwnerNumber[celli]++;
1722               
1723               if(myPrimalMesh->segmentNumber % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && myPrimalMesh->segmentNumber != 0)
1724                 release_print( "%-40s :  (%ld)\n","SEGMENTS", myPrimalMesh->segmentNumber);
1725               
1726               for (connectivity_int cellAi=0;cellAi<myPrimalMesh->cellToCellsNumbers[celli];cellAi++)
1727                 {      
1728                   connectivity_int selectedCell = myPrimalMesh->cellToCells[celli][cellAi];
1729                   
1730                   
1731                   if(celli!=selectedCell) {
1732                       
1733                       //                      affiche("selectedCell %ld:%ld\n",celli,selectedCell);
1734                       
1735                       for (connectivity_int segmentBi=0;segmentBi<HEXAHEDRON_SEGMENTS;segmentBi++)
1736                         { 
1737                           if(cellSegmentTest[selectedCell][segmentBi] ==0)
1738                             {
1739                               
1740                               connectivity_int localPointB1 = hexaedron_localSegmentList[segmentBi][PT1];
1741                               connectivity_int localPointB2 = hexaedron_localSegmentList[segmentBi][PT2];
1742                               
1743                               if(
1744                                  myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
1745                                  &&
1746                                  myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
1747                                  )
1748                                 {
1749                                   //                                  affiche("A: %ld %ld\n",segmentAi,segmentBi);
1750                                   cellSegmentTest[selectedCell][segmentBi] = 1;
1751                                   
1752                                 }
1753                               else 
1754                                 if(
1755                                    myPrimalMesh->cellToVertex[celli][localPointA1] == myPrimalMesh->cellToVertex[selectedCell][localPointB2] 
1756                                    &&
1757                                    myPrimalMesh->cellToVertex[celli][localPointA2] == myPrimalMesh->cellToVertex[selectedCell][localPointB1] 
1758                                    )
1759                                   {
1760                                     //                                    affiche("B: %ld %ld\n",segmentAi,segmentBi);
1761                                     cellSegmentTest[selectedCell][segmentBi] = -1;
1762                                   }
1763                               
1764                               
1765                               if(cellSegmentTest[selectedCell][segmentBi] !=0)
1766                                 {
1767                                   
1768                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell]  = (connectivity_int * ) realloc( myPrimalMesh->cellToSegmentNeighbour[selectedCell], sizeof (connectivity_int)*(myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]+1));
1769                                   
1770                                   myPrimalMesh->cellToSegmentNeighbour[selectedCell][myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]] = myPrimalMesh->segmentNumber -1;
1771                                   
1772                                   myPrimalMesh->cellToSegmentNeighbourNumber[selectedCell]++;
1773                                   
1774                                 }
1775                               
1776                             }
1777                           
1778                         }
1779                     }  
1780                   else 
1781                     {
1782                       //                      affiche("Noeud voisin \n");
1783                     }
1784                   
1785                 }
1786             }
1787         }
1788       
1789     }
1790   
1791 #ifdef DEBUG  
1792   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1793     {
1794       debug_print( "%-40s : %ld ( ","cellToSegmentOwner", celli);
1795       
1796       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli];segmenti++)
1797         {
1798           debug_print( "%ld " , myPrimalMesh->cellToSegmentOwner[celli][segmenti] );
1799         }
1800       debug_print(")\n");
1801     }  
1802   
1803   
1804   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1805     {
1806       debug_print( "%-40s : %ld ( ","cellToSegmentNeighbour", celli);
1807       
1808       for(connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli];segmenti++)
1809         {
1810           debug_print( "%ld " , myPrimalMesh->cellToSegmentNeighbour[celli][segmenti] );
1811         }
1812       debug_print(")\n");
1813     }  
1814   
1815   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1816     {
1817       debug_print( "%-40s : %ld ( ","myPrimalMesh->segments", segmentAi);
1818       
1819       for(connectivity_int vertexAi=0;vertexAi<SEGMENTVERTEX;vertexAi++)
1820         {
1821           debug_print( "%ld " , myPrimalMesh->segments[segmentAi][vertexAi] );
1822         }
1823       debug_print(")\n");
1824     }
1825   
1826 #endif
1827   
1828   free(copyFace);
1829   free(vertexList);
1830   free(vertexListInverted);
1831   
1832   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
1833     free(cellSegmentTest[celli]);
1834   free(cellSegmentTest);  
1835   
1836   
1837   
1838   
1839   endFunction(__FUNCTION__, t);
1840   
1845 void setHexahedreSegmentsCentres(struct primalMesh * myPrimalMesh)
1847   
1848   clock_t t=startFunction(__FUNCTION__);
1849   affiche("\n");
1850   
1851   myPrimalMesh->segmentCentres = (dataType ** ) malloc(sizeof (dataType*)*myPrimalMesh->segmentNumber);
1852   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1853     {
1854       myPrimalMesh->segmentCentres[segmenti]  = (dataType * ) malloc(sizeof (dataType)*DIM3D);
1855     }
1856   
1857   
1858   for (connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
1859     {
1860       
1861       if(segmenti % 100000 == 0 && segmenti != 0)
1862         release_print( "%-40s :  (%ld)\n","SEGMENT", segmenti);
1863       
1864       myPrimalMesh->segmentCentres[segmenti][0] = (
1865             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][0]
1866           +
1867           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][0]
1868           )/2.0;
1869       myPrimalMesh->segmentCentres[segmenti][1] = (
1870             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][1]
1871           +
1872           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][1]
1873           )/2.0;
1874       myPrimalMesh->segmentCentres[segmenti][2] = (
1875             myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT1] ][2]
1876           +
1877           myPrimalMesh->vertex[ myPrimalMesh->segments[segmenti][PT2] ][2]
1878           )/2.0;
1879     }
1880   
1881 #ifdef DEBUG
1882   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1883     {
1884       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentCentres", segmentAi);
1885       
1886       for(connectivity_int vertexAi=0;vertexAi<DIM3D;vertexAi++)
1887         {
1888           debug_print( "%lf " , myPrimalMesh->segmentCentres[segmentAi][vertexAi] );
1889         }
1890       debug_print(")\n");
1891     }
1892   
1893 #endif  
1894   
1895   endFunction(__FUNCTION__, t);
1896   
1901 void setHexahedreSegments2(struct primalMesh * myPrimalMesh)
1903   
1904   clock_t t=startFunction(__FUNCTION__);
1905   affiche("\n");
1906   affiche("%s : Nouvelle fonction setHexahedreSegments pour réduire le temps de calcul\n",__FUNCTION__);
1907   
1908   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1909   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1910   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
1911   
1912   connectivity_int testSegment = 0;
1913   
1914   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->cellNumber*6));
1915   for (connectivity_int segmentAi=0;segmentAi<(myPrimalMesh->cellNumber*6);segmentAi++)
1916     myPrimalMesh->segments[segmentAi]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1917   
1918   for (connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
1919     //  for (long faceAi=(long)(myPrimalMesh->faceNumber-1);faceAi>=0;faceAi--)
1920     {
1921       
1922       if(faceAi % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && faceAi != 0)
1923         release_print( "%-40s :  (%ld)\n","FACE", faceAi);
1924       
1925       
1926       memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
1927       copyFace[QUAD] = copyFace[0];
1928       
1929       for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
1930         {
1931           
1932           vertexList[0] = copyFace[vertexAi];
1933           vertexList[1] = copyFace[(vertexAi+1) % QUAD];
1934           
1935           vertexListInverted[0] = vertexList[1];
1936           vertexListInverted[1] = vertexList[0];
1937           
1938           testSegment = 0;
1939           
1940           /*
1941            * AMELIORATION du temps de calcul :
1942            * Augmenter la vitesse de recherche
1943            * en incluant une fonction faceToCells puis cellToSegments pour réduire l'impact de la recherche
1944            */
1945           for (connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber  && testSegment == 0;segmentAi++)
1946             {
1947               
1948               if( vertexList[0] == myPrimalMesh->segments[segmentAi][0])
1949                 {
1950                   if( vertexList[1] == myPrimalMesh->segments[segmentAi][1])
1951                     {
1952                       testSegment = 1;
1953                     }
1954                 }
1955               else
1956                 if( vertexListInverted[0] == myPrimalMesh->segments[segmentAi][0])
1957                   {
1958                     if( vertexListInverted[1] == myPrimalMesh->segments[segmentAi][1])
1959                       
1960                       {
1961                         testSegment = 2;
1962                       }
1963                   }
1964             }
1965           
1966           if(testSegment == 0)
1967             {
1968               if(myPrimalMesh->cellNumber*HEXAHEDRON_FACES<=myPrimalMesh->segmentNumber)
1969                 {
1970                   myPrimalMesh->segments  = (connectivity_int ** ) realloc(myPrimalMesh->segments, sizeof (connectivity_int*)*(myPrimalMesh->segmentNumber+1));
1971                   myPrimalMesh->segments[myPrimalMesh->segmentNumber]  = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
1972                 }
1973               memcpy(myPrimalMesh->segments[myPrimalMesh->segmentNumber], vertexList, sizeof(connectivity_int)*SEGMENTVERTEX);
1974               
1975               myPrimalMesh->segmentNumber++;
1976             }
1977           
1978         }
1979     }
1980   
1981 #ifdef DEBUG
1982   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
1983     {
1984       debug_print( "%-40s : %ld ( ","myPrimalMesh->segments", segmentAi);
1985       
1986       for(connectivity_int vertexAi=0;vertexAi<SEGMENTVERTEX;vertexAi++)
1987         {
1988           debug_print( "%ld " , myPrimalMesh->segments[segmentAi][vertexAi] );
1989         }
1990       debug_print(")\n");
1991     }
1992   
1993 #endif
1994   
1995   free(copyFace);
1996   free(vertexList);
1997   free(vertexListInverted);
1998   
1999   
2000   endFunction(__FUNCTION__, t);
2001   
2007 void setHexahedreVertexToSegments(struct primalMesh * myPrimalMesh)
2009   
2010   clock_t t=startFunction(__FUNCTION__);
2011   affiche("\n");
2012   
2013   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2014   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2015   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
2016   
2017   //connectivity_int testSegment = 0;
2018   
2019   myPrimalMesh->vertexToSegmentOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
2020   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
2021       myPrimalMesh->vertexToSegmentOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2022     }
2023   
2024   myPrimalMesh->vertexToSegmentOwnerNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->vertexNumber);
2025   memset(myPrimalMesh->vertexToSegmentOwnerNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_int));
2026   
2027   myPrimalMesh->vertexToSegmentNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->vertexNumber);
2028   for (connectivity_int ii=0;ii<myPrimalMesh->vertexNumber;ii++) {
2029       myPrimalMesh->vertexToSegmentNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2030     }
2031   
2032   myPrimalMesh->vertexToSegmentNeighbourNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->vertexNumber);
2033   memset(myPrimalMesh->vertexToSegmentNeighbourNumber, 0, myPrimalMesh->vertexNumber*sizeof(connectivity_int));
2034   
2035   
2036   connectivity_int pointOwner = 0;
2037   connectivity_int pointNeighbour = 0;
2038   
2039   for (connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
2040     {
2041       if(segmentAi % (int)floor(myPrimalMesh->cellNumber*10/NBAFFICHE) == 0 && segmentAi!=0)
2042         release_print( "%-40s :  (%ld)\n","SEGMENTS", segmentAi);
2043       
2044       pointOwner = myPrimalMesh->segments[segmentAi][1];
2045       
2046       myPrimalMesh->vertexToSegmentOwner[pointOwner] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentOwner[pointOwner],(myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]+1)*sizeof (connectivity_int));
2047       myPrimalMesh->vertexToSegmentOwner[pointOwner][myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]] = segmentAi;
2048       myPrimalMesh->vertexToSegmentOwnerNumber[pointOwner]++;
2049       
2050       pointNeighbour = myPrimalMesh->segments[segmentAi][0];
2051       
2052       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour] = (connectivity_int *) realloc(myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour],(myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]+1)*sizeof (connectivity_int));
2053       myPrimalMesh->vertexToSegmentNeighbour[pointNeighbour][myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]] = segmentAi;
2054       myPrimalMesh->vertexToSegmentNeighbourNumber[pointNeighbour]++;
2055       
2056     }
2057   
2058 #ifdef DEBUG
2059   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
2060     {
2061       
2062       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentOwner", vertexAi);
2063       
2064       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentOwnerNumber[vertexAi];segmentAi++)
2065         {
2066           debug_print( "%ld " , myPrimalMesh->vertexToSegmentOwner[vertexAi][segmentAi] );
2067         }
2068       
2069       debug_print(")\n");
2070       
2071     }
2072   
2073   for(connectivity_int vertexAi=0;vertexAi<myPrimalMesh->vertexNumber;vertexAi++)
2074     {
2075       
2076       debug_print( "%-40s : %ld ( ","myPrimalMesh->vertexToSegmentNeighbour", vertexAi);
2077       
2078       for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->vertexToSegmentNeighbourNumber[vertexAi];segmentAi++)
2079         {
2080           debug_print( "%ld " , myPrimalMesh->vertexToSegmentNeighbour[vertexAi][segmentAi] );
2081         }
2082       
2083       debug_print(")\n");
2084       
2085     }
2086 #endif
2087   
2088   free(copyFace);
2089   free(vertexList);
2090   free(vertexListInverted);
2091   
2092   
2093   endFunction(__FUNCTION__, t);
2094   
2098 void setHexahedreSegmentToFaces(struct primalMesh * myPrimalMesh)
2100   
2101   clock_t t=startFunction(__FUNCTION__);
2102   affiche("\n");
2103   
2104   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2105   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2106   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
2107   connectivity_int segmentAi = 0;
2108   connectivity_int testSegment=0;
2109   connectivity_int faceAi =0;
2110   
2111   //connectivity_int testSegment = 0;
2112   
2113   //  affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
2114   
2115   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
2116   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
2117       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2118     }
2119   
2120   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
2121   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_int));
2122   
2123   
2124   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
2125   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
2126       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2127     }
2128   
2129   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
2130   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_int));
2131   
2132   
2133   for (connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
2134     {
2135       if(celli % 100000 == 0 && celli != 0)
2136         release_print( "%-40s :  (%ld)\n","CELL", celli);
2137       
2138       
2139       
2140       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesOwnerNumber[celli];facei++)
2141         {
2142           faceAi = myPrimalMesh->cellToFacesOwner[celli][facei];
2143           
2144           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
2145           copyFace[QUAD] = copyFace[0];  
2146           
2147           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
2148             {
2149               testSegment=0;
2150               
2151               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
2152                 {
2153                   
2154                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
2155                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
2156                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
2157                   
2158                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2159                     {
2160                       
2161                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
2162                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
2163                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
2164                       testSegment=1;
2165                     }
2166                   else
2167                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2168                       {
2169                         
2170                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
2171                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
2172                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
2173                         testSegment=2;
2174                       }
2175                   
2176                 }
2177               
2178               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
2179                 {
2180                   
2181                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
2182                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
2183                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
2184                   
2185                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2186                     {
2187                       
2188                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
2189                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
2190                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
2191                       testSegment=1;
2192                     }
2193                   else
2194                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2195                       {
2196                         
2197                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
2198                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
2199                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
2200                         testSegment=2;
2201                       }
2202                   
2203                   
2204                   
2205                 }
2206             }
2207         }
2208       
2209       for (connectivity_int facei=0;facei<myPrimalMesh->cellToFacesNeighbourNumber[celli];facei++)
2210         {
2211           faceAi = myPrimalMesh->cellToFacesNeighbour[celli][facei];
2212           
2213           memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
2214           copyFace[QUAD] = copyFace[0];  
2215           
2216           for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
2217             {
2218               testSegment=0;
2219               
2220               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentOwnerNumber[celli] && testSegment==0;segmenti++)
2221                 {
2222                   
2223                   segmentAi = myPrimalMesh->cellToSegmentOwner[celli][segmenti];
2224                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
2225                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
2226                   
2227                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2228                     {
2229                       
2230                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
2231                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
2232                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
2233                       testSegment=1;
2234                     }
2235                   else
2236                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2237                       {
2238                         
2239                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
2240                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
2241                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
2242                         testSegment=2;
2243                       }
2244                 }
2245               
2246               for (connectivity_int segmenti=0;segmenti<myPrimalMesh->cellToSegmentNeighbourNumber[celli] && testSegment==0;segmenti++)
2247                 {
2248                   segmentAi = myPrimalMesh->cellToSegmentNeighbour[celli][segmenti];
2249                   vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
2250                   vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
2251                   
2252                   if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2253                     {
2254                       
2255                       myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
2256                       myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
2257                       myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
2258                       testSegment=1;
2259                     }
2260                   else
2261                     if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2262                       {
2263                         
2264                         myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
2265                         myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
2266                         myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
2267                         testSegment=2;
2268                       }
2269                   
2270                   
2271                 }
2272             }
2273         }
2274       
2275       
2276       
2277       
2278       
2279     }
2280   
2281 #ifdef DEBUG
2282   
2283   
2284   
2285   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
2286     {
2287       
2288       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
2289       
2290       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
2291         {
2292           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
2293         }
2294       
2295       debug_print(")\n");
2296       
2297     }
2298   
2299   
2300   
2301   
2302   
2303   
2304   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
2305     {
2306       
2307       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
2308       
2309       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
2310         {
2311           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
2312         }
2313       
2314       debug_print(")\n");
2315       
2316     }
2317   
2318 #endif
2319   
2320   free(copyFace);
2321   free(vertexList);
2322   free(vertexListInverted);
2323   
2324   
2325   endFunction(__FUNCTION__, t);
2326   
2331 void setHexahedreSegmentToFaces2(struct primalMesh * myPrimalMesh)
2333   
2334   clock_t t=startFunction(__FUNCTION__);
2335   affiche("\n");
2336   
2337   connectivity_int * vertexList         = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2338   connectivity_int * vertexListInverted = (connectivity_int * ) malloc(sizeof (connectivity_int)*SEGMENTVERTEX);
2339   connectivity_int * copyFace           = (connectivity_int * ) malloc(sizeof (connectivity_int)*(QUAD+1));
2340   
2341   //connectivity_int testSegment = 0;
2342   
2343   affiche("%s : REDUIRE LE TEMPS DE CETTE FONCTION EN DEFINISSANT DE NOUVELLES MATRICES avec FACETOCELL puis CELLTOSEGMENTS \n",__FUNCTION__);
2344   
2345   myPrimalMesh->segmentToFaceOwner  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
2346   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
2347       myPrimalMesh->segmentToFaceOwner[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2348     }
2349   
2350   myPrimalMesh->segmentToFaceOwnerNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
2351   memset(myPrimalMesh->segmentToFaceOwnerNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_int));
2352   
2353   
2354   myPrimalMesh->segmentToFaceNeighbour  = (connectivity_int ** ) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
2355   for (connectivity_int ii=0;ii<myPrimalMesh->segmentNumber;ii++) {
2356       myPrimalMesh->segmentToFaceNeighbour[ii]  =  NULL;//(connectivity_int * ) malloc(sizeof (connectivity_int*)*0);
2357     }
2358   
2359   myPrimalMesh->segmentToFaceNeighbourNumber  = (connectivity_int * ) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
2360   memset(myPrimalMesh->segmentToFaceNeighbourNumber, 0, myPrimalMesh->segmentNumber*sizeof(connectivity_int));
2361   
2362   
2363   for (connectivity_int faceAi=0;faceAi<myPrimalMesh->faceNumber;faceAi++)
2364     {
2365       if(faceAi % ((int)(floor(myPrimalMesh->cellNumber*10/NBAFFICHE))) == 0 && faceAi!=0)
2366         release_print( "%-40s :  (%ld)\n","FACES", faceAi);
2367       
2368       memcpy(copyFace,myPrimalMesh->faces[faceAi],sizeof (connectivity_int)*(QUAD));
2369       copyFace[QUAD] = copyFace[0];
2370       
2371       for (connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
2372         {
2373           connectivity_int testSegment=0;
2374           /*
2375            * AMELIORATION du temps de calcul :
2376            * Augmenter la vitesse de recherche
2377            * en incluant une fonction faceToCells puis cellToSegments pour réduire l'impact de la recherche
2378            * de segmentToFaceOwner et segmentToFaceNeighbour
2379            */
2380           for (connectivity_int segmentAi=0 ; segmentAi<myPrimalMesh->segmentNumber && testSegment==0 ;segmentAi++)
2381             {
2382               
2383               vertexListInverted[0]=myPrimalMesh->segments[segmentAi][1];
2384               vertexListInverted[1]=myPrimalMesh->segments[segmentAi][0];
2385               
2386               if(memcmp(&(copyFace[vertexAi]),(myPrimalMesh->segments[segmentAi]),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2387                 {
2388                   
2389                   myPrimalMesh->segmentToFaceOwner[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceOwner[segmentAi], (myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]+1) * sizeof (connectivity_int));
2390                   myPrimalMesh->segmentToFaceOwner[segmentAi][myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]] = faceAi;
2391                   myPrimalMesh->segmentToFaceOwnerNumber[segmentAi]++;
2392                   testSegment=1;
2393                 }
2394               else
2395                 if(memcmp(&(copyFace[vertexAi]),(vertexListInverted),sizeof (connectivity_int)*SEGMENTVERTEX)==0)
2396                   {
2397                     
2398                     myPrimalMesh->segmentToFaceNeighbour[segmentAi] = (connectivity_int *) realloc(myPrimalMesh->segmentToFaceNeighbour[segmentAi], (myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]+1)*sizeof (connectivity_int));
2399                     myPrimalMesh->segmentToFaceNeighbour[segmentAi][myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]] = faceAi;
2400                     myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi]++;
2401                     testSegment=2;
2402                   }
2403             }
2404         }
2405     }
2406   
2407   
2408   
2409 #ifdef DEBUG
2410   
2411   
2412   
2413   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
2414     {
2415       
2416       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceOwner", segmentAi);
2417       
2418       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceOwnerNumber[segmentAi];faceAi++)
2419         {
2420           debug_print( "%ld " , myPrimalMesh->segmentToFaceOwner[segmentAi][faceAi] );
2421         }
2422       
2423       debug_print(")\n");
2424       
2425     }
2426   
2427   
2428   
2429   
2430   
2431   
2432   for(connectivity_int segmentAi=0;segmentAi<myPrimalMesh->segmentNumber;segmentAi++)
2433     {
2434       
2435       debug_print( "%-40s : %ld ( ","myPrimalMesh->segmentToFaceNeighbour", segmentAi);
2436       
2437       for(connectivity_int faceAi=0;faceAi<myPrimalMesh->segmentToFaceNeighbourNumber[segmentAi];faceAi++)
2438         {
2439           debug_print( "%ld " , myPrimalMesh->segmentToFaceNeighbour[segmentAi][faceAi] );
2440         }
2441       
2442       debug_print(")\n");
2443       
2444     }
2445   
2446 #endif
2447   
2448   free(copyFace);
2449   free(vertexList);
2450   free(vertexListInverted);
2451   
2452   
2453   endFunction(__FUNCTION__, t);
2454   
2461 void setDualFaces(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
2463   
2464   clock_t t=startFunction(__FUNCTION__);
2465   affiche("\n");
2466   
2467   connectivity_int selectedCell = 0;
2468   connectivity_int selectedVertex = 0;
2469   connectivity_int * selectedSegments = 0;
2470   connectivity_int selectedSegmentNumber = 0;
2471   connectivity_int * selectedFaces;
2472   connectivity_int   selectedFacesNumber = 0;
2473   
2474   myDualMesh->internalDualFaces  = (dataType *** ) malloc(sizeof (dataType**)*myDualMesh->internalDualFacesNumber);
2475   
2476   myDualMesh->segmentToInternalDualFace = (connectivity_int **) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
2477   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
2478     {
2479       myDualMesh->segmentToInternalDualFace[segmenti] = NULL;
2480     }  
2481   myDualMesh->segmentToInternalDualFaceNumber = (connectivity_int *) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
2482   memset(myDualMesh->segmentToInternalDualFaceNumber,0,myPrimalMesh->segmentNumber*sizeof (connectivity_int) );
2483   
2484   
2485   connectivity_int * output1;
2486   connectivity_int lenOutput1;
2487   
2488   
2489   connectivity_int * output2;
2490   connectivity_int lenOutput2;
2491   
2492   
2493   connectivity_int * output3;
2494   connectivity_int lenOutput3;
2495   
2496   output1  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
2497   output2  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
2498   output3  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
2499   
2500   
2501   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
2502     {
2503       
2504       selectedCell = celli;
2505       if(celli % 100000== 0 && celli != 0)
2506         release_print( "%-40s :  (%ld)\n","CELL", celli);
2507       
2508       unionTwoIntegerSet(&selectedFaces, &selectedFacesNumber,
2509                          (myPrimalMesh->cellToFacesOwner[celli]) , myPrimalMesh->cellToFacesOwnerNumber[celli] , 
2510                          (myPrimalMesh->cellToFacesNeighbour[celli]) , myPrimalMesh->cellToFacesNeighbourNumber[celli]);
2511       
2512       
2513 #ifdef DEBUG      
2514       debug_print( " ( ");
2515       
2516       for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
2517         {
2518           debug_print( "%ld " , selectedFaces[vertexi] );
2519         }
2520       debug_print(" )");
2521       debug_print(" \n");
2522 #endif  
2523       
2524       unionTwoIntegerSet(&selectedSegments, &selectedSegmentNumber,
2525                          (myPrimalMesh->cellToSegmentOwner[celli]) , myPrimalMesh->cellToSegmentOwnerNumber[celli] , 
2526                          (myPrimalMesh->cellToSegmentNeighbour[celli]) , myPrimalMesh->cellToSegmentNeighbourNumber[celli]);
2527       
2528       
2529       
2530       for(connectivity_int segmenti=0;segmenti<selectedSegmentNumber;segmenti++)
2531         {
2532           // select segments attached to one cell          
2533           
2534           
2535           //          for(connectivity_int facei=0;facei<selectedFacesNumber;facei++)
2536           //            {          
2537           // select faces attached to one segment
2538           // We have to intersect previous set and this one
2539           
2540           intersectTwoIntegerSet(&output1,&lenOutput1,
2541                                  selectedFaces,selectedFacesNumber,
2542                                  myPrimalMesh->segmentToFaceOwner[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceOwnerNumber[selectedSegments[segmenti]]
2543               );
2544           
2545           intersectTwoIntegerSet(&output2,&lenOutput2,
2546                                  selectedFaces,selectedFacesNumber,
2547                                  myPrimalMesh->segmentToFaceNeighbour[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceNeighbourNumber[selectedSegments[segmenti]]
2548               );
2549           
2550           unionTwoIntegerSet(&output3, &lenOutput3,
2551                              output1 , lenOutput1 , 
2552                              output2 , lenOutput2);
2553           
2554           
2555 #ifdef DEBUG
2556           
2557           debug_print( "Intersect %ld ( ",celli);
2558           
2559           for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
2560             {
2561               debug_print( "%ld " , selectedFaces[vertexi] );
2562             }
2563           debug_print(" ) ");
2564           debug_print( "( ");
2565           
2566           for(connectivity_int vertexi=0;vertexi<selectedSegmentNumber;vertexi++)
2567             {
2568               debug_print( "%ld " , selectedSegments[vertexi] );
2569             }
2570           debug_print(" ) ");
2571           debug_print( "( ");
2572           
2573           for(connectivity_int vertexi=0;vertexi<lenOutput1;vertexi++)
2574             {
2575               debug_print( "%ld " ,output1[vertexi] );
2576             }
2577           for(connectivity_int vertexi=0;vertexi<lenOutput2;vertexi++)
2578             {
2579               debug_print( "%ld " ,output2[vertexi] );
2580             }
2581           debug_print(" )");
2582           debug_print(" \n");
2583           
2584 #endif              
2585           
2586           if(lenOutput3==2) {
2587               //selectedSegments[segmenti]
2588               myDualMesh->internalDualFaces = (dataType ***) realloc(myDualMesh->internalDualFaces, (myDualMesh->internalDualFacesNumber+1) * sizeof (dataType**));
2589               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber] = (dataType **) malloc( HEXAHEDRON_DUAL_FACE_POINT * sizeof (dataType*));   
2590               //              memset(myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber], 0, HEXAHEDRON_DUAL_FACE_POINT*sizeof(connectivity_int));
2591               
2592               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT1] = (myPrimalMesh->faceCentres[output3[0]]); 
2593               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT2] = (myPrimalMesh->volumeCentroid[celli]); 
2594               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT3] = (myPrimalMesh->faceCentres[output3[1]]); 
2595               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT4] = (myPrimalMesh->segmentCentres[selectedSegments[segmenti]]); 
2596               
2597               
2598               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]]= (connectivity_int*) realloc(myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]] ,sizeof (connectivity_int)*(myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]+1));
2599               
2600               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]][ myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]] ] = myDualMesh->internalDualFacesNumber;
2601               myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]++;
2602               
2603               myDualMesh->internalDualFacesNumber++;              
2604             }
2605           //              else 
2606           //                {
2607           //                  affiche("Error on setDualFaces definition\n");
2608           //                  exit(1);
2609           //                }
2610         }
2611       //        }
2612     }
2613   
2614 #ifdef DEBUG
2615   
2616   
2617   
2618   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
2619     {
2620       
2621       debug_print( "%-40s : %ld","myDualMesh->internalDualFaces", faceAi);
2622       
2623       for(connectivity_int pointi=0;pointi<HEXAHEDRON_DUAL_FACE_POINT;pointi++)
2624         {
2625           debug_print( " ( ");
2626           
2627           for(connectivity_int vertexi=0;vertexi<DIM3D;vertexi++)
2628             {
2629               debug_print( "%lf " , myDualMesh->internalDualFaces[faceAi][pointi][vertexi] );
2630             }
2631           debug_print(" )");
2632         }
2633       debug_print(" \n");
2634       
2635     }
2636   
2637   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
2638     {
2639       
2640       debug_print( "%-40s : %ld","myDualMesh->segmentToInternalDualFace", segmenti);
2641       
2642       debug_print( " ( ");
2643       
2644       
2645       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
2646         {
2647           debug_print( "%ld " , myDualMesh->segmentToInternalDualFace[segmenti][facei] );
2648         }
2649       
2650       debug_print(" )");
2651       debug_print(" \n");
2652       
2653     }
2654   
2655   
2656   
2657   
2658   
2659   
2660 #endif
2661   
2662   free(output1);
2663   free(output2);
2664   free(output3);
2665   
2666   
2667   endFunction(__FUNCTION__, t);
2668   
2673 void setDualFacesCentres(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
2675   
2676   clock_t t=startFunction(__FUNCTION__);
2677   affiche("\n");
2678   
2679   connectivity_int actualPoint;
2680   connectivity_int nextPoint;
2681   
2682   dataType *fC;
2684   dataType  * result1 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
2685   dataType  * result2 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
2686   dataType    result3 = 0.0;
2687   dataType    result4 = 0.0;
2688   connectivity_int  face  =0.0;  
2689   connectivity_int * segmentVector=NULL;
2690   
2691   fC =(dataType *) calloc(DIM3D, sizeof(dataType));
2692   
2693   
2694   myDualMesh->internalDualFaceCentres  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
2695   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
2696     {
2697       myDualMesh->internalDualFaceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
2698       memset(myDualMesh->internalDualFaceCentres[i], 0.0, DIM3D*sizeof(dataType));
2699     }
2700   
2701   
2702   
2703   myDualMesh->internalDualFaceArea  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
2704   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
2705     {
2706       myDualMesh->internalDualFaceArea[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
2707       memset(myDualMesh->internalDualFaceArea[i], 0.0, DIM3D*sizeof(dataType));
2708     }
2709   
2710   
2711   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
2712     {
2713       if(faceAi % 1000000== 0 && faceAi != 0)
2714         release_print( "%-40s :  (%ld)\n","FACE", faceAi);
2715       
2716       
2717       zeroVector(&fC);
2718       
2719       for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
2720         {
2721           actualPoint = vertexAi;
2722           
2723           sumTwoVectors(&fC,fC,myDualMesh->internalDualFaces[faceAi][actualPoint]);
2724           
2725         }
2726       
2727       
2728       scalarDivVector(&(myDualMesh->internalDualFaceCentres[faceAi]), QUAD,fC);
2729       
2730 #ifdef DEBUG      
2731       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres Esti", faceAi, QUAD);
2732       
2733       for(connectivity_int ii=0;ii<DIM3D;ii++)
2734         {
2735           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
2736         }
2737       debug_print(")\n");
2738 #endif
2739       if(QUAD==3) {
2740           
2741           // Triangle
2742           affiche("NON IMPLEMENTED ...");
2743           exit(1);
2744           
2745         }
2746       else
2747         { //https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
2748           dataType *sumN;
2749           dataType sumA       = ZEROSCALAR;
2750           dataType *sumAc;
2751           
2752           dataType *c;
2753           dataType a          = ZEROSCALAR;
2754           dataType * n       ;
2755           
2756           dataType *result1;
2757           dataType *result2;
2758           
2759           
2760           result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
2761           sumN =(dataType *) calloc(DIM3D, sizeof(dataType));
2762           sumAc =(dataType *) calloc(DIM3D, sizeof(dataType));
2763           c =(dataType *) calloc(DIM3D, sizeof(dataType));
2764           n =(dataType *) calloc(DIM3D, sizeof(dataType));
2765           result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
2766           
2767           
2768           
2769           for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++) 
2770             {
2771               
2772               actualPoint = vertexAi;
2773               nextPoint = (vertexAi + 1) % QUAD;
2774               
2775               sumThreeVectors(&c, (myDualMesh->internalDualFaces[faceAi][actualPoint]), (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaceCentres[faceAi]) );
2776               
2777               subTwoVectors(&result1, (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
2778               subTwoVectors(&result2, (myDualMesh->internalDualFaceCentres[faceAi]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
2779               
2780               crossProduct(&n, result1, result2);
2781               mag(&(a), n);
2782               
2783               sumTwoVectors(&sumN, sumN, n);
2784               
2785               sumA+=a;
2786               scalarDotVector(&result1, a, c);
2787               sumTwoVectors(&sumAc, sumAc, result1);
2788             }
2789           
2790           if (sumA < SMALL)
2791             {
2792               zeroVector( &(myDualMesh->internalDualFaceArea[faceAi]) );
2793             }
2794           else
2795             {
2796               scalarDotVector(&(myDualMesh->internalDualFaceCentres[faceAi]), 1.0/(3.0*sumA), sumAc); //correct faceCentres after estimating them
2797               scalarDotVector(&(myDualMesh->internalDualFaceArea[faceAi]), 0.5, sumN);
2798             }
2799           
2800           
2801           free(sumN);
2802           free(sumAc);
2803           free(c);
2804           free(n)       ;
2805           free(result1);
2806           free(result2);
2807           
2808         }
2809       
2810     }
2811   
2812   // correct direction of seg and area vector of dual mesh
2813   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
2814     {    
2815       segmentVector = myPrimalMesh->segments[segmenti];
2816       
2818       subTwoVectors(&result1, myPrimalMesh->vertex[segmentVector[PT2]], myPrimalMesh->vertex[segmentVector[PT1]]);      
2819       
2820       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
2821         {
2822           face  = myDualMesh->segmentToInternalDualFace[segmenti][facei];
2823           
2824           dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
2825           if(result4<0.0)
2826             {
2827               for(connectivity_int pointi=0;pointi<DIM3D;pointi++)
2828                 {
2829                   myDualMesh->internalDualFaceArea[face][pointi] *= -1.0;
2830                 }              
2831             }
2832         }
2833       
2834       
2835 #ifdef DEBUG      
2836       
2837       dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
2839       debug_print( "%-40s : %ld ( ","VECTOR", segmenti);
2840       
2841       for(connectivity_int ii=0;ii<DIM3D;ii++)
2842         {
2843           debug_print( "%lf " , result1[ii]);
2844         }
2845       
2846       debug_print(") ");     
2847       
2848       debug_print( "%lf ", result4);
2849       
2850       debug_print("\n");      
2851 #endif      
2852     }
2853   
2854   
2855   
2856   
2857   
2858   
2859   
2860 #ifdef DEBUG     
2861   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
2862     {          
2863       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres", faceAi, QUAD);
2864       
2865       for(connectivity_int ii=0;ii<DIM3D;ii++)
2866         {
2867           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
2868         }
2869       
2870       debug_print(")\n");
2871       
2872       debug_print( "%-40s : %ld ( ","myDualMesh->internalDualFaceArea", faceAi);
2873       
2874       for(connectivity_int ii=0;ii<DIM3D;ii++)
2875         {
2876           debug_print( "%lf " , myDualMesh->internalDualFaceArea[faceAi][ii]);
2877         }
2878       
2879       debug_print(")\n");
2880     }     
2881 #endif          
2882   
2883  
2884   
2885   free(result1);
2886   free(result2);
2887   free(fC);    
2888   endFunction(__FUNCTION__, t);
2889   
2894 void setHexahedreSortFaceToVertex(struct primalMesh * myPrimalMesh)
2896   //////////////////// TRIE LES NOEUDS
2897   
2898   clock_t t=startFunction(__FUNCTION__);
2899   
2900   //    connectivity_int actualPoint;
2901   //    connectivity_int nextPoint;
2902   
2903   // on trie les noeund dans les faces
2904   for(connectivity_int k=0;k<myPrimalMesh->faceNumber;k++)
2905     {
2906       connectivity_int              position=0;
2907       
2908       
2909       for(connectivity_int ii=1;ii<myPrimalMesh->faceToVertexNumber[k];ii++)
2910         {
2911           if(myPrimalMesh->faceToVertex[k][position]>myPrimalMesh->faceToVertex[k][ii])
2912             {
2913               position=ii;
2914             }
2915         }
2916       
2917       for(connectivity_int ii=0;ii<position;ii++)
2918         rotate(myPrimalMesh->faceToVertex[k], myPrimalMesh->faceToVertexNumber[k]);
2919       
2920       
2921     }
2922   
2923   
2924   
2925   
2926   //        for(conctivity_int k=0;k<myPrimalMesh->faceNumber;k++)
2927   //        {
2928   //            conctivity_int it_surface=0;
2929   //            int             bool=0;
2930   //            do
2931   //            {
2932   //                for(conctivity_int ii=0;ii<myPrimalMesh->faceToVertexNumber[k];ii++)
2933   //                {
2934   //                    //     ((ii + 1) % myPrimalMesh->faceToVertexNumber[k]) ;
2935   
2936   //                    if(myPrimalMesh->faceToVertex[k][ii] == myPrimalMesh->faceToVertex[it_surface][0])
2937   //                    {
2938   //                        bool = 1;
2939   //                        debug_print("TROUVE : (%ld) myPrimalMesh->faceToVertex[%ld][%ld] = %ld\n", it_surface,k,ii,myPrimalMesh->faceToVertex[k][ii])
2940   //                    }
2941   //                }
2942   //                it_surface++;
2943   //            }
2944   //            while( (it_surface<myPrimalMesh->faceNumber) && (bool == 0) );
2945   //        }
2946   endFunction(__FUNCTION__, t);
2947   
2951 void setHexahedreAllocVertexToSegment(struct primalMesh * myPrimalMesh)
2953   
2954   
2955   clock_t t=startFunction(__FUNCTION__);
2956   
2957   
2958   
2959   myPrimalMesh->vertexToSegments = (connectivity_int **) malloc(myPrimalMesh->vertexNumber * sizeof(connectivity_int*));
2960   for(connectivity_int i=0; i<myPrimalMesh->vertexNumber; i++)
2961     {
2962       myPrimalMesh->vertexToSegments[i] = (connectivity_int *) malloc( myPrimalMesh->vertexToSegmentNumber[i] * sizeof(connectivity_int));
2963       memset(myPrimalMesh->vertexToSegments[i], 0.0, QUAD*sizeof(connectivity_int));
2964     }
2965   
2966   
2967 #ifdef DEBUG  
2968   for(connectivity_int i=0;i<myPrimalMesh->vertexNumber;i++)
2969     {
2970