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 "ProData.h"
9 #include "DofData.h"
10 #include "GeoData.h"
11 #include "Get_DofOfElement.h"
12 #include "Cal_Quantity.h"
13 #include "Pos_Print.h"
14 #include "Pos_Format.h"
15 #include "ListUtils.h"
16 #include "Message.h"
17 #include "OS.h"
18 #if defined(HAVE_GMSH)
19 #include <gmsh/GmshGlobal.h>
20 #include <gmsh/MVertex.h>
21 #include <gmsh/GModel.h>
22 #include <gmsh/PView.h>
23 #include <gmsh/PViewData.h>
24 #endif
25 #include "MallocUtils.h"
26 #include "SolvingAnalyse.h"
27 #if defined(HAVE_GSL)
28 #include <gsl/gsl_errno.h>
29 #include <gsl/gsl_spline.h>
30 #endif
31
32 #define TWO_PI 6.2831853071795865
33
34 extern struct Problem Problem_S ;
35 extern struct CurrentData Current ;
36 extern int Flag_BIN, Flag_GMSH_VERSION ;
37 extern char *Name_Path ;
38
39 FILE *PostStream = stdout;
40 char PostFileName[256];
41
42 /* ------------------------------------------------------------------------ */
43 /* P o s _ F e m F o r m u l a t i o n */
44 /* ------------------------------------------------------------------------ */
45
Pos_FemFormulation(struct Formulation * Formulation_P,struct PostQuantity * NCPQ_P,struct PostQuantity * CPQ_P,int Order,struct PostSubOperation * PostSubOperation_P)46 void Pos_FemFormulation(struct Formulation *Formulation_P,
47 struct PostQuantity *NCPQ_P,
48 struct PostQuantity *CPQ_P,
49 int Order,
50 struct PostSubOperation *PostSubOperation_P)
51 {
52 struct Element Element ;
53 struct DefineQuantity *DefineQuantity_P0 ;
54 struct QuantityStorage *QuantityStorage_P0, QuantityStorage ;
55
56 List_T *QuantityStorage_L ;
57 int i ;
58
59 Get_InitDofOfElement(&Element) ;
60
61 DefineQuantity_P0 = (struct DefineQuantity*)
62 List_Pointer(Formulation_P->DefineQuantity, 0) ;
63 QuantityStorage_L = List_Create(List_Nbr(Formulation_P->DefineQuantity), 1,
64 sizeof (struct QuantityStorage) ) ;
65
66 for(i = 0 ; i < List_Nbr(Formulation_P->DefineQuantity) ; i++) {
67 QuantityStorage.DefineQuantity = DefineQuantity_P0 + i ;
68
69 if(QuantityStorage.DefineQuantity->Type == INTEGRALQUANTITY &&
70 QuantityStorage.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){
71 QuantityStorage.TypeQuantity = VECTOR ; /* on ne sait pas... */
72 }
73 else{
74 QuantityStorage.TypeQuantity =
75 ((struct FunctionSpace *)
76 List_Pointer(Problem_S.FunctionSpace,
77 (DefineQuantity_P0+i)->FunctionSpaceIndex))->Type ;
78 }
79
80 QuantityStorage.NumLastElementForFunctionSpace = 0 ;
81 List_Add(QuantityStorage_L, &QuantityStorage) ;
82 }
83
84 QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L, 0) ;
85
86 switch (PostSubOperation_P->Type) {
87
88 case POP_PRINT :
89 switch (PostSubOperation_P->SubType) {
90 case PRINT_ONREGION :
91 Pos_PrintOnRegion(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
92 QuantityStorage_P0, PostSubOperation_P) ;
93 break ;
94 case PRINT_ONELEMENTSOF :
95 case PRINT_ONGRID :
96 Pos_PrintOnElementsOf(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
97 QuantityStorage_P0, PostSubOperation_P) ;
98 break ;
99 case PRINT_ONSECTION_1D :
100 case PRINT_ONSECTION_2D :
101 Pos_PrintOnSection(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
102 QuantityStorage_P0, PostSubOperation_P) ;
103 break ;
104 case PRINT_ONGRID_0D :
105 case PRINT_ONGRID_1D :
106 case PRINT_ONGRID_2D :
107 case PRINT_ONGRID_3D :
108 case PRINT_ONGRID_PARAM :
109 Pos_PrintOnGrid(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
110 QuantityStorage_P0, PostSubOperation_P) ;
111 break ;
112 case PRINT_WITHARGUMENT :
113 Pos_PrintWithArgument(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
114 QuantityStorage_P0, PostSubOperation_P) ;
115 break ;
116 default :
117 Message::Error("Unknown Operation type for Print");
118 break;
119 }
120 break ;
121
122 case POP_EXPRESSION :
123 Pos_PrintExpression(PostSubOperation_P);
124 break;
125
126 case POP_GROUP :
127 Pos_PrintGroup(PostSubOperation_P);
128 break;
129
130 default :
131 Message::Error("Unknown PostSubOperation type") ;
132 break;
133 }
134
135 List_Delete(QuantityStorage_L);
136 }
137
138 /* ------------------------------------------------------------------------ */
139 /* P o s _ I n i t T i m e S t e p s */
140 /* ------------------------------------------------------------------------ */
141
Pos_InitTimeSteps(struct PostSubOperation * PostSubOperation_P)142 int Pos_InitTimeSteps(struct PostSubOperation *PostSubOperation_P)
143 {
144 int iTime, NbTimeStep;
145 double TOL = 1.e-15;
146
147 // last time step only
148 if(PostSubOperation_P->LastTimeStepOnly ||
149 PostSubOperation_P->AppendTimeStepToFileName){
150 iTime = List_Nbr(Current.DofData->Solutions) - 1;
151 List_Reset(PostSubOperation_P->TimeStep_L);
152 List_Add(PostSubOperation_P->TimeStep_L, &iTime);
153 return 1;
154 }
155
156 // specific time values or time interval
157 if(PostSubOperation_P->TimeInterval_Flag ||
158 List_Nbr(PostSubOperation_P->TimeValue_L) ||
159 List_Nbr(PostSubOperation_P->TimeImagValue_L)){
160 List_Reset(PostSubOperation_P->TimeStep_L);
161 for(int i = 0; i < List_Nbr(Current.DofData->Solutions); i++){
162 Solution *s = (struct Solution*)List_Pointer(Current.DofData->Solutions, i);
163 int step = s->TimeStep;
164 double time = s->Time, timeImag = s->TimeImag;
165 if(PostSubOperation_P->TimeInterval_Flag){
166 if((time >= PostSubOperation_P->TimeInterval[0]-TOL) &&
167 (time <= PostSubOperation_P->TimeInterval[1]+TOL)){
168 List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
169 }
170 }
171 else{
172 for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeValue_L); j++){
173 double t;
174 List_Read(PostSubOperation_P->TimeValue_L, j, &t);
175 if(fabs(t - time) < TOL){
176 List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
177 }
178 }
179 for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeImagValue_L); j++){
180 double t;
181 List_Read(PostSubOperation_P->TimeImagValue_L, j, &t);
182 if(fabs(t - timeImag) < TOL)
183 List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
184 }
185 }
186 }
187 NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L);
188 if(NbTimeStep) return NbTimeStep;
189 }
190
191 // specific time steps
192 NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L);
193
194 if(!NbTimeStep || !PostSubOperation_P->FrozenTimeStepList){
195 NbTimeStep = List_Nbr(Current.DofData->Solutions);
196 List_Reset(PostSubOperation_P->TimeStep_L);
197 for(iTime = 0 ; iTime < NbTimeStep ; iTime++)
198 List_Add(PostSubOperation_P->TimeStep_L, &iTime);
199 }
200
201 return NbTimeStep;
202 }
203
204 /* ------------------------------------------------------------------------ */
205 /* P o s _ I n i t A l l S o l u t i o n s */
206 /* ------------------------------------------------------------------------ */
207
Pos_InitAllSolutions(List_T * TimeStep_L,int Index_TimeStep)208 void Pos_InitAllSolutions(List_T * TimeStep_L, int Index_TimeStep)
209 {
210 int TimeStepIndex, k, Num_Solution ;
211
212 List_Read(TimeStep_L, Index_TimeStep, &TimeStepIndex) ;
213
214 for(k = 0 ; k < Current.NbrSystem ; k++)
215 if( (Num_Solution = std::min(List_Nbr((Current.DofData_P0+k)->Solutions)-1,
216 TimeStepIndex)) >=0 )
217 (Current.DofData_P0+k)->CurrentSolution = (struct Solution*)
218 List_Pointer((Current.DofData_P0+k)->Solutions, Num_Solution) ;
219
220 if(TimeStepIndex >= 0 && TimeStepIndex < List_Nbr(Current.DofData->Solutions)){
221 Solution *Solution_P = ((struct Solution*)List_Pointer
222 (Current.DofData->Solutions, TimeStepIndex));
223 Current.TimeStep = Solution_P->TimeStep ;
224 Current.Time = Solution_P->Time ;
225 Current.TimeImag = Solution_P->TimeImag ;
226 }
227 else{ // Warning: this can be wrong
228 Current.TimeStep = TimeStepIndex;
229 if(Current.DofData->CurrentSolution){
230 Current.Time = Current.DofData->CurrentSolution->Time;
231 Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
232 }
233 }
234 }
235
236 /* ------------------------------------------------------------------------ */
237 /* P o s _ R e s a m p l e T i m e */
238 /* ------------------------------------------------------------------------ */
239
240 #if !defined(HAVE_GSL)
241
Pos_ResampleTime(struct PostOperation * PostOperation_P)242 void Pos_ResampleTime(struct PostOperation *PostOperation_P)
243 {
244 Message::Error("ResampleTime requires the GSL");
245 }
246
247 #else
248
Pos_ResampleTime(struct PostOperation * PostOperation_P)249 void Pos_ResampleTime(struct PostOperation *PostOperation_P)
250 {
251 double ResampleTimeStart, ResampleTimeStop, ResampleTimeStep;
252 double OriginalStopTime, *OriginalTime_P, *OriginalValueR_P, *OriginalValueI_P;
253 double InterpValueRe, InterpValueIm;
254 int OriginalNbrOfSolutions, NewNbrOfSolutions, xLength;
255 Solution *Solution_P, Solution_S;
256 List_T *NewSolutions_L;
257
258 ResampleTimeStart = PostOperation_P->ResampleTimeStart;
259 ResampleTimeStop = PostOperation_P->ResampleTimeStop;
260 ResampleTimeStep = PostOperation_P->ResampleTimeStep;
261 OriginalNbrOfSolutions = List_Nbr(Current.DofData->Solutions);
262
263 OriginalTime_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
264 OriginalValueR_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
265 if (gSCALAR_SIZE == 2)
266 OriginalValueI_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
267 else
268 OriginalValueI_P = NULL;
269
270 Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions,
271 OriginalNbrOfSolutions-1);
272 OriginalStopTime = Solution_P->Time;
273 ResampleTimeStop = (OriginalStopTime < ResampleTimeStop) ? OriginalStopTime :
274 ResampleTimeStop;
275 for (int i=0; i<OriginalNbrOfSolutions; i++) {
276 Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, i);
277 if (!Solution_P->SolutionExist)
278 Message::Error("Empty solution(s) found");
279 OriginalTime_P[i] = Solution_P->Time;
280 }
281
282 LinAlg_GetVectorSize(&Solution_P->x, &xLength);
283
284 NewNbrOfSolutions = floor((ResampleTimeStop-ResampleTimeStart) /
285 ResampleTimeStep) + 1;
286 if (NewNbrOfSolutions < 1)
287 Message::Error("Invalid ResampleTime settings - t_start: %.6g t_stop: %.6g "
288 "t_sample: %.6g", ResampleTimeStart, ResampleTimeStop,
289 ResampleTimeStep);
290 NewSolutions_L = List_Create(NewNbrOfSolutions, 1, sizeof(Solution));
291 for (int i=0; i<NewNbrOfSolutions; i++) {
292 // Create new Solutions list
293 Solution_S.TimeStep = i ;
294 Solution_S.Time = ResampleTimeStart + i * ResampleTimeStep;
295 Solution_S.TimeImag = 0.0;
296 Solution_S.SolutionExist = 1 ;
297 LinAlg_CreateVector(&Solution_S.x, &Current.DofData->Solver, xLength);
298 List_Add(NewSolutions_L, &Solution_S);
299 }
300
301 for (int i=0; i<xLength; i++) {
302 if (gSCALAR_SIZE == 1) {
303 for (int j=0; j<OriginalNbrOfSolutions; j++) {
304 Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, j);
305 LinAlg_GetDoubleInVector(&OriginalValueR_P[j], &Solution_P->x, i);
306 }
307 gsl_interp_accel *acc = gsl_interp_accel_alloc ();
308 gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOfSolutions);
309 gsl_spline_init (spline, OriginalTime_P, OriginalValueR_P, OriginalNbrOfSolutions);
310
311 for (int j=0; j<NewNbrOfSolutions; j++) {
312 Solution_P = (struct Solution*)List_Pointer(NewSolutions_L, j);
313 InterpValueRe = gsl_spline_eval (spline, Solution_P->Time, acc);
314 LinAlg_SetDoubleInVector(InterpValueRe, &Solution_P->x, i);
315 }
316
317 gsl_spline_free (spline);
318 gsl_interp_accel_free (acc);
319 }
320 if (gSCALAR_SIZE == 2) {
321 for (int j=0; j<OriginalNbrOfSolutions; j++) {
322 Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, j);
323 LinAlg_GetComplexInVector(&OriginalValueR_P[j], &OriginalValueI_P[j], &Solution_P->x, i, -1);
324 }
325 gsl_interp_accel *accRe = gsl_interp_accel_alloc ();
326 gsl_interp_accel *accIm = gsl_interp_accel_alloc ();
327 gsl_spline *splineRe = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOfSolutions);
328 gsl_spline *splineIm = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOfSolutions);
329 gsl_spline_init (splineRe, OriginalTime_P, OriginalValueR_P, OriginalNbrOfSolutions);
330 gsl_spline_init (splineIm, OriginalTime_P, OriginalValueI_P, OriginalNbrOfSolutions);
331
332 for (int j=0; j<NewNbrOfSolutions; j++) {
333 Solution_P = (struct Solution*)List_Pointer(NewSolutions_L, j);
334 InterpValueRe = gsl_spline_eval (splineRe, Solution_P->Time, accRe);
335 InterpValueIm = gsl_spline_eval (splineIm, Solution_P->Time, accIm);
336 LinAlg_SetComplexInVector(InterpValueRe, InterpValueIm, &Solution_P->x, i, -1);
337 }
338
339 gsl_spline_free (splineRe);
340 gsl_spline_free (splineIm);
341 gsl_interp_accel_free (accRe);
342 gsl_interp_accel_free (accIm);
343 }
344 }
345 Current.DofData->Solutions = NewSolutions_L;
346 Current.DofData->CurrentSolution = (struct Solution*)
347 List_Pointer(NewSolutions_L, List_Nbr(NewSolutions_L)-1) ;
348 for (int j=0; j<NewNbrOfSolutions; j++) {
349 Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, j);
350 Solution_P->TimeFunctionValues = Get_TimeFunctionValues(Current.DofData) ;
351 }
352
353 Free(OriginalTime_P);
354 Free(OriginalValueR_P);
355 Free(OriginalValueI_P);
356 }
357 #endif
358
359 /* ------------------------------------------------------------------------ */
360 /* P o s _ F o r m u l a t i o n */
361 /* ------------------------------------------------------------------------ */
362
Pos_Formulation(struct Formulation * Formulation_P,struct PostProcessing * PostProcessing_P,struct PostSubOperation * PostSubOperation_P)363 void Pos_Formulation(struct Formulation *Formulation_P,
364 struct PostProcessing *PostProcessing_P,
365 struct PostSubOperation *PostSubOperation_P)
366 {
367 struct PostQuantity *NCPQ_P = NULL, *CPQ_P = NULL ;
368 double Pulsation ;
369 int i, Order = 0 ;
370
371 if(PostSubOperation_P->Type == POP_MERGE){
372 Message::SendMergeFileRequest(PostSubOperation_P->FileOut);
373 return;
374 }
375
376 if(PostSubOperation_P->Type == POP_DELETEFILE){
377 Message::Info("DeleteFile[%s]", PostSubOperation_P->FileOut) ;
378 RemoveFile(PostSubOperation_P->FileOut);
379 return;
380 }
381
382 if(PostSubOperation_P->Type == POP_CREATEDIR){
383 Message::Info("CreateDir[%s]", PostSubOperation_P->FileOut) ;
384 CreateDirs(PostSubOperation_P->FileOut);
385 return;
386 }
387
388 if(PostSubOperation_P->FileOut){
389 strcpy(PostFileName, Fix_RelativePath(PostSubOperation_P->FileOut,
390 Name_Path).c_str());
391
392 if(PostSubOperation_P->AppendExpressionToFileName >= 0) {
393 struct Value Value ;
394 Get_ValueOfExpressionByIndex(PostSubOperation_P->AppendExpressionToFileName,
395 NULL, 0., 0., 0., &Value) ;
396 char AddExt[100];
397 if(PostSubOperation_P->AppendExpressionFormat)
398 sprintf(AddExt, PostSubOperation_P->AppendExpressionFormat, Value.Val[0]);
399 else
400 sprintf(AddExt, "%.16g", Value.Val[0]);
401 strcat(PostFileName, AddExt);
402 }
403
404 if(PostSubOperation_P->AppendTimeStepToFileName) {
405 char AddExt[100] ;
406 sprintf(AddExt, "_%03d", (PostSubOperation_P->OverrideTimeStepValue >= 0) ?
407 PostSubOperation_P->OverrideTimeStepValue : (int)Current.TimeStep) ;
408 strcat(PostFileName, AddExt);
409 }
410
411 if(PostSubOperation_P->AppendStringToFileName) {
412 strcat(PostFileName, PostSubOperation_P->AppendStringToFileName);
413 }
414
415 if(!strlen(PostFileName) ||
416 (Message::GetIsCommWorld() && Message::GetCommRank())){
417 // in parallel mode (SetCommWorld), only rank 0 prints output
418 PostStream = NULL ;
419 }
420 else if(!PostSubOperation_P->CatFile) {
421 if((PostStream = FOpen(PostFileName, Flag_BIN ? "wb" : "w")))
422 Message::Direct(4, " > '%s'", PostFileName) ;
423 else{
424 Message::Error("Unable to open file '%s'", PostFileName) ;
425 PostStream = stdout ;
426 }
427 }
428 else {
429 if((PostStream = FOpen(PostFileName, Flag_BIN ?
430 (PostSubOperation_P->Format == FORMAT_NXUNV ? "r+b" : "ab") :
431 "a")))
432 Message::Direct(4, " >> '%s'", PostFileName) ;
433 else{
434 Message::Error("Unable to open file '%s'", PostFileName) ;
435 PostStream = stdout ;
436 }
437 }
438 }
439 else{
440 PostStream = stdout ;
441 }
442
443 // force Gmsh version 1 for anything else than OnElementsOf, or if we store in
444 // memory (which requires old-style list ordering)
445 int oldVersion = Flag_GMSH_VERSION;
446 if(PostSubOperation_P->Type != POP_PRINT ||
447 PostSubOperation_P->SubType != PRINT_ONELEMENTSOF ||
448 PostSubOperation_P->Depth != 1 ||
449 PostSubOperation_P->StoreInField >= 0)
450 Flag_GMSH_VERSION = 1;
451
452 if(PostSubOperation_P->StoreInField >= 0 &&
453 PostSubOperation_P->Format != FORMAT_GMSH)
454 Message::Warning("StoreInField only available with Gmsh output format");
455
456 if(PostSubOperation_P->StoreInMeshBasedField >= 0){
457 Flag_GMSH_VERSION = 2;
458 if(PostSubOperation_P->SubType != PRINT_ONELEMENTSOF ||
459 PostSubOperation_P->Depth != 1)
460 Message::Error("StoreInMeshBasedField not compatible with selected options");
461 }
462
463 if(PostStream && PostSubOperation_P->CatFile == 2) fprintf(PostStream, "\n\n") ;
464 /* two blanks lines for -index in gnuplot */
465
466 Format_PostFormat(PostSubOperation_P) ;
467
468 if(PostSubOperation_P->PostQuantityIndex[0] >= 0) {
469 if(PostSubOperation_P->PostQuantitySupport[0] < 0) { /* Noncumulative */
470 NCPQ_P =
471 (struct PostQuantity *) List_Pointer(PostProcessing_P->PostQuantity,
472 PostSubOperation_P->PostQuantityIndex[0]) ;
473 CPQ_P =
474 (PostSubOperation_P->PostQuantityIndex[1] >= 0) ?
475 (struct PostQuantity *)List_Pointer(PostProcessing_P->PostQuantity,
476 PostSubOperation_P->PostQuantityIndex[1]) :
477 NULL ;
478 Order = 1 ;
479 }
480 else {
481 CPQ_P =
482 (struct PostQuantity *) List_Pointer(PostProcessing_P->PostQuantity,
483 PostSubOperation_P->PostQuantityIndex[0]) ;
484 NCPQ_P =
485 (PostSubOperation_P->PostQuantityIndex[1] >= 0) ?
486 (struct PostQuantity *)List_Pointer(PostProcessing_P->PostQuantity,
487 PostSubOperation_P->PostQuantityIndex[1]) :
488 NULL ;
489 Order = 0 ;
490 }
491 }
492
493 if(List_Nbr(PostSubOperation_P->Frequency_L)){
494 if(List_Nbr(PostSubOperation_P->Frequency_L) > List_Nbr(Current.DofData->Pulsation))
495 Message::Error("Too many frequencies specified in PostOperation");
496 else{
497 for(i = 0 ; i < List_Nbr(PostSubOperation_P->Frequency_L) ; i++){
498 Pulsation = *((double *)List_Pointer(PostSubOperation_P->Frequency_L, i)) * TWO_PI ;
499 List_Write(Current.DofData->Pulsation, i, &Pulsation) ;
500 }
501 }
502 }
503
504 switch (Formulation_P->Type) {
505
506 case FEMEQUATION :
507 Pos_FemFormulation(Formulation_P, NCPQ_P, CPQ_P, Order, PostSubOperation_P) ;
508 break ;
509
510 case GLOBALEQUATION :
511 break ;
512
513 default :
514 Message::Error("Unknown Type for Formulation (%s)", Formulation_P->Name) ;
515 break;
516
517 }
518
519 Flag_GMSH_VERSION = oldVersion;
520
521 if(PostStream && PostSubOperation_P->FileOut){
522 fclose(PostStream) ;
523
524 if(PostSubOperation_P->SendToServer == NULL ||
525 strcmp(PostSubOperation_P->SendToServer, "No")){
526 if(PostSubOperation_P->Format == FORMAT_GMSH_PARSED ||
527 PostSubOperation_P->Format == FORMAT_GMSH){
528 // send merge request
529 Message::SendMergeFileRequest(PostFileName);
530 }
531 // Add link to file
532 Message::AddOnelabStringChoice(Message::GetOnelabClientName() + "/{Output files",
533 "file", PostFileName, true, true);
534 }
535
536 /* NewCoordinates print option: write a new mesh */
537 if(PostSubOperation_P->NewCoordinates){
538
539 #if defined(HAVE_GMSH)
540
541 GmshMergeFile(std::string(PostFileName));
542 int iview = PView::list.size() - 1;
543 PViewData *data = PView::list[iview]->getData();
544
545 GModel* m = new GModel();
546 m->readMSH(std::string(Current.GeoData->Name));
547
548 std::vector<GEntity*> entities;
549 m->getEntities(entities);
550 std::map<MVertex*, std::vector<double>, MVertexPtrLessThan> newcoords;
551 for(unsigned int i = 0; i < entities.size(); i++) {
552 for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) {
553 MVertex* v = entities[i]->mesh_vertices[j];
554 std::vector<double> xyz(3);
555 if(!data->searchVector(v->x(), v->y(), v->z(), &xyz[0]))
556 Message::Error("Did not find new coordinate Vector at point (%g,%g,%g) "
557 "from file %s", v->x(), v->y(), v->z(), PostFileName);
558 newcoords[v] = xyz;
559 }
560 }
561
562 for(std::map<MVertex*, std::vector<double>, MVertexPtrLessThan>::iterator
563 it = newcoords.begin(); it != newcoords.end(); it++) {
564 it->first->x() = it->second[0];
565 it->first->y() = it->second[1];
566 it->first->z() = it->second[2];
567 }
568
569 char NewCoordsFileName[256];
570 strcpy(NewCoordsFileName, Fix_RelativePath(PostSubOperation_P->NewCoordinatesFile,
571 Name_Path).c_str());
572 m->writeMSH(NewCoordsFileName);
573 Message::Info("Wrote new coordinates in file %s", NewCoordsFileName);
574 delete m;
575 delete PView::list[iview];
576 PView::list.pop_back();
577
578 #else
579 Message::Error("You need to compile GetDP with Gmsh support to use 'NewCoordinates'");
580 #endif
581
582 }
583
584 }
585
586 }
587