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

eac1835b2e8b44d45ca63d2860ff98d7942e4244
1 #include "mmd.h"
2 #include "dual.h"
3 #include "tensor.h"
4 #include "listop.h"
9 void setDualFaces(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
10 {
11   
12   clock_t t=startFunction(__FUNCTION__);
13   affiche("\n");
14   
15   connectivity_int selectedCell = 0;
16   connectivity_int selectedVertex = 0;
17   connectivity_int * selectedSegments = NULL;
18   connectivity_int selectedSegmentNumber = 0;
19   connectivity_int * selectedFaces = NULL;
20   connectivity_int   selectedFacesNumber = 0;
21   
22   myDualMesh->internalDualFaces  = (dataType *** ) malloc(sizeof (dataType**)*myDualMesh->internalDualFacesNumber);
23   
24   myDualMesh->segmentToInternalDualFace = (connectivity_int **) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
25   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
26     {
27       myDualMesh->segmentToInternalDualFace[segmenti] = NULL;
28     }  
29   myDualMesh->segmentToInternalDualFaceNumber = (connectivity_int *) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
30   memset(myDualMesh->segmentToInternalDualFaceNumber,0,myPrimalMesh->segmentNumber*sizeof (connectivity_int) );
31   
32   
33   connectivity_int * output1;
34   connectivity_int lenOutput1;
35   
36   
37   connectivity_int * output2;
38   connectivity_int lenOutput2;
39   
40   
41   connectivity_int * output3;
42   connectivity_int lenOutput3;
43   
44   output1  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
45   output2  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
46   output3  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
47   selectedFaces  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
48   selectedSegments = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
49   
50   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
51     {
52       
53       selectedCell = celli;
54       if(celli % 100000== 0 && celli != 0)
55         release_print( "%-40s :  (%ld)\n","CELL", celli);
56       
57       unionTwoIntegerSet(&selectedFaces, &selectedFacesNumber,
58                          (myPrimalMesh->cellToFacesOwner[celli]) , myPrimalMesh->cellToFacesOwnerNumber[celli] , 
59                          (myPrimalMesh->cellToFacesNeighbour[celli]) , myPrimalMesh->cellToFacesNeighbourNumber[celli]);
60       
61       
62 #ifdef DEBUG      
63       debug_print( " ( ");
64       
65       for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
66         {
67           debug_print( "%ld " , selectedFaces[vertexi] );
68         }
69       debug_print(" )");
70       debug_print(" \n");
71 #endif  
72       
73       unionTwoIntegerSet(&selectedSegments, &selectedSegmentNumber,
74                          (myPrimalMesh->cellToSegmentOwner[celli]) , myPrimalMesh->cellToSegmentOwnerNumber[celli] , 
75                          (myPrimalMesh->cellToSegmentNeighbour[celli]) , myPrimalMesh->cellToSegmentNeighbourNumber[celli]);
76       
77       
78       
79       for(connectivity_int segmenti=0;segmenti<selectedSegmentNumber;segmenti++)
80         {
81           // select segments attached to one cell          
82           
83           
84           //          for(connectivity_int facei=0;facei<selectedFacesNumber;facei++)
85           //            {          
86           // select faces attached to one segment
87           // We have to intersect previous set and this one
88           
89           intersectTwoIntegerSet(&output1,&lenOutput1,
90                                  selectedFaces,selectedFacesNumber,
91                                  myPrimalMesh->segmentToFaceOwner[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceOwnerNumber[selectedSegments[segmenti]]
92               );
93           
94           intersectTwoIntegerSet(&output2,&lenOutput2,
95                                  selectedFaces,selectedFacesNumber,
96                                  myPrimalMesh->segmentToFaceNeighbour[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceNeighbourNumber[selectedSegments[segmenti]]
97               );
98           
99           unionTwoIntegerSet(&output3, &lenOutput3,
100                              output1 , lenOutput1 , 
101                              output2 , lenOutput2);
102           
103           
104 #ifdef DEBUG
105           
106           debug_print( "Intersect %ld ( ",celli);
107           
108           for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
109             {
110               debug_print( "%ld " , selectedFaces[vertexi] );
111             }
112           debug_print(" ) ");
113           debug_print( "( ");
114           
115           for(connectivity_int vertexi=0;vertexi<selectedSegmentNumber;vertexi++)
116             {
117               debug_print( "%ld " , selectedSegments[vertexi] );
118             }
119           debug_print(" ) ");
120           debug_print( "( ");
121           
122           for(connectivity_int vertexi=0;vertexi<lenOutput1;vertexi++)
123             {
124               debug_print( "%ld " ,output1[vertexi] );
125             }
126           for(connectivity_int vertexi=0;vertexi<lenOutput2;vertexi++)
127             {
128               debug_print( "%ld " ,output2[vertexi] );
129             }
130           debug_print(" )");
131           debug_print(" \n");
132           
133 #endif              
134           
135           if(lenOutput3==2) {
136               //selectedSegments[segmenti]
137               myDualMesh->internalDualFaces = (dataType ***) realloc(myDualMesh->internalDualFaces, (myDualMesh->internalDualFacesNumber+1) * sizeof (dataType**));
138               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber] = (dataType **) malloc( HEXAHEDRON_DUAL_FACE_POINT * sizeof (dataType*));   
139               //              memset(myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber], 0, HEXAHEDRON_DUAL_FACE_POINT*sizeof(connectivity_int));
140               
141               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT1] = (myPrimalMesh->faceCentres[output3[0]]); 
142               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT2] = (myPrimalMesh->volumeCentroid[celli]); 
143               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT3] = (myPrimalMesh->faceCentres[output3[1]]); 
144               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT4] = (myPrimalMesh->segmentCentres[selectedSegments[segmenti]]); 
145               
146               
147               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]]= (connectivity_int*) realloc(myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]] ,sizeof (connectivity_int)*(myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]+1));
148               
149               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]][ myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]] ] = myDualMesh->internalDualFacesNumber;
150               myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]++;
151               
152               myDualMesh->internalDualFacesNumber++;              
153             }
154           //              else 
155           //                {
156           //                  affiche("Error on setDualFaces definition\n");
157           //                  exit(1);
158           //                }
159         }
160       //        }
161     }
162   
163 #ifdef DEBUG
164   
165   
166   
167   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
168     {
169       
170       debug_print( "%-40s : %ld","myDualMesh->internalDualFaces", faceAi);
171       
172       for(connectivity_int pointi=0;pointi<HEXAHEDRON_DUAL_FACE_POINT;pointi++)
173         {
174           debug_print( " ( ");
175           
176           for(connectivity_int vertexi=0;vertexi<DIM3D;vertexi++)
177             {
178               debug_print( "%lf " , myDualMesh->internalDualFaces[faceAi][pointi][vertexi] );
179             }
180           debug_print(" )");
181         }
182       debug_print(" \n");
183       
184     }
185   
186   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
187     {
188       
189       debug_print( "%-40s : %ld","myDualMesh->segmentToInternalDualFace", segmenti);
190       
191       debug_print( " ( ");
192       
193       
194       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
195         {
196           debug_print( "%ld " , myDualMesh->segmentToInternalDualFace[segmenti][facei] );
197         }
198       
199       debug_print(" )");
200       debug_print(" \n");
201       
202     }
203   
204   
205   
206   
207   
208   
209 #endif
210   
211   free(selectedFaces);
212   free(selectedSegments);
213   free(output1);
214   free(output2);
215   free(output3);
216   
217   
218   endFunction(__FUNCTION__, t);
219   
224 void setDualFacesCentres(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
226   
227   clock_t t=startFunction(__FUNCTION__);
228   affiche("\n");
229   
230   connectivity_int actualPoint;
231   connectivity_int nextPoint;
232   
233   dataType *fC;
235   dataType  * result1 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
236   dataType  * result2 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
237   dataType    result3 = 0.0;
238   dataType    result4 = 0.0;
239   connectivity_int  face  =0.0;  
240   connectivity_int * segmentVector=NULL;
241   
242   fC =(dataType *) calloc(DIM3D, sizeof(dataType));
243   
244   
245   myDualMesh->internalDualFaceCentres  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
246   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
247     {
248       myDualMesh->internalDualFaceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
249       memset(myDualMesh->internalDualFaceCentres[i], 0.0, DIM3D*sizeof(dataType));
250     }
251   
252   
253   
254   myDualMesh->internalDualFaceArea  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
255   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
256     {
257       myDualMesh->internalDualFaceArea[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
258       memset(myDualMesh->internalDualFaceArea[i], 0.0, DIM3D*sizeof(dataType));
259     }
260   
261   
262   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
263     {
264       if(faceAi % 1000000== 0 && faceAi != 0)
265         release_print( "%-40s :  (%ld)\n","FACE", faceAi);
266       
267       
268       zeroVector(&fC);
269       
270       for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
271         {
272           actualPoint = vertexAi;
273           
274           sumTwoVectors(&fC,fC,myDualMesh->internalDualFaces[faceAi][actualPoint]);
275           
276         }
277       
278       
279       scalarDivVector(&(myDualMesh->internalDualFaceCentres[faceAi]), QUAD,fC);
280       
281 #ifdef DEBUG      
282       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres Esti", faceAi, QUAD);
283       
284       for(connectivity_int ii=0;ii<DIM3D;ii++)
285         {
286           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
287         }
288       debug_print(")\n");
289 #endif
290       if(QUAD==3) {
291           
292           // Triangle
293           affiche("NON IMPLEMENTED ...");
294           exit(1);
295           
296         }
297       else
298         { //https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
299           dataType *sumN;
300           dataType sumA       = ZEROSCALAR;
301           dataType *sumAc;
302           
303           dataType *c;
304           dataType a          = ZEROSCALAR;
305           dataType * n       ;
306           
307           dataType *result1;
308           dataType *result2;
309           
310           
311           result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
312           sumN =(dataType *) calloc(DIM3D, sizeof(dataType));
313           sumAc =(dataType *) calloc(DIM3D, sizeof(dataType));
314           c =(dataType *) calloc(DIM3D, sizeof(dataType));
315           n =(dataType *) calloc(DIM3D, sizeof(dataType));
316           result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
317           
318           
319           
320           for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++) 
321             {
322               
323               actualPoint = vertexAi;
324               nextPoint = (vertexAi + 1) % QUAD;
325               
326               sumThreeVectors(&c, (myDualMesh->internalDualFaces[faceAi][actualPoint]), (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaceCentres[faceAi]) );
327               
328               subTwoVectors(&result1, (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
329               subTwoVectors(&result2, (myDualMesh->internalDualFaceCentres[faceAi]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
330               
331               crossProduct(&n, result1, result2);
332               mag(&(a), n);
333               
334               sumTwoVectors(&sumN, sumN, n);
335               
336               sumA+=a;
337               scalarDotVector(&result1, a, c);
338               sumTwoVectors(&sumAc, sumAc, result1);
339             }
340           
341           if (sumA < SMALL)
342             {
343               zeroVector( &(myDualMesh->internalDualFaceArea[faceAi]) );
344             }
345           else
346             {
347               scalarDotVector(&(myDualMesh->internalDualFaceCentres[faceAi]), 1.0/(3.0*sumA), sumAc); //correct faceCentres after estimating them
348               scalarDotVector(&(myDualMesh->internalDualFaceArea[faceAi]), 0.5, sumN);
349             }
350           
351           
352           free(sumN);
353           free(sumAc);
354           free(c);
355           free(n)       ;
356           free(result1);
357           free(result2);
358           
359         }
360       
361     }
362   
363   // correct direction of seg and area vector of dual mesh
364   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
365     {    
366       segmentVector = myPrimalMesh->segments[segmenti];
367       
369       subTwoVectors(&result1, myPrimalMesh->vertex[segmentVector[PT2]], myPrimalMesh->vertex[segmentVector[PT1]]);      
370       
371       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
372         {
373           face  = myDualMesh->segmentToInternalDualFace[segmenti][facei];
374           
375           dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
376           if(result4<0.0)
377             {
378               for(connectivity_int pointi=0;pointi<DIM3D;pointi++)
379                 {
380                   myDualMesh->internalDualFaceArea[face][pointi] *= -1.0;
381                 }              
382             }
383         }
384       
385       
386 #ifdef DEBUG      
387       
388       dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
390       debug_print( "%-40s : %ld ( ","VECTOR", segmenti);
391       
392       for(connectivity_int ii=0;ii<DIM3D;ii++)
393         {
394           debug_print( "%lf " , result1[ii]);
395         }
396       
397       debug_print(") ");     
398       
399       debug_print( "%lf ", result4);
400       
401       debug_print("\n");      
402 #endif      
403     }
404   
405   
406   
407   
408   
409   
410   
411 #ifdef DEBUG     
412   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
413     {          
414       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres", faceAi, QUAD);
415       
416       for(connectivity_int ii=0;ii<DIM3D;ii++)
417         {
418           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
419         }
420       
421       debug_print(")\n");
422       
423       debug_print( "%-40s : %ld ( ","myDualMesh->internalDualFaceArea", faceAi);
424       
425       for(connectivity_int ii=0;ii<DIM3D;ii++)
426         {
427           debug_print( "%lf " , myDualMesh->internalDualFaceArea[faceAi][ii]);
428         }
429       
430       debug_print(")\n");
431     }     
432 #endif          
433   
434  
435   
436   free(result1);
437   free(result2);
438   free(fC);    
439   endFunction(__FUNCTION__, t);
440