"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/Pos_Formulation.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

Pos_Formulation.cpp  (getdp-3.4.0-source.tgz):Pos_Formulation.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include "ProData.h" #include "ProData.h"
#include "DofData.h" #include "DofData.h"
#include "GeoData.h" #include "GeoData.h"
#include "Get_DofOfElement.h" #include "Get_DofOfElement.h"
skipping to change at line 32 skipping to change at line 32
#include <gmsh/PView.h> #include <gmsh/PView.h>
#include <gmsh/PViewData.h> #include <gmsh/PViewData.h>
#endif #endif
#include "MallocUtils.h" #include "MallocUtils.h"
#include "SolvingAnalyse.h" #include "SolvingAnalyse.h"
#if defined(HAVE_GSL) #if defined(HAVE_GSL)
#include <gsl/gsl_errno.h> #include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h> #include <gsl/gsl_spline.h>
#endif #endif
#define TWO_PI 6.2831853071795865 #define TWO_PI 6.2831853071795865
extern struct Problem Problem_S ; extern struct Problem Problem_S;
extern struct CurrentData Current ; extern struct CurrentData Current;
extern int Flag_BIN, Flag_GMSH_VERSION ; extern int Flag_BIN, Flag_GMSH_VERSION;
extern char *Name_Path ; extern char *Name_Path;
FILE *PostStream = stdout; FILE *PostStream = stdout;
char PostFileName[256]; char PostFileName[256];
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* P o s _ F e m F o r m u l a t i o n */ /* P o s _ F e m F o r m u l a t i o n */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Pos_FemFormulation(struct Formulation *Formulation_P, void Pos_FemFormulation(struct Formulation *Formulation_P,
struct PostQuantity *NCPQ_P, struct PostQuantity *NCPQ_P, struct PostQuantity *CPQ_P,
struct PostQuantity *CPQ_P, int Order, struct PostSubOperation *PostSubOperation_P)
int Order,
struct PostSubOperation *PostSubOperation_P)
{ {
struct Element Element ; struct Element Element;
struct DefineQuantity *DefineQuantity_P0 ; struct DefineQuantity *DefineQuantity_P0;
struct QuantityStorage *QuantityStorage_P0, QuantityStorage ; struct QuantityStorage *QuantityStorage_P0, QuantityStorage;
List_T *QuantityStorage_L ; List_T *QuantityStorage_L;
int i ; int i;
Get_InitDofOfElement(&Element) ; Get_InitDofOfElement(&Element);
DefineQuantity_P0 = (struct DefineQuantity*) DefineQuantity_P0 =
List_Pointer(Formulation_P->DefineQuantity, 0) ; (struct DefineQuantity *)List_Pointer(Formulation_P->DefineQuantity, 0);
QuantityStorage_L = List_Create(List_Nbr(Formulation_P->DefineQuantity), 1, QuantityStorage_L = List_Create(List_Nbr(Formulation_P->DefineQuantity), 1,
sizeof (struct QuantityStorage) ) ; sizeof(struct QuantityStorage));
for(i = 0 ; i < List_Nbr(Formulation_P->DefineQuantity) ; i++) { for(i = 0; i < List_Nbr(Formulation_P->DefineQuantity); i++) {
QuantityStorage.DefineQuantity = DefineQuantity_P0 + i ; QuantityStorage.DefineQuantity = DefineQuantity_P0 + i;
if(QuantityStorage.DefineQuantity->Type == INTEGRALQUANTITY && if(QuantityStorage.DefineQuantity->Type == INTEGRALQUANTITY &&
QuantityStorage.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < QuantityStorage.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof <
0){ 0) {
QuantityStorage.TypeQuantity = VECTOR ; /* on ne sait pas... */ QuantityStorage.TypeQuantity = VECTOR; /* on ne sait pas... */
} }
else{ else {
QuantityStorage.TypeQuantity = QuantityStorage.TypeQuantity =
((struct FunctionSpace *) ((struct FunctionSpace *)List_Pointer(
List_Pointer(Problem_S.FunctionSpace, Problem_S.FunctionSpace,
(DefineQuantity_P0+i)->FunctionSpaceIndex))->Type ; (DefineQuantity_P0 + i)->FunctionSpaceIndex))
->Type;
} }
QuantityStorage.NumLastElementForFunctionSpace = 0 ; QuantityStorage.NumLastElementForFunctionSpace = 0;
List_Add(QuantityStorage_L, &QuantityStorage) ; List_Add(QuantityStorage_L, &QuantityStorage);
} }
QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L, QuantityStorage_P0 =
0) ; (struct QuantityStorage *)List_Pointer(QuantityStorage_L, 0);
switch (PostSubOperation_P->Type) {
case POP_PRINT : switch(PostSubOperation_P->Type) {
switch (PostSubOperation_P->SubType) { case POP_PRINT:
case PRINT_ONREGION : switch(PostSubOperation_P->SubType) {
case PRINT_ONREGION:
Pos_PrintOnRegion(NCPQ_P, CPQ_P, Order, DefineQuantity_P0, Pos_PrintOnRegion(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
QuantityStorage_P0, PostSubOperation_P) ; QuantityStorage_P0, PostSubOperation_P);
break ; break;
case PRINT_ONELEMENTSOF : case PRINT_ONELEMENTSOF:
case PRINT_ONGRID : case PRINT_ONGRID:
Pos_PrintOnElementsOf(NCPQ_P, CPQ_P, Order, DefineQuantity_P0, Pos_PrintOnElementsOf(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
QuantityStorage_P0, PostSubOperation_P) ; QuantityStorage_P0, PostSubOperation_P);
break ; break;
case PRINT_ONSECTION_1D : case PRINT_ONSECTION_1D:
case PRINT_ONSECTION_2D : case PRINT_ONSECTION_2D:
Pos_PrintOnSection(NCPQ_P, CPQ_P, Order, DefineQuantity_P0, Pos_PrintOnSection(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
QuantityStorage_P0, PostSubOperation_P) ; QuantityStorage_P0, PostSubOperation_P);
break ; break;
case PRINT_ONGRID_0D : case PRINT_ONGRID_0D:
case PRINT_ONGRID_1D : case PRINT_ONGRID_1D:
case PRINT_ONGRID_2D : case PRINT_ONGRID_2D:
case PRINT_ONGRID_3D : case PRINT_ONGRID_3D:
case PRINT_ONGRID_PARAM : case PRINT_ONGRID_PARAM:
Pos_PrintOnGrid(NCPQ_P, CPQ_P, Order, DefineQuantity_P0, Pos_PrintOnGrid(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
QuantityStorage_P0, PostSubOperation_P) ; QuantityStorage_P0, PostSubOperation_P);
break ; break;
case PRINT_WITHARGUMENT : case PRINT_WITHARGUMENT:
Pos_PrintWithArgument(NCPQ_P, CPQ_P, Order, DefineQuantity_P0, Pos_PrintWithArgument(NCPQ_P, CPQ_P, Order, DefineQuantity_P0,
QuantityStorage_P0, PostSubOperation_P) ; QuantityStorage_P0, PostSubOperation_P);
break ;
default :
Message::Error("Unknown Operation type for Print");
break; break;
default: Message::Error("Unknown Operation type for Print"); break;
} }
break ;
case POP_EXPRESSION :
Pos_PrintExpression(PostSubOperation_P);
break; break;
case POP_GROUP : case POP_EXPRESSION: Pos_PrintExpression(PostSubOperation_P); break;
Pos_PrintGroup(PostSubOperation_P);
break;
default : case POP_GROUP: Pos_PrintGroup(PostSubOperation_P); break;
Message::Error("Unknown PostSubOperation type") ;
break; default: Message::Error("Unknown PostSubOperation type"); break;
} }
List_Delete(QuantityStorage_L); List_Delete(QuantityStorage_L);
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* P o s _ I n i t T i m e S t e p s */ /* P o s _ I n i t T i m e S t e p s */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
int Pos_InitTimeSteps(struct PostSubOperation *PostSubOperation_P) int Pos_InitTimeSteps(struct PostSubOperation *PostSubOperation_P)
{ {
int iTime, NbTimeStep; int iTime, NbTimeStep;
double TOL = 1.e-15; double TOL = 1.e-15;
// last time step only // last time step only
if(PostSubOperation_P->LastTimeStepOnly || if(PostSubOperation_P->LastTimeStepOnly ||
PostSubOperation_P->AppendTimeStepToFileName){ PostSubOperation_P->AppendTimeStepToFileName) {
iTime = List_Nbr(Current.DofData->Solutions) - 1; iTime = List_Nbr(Current.DofData->Solutions) - 1;
List_Reset(PostSubOperation_P->TimeStep_L); List_Reset(PostSubOperation_P->TimeStep_L);
List_Add(PostSubOperation_P->TimeStep_L, &iTime); List_Add(PostSubOperation_P->TimeStep_L, &iTime);
return 1; return 1;
} }
// specific time values or time interval // specific time values or time interval
if(PostSubOperation_P->TimeInterval_Flag || if(PostSubOperation_P->TimeInterval_Flag ||
List_Nbr(PostSubOperation_P->TimeValue_L) || List_Nbr(PostSubOperation_P->TimeValue_L) ||
List_Nbr(PostSubOperation_P->TimeImagValue_L)){ List_Nbr(PostSubOperation_P->TimeImagValue_L)) {
List_Reset(PostSubOperation_P->TimeStep_L); List_Reset(PostSubOperation_P->TimeStep_L);
for(int i = 0; i < List_Nbr(Current.DofData->Solutions); i++){ for(int i = 0; i < List_Nbr(Current.DofData->Solutions); i++) {
Solution *s = (struct Solution*)List_Pointer(Current.DofData->Solutions, i Solution *s =
); (struct Solution *)List_Pointer(Current.DofData->Solutions, i);
int step = s->TimeStep; int step = s->TimeStep;
double time = s->Time, timeImag = s->TimeImag; double time = s->Time, timeImag = s->TimeImag;
if(PostSubOperation_P->TimeInterval_Flag){ if(PostSubOperation_P->TimeInterval_Flag) {
if((time >= PostSubOperation_P->TimeInterval[0]-TOL) && if((time >= PostSubOperation_P->TimeInterval[0] - TOL) &&
(time <= PostSubOperation_P->TimeInterval[1]+TOL)){ (time <= PostSubOperation_P->TimeInterval[1] + TOL)) {
List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int); List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
} }
} }
else{ else {
for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeValue_L); j++){ for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeValue_L); j++) {
double t; double t;
List_Read(PostSubOperation_P->TimeValue_L, j, &t); List_Read(PostSubOperation_P->TimeValue_L, j, &t);
if(fabs(t - time) < TOL){ if(fabs(t - time) < TOL) {
List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int); List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
} }
} }
for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeImagValue_L); j++){ for(int j = 0; j < List_Nbr(PostSubOperation_P->TimeImagValue_L); j++) {
double t; double t;
List_Read(PostSubOperation_P->TimeImagValue_L, j, &t); List_Read(PostSubOperation_P->TimeImagValue_L, j, &t);
if(fabs(t - timeImag) < TOL) if(fabs(t - timeImag) < TOL)
List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int); List_Insert(PostSubOperation_P->TimeStep_L, &step, fcmp_int);
} }
} }
} }
NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L); NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L);
if(NbTimeStep) return NbTimeStep; if(NbTimeStep) return NbTimeStep;
} }
// specific time steps // specific time steps
NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L); NbTimeStep = List_Nbr(PostSubOperation_P->TimeStep_L);
if(!NbTimeStep || !PostSubOperation_P->FrozenTimeStepList){ if(!NbTimeStep || !PostSubOperation_P->FrozenTimeStepList) {
NbTimeStep = List_Nbr(Current.DofData->Solutions); NbTimeStep = List_Nbr(Current.DofData->Solutions);
List_Reset(PostSubOperation_P->TimeStep_L); List_Reset(PostSubOperation_P->TimeStep_L);
for(iTime = 0 ; iTime < NbTimeStep ; iTime++) for(iTime = 0; iTime < NbTimeStep; iTime++)
List_Add(PostSubOperation_P->TimeStep_L, &iTime); List_Add(PostSubOperation_P->TimeStep_L, &iTime);
} }
return NbTimeStep; return NbTimeStep;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* P o s _ I n i t A l l S o l u t i o n s */ /* P o s _ I n i t A l l S o l u t i o n s */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Pos_InitAllSolutions(List_T * TimeStep_L, int Index_TimeStep) void Pos_InitAllSolutions(List_T *TimeStep_L, int Index_TimeStep)
{ {
int TimeStepIndex, k, Num_Solution ; int TimeStepIndex, k, Num_Solution;
List_Read(TimeStep_L, Index_TimeStep, &TimeStepIndex) ; List_Read(TimeStep_L, Index_TimeStep, &TimeStepIndex);
for(k = 0 ; k < Current.NbrSystem ; k++) for(k = 0; k < Current.NbrSystem; k++)
if( (Num_Solution = std::min(List_Nbr((Current.DofData_P0+k)->Solutions)-1, if((Num_Solution =
TimeStepIndex)) >=0 ) std::min(List_Nbr((Current.DofData_P0 + k)->Solutions) - 1,
(Current.DofData_P0+k)->CurrentSolution = (struct Solution*) TimeStepIndex)) >= 0)
List_Pointer((Current.DofData_P0+k)->Solutions, Num_Solution) ; (Current.DofData_P0 + k)->CurrentSolution =
(struct Solution *)List_Pointer((Current.DofData_P0 + k)->Solutions,
if(TimeStepIndex >= 0 && TimeStepIndex < List_Nbr(Current.DofData->Solutions)) Num_Solution);
{
Solution *Solution_P = ((struct Solution*)List_Pointer if(TimeStepIndex >= 0 &&
(Current.DofData->Solutions, TimeStepIndex)); TimeStepIndex < List_Nbr(Current.DofData->Solutions)) {
Current.TimeStep = Solution_P->TimeStep ; Solution *Solution_P = ((struct Solution *)List_Pointer(
Current.Time = Solution_P->Time ; Current.DofData->Solutions, TimeStepIndex));
Current.TimeImag = Solution_P->TimeImag ; Current.TimeStep = Solution_P->TimeStep;
Current.Time = Solution_P->Time;
Current.TimeImag = Solution_P->TimeImag;
} }
else{ // Warning: this can be wrong else { // Warning: this can be wrong
Current.TimeStep = TimeStepIndex; Current.TimeStep = TimeStepIndex;
if(Current.DofData->CurrentSolution){ if(Current.DofData->CurrentSolution) {
Current.Time = Current.DofData->CurrentSolution->Time; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
} }
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* P o s _ R e s a m p l e T i m e */ /* P o s _ R e s a m p l e T i m e */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
skipping to change at line 251 skipping to change at line 247
void Pos_ResampleTime(struct PostOperation *PostOperation_P) void Pos_ResampleTime(struct PostOperation *PostOperation_P)
{ {
Message::Error("ResampleTime requires the GSL"); Message::Error("ResampleTime requires the GSL");
} }
#else #else
void Pos_ResampleTime(struct PostOperation *PostOperation_P) void Pos_ResampleTime(struct PostOperation *PostOperation_P)
{ {
double ResampleTimeStart, ResampleTimeStop, ResampleTimeStep; double ResampleTimeStart, ResampleTimeStop, ResampleTimeStep;
double OriginalStopTime, *OriginalTime_P, *OriginalValueR_P, *OriginalValueI double OriginalStopTime, *OriginalTime_P, *OriginalValueR_P,
_P; *OriginalValueI_P;
double InterpValueRe, InterpValueIm; double InterpValueRe, InterpValueIm;
int OriginalNbrOfSolutions, NewNbrOfSolutions, xLength; int OriginalNbrOfSolutions, NewNbrOfSolutions, xLength;
Solution *Solution_P, Solution_S; Solution *Solution_P, Solution_S;
List_T *NewSolutions_L; List_T *NewSolutions_L;
ResampleTimeStart = PostOperation_P->ResampleTimeStart; ResampleTimeStart = PostOperation_P->ResampleTimeStart;
ResampleTimeStop = PostOperation_P->ResampleTimeStop; ResampleTimeStop = PostOperation_P->ResampleTimeStop;
ResampleTimeStep = PostOperation_P->ResampleTimeStep; ResampleTimeStep = PostOperation_P->ResampleTimeStep;
OriginalNbrOfSolutions = List_Nbr(Current.DofData->Solutions); OriginalNbrOfSolutions = List_Nbr(Current.DofData->Solutions);
OriginalTime_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double)); OriginalTime_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
OriginalValueR_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double)); OriginalValueR_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
if (gSCALAR_SIZE == 2) if(gSCALAR_SIZE == 2)
OriginalValueI_P = (double *)Malloc(OriginalNbrOfSolutions * sizeof(double)) OriginalValueI_P =
; (double *)Malloc(OriginalNbrOfSolutions * sizeof(double));
else else
OriginalValueI_P = NULL; OriginalValueI_P = NULL;
Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, Solution_P = (struct Solution *)List_Pointer(Current.DofData->Solutions,
OriginalNbrOfSolutions-1); OriginalNbrOfSolutions - 1);
OriginalStopTime = Solution_P->Time; OriginalStopTime = Solution_P->Time;
ResampleTimeStop = (OriginalStopTime < ResampleTimeStop) ? OriginalStopTime : ResampleTimeStop =
ResampleTimeStop; (OriginalStopTime < ResampleTimeStop) ? OriginalStopTime : ResampleTimeStop;
for (int i=0; i<OriginalNbrOfSolutions; i++) { for(int i = 0; i < OriginalNbrOfSolutions; i++) {
Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, i); Solution_P = (struct Solution *)List_Pointer(Current.DofData->Solutions, i);
if (!Solution_P->SolutionExist) if(!Solution_P->SolutionExist) Message::Error("Empty solution(s) found");
Message::Error("Empty solution(s) found");
OriginalTime_P[i] = Solution_P->Time; OriginalTime_P[i] = Solution_P->Time;
} }
LinAlg_GetVectorSize(&Solution_P->x, &xLength); LinAlg_GetVectorSize(&Solution_P->x, &xLength);
NewNbrOfSolutions = floor((ResampleTimeStop-ResampleTimeStart) / NewNbrOfSolutions =
ResampleTimeStep) + 1; floor((ResampleTimeStop - ResampleTimeStart) / ResampleTimeStep) + 1;
if (NewNbrOfSolutions < 1) if(NewNbrOfSolutions < 1)
Message::Error("Invalid ResampleTime settings - t_start: %.6g t_stop: %.6g Message::Error(
" "Invalid ResampleTime settings - t_start: %.6g t_stop: %.6g "
"t_sample: %.6g", ResampleTimeStart, ResampleTimeStop, "t_sample: %.6g",
ResampleTimeStep); ResampleTimeStart, ResampleTimeStop, ResampleTimeStep);
NewSolutions_L = List_Create(NewNbrOfSolutions, 1, sizeof(Solution)); NewSolutions_L = List_Create(NewNbrOfSolutions, 1, sizeof(Solution));
for (int i=0; i<NewNbrOfSolutions; i++) { for(int i = 0; i < NewNbrOfSolutions; i++) {
// Create new Solutions list // Create new Solutions list
Solution_S.TimeStep = i ; Solution_S.TimeStep = i;
Solution_S.Time = ResampleTimeStart + i * ResampleTimeStep; Solution_S.Time = ResampleTimeStart + i * ResampleTimeStep;
Solution_S.TimeImag = 0.0; Solution_S.TimeImag = 0.0;
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &Current.DofData->Solver, xLength); LinAlg_CreateVector(&Solution_S.x, &Current.DofData->Solver, xLength);
List_Add(NewSolutions_L, &Solution_S); List_Add(NewSolutions_L, &Solution_S);
} }
for (int i=0; i<xLength; i++) { for(int i = 0; i < xLength; i++) {
if (gSCALAR_SIZE == 1) { if(gSCALAR_SIZE == 1) {
for (int j=0; j<OriginalNbrOfSolutions; j++) { for(int j = 0; j < OriginalNbrOfSolutions; j++) {
Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, Solution_P =
j); (struct Solution *)List_Pointer(Current.DofData->Solutions, j);
LinAlg_GetDoubleInVector(&OriginalValueR_P[j], &Solution_P->x, i); LinAlg_GetDoubleInVector(&OriginalValueR_P[j], &Solution_P->x, i);
} }
gsl_interp_accel *acc = gsl_interp_accel_alloc (); gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOfSo gsl_spline *spline =
lutions); gsl_spline_alloc(gsl_interp_cspline, OriginalNbrOfSolutions);
gsl_spline_init (spline, OriginalTime_P, OriginalValueR_P, OriginalNbrOfSo gsl_spline_init(spline, OriginalTime_P, OriginalValueR_P,
lutions); OriginalNbrOfSolutions);
for (int j=0; j<NewNbrOfSolutions; j++) { for(int j = 0; j < NewNbrOfSolutions; j++) {
Solution_P = (struct Solution*)List_Pointer(NewSolutions_L, j); Solution_P = (struct Solution *)List_Pointer(NewSolutions_L, j);
InterpValueRe = gsl_spline_eval (spline, Solution_P->Time, acc); InterpValueRe = gsl_spline_eval(spline, Solution_P->Time, acc);
LinAlg_SetDoubleInVector(InterpValueRe, &Solution_P->x, i); LinAlg_SetDoubleInVector(InterpValueRe, &Solution_P->x, i);
} }
gsl_spline_free (spline); gsl_spline_free(spline);
gsl_interp_accel_free (acc); gsl_interp_accel_free(acc);
} }
if (gSCALAR_SIZE == 2) { if(gSCALAR_SIZE == 2) {
for (int j=0; j<OriginalNbrOfSolutions; j++) { for(int j = 0; j < OriginalNbrOfSolutions; j++) {
Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, Solution_P =
j); (struct Solution *)List_Pointer(Current.DofData->Solutions, j);
LinAlg_GetComplexInVector(&OriginalValueR_P[j], &OriginalValueI_P[j], &S LinAlg_GetComplexInVector(&OriginalValueR_P[j], &OriginalValueI_P[j],
olution_P->x, i, -1); &Solution_P->x, i, -1);
} }
gsl_interp_accel *accRe = gsl_interp_accel_alloc (); gsl_interp_accel *accRe = gsl_interp_accel_alloc();
gsl_interp_accel *accIm = gsl_interp_accel_alloc (); gsl_interp_accel *accIm = gsl_interp_accel_alloc();
gsl_spline *splineRe = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOf gsl_spline *splineRe =
Solutions); gsl_spline_alloc(gsl_interp_cspline, OriginalNbrOfSolutions);
gsl_spline *splineIm = gsl_spline_alloc (gsl_interp_cspline, OriginalNbrOf gsl_spline *splineIm =
Solutions); gsl_spline_alloc(gsl_interp_cspline, OriginalNbrOfSolutions);
gsl_spline_init (splineRe, OriginalTime_P, OriginalValueR_P, OriginalNbrOf gsl_spline_init(splineRe, OriginalTime_P, OriginalValueR_P,
Solutions); OriginalNbrOfSolutions);
gsl_spline_init (splineIm, OriginalTime_P, OriginalValueI_P, OriginalNbrOf gsl_spline_init(splineIm, OriginalTime_P, OriginalValueI_P,
Solutions); OriginalNbrOfSolutions);
for (int j=0; j<NewNbrOfSolutions; j++) { for(int j = 0; j < NewNbrOfSolutions; j++) {
Solution_P = (struct Solution*)List_Pointer(NewSolutions_L, j); Solution_P = (struct Solution *)List_Pointer(NewSolutions_L, j);
InterpValueRe = gsl_spline_eval (splineRe, Solution_P->Time, accRe); InterpValueRe = gsl_spline_eval(splineRe, Solution_P->Time, accRe);
InterpValueIm = gsl_spline_eval (splineIm, Solution_P->Time, accIm); InterpValueIm = gsl_spline_eval(splineIm, Solution_P->Time, accIm);
LinAlg_SetComplexInVector(InterpValueRe, InterpValueIm, &Solution_P->x, LinAlg_SetComplexInVector(InterpValueRe, InterpValueIm, &Solution_P->x,
i, -1); i, -1);
} }
gsl_spline_free (splineRe); gsl_spline_free(splineRe);
gsl_spline_free (splineIm); gsl_spline_free(splineIm);
gsl_interp_accel_free (accRe); gsl_interp_accel_free(accRe);
gsl_interp_accel_free (accIm); gsl_interp_accel_free(accIm);
} }
} }
Current.DofData->Solutions = NewSolutions_L; Current.DofData->Solutions = NewSolutions_L;
Current.DofData->CurrentSolution = (struct Solution*) Current.DofData->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(NewSolutions_L, List_Nbr(NewSolutions_L)-1) ; NewSolutions_L, List_Nbr(NewSolutions_L) - 1);
for (int j=0; j<NewNbrOfSolutions; j++) { for(int j = 0; j < NewNbrOfSolutions; j++) {
Solution_P = (struct Solution*)List_Pointer(Current.DofData->Solutions, j); Solution_P = (struct Solution *)List_Pointer(Current.DofData->Solutions, j);
Solution_P->TimeFunctionValues = Get_TimeFunctionValues(Current.DofData) ; Solution_P->TimeFunctionValues = Get_TimeFunctionValues(Current.DofData);
} }
Free(OriginalTime_P); Free(OriginalTime_P);
Free(OriginalValueR_P); Free(OriginalValueR_P);
Free(OriginalValueI_P); Free(OriginalValueI_P);
} }
#endif #endif
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* P o s _ F o r m u l a t i o n */ /* P o s _ F o r m u l a t i o n */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Pos_Formulation(struct Formulation *Formulation_P, void Pos_Formulation(struct Formulation *Formulation_P,
struct PostProcessing *PostProcessing_P, struct PostProcessing *PostProcessing_P,
struct PostSubOperation *PostSubOperation_P) struct PostSubOperation *PostSubOperation_P)
{ {
struct PostQuantity *NCPQ_P = NULL, *CPQ_P = NULL ; struct PostQuantity *NCPQ_P = NULL, *CPQ_P = NULL;
double Pulsation ; double Pulsation;
int i, Order = 0 ; int i, Order = 0;
if(PostSubOperation_P->Type == POP_MERGE){ if(PostSubOperation_P->Type == POP_MERGE) {
Message::SendMergeFileRequest(PostSubOperation_P->FileOut); Message::SendMergeFileRequest(PostSubOperation_P->FileOut);
return; return;
} }
if(PostSubOperation_P->Type == POP_DELETEFILE){ if(PostSubOperation_P->Type == POP_DELETEFILE) {
Message::Info("DeleteFile[%s]", PostSubOperation_P->FileOut) ; Message::Info("DeleteFile[%s]", PostSubOperation_P->FileOut);
RemoveFile(PostSubOperation_P->FileOut); RemoveFile(PostSubOperation_P->FileOut);
return; return;
} }
if(PostSubOperation_P->Type == POP_CREATEDIR){ if(PostSubOperation_P->Type == POP_CREATEDIR) {
Message::Info("CreateDir[%s]", PostSubOperation_P->FileOut) ; Message::Info("CreateDir[%s]", PostSubOperation_P->FileOut);
CreateDirs(PostSubOperation_P->FileOut); CreateDirs(PostSubOperation_P->FileOut);
return; return;
} }
if(PostSubOperation_P->FileOut){ if(PostSubOperation_P->FileOut) {
strcpy(PostFileName, Fix_RelativePath(PostSubOperation_P->FileOut, strcpy(PostFileName,
Name_Path).c_str()); Fix_RelativePath(PostSubOperation_P->FileOut, Name_Path).c_str());
if(PostSubOperation_P->AppendExpressionToFileName >= 0) { if(PostSubOperation_P->AppendExpressionToFileName >= 0) {
struct Value Value ; struct Value Value;
Get_ValueOfExpressionByIndex(PostSubOperation_P->AppendExpressionToFileNam Get_ValueOfExpressionByIndex(
e, PostSubOperation_P->AppendExpressionToFileName, NULL, 0., 0., 0.,
NULL, 0., 0., 0., &Value) ; &Value);
char AddExt[100]; char AddExt[100];
if(PostSubOperation_P->AppendExpressionFormat) if(PostSubOperation_P->AppendExpressionFormat)
sprintf(AddExt, PostSubOperation_P->AppendExpressionFormat, Value.Val[0] sprintf(AddExt, PostSubOperation_P->AppendExpressionFormat,
); Value.Val[0]);
else else
sprintf(AddExt, "%.16g", Value.Val[0]); sprintf(AddExt, "%.16g", Value.Val[0]);
strcat(PostFileName, AddExt); strcat(PostFileName, AddExt);
} }
if(PostSubOperation_P->AppendTimeStepToFileName) { if(PostSubOperation_P->AppendTimeStepToFileName) {
char AddExt[100] ; char AddExt[100];
sprintf(AddExt, "_%03d", (PostSubOperation_P->OverrideTimeStepValue >= 0) sprintf(AddExt, "_%03d",
? (PostSubOperation_P->OverrideTimeStepValue >= 0) ?
PostSubOperation_P->OverrideTimeStepValue : (int)Current.TimeStep) PostSubOperation_P->OverrideTimeStepValue :
; (int)Current.TimeStep);
strcat(PostFileName, AddExt); strcat(PostFileName, AddExt);
} }
if(PostSubOperation_P->AppendStringToFileName) { if(PostSubOperation_P->AppendStringToFileName) {
strcat(PostFileName, PostSubOperation_P->AppendStringToFileName); strcat(PostFileName, PostSubOperation_P->AppendStringToFileName);
} }
if(!strlen(PostFileName) || if(!strlen(PostFileName) ||
(Message::GetIsCommWorld() && Message::GetCommRank())){ (Message::GetIsCommWorld() && Message::GetCommRank())) {
// in parallel mode (SetCommWorld), only rank 0 prints output // in parallel mode (SetCommWorld), only rank 0 prints output
PostStream = NULL ; PostStream = NULL;
} }
else if(!PostSubOperation_P->CatFile) { else if(!PostSubOperation_P->CatFile) {
if((PostStream = FOpen(PostFileName, Flag_BIN ? "wb" : "w"))) if((PostStream = FOpen(PostFileName, Flag_BIN ? "wb" : "w")))
Message::Direct(4, " > '%s'", PostFileName) ; Message::Direct(4, " > '%s'", PostFileName);
else{ else {
Message::Error("Unable to open file '%s'", PostFileName) ; Message::Error("Unable to open file '%s'", PostFileName);
PostStream = stdout ; PostStream = stdout;
} }
} }
else { else {
if((PostStream = FOpen(PostFileName, Flag_BIN ? if((PostStream = FOpen(
(PostSubOperation_P->Format == FORMAT_NXUNV ? "r+b" PostFileName,
: "ab") : Flag_BIN ?
"a"))) (PostSubOperation_P->Format == FORMAT_NXUNV ? "r+b" : "ab") :
Message::Direct(4, " >> '%s'", PostFileName) ; "a")))
else{ Message::Direct(4, " >> '%s'", PostFileName);
Message::Error("Unable to open file '%s'", PostFileName) ; else {
PostStream = stdout ; Message::Error("Unable to open file '%s'", PostFileName);
PostStream = stdout;
} }
} }
} }
else{ else {
PostStream = stdout ; PostStream = stdout;
} }
// force Gmsh version 1 for anything else than OnElementsOf, or if we store in // force Gmsh version 1 for anything else than OnElementsOf, or if we store in
// memory (which requires old-style list ordering) // memory (which requires old-style list ordering)
int oldVersion = Flag_GMSH_VERSION; int oldVersion = Flag_GMSH_VERSION;
if(PostSubOperation_P->Type != POP_PRINT || if(PostSubOperation_P->Type != POP_PRINT ||
PostSubOperation_P->SubType != PRINT_ONELEMENTSOF || PostSubOperation_P->SubType != PRINT_ONELEMENTSOF ||
PostSubOperation_P->Depth != 1 || PostSubOperation_P->Depth != 1 || PostSubOperation_P->StoreInField >= 0)
PostSubOperation_P->StoreInField >= 0)
Flag_GMSH_VERSION = 1; Flag_GMSH_VERSION = 1;
if(PostSubOperation_P->StoreInField >= 0 && if(PostSubOperation_P->StoreInField >= 0 &&
PostSubOperation_P->Format != FORMAT_GMSH) PostSubOperation_P->Format != FORMAT_GMSH)
Message::Warning("StoreInField only available with Gmsh output format"); Message::Warning("StoreInField only available with Gmsh output format");
if(PostSubOperation_P->StoreInMeshBasedField >= 0){ if(PostSubOperation_P->StoreInMeshBasedField >= 0) {
Flag_GMSH_VERSION = 2; Flag_GMSH_VERSION = 2;
if(PostSubOperation_P->SubType != PRINT_ONELEMENTSOF || if(PostSubOperation_P->SubType != PRINT_ONELEMENTSOF ||
PostSubOperation_P->Depth != 1) PostSubOperation_P->Depth != 1)
Message::Error("StoreInMeshBasedField not compatible with selected options Message::Error(
"); "StoreInMeshBasedField not compatible with selected options");
} }
if(PostStream && PostSubOperation_P->CatFile == 2) fprintf(PostStream, "\n\n" if(PostStream && PostSubOperation_P->CatFile == 2)
) ; fprintf(PostStream, "\n\n");
/* two blanks lines for -index in gnuplot */ /* two blanks lines for -index in gnuplot */
Format_PostFormat(PostSubOperation_P) ; Format_PostFormat(PostSubOperation_P);
if(PostSubOperation_P->PostQuantityIndex[0] >= 0) { if(PostSubOperation_P->PostQuantityIndex[0] >= 0) {
if(PostSubOperation_P->PostQuantitySupport[0] < 0) { /* Noncumulative */ if(PostSubOperation_P->PostQuantitySupport[0] < 0) { /* Noncumulative */
NCPQ_P = NCPQ_P = (struct PostQuantity *)List_Pointer(
(struct PostQuantity *) List_Pointer(PostProcessing_P->PostQuantity, PostProcessing_P->PostQuantity,
PostSubOperation_P->PostQuantityInde PostSubOperation_P->PostQuantityIndex[0]);
x[0]) ; CPQ_P = (PostSubOperation_P->PostQuantityIndex[1] >= 0) ?
CPQ_P = (struct PostQuantity *)List_Pointer(
(PostSubOperation_P->PostQuantityIndex[1] >= 0) ? PostProcessing_P->PostQuantity,
(struct PostQuantity *)List_Pointer(PostProcessing_P->PostQuantity, PostSubOperation_P->PostQuantityIndex[1]) :
PostSubOperation_P->PostQuantityIndex NULL;
[1]) : Order = 1;
NULL ;
Order = 1 ;
} }
else { else {
CPQ_P = CPQ_P = (struct PostQuantity *)List_Pointer(
(struct PostQuantity *) List_Pointer(PostProcessing_P->PostQuantity, PostProcessing_P->PostQuantity,
PostSubOperation_P->PostQuantityInde PostSubOperation_P->PostQuantityIndex[0]);
x[0]) ; NCPQ_P = (PostSubOperation_P->PostQuantityIndex[1] >= 0) ?
NCPQ_P = (struct PostQuantity *)List_Pointer(
(PostSubOperation_P->PostQuantityIndex[1] >= 0) ? PostProcessing_P->PostQuantity,
(struct PostQuantity *)List_Pointer(PostProcessing_P->PostQuantity, PostSubOperation_P->PostQuantityIndex[1]) :
PostSubOperation_P->PostQuantityIndex NULL;
[1]) : Order = 0;
NULL ;
Order = 0 ;
} }
} }
if(List_Nbr(PostSubOperation_P->Frequency_L)){ if(List_Nbr(PostSubOperation_P->Frequency_L)) {
if(List_Nbr(PostSubOperation_P->Frequency_L) > List_Nbr(Current.DofData->Pul if(List_Nbr(PostSubOperation_P->Frequency_L) >
sation)) List_Nbr(Current.DofData->Pulsation))
Message::Error("Too many frequencies specified in PostOperation"); Message::Error("Too many frequencies specified in PostOperation");
else{ else {
for(i = 0 ; i < List_Nbr(PostSubOperation_P->Frequency_L) ; i++){ for(i = 0; i < List_Nbr(PostSubOperation_P->Frequency_L); i++) {
Pulsation = *((double *)List_Pointer(PostSubOperation_P->Frequency_L, i) Pulsation =
) * TWO_PI ; *((double *)List_Pointer(PostSubOperation_P->Frequency_L, i)) *
List_Write(Current.DofData->Pulsation, i, &Pulsation) ; TWO_PI;
List_Write(Current.DofData->Pulsation, i, &Pulsation);
} }
} }
} }
switch (Formulation_P->Type) { switch(Formulation_P->Type) {
case FEMEQUATION:
case FEMEQUATION : Pos_FemFormulation(Formulation_P, NCPQ_P, CPQ_P, Order, PostSubOperation_P);
Pos_FemFormulation(Formulation_P, NCPQ_P, CPQ_P, Order, PostSubOperation_P) break;
;
break ;
case GLOBALEQUATION : case GLOBALEQUATION: break;
break ;
default : default:
Message::Error("Unknown Type for Formulation (%s)", Formulation_P->Name) ; Message::Error("Unknown Type for Formulation (%s)", Formulation_P->Name);
break; break;
} }
Flag_GMSH_VERSION = oldVersion; Flag_GMSH_VERSION = oldVersion;
if(PostStream && PostSubOperation_P->FileOut){ if(PostStream && PostSubOperation_P->FileOut) {
fclose(PostStream) ; fclose(PostStream);
if(PostSubOperation_P->SendToServer == NULL || if(PostSubOperation_P->SendToServer == NULL ||
strcmp(PostSubOperation_P->SendToServer, "No")){ strcmp(PostSubOperation_P->SendToServer, "No")) {
if(PostSubOperation_P->Format == FORMAT_GMSH_PARSED || if(PostSubOperation_P->Format == FORMAT_GMSH_PARSED ||
PostSubOperation_P->Format == FORMAT_GMSH){ PostSubOperation_P->Format == FORMAT_GMSH) {
// send merge request // send merge request
Message::SendMergeFileRequest(PostFileName); Message::SendMergeFileRequest(PostFileName);
} }
// Add link to file // Add link to file
Message::AddOnelabStringChoice(Message::GetOnelabClientName() + "/{Output Message::AddOnelabStringChoice(Message::GetOnelabClientName() +
files", "/{Output files",
"file", PostFileName, true, true); "file", PostFileName, true, true);
} }
/* NewCoordinates print option: write a new mesh */ /* NewCoordinates print option: write a new mesh */
if(PostSubOperation_P->NewCoordinates){ if(PostSubOperation_P->NewCoordinates) {
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
GmshMergeFile(std::string(PostFileName)); GmshMergeFile(std::string(PostFileName));
int iview = PView::list.size() - 1; int iview = PView::list.size() - 1;
PViewData *data = PView::list[iview]->getData(); PViewData *data = PView::list[iview]->getData();
GModel* m = new GModel(); GModel *m = new GModel();
m->readMSH(std::string(Current.GeoData->Name)); m->readMSH(std::string(Current.GeoData->Name));
std::vector<GEntity*> entities; std::vector<GEntity *> entities;
m->getEntities(entities); m->getEntities(entities);
std::map<MVertex*, std::vector<double>, MVertexPtrLessThan> newcoords; std::map<MVertex *, std::vector<double>, MVertexPtrLessThan> newcoords;
for(unsigned int i = 0; i < entities.size(); i++) { for(unsigned int i = 0; i < entities.size(); i++) {
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) { for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) {
MVertex* v = entities[i]->mesh_vertices[j]; MVertex *v = entities[i]->mesh_vertices[j];
std::vector<double> xyz(3); std::vector<double> xyz(3);
if(!data->searchVector(v->x(), v->y(), v->z(), &xyz[0])) if(!data->searchVector(v->x(), v->y(), v->z(), &xyz[0]))
Message::Error("Did not find new coordinate Vector at point (%g,%g,% Message::Error(
g) " "Did not find new coordinate Vector at point (%g,%g,%g) "
"from file %s", v->x(), v->y(), v->z(), PostFileName) "from file %s",
; v->x(), v->y(), v->z(), PostFileName);
newcoords[v] = xyz; newcoords[v] = xyz;
} }
} }
for(std::map<MVertex*, std::vector<double>, MVertexPtrLessThan>::iterator for(std::map<MVertex *, std::vector<double>, MVertexPtrLessThan>::iterator
it = newcoords.begin(); it != newcoords.end(); it++) { it = newcoords.begin();
it != newcoords.end(); it++) {
it->first->x() = it->second[0]; it->first->x() = it->second[0];
it->first->y() = it->second[1]; it->first->y() = it->second[1];
it->first->z() = it->second[2]; it->first->z() = it->second[2];
} }
char NewCoordsFileName[256]; char NewCoordsFileName[256];
strcpy(NewCoordsFileName, Fix_RelativePath(PostSubOperation_P->NewCoordina strcpy(NewCoordsFileName,
tesFile, Fix_RelativePath(PostSubOperation_P->NewCoordinatesFile, Name_Path)
Name_Path).c_str()); .c_str());
m->writeMSH(NewCoordsFileName); m->writeMSH(NewCoordsFileName);
Message::Info("Wrote new coordinates in file %s", NewCoordsFileName); Message::Info("Wrote new coordinates in file %s", NewCoordsFileName);
delete m; delete m;
delete PView::list[iview]; delete PView::list[iview];
PView::list.pop_back(); PView::list.pop_back();
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'NewCoord Message::Error(
inates'"); "You need to compile GetDP with Gmsh support to use 'NewCoordinates'");
#endif #endif
} }
} }
} }
 End of changes. 94 change blocks. 
301 lines changed or deleted 282 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)