1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 #include <string.h>
7 #include <math.h>
8 #include <vector>
9 #include "GetDPConfig.h"
10 #include "GeoData.h"
11 #include "ProData.h"
12 #include "ProDefine.h"
13 #include "Pos_Search.h"
14 #include "MallocUtils.h"
15 #include "Message.h"
16 #include "OS.h"
17
18 #if defined(HAVE_GMSH)
19 #include "gmsh.h"
20 #endif
21
22 #define SQU(a) ((a)*(a))
23
24 extern double Flag_ORDER, Flag_MSH_SCALING ;
25
26 FILE * File_GEO ;
27 char *name;
28
29 struct GeoData * CurrentGeoData ;
30
swapBytes(char * array,int size,int n)31 static void swapBytes(char *array, int size, int n)
32 {
33 int i, c;
34
35 char *x = (char *)Malloc(size * sizeof(char)) ;
36 for(i = 0; i < n; i++) {
37 char *a = &array[i * size];
38 memcpy(x, a, size);
39 for(c = 0; c < size; c++)
40 a[size - 1 - c] = x[c];
41 }
42 Free(x);
43 }
44
45 /* ------------------------------------------------------------------------ */
46 /* G e o _ A d d G e o D a t a */
47 /* ------------------------------------------------------------------------ */
48
Geo_AddGeoData(List_T * GeoData_L,char * Name_MshFile,char * Name_DefaultMshFile,char * Name_AdaptFile,char * Name_DefaultAdaptFile)49 int Geo_AddGeoData(List_T * GeoData_L,
50 char * Name_MshFile, char * Name_DefaultMshFile,
51 char * Name_AdaptFile, char * Name_DefaultAdaptFile)
52 {
53 struct GeoData GeoData_S ;
54 int i ;
55
56 if (!Name_MshFile) Name_MshFile = Name_DefaultMshFile ;
57
58 if ((i = List_ISearchSeq(GeoData_L, Name_MshFile, fcmp_GeoData_Name)) < 0) {
59 Message::Info("Loading Geometric data '%s'", Name_MshFile) ;
60 i = List_Nbr(GeoData_L) ;
61 Geo_InitGeoData(&GeoData_S, i, Name_MshFile) ;
62 Geo_OpenFile(Name_MshFile, "rb") ;
63 Geo_ReadFile(&GeoData_S) ;
64 Geo_CloseFile() ;
65
66 if (!Name_AdaptFile) Name_AdaptFile = Name_DefaultAdaptFile ;
67 if (Name_AdaptFile) {
68 Message::Info("Loading Adaptation data '%s'", Name_AdaptFile) ;
69 Geo_OpenFile(Name_AdaptFile, "r") ;
70 Geo_SetCurrentGeoData(&GeoData_S) ;
71 Geo_ReadFileAdapt(&GeoData_S) ;
72 Geo_CloseFile() ;
73 }
74 List_Add(GeoData_L, &GeoData_S) ;
75 }
76
77 return(i) ;
78 }
79
fcmp_GeoData_Name(const void * a,const void * b)80 int fcmp_GeoData_Name(const void * a, const void * b)
81 {
82 return ( strcmp((char *)a, ((struct GeoData *)b)->Name ) ) ;
83 }
84
85 /* ------------------------------------------------------------------------ */
86 /* G e o _ I n i t G e o D a t a */
87 /* ------------------------------------------------------------------------ */
88
Geo_InitGeoData(struct GeoData * GeoData_P,int Num,char * Name)89 void Geo_InitGeoData(struct GeoData * GeoData_P, int Num, char * Name)
90 {
91 GeoData_P->Num = Num ;
92 GeoData_P->Name = Name ;
93
94 GeoData_P->Nodes = NULL ;
95 GeoData_P->Elements = NULL ;
96
97 GeoData_P->NbrElementsWithEdges = GeoData_P->NbrElementsWithFacets = 0 ;
98 GeoData_P->NumCurrentEdge = GeoData_P->NumCurrentFacet = 0 ;
99 GeoData_P->EdgesXNodes =
100 Tree_Create(sizeof(struct Entity2XEntity1), fcmp_E2XE1) ;
101 GeoData_P->FacetsXEdges =
102 Tree_Create(sizeof(struct Entity2XEntity1), fcmp_E2XE1) ;
103
104 GeoData_P->NodesXElements = NULL ;
105 GeoData_P->Normals =
106 Tree_Create(sizeof(struct EntityXVector), fcmp_EXVector) ;
107
108 GeoData_P->GroupForPRE = NULL ;
109
110 GeoData_P->Grid.Init = 0 ;
111
112 GeoData_P->H = GeoData_P->P = NULL ;
113
114 GeoData_P->PeriodicNodes = NULL ;
115 }
116
117 /* ------------------------------------------------------------------------ */
118 /* G e o _ F r e e G e o D a t a */
119 /* ------------------------------------------------------------------------ */
120
Geo_FreeGeoData(struct GeoData * GeoData_P)121 void Geo_FreeGeoData(struct GeoData * GeoData_P)
122 {
123 Message::Debug("Freeing GeoData %d", GeoData_P->Num);
124 List_Delete(GeoData_P->Nodes);
125 if(GeoData_P->Elements){
126 for(int i = 0; i < List_Nbr(GeoData_P->Elements); i++){
127 Geo_Element *e = (Geo_Element*)List_Pointer(GeoData_P->Elements, i);
128 Free(e->NumNodes);
129 Free(e->NumEdges);
130 Free(e->NumFacets);
131 }
132 List_Delete(GeoData_P->Elements);
133 }
134 if(GeoData_P->EdgesXNodes){
135 Tree_Action(GeoData_P->EdgesXNodes, free_E2XE1);
136 Tree_Delete(GeoData_P->EdgesXNodes);
137 }
138 if(GeoData_P->FacetsXEdges){
139 Tree_Action(GeoData_P->FacetsXEdges, free_E2XE1);
140 Tree_Delete(GeoData_P->FacetsXEdges);
141 }
142 if(GeoData_P->NodesXElements){
143 Tree_Action(GeoData_P->NodesXElements, free_E2XE1);
144 Tree_Delete(GeoData_P->NodesXElements);
145 }
146 Tree_Delete(GeoData_P->Normals);
147 List_Delete(GeoData_P->GroupForPRE);
148 Free_SearchGrid(&GeoData_P->Grid);
149 Free(GeoData_P->H);
150 Free(GeoData_P->P);
151 }
152
153 /* ------------------------------------------------------------------------ */
154 /* G e o _ S e t C u r r e n t G e o D a t a B a s e */
155 /* ------------------------------------------------------------------------ */
156
Geo_SetCurrentGeoData(struct GeoData * GeoData_P)157 void Geo_SetCurrentGeoData(struct GeoData * GeoData_P)
158 {
159 CurrentGeoData = GeoData_P ;
160 }
161
162 /* ------------------------------------------------------------------------ */
163 /* G e o _ O p e n F i l e */
164 /* ------------------------------------------------------------------------ */
165
Geo_OpenFile(char * Name,const char * Mode)166 void Geo_OpenFile(char * Name, const char * Mode)
167 {
168 File_GEO = FOpen(Name, Mode) ;
169 name = Name;
170
171 if (!File_GEO) Message::Error("Unable to open file '%s'", Name);
172 }
173
174 /* ------------------------------------------------------------------------ */
175 /* G e o _ C l o s e F i l e */
176 /* ------------------------------------------------------------------------ */
177
Geo_CloseFile(void)178 void Geo_CloseFile(void)
179 {
180 fclose(File_GEO) ;
181 }
182
183 /* ------------------------------------------------------------------------ */
184 /* G e o _ R e a d F i l e */
185 /* ------------------------------------------------------------------------ */
186
Geo_GetNbNodesPerElement(int Type)187 int Geo_GetNbNodesPerElement(int Type)
188 {
189 switch(Type){
190 case POINT_ELEMENT : return 1;
191
192 case LINE : return 2;
193 case TRIANGLE : return 3;
194 case QUADRANGLE : return 4;
195 case TETRAHEDRON : return 4;
196 case HEXAHEDRON : return 8;
197 case PRISM : return 6;
198 case PYRAMID : return 5;
199
200 case LINE_2 : return 3;
201 case TRIANGLE_2 : return 6;
202 case QUADRANGLE_2 : return 9;
203 case TETRAHEDRON_2 : return 10;
204 case HEXAHEDRON_2 : return 27;
205 case PRISM_2 : return 18;
206 case PYRAMID_2 : return 13;
207 case QUADRANGLE_2_8N: return 8;
208 case HEXAHEDRON_2_20N : return 20;
209 case PRISM_2_15N : return 15;
210 case PYRAMID_2_13N : return 13;
211
212 case LINE_3 : return 4;
213 case TRIANGLE_3 : return 10;
214 case QUADRANGLE_3 : return 16;
215 case TETRAHEDRON_3 : return 20;
216 case HEXAHEDRON_3 : return 64;
217 case PRISM_3 : return 40;
218 case PYRAMID_3 : return 30;
219
220 case LINE_4 : return 5;
221 case TRIANGLE_4 : return 15;
222 case QUADRANGLE_4 : return 25;
223 case TETRAHEDRON_4 : return 35;
224 case HEXAHEDRON_4 : return 125;
225 case PRISM_4 : return 75;
226 //case PYRAMID_4 : return 55;
227
228 default :
229 Message::Error("Unknown type of Element");
230 return -1;
231 }
232 }
233
Geo_GetDimOfElement(int Type)234 int Geo_GetDimOfElement(int Type)
235 {
236 switch(Type){
237 case POINT_ELEMENT : return 0;
238
239 case LINE : return 1;
240 case TRIANGLE : return 2;
241 case QUADRANGLE : return 2;
242 case TETRAHEDRON : return 3;
243 case HEXAHEDRON : return 3;
244 case PRISM : return 3;
245 case PYRAMID : return 3;
246
247 case LINE_2 : return 1;
248 case TRIANGLE_2 : return 2;
249 case QUADRANGLE_2 : return 2;
250 case TETRAHEDRON_2 : return 3;
251 case HEXAHEDRON_2 : return 3;
252 case PRISM_2 : return 3;
253 case PYRAMID_2 : return 3;
254 case QUADRANGLE_2_8N: return 2;
255 case HEXAHEDRON_2_20N : return 3;
256 case PRISM_2_15N : return 3;
257 case PYRAMID_2_13N : return 3;
258
259 case LINE_3 : return 1;
260 case TRIANGLE_3 : return 2;
261 case QUADRANGLE_3 : return 2;
262 case TETRAHEDRON_3 : return 3;
263 case HEXAHEDRON_3 : return 3;
264 case PRISM_3 : return 3;
265 case PYRAMID_3 : return 3;
266
267 case LINE_4 : return 1;
268 case TRIANGLE_4 : return 2;
269 case QUADRANGLE_4 : return 2;
270 case TETRAHEDRON_4 : return 3;
271 case HEXAHEDRON_4 : return 3;
272 case PRISM_4 : return 3;
273 //case PYRAMID_4 : return 3;
274 default :
275 Message::Error("Unknown type of Element");
276 return -1;
277 }
278 }
279
Geo_ReverseElement(Geo_Element * Geo_Element)280 void Geo_ReverseElement(Geo_Element *Geo_Element)
281 {
282 int tmp;
283 switch(Geo_Element->Type){
284 case POINT_ELEMENT :
285 break;
286 case LINE :
287 case LINE_2 :
288 tmp = Geo_Element->NumNodes[0];
289 Geo_Element->NumNodes[0] = Geo_Element->NumNodes[1];
290 Geo_Element->NumNodes[1] = tmp;
291 break;
292 case TRIANGLE :
293 tmp = Geo_Element->NumNodes[1];
294 Geo_Element->NumNodes[1] = Geo_Element->NumNodes[2];
295 Geo_Element->NumNodes[2] = tmp;
296 break;
297 case TRIANGLE_2 :
298 tmp = Geo_Element->NumNodes[1];
299 Geo_Element->NumNodes[1] = Geo_Element->NumNodes[2];
300 Geo_Element->NumNodes[2] = tmp;
301 tmp = Geo_Element->NumNodes[3+0];
302 Geo_Element->NumNodes[3] = Geo_Element->NumNodes[5];
303 Geo_Element->NumNodes[5] = tmp;
304 break;
305 case QUADRANGLE :
306 tmp = Geo_Element->NumNodes[1];
307 Geo_Element->NumNodes[1] = Geo_Element->NumNodes[3];
308 Geo_Element->NumNodes[3] = tmp;
309 break;
310 case QUADRANGLE_2 :
311 case QUADRANGLE_2_8N:
312 tmp = Geo_Element->NumNodes[1];
313 Geo_Element->NumNodes[1] = Geo_Element->NumNodes[3];
314 Geo_Element->NumNodes[3] = tmp;
315 tmp = Geo_Element->NumNodes[4+0];
316 Geo_Element->NumNodes[4] = Geo_Element->NumNodes[7];
317 Geo_Element->NumNodes[7] = tmp;
318 tmp = Geo_Element->NumNodes[5];
319 Geo_Element->NumNodes[5] = Geo_Element->NumNodes[6];
320 Geo_Element->NumNodes[6] = tmp;
321 break;
322 default :
323 Message::Error("Unknown type of element (%d) to reverse",
324 Geo_Element->Type);
325 break;
326 }
327 }
328
Geo_SaveMesh(struct GeoData * GeoData_P,List_T * InitialList,char * FileName)329 void Geo_SaveMesh(struct GeoData * GeoData_P, List_T * InitialList, char * FileName)
330 {
331 FILE * file;
332 struct Geo_Node Geo_Node ;
333 struct Geo_Node * Geo_Node_P ;
334 struct Geo_Element Geo_Element ;
335 struct GeoData GeoData ;
336 int fcmp_int(const void * a, const void * b);
337
338 GeoData.Nodes = List_Create(1000, 1000, sizeof(struct Geo_Node)) ;
339 GeoData.Elements = List_Create(1000, 1000, sizeof(struct Geo_Element)) ;
340
341 int maxdim = 0;
342 for (int i = 0 ; i < List_Nbr(GeoData_P->Elements) ; i++) {
343 List_Read(GeoData_P->Elements, i, &Geo_Element) ;
344 maxdim = std::max(maxdim, Geo_GetDimOfElement(Geo_Element.Type));
345 if (List_Search(InitialList, &Geo_Element.Region, fcmp_int) ) {
346 List_Add(GeoData.Elements, &Geo_Element) ;
347 for (int j = 0 ; j < Geo_Element.NbrNodes ; j++)
348 if (!List_Search(GeoData.Nodes,
349 Geo_Node_P = Geo_GetGeoNodeOfNum(Geo_Element.NumNodes[j]),
350 fcmp_Nod) )
351 List_Add(GeoData.Nodes, Geo_Node_P) ;
352 }
353 }
354
355 file = FOpen(FileName,"w");
356 if(file){
357 Message::Error("Could not open file '%s'", FileName);
358 return;
359 }
360 Message::Info("Saving mesh in file \"%s\" (%d nodes, %d elements)",
361 FileName, List_Nbr(GeoData.Nodes), List_Nbr(GeoData.Elements));
362 fprintf(file, "$NOD\n%d\n", List_Nbr(GeoData.Nodes));
363 for (int i = 0 ; i < List_Nbr(GeoData.Nodes) ; i++) {
364 List_Read(GeoData.Nodes, i, &Geo_Node) ;
365 fprintf(file, "%d %.16g %.16g %.16g\n",
366 Geo_Node.Num, Geo_Node.x, Geo_Node.y, Geo_Node.z) ;
367 }
368 fprintf(file, "$ENDNOD\n$ELM\n%d\n", List_Nbr(GeoData.Elements));
369 for (int i = 0 ; i < List_Nbr(GeoData.Elements) ; i++) {
370 List_Read(GeoData.Elements, i, &Geo_Element) ;
371 int Type = GetDP2Gmsh(Geo_Element.Type) ;
372 fprintf(file, "%d %d %d %d %d ",
373 Geo_Element.Num, Type, Geo_Element.Region,
374 Geo_Element.Region, Geo_Element.NbrNodes) ;
375 for (int j = 0 ; j < Geo_Element.NbrNodes ; j++)
376 fprintf(file, "%d ", Geo_Element.NumNodes[j]) ;
377 fprintf(file, "\n") ;
378 }
379 fprintf(file, "$ENDELM\n");
380 fclose(file);
381
382 List_Delete(GeoData.Nodes);
383 List_Delete(GeoData.Elements);
384 }
385
ExtractDoubleQuotedString(const char * str,int len)386 static std::string ExtractDoubleQuotedString(const char *str, int len)
387 {
388 char *c = strstr((char*)str, "\"");
389 if(!c) return "";
390 std::string ret;
391 for(int i = 1; i < len; i++) {
392 if(c[i] == '"' || c[i] == EOF || c[i] == '\n' || c[i] == '\r') break;
393 ret.push_back(c[i]);
394 }
395 return ret;
396 }
397
Geo_SnapNodes(struct GeoData * GeoData_P)398 static void Geo_SnapNodes(struct GeoData * GeoData_P)
399 {
400 // 2D meshes used for axisymmetric calculations should have the minimum
401 // x-coordinate *exactly* equal to 0: snap x coordinate of all vertices of 2D
402 // meshes with fabs(x) < 1e-13 to 0.
403 if(GeoData_P->Dimension != DIM_2D || fabs(GeoData_P->Xmin) > 1e-13) return;
404 int snaps = 0;
405 for(int i = 0; i < List_Nbr(GeoData_P->Nodes); i++) {
406 Geo_Node *n = (Geo_Node*)List_Pointer(GeoData_P->Nodes, i);
407 if(fabs(n->x) < 1e-13) {
408 n->x = 0.;
409 snaps++;
410 }
411 }
412 if(snaps)
413 Message::Info("Snapped x coordinate of %d node%s to 0", snaps,
414 (snaps > 1) ? "s" : "");
415 }
416
Geo_ReadFileWithGmsh(struct GeoData * GeoData_P)417 static void Geo_ReadFileWithGmsh(struct GeoData * GeoData_P)
418 {
419 #if defined(HAVE_GMSH)
420 gmsh::open(name);
421
422 /* N O D E S */
423
424 struct Geo_Node Geo_Node ;
425
426 std::vector<std::size_t> nodeTags;
427 std::vector<double> coord;
428 std::vector<double> parametricCoord;
429 gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord, -1, -1);
430
431 if(GeoData_P->Nodes == NULL)
432 GeoData_P->Nodes = List_Create(nodeTags.size(), 1000, sizeof(struct Geo_Node)) ;
433
434 for(unsigned int i = 0; i < nodeTags.size(); i++){
435 Geo_Node.Num = nodeTags[i];
436 Geo_Node.x = coord[3*i + 0];
437 Geo_Node.y = coord[3*i + 1];
438 Geo_Node.z = coord[3*i + 2];
439 if(Flag_MSH_SCALING != 1.0){
440 Geo_Node.x *= Flag_MSH_SCALING;
441 Geo_Node.y *= Flag_MSH_SCALING;
442 Geo_Node.z *= Flag_MSH_SCALING;
443 }
444 List_Add(GeoData_P->Nodes, &Geo_Node) ;
445
446 if(!i){
447 GeoData_P->Xmin = GeoData_P->Xmax = Geo_Node.x;
448 GeoData_P->Ymin = GeoData_P->Ymax = Geo_Node.y;
449 GeoData_P->Zmin = GeoData_P->Zmax = Geo_Node.z;
450 }
451 else{
452 GeoData_P->Xmin = std::min(GeoData_P->Xmin, Geo_Node.x);
453 GeoData_P->Xmax = std::max(GeoData_P->Xmax, Geo_Node.x);
454 GeoData_P->Ymin = std::min(GeoData_P->Ymin, Geo_Node.y);
455 GeoData_P->Ymax = std::max(GeoData_P->Ymax, Geo_Node.y);
456 GeoData_P->Zmin = std::min(GeoData_P->Zmin, Geo_Node.z);
457 GeoData_P->Zmax = std::max(GeoData_P->Zmax, Geo_Node.z);
458 }
459 }
460
461 if(GeoData_P->Xmin != GeoData_P->Xmax &&
462 GeoData_P->Ymin != GeoData_P->Ymax &&
463 GeoData_P->Zmin != GeoData_P->Zmax)
464 GeoData_P->Dimension = DIM_3D;
465 else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Ymin != GeoData_P->Ymax)
466 GeoData_P->Dimension = DIM_2D;
467 else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Zmin != GeoData_P->Zmax)
468 GeoData_P->Dimension = DIM_2D;
469 else if(GeoData_P->Ymin != GeoData_P->Ymax && GeoData_P->Zmin != GeoData_P->Zmax)
470 GeoData_P->Dimension = DIM_2D;
471 else if(GeoData_P->Xmin != GeoData_P->Xmax)
472 GeoData_P->Dimension = DIM_1D;
473 else if(GeoData_P->Ymin != GeoData_P->Ymax)
474 GeoData_P->Dimension = DIM_1D;
475 else if(GeoData_P->Zmin != GeoData_P->Zmax)
476 GeoData_P->Dimension = DIM_1D;
477 else
478 GeoData_P->Dimension = DIM_0D;
479
480
481 GeoData_P->CharacteristicLength =
482 sqrt(SQU(GeoData_P->Xmax - GeoData_P->Xmin) +
483 SQU(GeoData_P->Ymax - GeoData_P->Ymin) +
484 SQU(GeoData_P->Zmax - GeoData_P->Zmin));
485 if(!GeoData_P->CharacteristicLength)
486 GeoData_P->CharacteristicLength = 1.;
487
488 coord.clear();
489 parametricCoord.clear();
490
491 /* E L E M E N T S */
492
493 struct Geo_Element Geo_Element ;
494
495 std::vector<int> elementTypes;
496 std::vector< std::vector<std::size_t> > elementTags;
497 std::vector< std::vector<std::size_t> > elementNodeTags;
498 gmsh::model::mesh::getElements(elementTypes, elementTags, elementNodeTags, -1, -1);
499
500 std::size_t nbr = 0, maxTag = 0;
501 for(unsigned int i = 0; i < elementTypes.size(); i++){
502 nbr += elementTags[i].size();
503 for(unsigned int j = 0; j < elementTags[i].size(); j++)
504 maxTag = std::max(maxTag, elementTags[i][j]);
505 }
506
507 if (GeoData_P->Elements == NULL)
508 GeoData_P->Elements = List_Create(nbr, 1000, sizeof(struct Geo_Element)) ;
509
510 Geo_Element.NbrEdges = Geo_Element.NbrFacets = 0 ;
511 Geo_Element.NumEdges = Geo_Element.NumFacets = NULL ;
512
513 gmsh::vectorpair dimTags;
514 gmsh::model::getEntities(dimTags, -1);
515 for(unsigned int entity = 0; entity < dimTags.size(); entity++){
516 gmsh::model::mesh::getElements(elementTypes, elementTags, elementNodeTags,
517 dimTags[entity].first, dimTags[entity].second);
518 std::vector<int> physicalsTags;
519 gmsh::model::getPhysicalGroupsForEntity(dimTags[entity].first,
520 dimTags[entity].second, physicalsTags);
521
522 for(unsigned int phys = 0; phys < physicalsTags.size(); phys++){
523 for(unsigned int i = 0; i < elementTypes.size(); i++){
524 Geo_Element.Type = Gmsh2GetDP(elementTypes[i]) ;
525 Geo_Element.NbrNodes = elementNodeTags[i].size()/elementTags[i].size();
526 for(unsigned int j = 0; j < elementTags[i].size(); j++){
527 // if more than one physical group, create new elements (with new
528 // tags) for all additional groups - this is consistent with the
529 // behavior of the old MSH2 file format
530 Geo_Element.Num = (phys == 0) ? elementTags[i][j] : ++maxTag;
531 Geo_Element.Region = physicalsTags[phys];
532 Geo_Element.ElementaryRegion = dimTags[entity].second;
533 Geo_Element.NumNodes = (int *)Malloc(Geo_Element.NbrNodes * sizeof(int)) ;
534 for (int k = 0; k < Geo_Element.NbrNodes; k++)
535 Geo_Element.NumNodes[k] = elementNodeTags[i][Geo_Element.NbrNodes*j + k];
536 if(Geo_Element.Region < 0){
537 Geo_ReverseElement(&Geo_Element);
538 Geo_Element.Region = -Geo_Element.Region;
539 }
540 List_Add(GeoData_P->Elements, &Geo_Element) ;
541 }
542 }
543 }
544 }
545
546 List_Sort(GeoData_P->Elements, fcmp_Elm) ;
547
548 // Keep track of periodic nodes correspondance, if any
549 for(unsigned int entity = 0; entity < dimTags.size(); entity++){
550 int tagMaster;
551 std::vector<std::size_t> nodeTags, nodeTagsMaster;
552 std::vector<double> affineTransform;
553 gmsh::model::mesh::getPeriodicNodes(dimTags[entity].first, dimTags[entity].second,
554 tagMaster, nodeTags, nodeTagsMaster,
555 affineTransform);
556 if(!nodeTags.empty()){
557 if(!GeoData_P->PeriodicNodes)
558 GeoData_P->PeriodicNodes =
559 List_Create(nodeTags.size(), 1000, sizeof(TwoInt));
560 for(std::size_t i = 0; i < nodeTags.size(); i++){
561 TwoInt pair = {(int)nodeTags[i], (int)nodeTagsMaster[i]};
562 List_Add(GeoData_P->PeriodicNodes, &pair);
563 }
564 }
565 }
566
567 Geo_SnapNodes(GeoData_P);
568
569 #else
570 Message::Error("You need to compile GetDP with Gmsh support to open '%s'",
571 name);
572 #endif
573 }
574
Geo_ReadFile(struct GeoData * GeoData_P)575 void Geo_ReadFile(struct GeoData * GeoData_P)
576 {
577 int Nbr, i, j, Type, iDummy, Format, Size, NbTags ;
578 double Version = 1.0 ;
579 struct Geo_Node Geo_Node ;
580 struct Geo_Element Geo_Element ;
581 char String[256] = "" ;
582 int binary = 0, swap = 0;
583 std::map<int, std::vector<int> > entities[4];
584
585 while (1) {
586
587 do {
588 if(!fgets(String, sizeof(String), File_GEO) || feof(File_GEO))
589 break;
590 } while (String[0] != '$') ;
591
592 if (feof(File_GEO))
593 break ;
594
595 /* F O R M A T */
596
597 if(!strncmp(&String[1], "MeshFormat", 10)) {
598
599 if(!fgets(String, sizeof(String), File_GEO)) return;
600 if(sscanf(String, "%lf %d %d", &Version, &Format, &Size) != 3) return;
601 if(Version < 2.0){
602 Message::Error("Unsupported or unknown mesh file version (%g)", Version);
603 return;
604 }
605 else if(Version > 3.1){
606 Geo_ReadFileWithGmsh(GeoData_P);
607 return;
608 }
609
610 if(Format){
611 binary = 1;
612 Message::Info("Mesh is in binary format");
613 int one;
614 if(fread(&one, sizeof(int), 1, File_GEO) != 1) return;
615 if(one != 1){
616 swap = 1;
617 Message::Info("Swapping bytes from binary file");
618 }
619 }
620 }
621
622 /* P H Y S I C A L N A M E S */
623
624 else if(!strncmp(&String[1], "PhysicalNames", 13)) {
625
626 // GetDP does not currently use the name information
627
628 if(!fgets(String, sizeof(String), File_GEO)) return ;
629 int numNames;
630 if(sscanf(String, "%d", &numNames) != 1) return ;
631 for(int i = 0; i < numNames; i++) {
632 int dim = -1, num;
633 if(Version > 2.0){
634 if(fscanf(File_GEO, "%d", &dim) != 1) return ;
635 }
636 if(fscanf(File_GEO, "%d", &num) != 1) return ;
637 if(!fgets(String, sizeof(String), File_GEO)) return ;
638 std::string name = ExtractDoubleQuotedString(String, 256);
639 Message::Debug("Physical group %d (dim %d) has name %s",
640 num, dim, name.c_str());
641 }
642 }
643
644 /* E N T I T I E S */
645
646 else if(!strncmp(&String[1], "Entities", 8)) {
647 if(!fgets(String, sizeof(String), File_GEO)) return;
648 int num[4];
649 if(sscanf(String, "%d %d %d %d", &num[0], &num[1], &num[2], &num[3]) != 4) return;
650 for(int dim = 0; dim < 4; dim++) {
651 for(int j = 0; j < num[dim]; j++) {
652 int num;
653 if(fscanf(File_GEO, "%d", &num) != 1) return;
654 int nbound = 0;
655 if(dim > 0){
656 if(fscanf(File_GEO, "%d", &nbound) != 1) return;
657 for(int k = 0; k < nbound; k++){
658 int dummy;
659 if(fscanf(File_GEO, "%d", &dummy) != 1) return;
660 }
661 }
662 int nphys;
663 if(fscanf(File_GEO, "%d", &nphys) != 1) return;
664 std::vector<int> physicals(nphys);
665 for(int k = 0; k < nphys; k++){
666 if(fscanf(File_GEO, "%d", &physicals[k]) != 1) return;
667 }
668 entities[dim][num] = physicals;
669 if(nphys > 1){
670 Message::Error("GetDP does not support multiple physical groups per element:"
671 " elementary entity %d belongs to %d physical groups",
672 num, nphys);
673 return;
674 }
675 }
676 }
677 }
678
679 /* N O D E S */
680
681 else if (!strncmp(&String[1], "NOE", 3) ||
682 !strncmp(&String[1], "NOD", 3) ||
683 !strncmp(&String[1], "Nodes", 5) ||
684 !strncmp(&String[1], "ParametricNodes", 15)) {
685
686 bool parametric =
687 !strncmp(&String[1], "ParametricNodes", 15) || (Version >= 3.0);
688
689 if(!fgets(String, sizeof(String), File_GEO)) return;
690 if(sscanf(String, "%d", &Nbr) != 1) return;
691 Message::Debug("%d nodes", Nbr);
692
693 if (GeoData_P->Nodes == NULL)
694 GeoData_P->Nodes = List_Create(Nbr, 1000, sizeof(struct Geo_Node)) ;
695
696 for (i = 0 ; i < Nbr ; i++) {
697 if(!parametric){
698 if(!binary){
699 if(fscanf(File_GEO, "%d %lf %lf %lf",
700 &Geo_Node.Num, &Geo_Node.x, &Geo_Node.y, &Geo_Node.z) != 4)
701 return;
702 if(Flag_MSH_SCALING != 1.0){
703 Geo_Node.x *= Flag_MSH_SCALING;
704 Geo_Node.y *= Flag_MSH_SCALING;
705 Geo_Node.z *= Flag_MSH_SCALING;
706 }
707 }
708 else {
709 if(fread(&Geo_Node.Num, sizeof(int), 1, File_GEO) != 1) return;
710 if(swap) swapBytes((char*)&Geo_Node.Num, sizeof(int), 1);
711 double xyz[3];
712 if(fread(xyz, sizeof(double), 3, File_GEO) != 3) return;
713 if(swap) swapBytes((char*)xyz, sizeof(double), 3);
714 Geo_Node.x = xyz[0];
715 Geo_Node.y = xyz[1];
716 Geo_Node.z = xyz[2];
717 if(Flag_MSH_SCALING != 1.0){
718 Geo_Node.x *= Flag_MSH_SCALING;
719 Geo_Node.y *= Flag_MSH_SCALING;
720 Geo_Node.z *= Flag_MSH_SCALING;
721 }
722 }
723 }
724 else{
725 int dim = -1, entity;
726 if(!binary){
727 if(Version < 3.0){
728 if(fscanf(File_GEO, "%d %lf %lf %lf %d %d",
729 &Geo_Node.Num, &Geo_Node.x, &Geo_Node.y, &Geo_Node.z,
730 &dim, &entity) != 6) return;
731 if(Flag_MSH_SCALING != 1.0){
732 Geo_Node.x *= Flag_MSH_SCALING;
733 Geo_Node.y *= Flag_MSH_SCALING;
734 Geo_Node.z *= Flag_MSH_SCALING;
735 }
736 }
737 else{
738 if(fscanf(File_GEO, "%d %lf %lf %lf %d",
739 &Geo_Node.Num, &Geo_Node.x, &Geo_Node.y, &Geo_Node.z,
740 &entity) != 5) return;
741 if(Flag_MSH_SCALING != 1.0){
742 Geo_Node.x *= Flag_MSH_SCALING;
743 Geo_Node.y *= Flag_MSH_SCALING;
744 Geo_Node.z *= Flag_MSH_SCALING;
745 }
746 if(entity){
747 if(fscanf(File_GEO, "%d", &dim) != 1) return;
748 }
749 }
750 }
751 else{
752 if(fread(&Geo_Node.Num, sizeof(int), 1, File_GEO) != 1) return;
753 if(swap) swapBytes((char*)&Geo_Node.Num, sizeof(int), 1);
754 double xyz[3];
755 if(fread(xyz, sizeof(double), 3, File_GEO) != 3) return;
756 if(swap) swapBytes((char*)xyz, sizeof(double), 3);
757 if(Version < 3.0){
758 if(fread(&dim, sizeof(int), 1, File_GEO) != 1) return;
759 if(swap) swapBytes((char*)&dim, sizeof(int), 1);
760 if(fread(&entity, sizeof(int), 1, File_GEO) != 1) return;
761 if(swap) swapBytes((char*)&entity, sizeof(int), 1);
762 }
763 else{
764 if(fread(&entity, sizeof(int), 1, File_GEO) != 1) return;
765 if(swap) swapBytes((char*)&entity, sizeof(int), 1);
766 if(entity){
767 if(fread(&dim, sizeof(int), 1, File_GEO) != 1) return;
768 if(swap) swapBytes((char*)&dim, sizeof(int), 1);
769 }
770 }
771 Geo_Node.x = xyz[0];
772 Geo_Node.y = xyz[1];
773 Geo_Node.z = xyz[2];
774 if(Flag_MSH_SCALING != 1.0){
775 Geo_Node.x *= Flag_MSH_SCALING;
776 Geo_Node.y *= Flag_MSH_SCALING;
777 Geo_Node.z *= Flag_MSH_SCALING;
778 }
779 }
780 if(dim == 1 && (Version < 3.0 || entity)){
781 double u[1];
782 if(!binary){
783 if(fscanf(File_GEO, "%lf", &u[0]) != 1) return;
784 }
785 else{
786 if(fread(u, sizeof(double), 1, File_GEO) != 1) return;
787 if(swap) swapBytes((char*)u, sizeof(double), 1);
788 }
789 }
790 else if(dim == 2 && (Version < 3.0 || entity)){
791 double uv[2];
792 if(!binary){
793 if(fscanf(File_GEO, "%lf %lf", &uv[0], &uv[1]) != 2) return;
794 }
795 else{
796 if(fread(uv, sizeof(double), 2, File_GEO) != 2) return;
797 if(swap) swapBytes((char*)uv, sizeof(double), 2);
798 }
799 }
800 else if(dim == 3 && Version >= 3.0 && entity){
801 double uvw[3];
802 if(!binary){
803 if(fscanf(File_GEO, "%lf %lf %lf", &uvw[0], &uvw[1], &uvw[2]) != 3)
804 return;
805 }
806 else{
807 if(fread(uvw, sizeof(double), 3, File_GEO) != 2) return;
808 if(swap) swapBytes((char*)uvw, sizeof(double), 3);
809 }
810 }
811 }
812 List_Add(GeoData_P->Nodes, &Geo_Node) ;
813 if(!i){
814 GeoData_P->Xmin = GeoData_P->Xmax = Geo_Node.x;
815 GeoData_P->Ymin = GeoData_P->Ymax = Geo_Node.y;
816 GeoData_P->Zmin = GeoData_P->Zmax = Geo_Node.z;
817 }
818 else{
819 GeoData_P->Xmin = std::min(GeoData_P->Xmin, Geo_Node.x);
820 GeoData_P->Xmax = std::max(GeoData_P->Xmax, Geo_Node.x);
821 GeoData_P->Ymin = std::min(GeoData_P->Ymin, Geo_Node.y);
822 GeoData_P->Ymax = std::max(GeoData_P->Ymax, Geo_Node.y);
823 GeoData_P->Zmin = std::min(GeoData_P->Zmin, Geo_Node.z);
824 GeoData_P->Zmax = std::max(GeoData_P->Zmax, Geo_Node.z);
825 }
826 }
827
828 if(GeoData_P->Xmin != GeoData_P->Xmax &&
829 GeoData_P->Ymin != GeoData_P->Ymax &&
830 GeoData_P->Zmin != GeoData_P->Zmax)
831 GeoData_P->Dimension = DIM_3D;
832 else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Ymin != GeoData_P->Ymax)
833 GeoData_P->Dimension = DIM_2D;
834 else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Zmin != GeoData_P->Zmax)
835 GeoData_P->Dimension = DIM_2D;
836 else if(GeoData_P->Ymin != GeoData_P->Ymax && GeoData_P->Zmin != GeoData_P->Zmax)
837 GeoData_P->Dimension = DIM_2D;
838 else if(GeoData_P->Xmin != GeoData_P->Xmax)
839 GeoData_P->Dimension = DIM_1D;
840 else if(GeoData_P->Ymin != GeoData_P->Ymax)
841 GeoData_P->Dimension = DIM_1D;
842 else if(GeoData_P->Zmin != GeoData_P->Zmax)
843 GeoData_P->Dimension = DIM_1D;
844 else
845 GeoData_P->Dimension = DIM_0D;
846
847 GeoData_P->CharacteristicLength =
848 sqrt(SQU(GeoData_P->Xmax - GeoData_P->Xmin) +
849 SQU(GeoData_P->Ymax - GeoData_P->Ymin) +
850 SQU(GeoData_P->Zmax - GeoData_P->Zmin));
851 if(!GeoData_P->CharacteristicLength)
852 GeoData_P->CharacteristicLength = 1.;
853 }
854
855 /* E L E M E N T S */
856
857 else if (!strncmp(&String[1], "ELM", 3) ||
858 !strncmp(&String[1], "Elements", 8)) {
859
860 if(!fgets(String, sizeof(String), File_GEO)) return;
861 if(sscanf(String, "%d", &Nbr) != 1) return;
862 Message::Debug("%d elements", Nbr);
863
864 if (GeoData_P->Elements == NULL)
865 GeoData_P->Elements =
866 List_Create(Nbr, 1000, sizeof(struct Geo_Element)) ;
867
868 Geo_Element.NbrEdges = Geo_Element.NbrFacets = 0 ;
869 Geo_Element.NumEdges = Geo_Element.NumFacets = NULL ;
870
871 if (!binary){
872 for (i = 0 ; i < Nbr ; i++) {
873 if(Version == 1.0){
874 if(fscanf(File_GEO, "%d %d %d %d %d",
875 &Geo_Element.Num, &Type, &Geo_Element.Region,
876 &Geo_Element.ElementaryRegion, &Geo_Element.NbrNodes) != 5)
877 return;
878 Geo_Element.Type = Gmsh2GetDP(Type) ;
879 }
880 else if(Version < 3.0){
881 if(fscanf(File_GEO, "%d %d %d", &Geo_Element.Num, &Type, &NbTags) != 3)
882 return;
883 Geo_Element.Region = Geo_Element.ElementaryRegion = 1;
884 for(j = 0; j < NbTags; j++){
885 if(fscanf(File_GEO, "%d", &iDummy) != 1) return;
886 if(j == 0)
887 Geo_Element.Region = iDummy;
888 else if(j == 1)
889 Geo_Element.ElementaryRegion = iDummy;
890 /* ignore any other tags for now */
891 }
892 Geo_Element.Type = Gmsh2GetDP(Type) ;
893 Geo_Element.NbrNodes = Geo_GetNbNodesPerElement(Geo_Element.Type);
894 }
895 else{
896 int numData;
897 if(fscanf(File_GEO, "%d %d %d %d", &Geo_Element.Num, &Type,
898 &Geo_Element.ElementaryRegion, &numData) != 4) return;
899 Geo_Element.Type = Gmsh2GetDP(Type) ;
900 Geo_Element.NbrNodes = Geo_GetNbNodesPerElement(Geo_Element.Type);
901 std::vector<int> phys = entities[Geo_GetDimOfElement(Geo_Element.Type)]
902 [Geo_Element.ElementaryRegion];
903 if(phys.empty()){
904 Message::Error("No physical group provided for element %d",
905 Geo_Element.Num);
906 return;
907 }
908 else
909 Geo_Element.Region = phys[0];
910
911 /* ignore any other tags for now */
912 for(j = 0; j < numData - Geo_Element.NbrNodes; j++){
913 if(fscanf(File_GEO, "%d", &iDummy) != 1) return;
914 }
915 }
916
917 Geo_Element.NumNodes = (int *)Malloc(Geo_Element.NbrNodes * sizeof(int)) ;
918 for (j = 0 ; j < Geo_Element.NbrNodes ; j++){
919 if(fscanf(File_GEO, "%d", &Geo_Element.NumNodes[j]) != 1)
920 return;
921 }
922 if(Geo_Element.Region < 0){
923 Geo_ReverseElement(&Geo_Element);
924 Geo_Element.Region = -Geo_Element.Region;
925 }
926
927 List_Add(GeoData_P->Elements, &Geo_Element) ;
928 }
929 }
930 else {
931 if(Version < 3.0){
932 int numElementsPartial = 0;
933 while(numElementsPartial < Nbr){
934 int header[3];
935 if(fread(header, sizeof(int), 3, File_GEO) != 3) return;
936 if(swap) swapBytes((char*)header, sizeof(int), 3);
937 Type = header[0];
938 int numElms = header[1];
939 int numTags = header[2];
940 Geo_Element.Type = Gmsh2GetDP(Type) ;
941 Geo_Element.NbrNodes = Geo_GetNbNodesPerElement(Geo_Element.Type);
942 unsigned int n = 1 + numTags + Geo_Element.NbrNodes;
943 int *data = (int *)Malloc(n * sizeof(int)) ;
944 for(i = 0; i < numElms; i++) {
945 if(fread(data, sizeof(int), n, File_GEO) != n) return;
946 if(swap) swapBytes((char*)data, sizeof(int), n);
947 Geo_Element.Num = data[0];
948 Geo_Element.Region = (numTags > 0) ? data[1] : 0;
949 Geo_Element.ElementaryRegion = (numTags > 1) ? data[2] : 0;
950 Geo_Element.NumNodes = (int *)Malloc(Geo_Element.NbrNodes * sizeof(int)) ;
951 for (j = 0 ; j < Geo_Element.NbrNodes ; j++)
952 Geo_Element.NumNodes[j] = data[numTags + 1 + j] ;
953 List_Add(GeoData_P->Elements, &Geo_Element) ;
954 }
955 Free(data);
956 numElementsPartial += numElms;
957 }
958 }
959 else{
960 for (i = 0 ; i < Nbr ; i++) {
961 int numData;
962 if(fread(&Geo_Element.Num, sizeof(int), 1, File_GEO) != 1) return;
963 if(swap) swapBytes((char*)&Geo_Element.Num, sizeof(int), 1);
964 if(fread(&Type, sizeof(int), 1, File_GEO) != 1) return;
965 if(swap) swapBytes((char*)&Type, sizeof(int), 1);
966 if(fread(&Geo_Element.ElementaryRegion, sizeof(int), 1, File_GEO) != 1) return;
967 if(swap) swapBytes((char*)&Geo_Element.ElementaryRegion, sizeof(int), 1);
968 if(fread(&numData, sizeof(int), 1, File_GEO) != 1) return;
969 if(swap) swapBytes((char*)&numData, sizeof(int), 1);
970 std::vector<int> data;
971 if(numData > 0){
972 data.resize(numData);
973 if((int)fread(&data[0], sizeof(int), numData, File_GEO) != numData) return;
974 if(swap) swapBytes((char*)&data[0], sizeof(int), numData);
975 }
976 Geo_Element.Type = Gmsh2GetDP(Type) ;
977 Geo_Element.NbrNodes = Geo_GetNbNodesPerElement(Geo_Element.Type);
978 Geo_Element.NumNodes = (int *)Malloc(Geo_Element.NbrNodes * sizeof(int)) ;
979 if((int)data.size() >= Geo_Element.NbrNodes){
980 for (j = 0 ; j < Geo_Element.NbrNodes ; j++){
981 Geo_Element.NumNodes[j] = data[numData - Geo_Element.NbrNodes + j] ;
982 }
983 }
984 else{
985 Message::Error("Missing node tags in element %d",
986 Geo_Element.Num);
987 return;
988 }
989 std::vector<int> phys = entities[Geo_GetDimOfElement(Geo_Element.Type)]
990 [Geo_Element.ElementaryRegion];
991 if(phys.empty()){
992 Message::Error("No physical group provided for element %d",
993 Geo_Element.Num);
994 return;
995 }
996 else
997 Geo_Element.Region = phys[0];
998
999 if(Geo_Element.Region < 0){
1000 Geo_ReverseElement(&Geo_Element);
1001 Geo_Element.Region = -Geo_Element.Region;
1002 }
1003
1004 List_Add(GeoData_P->Elements, &Geo_Element) ;
1005 }
1006 }
1007 }
1008
1009 List_Sort(GeoData_P->Elements, fcmp_Elm) ;
1010
1011 }
1012
1013 do {
1014 if(!fgets(String, sizeof(String), File_GEO) || feof(File_GEO))
1015 break;
1016 } while (String[0] != '$') ;
1017
1018 } /* while 1 ... */
1019
1020
1021 Geo_SnapNodes(GeoData_P);
1022 }
1023
Geo_ReadFileAdapt(struct GeoData * GeoData_P)1024 void Geo_ReadFileAdapt(struct GeoData * GeoData_P)
1025 {
1026 struct Geo_Element Geo_Element, * Geo_Element_P ;
1027 int Nbr, i, Index_GeoElement ;
1028 double E, H, P, Max_Order = -1.0 ;
1029 char String[256] ;
1030
1031 Nbr = List_Nbr(GeoData_P->Elements) ;
1032
1033 if(!GeoData_P->H){
1034 GeoData_P->H = (double*)Malloc((Nbr+2)*sizeof(double)) ;
1035 for (i = 0 ; i < Nbr ; i++) GeoData_P->H[i+1] = -1.0 ;
1036 }
1037 if(!GeoData_P->P){
1038 GeoData_P->P = (double*)Malloc((Nbr+2)*sizeof(double)) ;
1039 for (i = 0 ; i < Nbr ; i++) GeoData_P->P[i+1] = -1.0 ;
1040 }
1041
1042 while (1) {
1043
1044 do {
1045 if(!fgets(String, sizeof(String), File_GEO) || feof(File_GEO))
1046 break ;
1047 } while (String[0] != '$') ;
1048
1049 if (feof(File_GEO)) break ;
1050
1051 if (!strncmp(&String[1], "Adapt", 5)) {
1052 if(fscanf(File_GEO, "%d", &Nbr) != 1){
1053 Message::Error("Could not read adapation data");
1054 break;
1055 }
1056 for (i = 0 ; i < Nbr ; i++) {
1057 if(fscanf(File_GEO, "%d %lf %lf %lf", &Geo_Element.Num, &E, &H, &P) != 4){
1058 Message::Error("Could not read adaptation data");
1059 break;
1060 }
1061 if(!(Geo_Element_P = (struct Geo_Element *)
1062 List_PQuery(GeoData_P->Elements, &Geo_Element, fcmp_Elm))){
1063 Message::Error("Element %d not found in database", Geo_Element.Num) ;
1064 break;
1065 }
1066 Index_GeoElement = Geo_GetGeoElementIndex(Geo_Element_P) ;
1067 GeoData_P->H[Index_GeoElement+1] = H ;
1068 GeoData_P->P[Index_GeoElement+1] = P ;
1069 if(P > Max_Order) Max_Order = P ;
1070 }
1071 }
1072
1073 do {
1074 if(!fgets(String, sizeof(String), File_GEO) || feof(File_GEO))
1075 break;
1076 } while (String[0] != '$') ;
1077
1078 } /* while 1 ... */
1079
1080 if(Flag_ORDER < 0) Flag_ORDER = Max_Order ;
1081
1082 Message::Info("Maximum interpolation order = %g", Flag_ORDER) ;
1083 }
1084
1085 /* ------------------------------------------------------------------------ */
1086 /* f c m p _ E l m & f c m p _ N o d */
1087 /* ------------------------------------------------------------------------ */
1088
fcmp_Elm(const void * a,const void * b)1089 int fcmp_Elm(const void * a, const void * b)
1090 {
1091 return ((struct Geo_Element *)a)->Num - ((struct Geo_Element *)b)->Num ;
1092 }
1093
fcmp_Nod(const void * a,const void * b)1094 int fcmp_Nod(const void * a, const void * b)
1095 {
1096 return ((struct Geo_Node *)a)->Num - ((struct Geo_Node *)b)->Num ;
1097 }
1098
1099 /* ------------------------------------------------------------------------ */
1100 /* G e o _ G e t N b r G e o E l e m e n t s */
1101 /* ------------------------------------------------------------------------ */
1102
Geo_GetNbrGeoElements(void)1103 int Geo_GetNbrGeoElements(void)
1104 {
1105 return(List_Nbr(CurrentGeoData->Elements)) ;
1106 }
1107
1108 /* ------------------------------------------------------------------------ */
1109 /* G e o _ G e t G e o E l e m e n t */
1110 /* ------------------------------------------------------------------------ */
1111
Geo_GetGeoElement(int Index_Element)1112 struct Geo_Element *Geo_GetGeoElement(int Index_Element)
1113 {
1114 return((struct Geo_Element *)List_Pointer(CurrentGeoData->Elements,
1115 Index_Element)) ;
1116 }
1117
1118 /* ------------------------------------------------------------------------ */
1119 /* G e o _ G e t G e o E l e m e n t I n d e x */
1120 /* ------------------------------------------------------------------------ */
1121
Geo_GetGeoElementIndex(struct Geo_Element * GeoElement)1122 int Geo_GetGeoElementIndex(struct Geo_Element * GeoElement)
1123 {
1124 return(GeoElement -
1125 (struct Geo_Element*)List_Pointer(CurrentGeoData->Elements, 0)) ;
1126 }
1127
1128 /* ------------------------------------------------------------------------ */
1129 /* G e o _ G e t G e o E l e m e n t O f N u m */
1130 /* ------------------------------------------------------------------------ */
1131
Geo_GetGeoElementOfNum(int Num_Element)1132 struct Geo_Element *Geo_GetGeoElementOfNum(int Num_Element)
1133 {
1134 struct Geo_Element elm ;
1135
1136 elm.Num = Num_Element ;
1137
1138 return((struct Geo_Element*)List_PQuery(CurrentGeoData->Elements,
1139 &elm, fcmp_Elm)) ;
1140 }
1141
1142 /* ------------------------------------------------------------------------ */
1143 /* G e o _ G e t N b r G e o N o d e s */
1144 /* ------------------------------------------------------------------------ */
1145
Geo_GetNbrGeoNodes(void)1146 int Geo_GetNbrGeoNodes(void)
1147 {
1148 return(List_Nbr(CurrentGeoData->Nodes)) ;
1149 }
1150
1151 /* ------------------------------------------------------------------------ */
1152 /* G e o _ G e t G e o N o d e */
1153 /* ------------------------------------------------------------------------ */
1154
Geo_GetGeoNode(int Index_Node)1155 struct Geo_Node *Geo_GetGeoNode(int Index_Node)
1156 {
1157 return((struct Geo_Node *)List_Pointer(CurrentGeoData->Nodes,
1158 Index_Node)) ;
1159 }
1160
1161 /* ------------------------------------------------------------------------ */
1162 /* G e o _ G e t G e o N o d e O f N u m */
1163 /* ------------------------------------------------------------------------ */
1164
Geo_GetGeoNodeOfNum(int Num_Node)1165 struct Geo_Node *Geo_GetGeoNodeOfNum(int Num_Node)
1166 {
1167 struct Geo_Node node ;
1168
1169 node.Num = Num_Node ;
1170
1171 return((struct Geo_Node*)List_PQuery(CurrentGeoData->Nodes,
1172 &node, fcmp_Nod)) ;
1173 }
1174
1175 /* ------------------------------------------------------------------------ */
1176 /* G e o _ G e t N o d e s C o o r d i n a t e s */
1177 /* ------------------------------------------------------------------------ */
1178
Geo_GetNodesCoordinates(int Nbr_Node,int * Num_Node,double * x,double * y,double * z)1179 void Geo_GetNodesCoordinates(int Nbr_Node, int * Num_Node,
1180 double * x, double * y, double * z)
1181 {
1182 int i ;
1183 struct Geo_Node Geo_Node, * Geo_Node_P ;
1184
1185 for (i = 0 ; i < Nbr_Node ; i++) {
1186 Geo_Node.Num = abs(Num_Node[i]) ;
1187
1188 if(!(Geo_Node_P = (struct Geo_Node*)
1189 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod))){
1190 Message::Error("Node %d does not exist", Geo_Node.Num) ;
1191 break;
1192 }
1193
1194 x[i] = Geo_Node_P->x ; y[i] = Geo_Node_P->y ; z[i] = Geo_Node_P->z ;
1195 }
1196 }
1197
1198 /* ------------------------------------------------------------------------ */
1199 /* G e o _ S e t N o d e s C o o r d i n a t e s */
1200 /* ------------------------------------------------------------------------ */
1201
Geo_SetNodesCoordinates(int Nbr_Node,int * Num_Node,double * x,double * y,double * z)1202 void Geo_SetNodesCoordinates(int Nbr_Node, int * Num_Node,
1203 double * x, double * y, double * z)
1204 {
1205 int i ;
1206 struct Geo_Node Geo_Node, * Geo_Node_P ;
1207
1208 for (i = 0 ; i < Nbr_Node ; i++) {
1209 Geo_Node.Num = abs(Num_Node[i]) ;
1210
1211 if(!(Geo_Node_P = (struct Geo_Node*)
1212 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod))){
1213 Message::Error("Node %d does not exist", Geo_Node.Num) ;
1214 break;
1215 }
1216
1217 Geo_Node_P->x = x[i] ; Geo_Node_P->y = y[i] ; Geo_Node_P->z = z[i] ;
1218 }
1219 }
1220
Geo_SetNodesCoordinatesX(int Nbr_Node,int * Num_Node,double * x)1221 void Geo_SetNodesCoordinatesX(int Nbr_Node, int * Num_Node,
1222 double *x)
1223 {
1224 int i ;
1225 struct Geo_Node Geo_Node, * Geo_Node_P ;
1226
1227 for (i = 0 ; i < Nbr_Node ; i++) {
1228 Geo_Node.Num = abs(Num_Node[i]) ;
1229
1230 if(!(Geo_Node_P = (struct Geo_Node*)
1231 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod))){
1232 Message::Error("Node %d does not exist", Geo_Node.Num) ;
1233 break;
1234 }
1235
1236 Geo_Node_P->x = x[i] ;
1237 }
1238 }
1239
Geo_SetNodesCoordinatesY(int Nbr_Node,int * Num_Node,double * y)1240 void Geo_SetNodesCoordinatesY(int Nbr_Node, int * Num_Node,
1241 double *y)
1242 {
1243 int i ;
1244 struct Geo_Node Geo_Node, * Geo_Node_P ;
1245
1246 for (i = 0 ; i < Nbr_Node ; i++) {
1247 Geo_Node.Num = abs(Num_Node[i]) ;
1248
1249 if(!(Geo_Node_P = (struct Geo_Node*)
1250 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod))){
1251 Message::Error("Node %d does not exist", Geo_Node.Num) ;
1252 break;
1253 }
1254
1255 Geo_Node_P->y = y[i] ;
1256 }
1257 }
1258
Geo_SetNodesCoordinatesZ(int Nbr_Node,int * Num_Node,double * z)1259 void Geo_SetNodesCoordinatesZ(int Nbr_Node, int * Num_Node,
1260 double *z)
1261 {
1262 int i ;
1263 struct Geo_Node Geo_Node, * Geo_Node_P ;
1264
1265 for (i = 0 ; i < Nbr_Node ; i++) {
1266 Geo_Node.Num = abs(Num_Node[i]) ;
1267
1268 if(!(Geo_Node_P = (struct Geo_Node*)
1269 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod))){
1270 Message::Error("Node %d does not exist", Geo_Node.Num) ;
1271 break;
1272 }
1273
1274 Geo_Node_P->z = z[i] ;
1275 }
1276 }
1277