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