The Higher Education and Research forge

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

SCM Repository

1 /*---------------------------------------------------------------------------*\
2     \o/\o/\o/
3     MMD
4     Version : 0.0.3
5     Web : https://github.com/alainbastide/MMD
6 -------------------------------------------------------------------------------
7 License
9     MMD is free software: you can redistribute it and/or modify it
10     under the terms of the GNU General Public License as published by the
11     Free Software Foundation, either version 3 of the License, or (at your
12     option) any later version.
14     MMD is distributed in the hope that it will be useful, but
15     WITHOUT ANY WARRANTY; without even the implied warranty of
16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17     General Public License for more details.
19     You should have received a copy of the GNU General Public License
20     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
22 Application
23     explain
25 Description
26     explain
28 Author
29     Alain Bastide, Université de La Réunion, FRANCE.  All rights reserved
31 \*---------------------------------------------------------------------------*/
36 #include "mmd.h"
37 #include "dual.h"
38 #include "tensor.h"
39 #include "listop.h"
44 void setDualFaces(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
45 {
46   
47   clock_t t=startFunction(__FUNCTION__);
48   affiche("\n");
49   
50   connectivity_int selectedCell = 0;
51   connectivity_int selectedVertex = 0;
52   connectivity_int * selectedSegments = NULL;
53   connectivity_int selectedSegmentNumber = 0;
54   connectivity_int * selectedFaces = NULL;
55   connectivity_int   selectedFacesNumber = 0;
56   
57   myDualMesh->internalDualFaces  = (dataType *** ) malloc(sizeof (dataType**)*myDualMesh->internalDualFacesNumber);
58   
59   myDualMesh->segmentToInternalDualFace = (connectivity_int **) malloc(sizeof (connectivity_int*)*myPrimalMesh->segmentNumber);
60   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
61     {
62       myDualMesh->segmentToInternalDualFace[segmenti] = NULL;
63     }  
64   myDualMesh->segmentToInternalDualFaceNumber = (connectivity_int *) malloc(sizeof (connectivity_int)*myPrimalMesh->segmentNumber);
65   memset(myDualMesh->segmentToInternalDualFaceNumber,0,myPrimalMesh->segmentNumber*sizeof (connectivity_int) );
66   
67   
68   connectivity_int * output1;
69   connectivity_int lenOutput1;
70   
71   
72   connectivity_int * output2;
73   connectivity_int lenOutput2;
74   
75   
76   connectivity_int * output3;
77   connectivity_int lenOutput3;
78   
79   output1  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
80   output2  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
81   output3  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
82   selectedFaces  = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
83   selectedSegments = (connectivity_int * ) malloc(sizeof (connectivity_int)*MAXSIZECONNECTIVITYVECTOR);
84   
85   for(connectivity_int celli=0;celli<myPrimalMesh->cellNumber;celli++)
86     {
87       
88       selectedCell = celli;
89       if(celli % 100000== 0 && celli != 0)
90         release_print( "%-40s :  (%ld)\n","CELL", celli);
91       
92       unionTwoIntegerSet(&selectedFaces, &selectedFacesNumber,
93                          (myPrimalMesh->cellToFacesOwner[celli]) , myPrimalMesh->cellToFacesOwnerNumber[celli] , 
94                          (myPrimalMesh->cellToFacesNeighbour[celli]) , myPrimalMesh->cellToFacesNeighbourNumber[celli]);
95       
96       
97 #ifdef DEBUG      
98       debug_print( " ( ");
99       
100       for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
101         {
102           debug_print( "%ld " , selectedFaces[vertexi] );
103         }
104       debug_print(" )");
105       debug_print(" \n");
106 #endif  
107       
108       unionTwoIntegerSet(&selectedSegments, &selectedSegmentNumber,
109                          (myPrimalMesh->cellToSegmentOwner[celli]) , myPrimalMesh->cellToSegmentOwnerNumber[celli] , 
110                          (myPrimalMesh->cellToSegmentNeighbour[celli]) , myPrimalMesh->cellToSegmentNeighbourNumber[celli]);
111       
112       
113       
114       for(connectivity_int segmenti=0;segmenti<selectedSegmentNumber;segmenti++)
115         {
116           // select segments attached to one cell          
117           
118           
119           //          for(connectivity_int facei=0;facei<selectedFacesNumber;facei++)
120           //            {          
121           // select faces attached to one segment
122           // We have to intersect previous set and this one
123           
124           intersectTwoIntegerSet(&output1,&lenOutput1,
125                                  selectedFaces,selectedFacesNumber,
126                                  myPrimalMesh->segmentToFaceOwner[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceOwnerNumber[selectedSegments[segmenti]]
127               );
128           
129           intersectTwoIntegerSet(&output2,&lenOutput2,
130                                  selectedFaces,selectedFacesNumber,
131                                  myPrimalMesh->segmentToFaceNeighbour[selectedSegments[segmenti]],myPrimalMesh->segmentToFaceNeighbourNumber[selectedSegments[segmenti]]
132               );
133           
134           unionTwoIntegerSet(&output3, &lenOutput3,
135                              output1 , lenOutput1 , 
136                              output2 , lenOutput2);
137           
138           
139 #ifdef DEBUG
140           
141           debug_print( "Intersect %ld ( ",celli);
142           
143           for(connectivity_int vertexi=0;vertexi<selectedFacesNumber;vertexi++)
144             {
145               debug_print( "%ld " , selectedFaces[vertexi] );
146             }
147           debug_print(" ) ");
148           debug_print( "( ");
149           
150           for(connectivity_int vertexi=0;vertexi<selectedSegmentNumber;vertexi++)
151             {
152               debug_print( "%ld " , selectedSegments[vertexi] );
153             }
154           debug_print(" ) ");
155           debug_print( "( ");
156           
157           for(connectivity_int vertexi=0;vertexi<lenOutput1;vertexi++)
158             {
159               debug_print( "%ld " ,output1[vertexi] );
160             }
161           for(connectivity_int vertexi=0;vertexi<lenOutput2;vertexi++)
162             {
163               debug_print( "%ld " ,output2[vertexi] );
164             }
165           debug_print(" )");
166           debug_print(" \n");
167           
168 #endif              
169           
170           if(lenOutput3==2) {
171               //selectedSegments[segmenti]
172               myDualMesh->internalDualFaces = (dataType ***) realloc(myDualMesh->internalDualFaces, (myDualMesh->internalDualFacesNumber+1) * sizeof (dataType**));
173               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber] = (dataType **) malloc( HEXAHEDRON_DUAL_FACE_POINT * sizeof (dataType*));   
174               //              memset(myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber], 0, HEXAHEDRON_DUAL_FACE_POINT*sizeof(connectivity_int));
175               
176               // HEXAHEDRE
177               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT1] = (myPrimalMesh->faceCentres[output3[0]]); 
178               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT2] = (myPrimalMesh->volumeCentroid[celli]); 
179               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT3] = (myPrimalMesh->faceCentres[output3[1]]); 
180               myDualMesh->internalDualFaces[myDualMesh->internalDualFacesNumber][PT4] = (myPrimalMesh->segmentCentres[selectedSegments[segmenti]]); 
181               
182               
183               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]]= (connectivity_int*) realloc(myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]] ,sizeof (connectivity_int)*(myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]+1));
184               
185               myDualMesh->segmentToInternalDualFace[selectedSegments[segmenti]][ myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]] ] = myDualMesh->internalDualFacesNumber;
186               myDualMesh->segmentToInternalDualFaceNumber[selectedSegments[segmenti]]++;
187               
188               myDualMesh->internalDualFacesNumber++;              
189             }
190           //              else 
191           //                {
192           //                  affiche("Error on setDualFaces definition\n");
193           //                  exit(1);
194           //                }
195         }
196       //        }
197     }
198   
199 #ifdef DEBUG
200   
201   
202   
203   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
204     {
205       
206       debug_print( "%-40s : %ld","myDualMesh->internalDualFaces", faceAi);
207       
208       for(connectivity_int pointi=0;pointi<HEXAHEDRON_DUAL_FACE_POINT;pointi++)
209         {
210           debug_print( " ( ");
211           
212           for(connectivity_int vertexi=0;vertexi<DIM3D;vertexi++)
213             {
214               debug_print( "%lf " , myDualMesh->internalDualFaces[faceAi][pointi][vertexi] );
215             }
216           debug_print(" )");
217         }
218       debug_print(" \n");
219       
220     }
221   
222   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
223     {
224       
225       debug_print( "%-40s : %ld","myDualMesh->segmentToInternalDualFace", segmenti);
226       
227       debug_print( " ( ");
228       
229       
230       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
231         {
232           debug_print( "%ld " , myDualMesh->segmentToInternalDualFace[segmenti][facei] );
233         }
234       
235       debug_print(" )");
236       debug_print(" \n");
237       
238     }
239   
240   
241   
242   
243   
244   
245 #endif
246   
247   free(selectedFaces);
248   free(selectedSegments);
249   free(output1);
250   free(output2);
251   free(output3);
252   
253   
254   endFunction(__FUNCTION__, t);
255   
260 void setDualFacesCentres(struct primalMesh * myPrimalMesh, struct dualMesh * myDualMesh)
262   
263   clock_t t=startFunction(__FUNCTION__);
264   affiche("\n");
265   
266   connectivity_int actualPoint;
267   connectivity_int nextPoint;
268   
269   dataType *fC;
271   dataType  * result1 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
272   dataType  * result2 = (dataType * ) malloc(sizeof (dataType)*DIM3D);
273   dataType    result3 = 0.0;
274   dataType    result4 = 0.0;
275   connectivity_int  face  =0.0;  
276   connectivity_int * segmentVector=NULL;
277   
278   fC =(dataType *) calloc(DIM3D, sizeof(dataType));
279   
280   
281   myDualMesh->internalDualFaceCentres  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
282   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
283     {
284       myDualMesh->internalDualFaceCentres[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
285       memset(myDualMesh->internalDualFaceCentres[i], 0.0, DIM3D*sizeof(dataType));
286     }
287   
288   
289   
290   myDualMesh->internalDualFaceArea  = (dataType ** ) malloc(sizeof (dataType*)*myDualMesh->internalDualFacesNumber);
291   for(connectivity_int i=0; i<myDualMesh->internalDualFacesNumber; i++)
292     {
293       myDualMesh->internalDualFaceArea[i] = (dataType *) malloc(DIM3D * sizeof(dataType));
294       memset(myDualMesh->internalDualFaceArea[i], 0.0, DIM3D*sizeof(dataType));
295     }
296   
297   
298   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
299     {
300       if(faceAi % 1000000== 0 && faceAi != 0)
301         release_print( "%-40s :  (%ld)\n","FACE", faceAi);
302       
303       
304       zeroVector(&fC);
305       
306       for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++)
307         {
308           actualPoint = vertexAi;
309           
310           sumTwoVectors(&fC,fC,myDualMesh->internalDualFaces[faceAi][actualPoint]);
311           
312         }
313       
314       
315       scalarDivVector(&(myDualMesh->internalDualFaceCentres[faceAi]), QUAD,fC);
316       
317 #ifdef DEBUG      
318       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres Esti", faceAi, QUAD);
319       
320       for(connectivity_int ii=0;ii<DIM3D;ii++)
321         {
322           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
323         }
324       debug_print(")\n");
325 #endif
326       if(QUAD==3) {
327           
328           // Triangle
329           affiche("NON IMPLEMENTED ...");
330           exit(1);
331           
332         }
333       else
334         { //https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
335           dataType *sumN;
336           dataType sumA       = ZEROSCALAR;
337           dataType *sumAc;
338           
339           dataType *c;
340           dataType a          = ZEROSCALAR;
341           dataType * n       ;
342           
343           dataType *result1;
344           dataType *result2;
345           
346           
347           result1 =(dataType *) calloc(DIM3D, sizeof(dataType));
348           sumN =(dataType *) calloc(DIM3D, sizeof(dataType));
349           sumAc =(dataType *) calloc(DIM3D, sizeof(dataType));
350           c =(dataType *) calloc(DIM3D, sizeof(dataType));
351           n =(dataType *) calloc(DIM3D, sizeof(dataType));
352           result2 =(dataType *) calloc(DIM3D, sizeof(dataType));
353           
354           
355           
356           for(connectivity_int vertexAi=0;vertexAi<QUAD;vertexAi++) 
357             {
358               
359               actualPoint = vertexAi;
360               nextPoint = (vertexAi + 1) % QUAD;
361               
362               sumThreeVectors(&c, (myDualMesh->internalDualFaces[faceAi][actualPoint]), (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaceCentres[faceAi]) );
363               
364               subTwoVectors(&result1, (myDualMesh->internalDualFaces[faceAi][nextPoint]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
365               subTwoVectors(&result2, (myDualMesh->internalDualFaceCentres[faceAi]), (myDualMesh->internalDualFaces[faceAi][actualPoint]));
366               
367               crossProduct(&n, result1, result2);
368               mag(&(a), n);
369               
370               sumTwoVectors(&sumN, sumN, n);
371               
372               sumA+=a;
373               scalarDotVector(&result1, a, c);
374               sumTwoVectors(&sumAc, sumAc, result1);
375             }
376           
377           if (sumA < SMALL)
378             {
379               zeroVector( &(myDualMesh->internalDualFaceArea[faceAi]) );
380             }
381           else
382             {
383               scalarDotVector(&(myDualMesh->internalDualFaceCentres[faceAi]), 1.0/(3.0*sumA), sumAc); //correct faceCentres after estimating them
384               scalarDotVector(&(myDualMesh->internalDualFaceArea[faceAi]), 0.5, sumN);
385             }
386           
387           
388           free(sumN);
389           free(sumAc);
390           free(c);
391           free(n)       ;
392           free(result1);
393           free(result2);
394           
395         }
396       
397     }
398   
399   // correct direction of seg and area vector of dual mesh
400   for(connectivity_int segmenti=0;segmenti<myPrimalMesh->segmentNumber;segmenti++)
401     {    
402       segmentVector = myPrimalMesh->segments[segmenti];
403       
405       subTwoVectors(&result1, myPrimalMesh->vertex[segmentVector[PT2]], myPrimalMesh->vertex[segmentVector[PT1]]);      
406       
407       for(connectivity_int facei=0;facei<myDualMesh->segmentToInternalDualFaceNumber[segmenti];facei++)
408         {
409           face  = myDualMesh->segmentToInternalDualFace[segmenti][facei];
410           
411           dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
412           if(result4<0.0)
413             {
414               for(connectivity_int pointi=0;pointi<DIM3D;pointi++)
415                 {
416                   myDualMesh->internalDualFaceArea[face][pointi] *= -1.0;
417                 }              
418             }
419         }
420       
421       
422 #ifdef DEBUG      
423       
424       dotProduct(&result4, result1, myDualMesh->internalDualFaceArea[face]) ;     
426       debug_print( "%-40s : %ld ( ","VECTOR", segmenti);
427       
428       for(connectivity_int ii=0;ii<DIM3D;ii++)
429         {
430           debug_print( "%lf " , result1[ii]);
431         }
432       
433       debug_print(") ");     
434       
435       debug_print( "%lf ", result4);
436       
437       debug_print("\n");      
438 #endif      
439     }
440   
441   
442   
443   
444   
445   
446   
447 #ifdef DEBUG     
448   for(connectivity_int faceAi=0;faceAi<myDualMesh->internalDualFacesNumber;faceAi++)
449     {          
450       debug_print( "%-40s : %ld %d ( ","myDualMesh->internalDualFaceCentres", faceAi, QUAD);
451       
452       for(connectivity_int ii=0;ii<DIM3D;ii++)
453         {
454           debug_print( "%lf " , myDualMesh->internalDualFaceCentres[faceAi][ii]);
455         }
456       
457       debug_print(")\n");
458       
459       debug_print( "%-40s : %ld ( ","myDualMesh->internalDualFaceArea", faceAi);
460       
461       for(connectivity_int ii=0;ii<DIM3D;ii++)
462         {
463           debug_print( "%lf " , myDualMesh->internalDualFaceArea[faceAi][ii]);
464         }
465       
466       debug_print(")\n");
467     }     
468 #endif          
469   
470  
471   
472   free(result1);
473   free(result2);
474   free(fC);    
475   endFunction(__FUNCTION__, t);
476