"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/SolvingOperations.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 ]

SolvingOperations.cpp  (getdp-3.4.0-source.tgz):SolvingOperations.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.
// //
// Contributor(s): // Contributor(s):
// Johan Gyselinck // Johan Gyselinck
// Ruth Sabariego // Ruth Sabariego
// //
#include <string.h> #include <string.h>
skipping to change at line 29 skipping to change at line 29
#include "Cal_Quantity.h" #include "Cal_Quantity.h"
#include "Cal_Value.h" #include "Cal_Value.h"
#include "MovingBand2D.h" #include "MovingBand2D.h"
#include "EigenSolve.h" #include "EigenSolve.h"
#include "Treatment_Formulation.h" #include "Treatment_Formulation.h"
#include "SolvingAnalyse.h" #include "SolvingAnalyse.h"
#include "SolvingOperations.h" #include "SolvingOperations.h"
#include "MallocUtils.h" #include "MallocUtils.h"
#include "OS.h" #include "OS.h"
#include "Message.h" #include "Message.h"
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
#include <gmsh/GmshGlobal.h> #include <gmsh/GmshGlobal.h>
#include <gmsh/PView.h> #include <gmsh/PView.h>
#endif #endif
#define TWO_PI 6.2831853071795865 #define TWO_PI 6.2831853071795865
#define GOLDENRATIO 1.6180339887498948482 #define GOLDENRATIO 1.6180339887498948482
// for performance tests extern struct Problem Problem_S;
#if !defined(WIN32) extern struct CurrentData Current;
//#define TIMER
#endif
extern struct Problem Problem_S ;
extern struct CurrentData Current ;
extern int TreatmentStatus ; extern int TreatmentStatus;
extern int Flag_POS ; extern int Flag_POS;
extern int Flag_RESTART ; extern int Flag_RESTART;
extern int Flag_BIN, Flag_SPLIT ; extern int Flag_BIN, Flag_SPLIT;
extern char *Name_Generic, *Name_Path ; extern char *Name_Generic, *Name_Path;
extern char *Name_MshFile, *Name_ResFile[NBR_MAX_RES] ; extern char *Name_MshFile, *Name_ResFile[NBR_MAX_RES];
extern List_T *GeoData_L ; extern List_T *GeoData_L;
static int Flag_IterativeLoop = 0 ; /* Attention: phase de test */ static int Flag_IterativeLoop = 0; /* Attention: phase de test */
struct Group * Generate_Group = NULL; struct Group *Generate_Group = NULL;
static int Flag_Break = 0; static int Flag_Break = 0;
// For adaptive time stepper (ugly, I know...) // For adaptive time stepper
int Flag_IterativeLoopConverged = 1; int Flag_IterativeLoopConverged = 1;
int Flag_IterativeLoopN = 0; int Flag_IterativeLoopN = 0;
// For IterativeTimeReduction (ugly also...) // For IterativeTimeReduction (ugly also...)
int Flag_NextThetaFixed = 0 ; int Flag_NextThetaFixed = 0;
// For Update // For Update
int Init_Update = 0 ; int Init_Update = 0;
// For Johan's multi-harmonic stuff (even uglier :-)
int Flag_RHS = 0, *DummyDof ;
double **MH_Moving_Matrix = NULL ;
int MHMoving_assemblyType = 0 ; // For multi-harmonic case
int Flag_AddMHMoving = 0 ; //one more :-) int Flag_RHS = 0, *DummyDof;
double **MH_Moving_Matrix = NULL;
Tree_T * DofTree_MH_moving ; int MHMoving_assemblyType = 0;
int Flag_AddMHMoving = 0; // one more :-)
Tree_T *DofTree_MH_moving;
// For extrapolation
int Flag_ExtrapolationOrder = 0; int Flag_ExtrapolationOrder = 0;
// For Cal_SolutionErrorRatio...
// Changed in OPERATION_ITERATIVELOOP
// double reltol = 1e-7;
// double abstol = 1e-5;
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* F r e e _ U n u s e d S o l u t i o n s */ /* F r e e _ U n u s e d S o l u t i o n s */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Free_UnusedSolutions(struct DofData * DofData_P) void Free_UnusedSolutions(struct DofData *DofData_P)
{ {
struct Solution * Solution_P ; struct Solution *Solution_P;
int index = -1; int index = -1;
// We store 1 solution too much (to allow for an imbricated iterative loop) // We store 1 solution too much (to allow for an imbricated iterative loop)
if(!Flag_POS){ if(!Flag_POS) {
switch (Current.TypeTime) { switch(Current.TypeTime) {
case TIME_THETA : case TIME_THETA:
index = List_Nbr(DofData_P->Solutions)-4 ; index = List_Nbr(DofData_P->Solutions) - 4;
// For TimeLoopAdaptive (Trapezoidal) we need 3 past solutions for the pre // For TimeLoopAdaptive (Trapezoidal) we need 3 past solutions for the
dictor // predictor
index = Message::GetOperatingInTimeLoopAdaptive() ? index - 1 : index; index = Message::GetOperatingInTimeLoopAdaptive() ? index - 1 : index;
break; break;
case TIME_GEAR : case TIME_GEAR:
// With -9 we store 7 past solutions (for Gear_6) // With -9 we store 7 past solutions (for Gear_6)
index = List_Nbr(DofData_P->Solutions)-9 ; index = List_Nbr(DofData_P->Solutions) - 9;
break;
case TIME_NEWMARK :
index = List_Nbr(DofData_P->Solutions)-4 ;
break; break;
case TIME_NEWMARK: index = List_Nbr(DofData_P->Solutions) - 4; break;
default: default:
// FIXME: when doing a handmade loop in the pro file - we (@julien.dular // FIXME: when doing a handmade loop in the pro file - we (@julien.dular
// :-) should clearly introduce a parameter for this. // :-) should clearly introduce a parameter for this.
index = List_Nbr(DofData_P->Solutions) - (Flag_ExtrapolationOrder + 10); index = List_Nbr(DofData_P->Solutions) - (Flag_ExtrapolationOrder + 10);
break; break;
} }
if(index >= 0){ if(index >= 0) {
Solution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, index); Solution_P = (struct Solution *)List_Pointer(DofData_P->Solutions, index);
if(Solution_P->SolutionExist){ if(Solution_P->SolutionExist) {
Message::Info("Freeing Solution %d", index); Message::Info("Freeing Solution %d", index);
LinAlg_DestroyVector(&Solution_P->x); LinAlg_DestroyVector(&Solution_P->x);
Free(Solution_P->TimeFunctionValues) ; Free(Solution_P->TimeFunctionValues);
Solution_P->SolutionExist = 0 ; Solution_P->SolutionExist = 0;
} }
} }
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* I n i t _ S y s t e m D a t a */ /* I n i t _ S y s t e m D a t a */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Init_SystemData(struct DofData * DofData_P, int Flag_Jac) void Init_SystemData(struct DofData *DofData_P, int Flag_Jac)
{ {
if (DofData_P->Flag_Init[0] < 1) { if(DofData_P->Flag_Init[0] < 1) {
DofData_P->Flag_Init[0] = 1 ; DofData_P->Flag_Init[0] = 1;
LinAlg_CreateSolver(&DofData_P->Solver, DofData_P->SolverDataFileName) ; LinAlg_CreateSolver(&DofData_P->Solver, DofData_P->SolverDataFileName);
LinAlg_CreateMatrix(&DofData_P->A, &DofData_P->Solver, LinAlg_CreateMatrix(&DofData_P->A, &DofData_P->Solver, DofData_P->NbrDof,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->b, &DofData_P->Solver, DofData_P->NbrDof) ; LinAlg_CreateVector(&DofData_P->b, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->res, &DofData_P->Solver, DofData_P->NbrDof) LinAlg_CreateVector(&DofData_P->res, &DofData_P->Solver, DofData_P->NbrDof);
; LinAlg_CreateVector(&DofData_P->dx, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->dx, &DofData_P->Solver, DofData_P->NbrDof) ;
} }
/* GenerateOnly: Taking advantage of the invariant parts of the matrix in /* GenerateOnly: Taking advantage of the invariant parts of the matrix in
every time-step */ every time-step */
if(DofData_P->Flag_InitOnly[0] == 1){ if(DofData_P->Flag_InitOnly[0] == 1) {
DofData_P->Flag_InitOnly[0] = 2; DofData_P->Flag_InitOnly[0] = 2;
Message::Info("Initializing System {A1,b1}"); Message::Info("Initializing System {A1,b1}");
LinAlg_CreateMatrix(&DofData_P->A1, &DofData_P->Solver, LinAlg_CreateMatrix(&DofData_P->A1, &DofData_P->Solver, DofData_P->NbrDof,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->b1, &DofData_P->Solver, DofData_P->NbrDof) ; LinAlg_CreateVector(&DofData_P->b1, &DofData_P->Solver, DofData_P->NbrDof);
} }
if(DofData_P->Flag_InitOnly[1] == 1){ if(DofData_P->Flag_InitOnly[1] == 1) {
DofData_P->Flag_InitOnly[1] = 2; DofData_P->Flag_InitOnly[1] = 2;
Message::Info("Initializing System {A2,b2}"); Message::Info("Initializing System {A2,b2}");
LinAlg_CreateMatrix(&DofData_P->A2, &DofData_P->Solver, LinAlg_CreateMatrix(&DofData_P->A2, &DofData_P->Solver, DofData_P->NbrDof,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->b2, &DofData_P->Solver, DofData_P->NbrDof) ; LinAlg_CreateVector(&DofData_P->b2, &DofData_P->Solver, DofData_P->NbrDof);
} }
if(DofData_P->Flag_InitOnly[2] == 1){ if(DofData_P->Flag_InitOnly[2] == 1) {
DofData_P->Flag_InitOnly[2] = 2; DofData_P->Flag_InitOnly[2] = 2;
Message::Info("Initializing System {A3,b3}"); Message::Info("Initializing System {A3,b3}");
LinAlg_CreateMatrix(&DofData_P->A3, &DofData_P->Solver, LinAlg_CreateMatrix(&DofData_P->A3, &DofData_P->Solver, DofData_P->NbrDof,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof);
LinAlg_CreateVector(&DofData_P->b3, &DofData_P->Solver, DofData_P->NbrDof) ; LinAlg_CreateVector(&DofData_P->b3, &DofData_P->Solver, DofData_P->NbrDof);
} }
if (DofData_P->Flag_Init[0] < 2 && Flag_Jac) { if(DofData_P->Flag_Init[0] < 2 && Flag_Jac) {
DofData_P->Flag_Init[0] = 2 ; DofData_P->Flag_Init[0] = 2;
LinAlg_CreateMatrix(&DofData_P->Jac, &DofData_P->Solver, LinAlg_CreateMatrix(&DofData_P->Jac, &DofData_P->Solver, DofData_P->NbrDof,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof);
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G e n e r a t e _ S y s t e m */ /* G e n e r a t e _ S y s t e m */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
static void ZeroMatrix(gMatrix *M, gSolver *S, int N) static void ZeroMatrix(gMatrix *M, gSolver *S, int N)
{ {
// We destroy and recreate the matrix to avoid filling-in the mask when // We destroy and recreate the matrix to avoid filling-in the mask when
// generating systems on meshes with changing topologies (remeshing, moving // generating systems on meshes with changing topologies (remeshing, moving
// band, ..., e.g. in time loops) or when constraints are updated. Using // band, ..., e.g. in time loops) or when constraints are updated. Using
// LinAlg_ZeroMatrix preserves the mask from iteration to iteration, which // LinAlg_ZeroMatrix preserves the mask from iteration to iteration, which
// increases memory every time we reassemble. // increases memory every time we reassemble.
LinAlg_DestroyMatrix(M); LinAlg_DestroyMatrix(M);
LinAlg_CreateMatrix(M, S, N, N); LinAlg_CreateMatrix(M, S, N, N);
} }
void ExtrapolatingPolynomial(int degreewanted, struct DofData * DofData_P, void ExtrapolatingPolynomial(int degreewanted, struct DofData *DofData_P,
struct Solution * Solution_S) struct Solution *Solution_S)
{ {
// Build Lagrange Interpolating Polynomial to predict current solution // Build Lagrange Interpolating Polynomial to predict current solution
// ... // ...
// Polynomial degree: // Polynomial degree:
// degree 0 ==> initialization at the previous Time Solution // degree 0 ==> initialization at the previous Time Solution
// degree 1 ==> linear extrapolation from the 2 previous Time Solutions // degree 1 ==> linear extrapolation from the 2 previous Time Solutions
// degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions // degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions
// ... // ...
// NB: The more data points that are used in the interpolation, // NB: The more data points that are used in the interpolation,
// the higher the degree of the resulting polynomial, // the higher the degree of the resulting polynomial,
// thus the greater oscillation it will exhibit between the data points, // thus the greater oscillation it will exhibit between the data points,
// and therefore the poorest the prediction may be ... // and therefore the poorest the prediction may be ...
double * xi; double *xi;
int * jvalid; int *jvalid;
double Pj=1., x = Solution_S->Time; double Pj = 1., x = Solution_S->Time;
//Vector to save the valid previous Time Solutions that will be interpolated // Vector to save the valid previous Time Solutions that will be interpolated
xi = new double[degreewanted+1]; xi = new double[degreewanted + 1];
jvalid = new int[degreewanted+1]; jvalid = new int[degreewanted + 1];
//Select the last Time Solution => degree 0 is achieved at least // Select the last Time Solution => degree 0 is achieved at least
xi[0]=((struct Solution *) xi[0] = ((struct Solution *)List_Pointer(DofData_P->Solutions,
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1))
List_Nbr(DofData_P->Solutions)-1))->Time; ->Time;
jvalid[0]=0; jvalid[0] = 0;
Message::Info("ExtrapolatingPolynomial: Using previous " Message::Info("ExtrapolatingPolynomial: Using previous "
"Theta Time = %g s (TimeStep %d)", "Theta Time = %g s (TimeStep %d)",
xi[0],Solution_S->TimeStep-1-jvalid[0]); xi[0], Solution_S->TimeStep - 1 - jvalid[0]);
int degree=0, j=1; int degree = 0, j = 1;
//Found other valid Time Solutions in order to reach // Found other valid Time Solutions in order to reach
//the desired polynomial degree [degreewanted] (if possible) // the desired polynomial degree [degreewanted] (if possible)
while(degree<degreewanted && j<=List_Nbr(DofData_P->Solutions)-1) { while(degree < degreewanted && j <= List_Nbr(DofData_P->Solutions) - 1) {
xi[degree+1]=((struct Solution *) xi[degree + 1] =
List_Pointer(DofData_P->Solutions, ((struct Solution *)List_Pointer(DofData_P->Solutions,
List_Nbr(DofData_P->Solutions)-1-j))->Time; List_Nbr(DofData_P->Solutions) - 1 - j))
if (xi[degree+1]<xi[degree]){ ->Time;
jvalid[degree+1]=j; if(xi[degree + 1] < xi[degree]) {
degree++; jvalid[degree + 1] = j;
Message::Info("ExtrapolatingPolynomial: Using previous " degree++;
"Theta Time = %g s (TimeStep %d)", Message::Info("ExtrapolatingPolynomial: Using previous "
xi[degree],Solution_S->TimeStep-1-jvalid[degree]); "Theta Time = %g s (TimeStep %d)",
} xi[degree], Solution_S->TimeStep - 1 - jvalid[degree]);
else{
Message::Info("ExtrapolatingPolynomial: Skipping previous "
"Theta Time = %g s (TimeStep %d) [Redundant]"
,xi[degree+1],Solution_S->TimeStep-1-j);
}
j++;
}
if (degree < degreewanted) {
Message::Info("ExtrapolatingPolynomial: "
"Impossible to build polynomial of degree %d "
"(%d usable previous solution%s found "
"while %d needed for degree %d)"
,degreewanted, degree+1
,((degree+1)==1)?"":"s"
,degreewanted+1,degreewanted) ;
} }
else {
Message::Info("ExtrapolatingPolynomial: Skipping previous "
"Theta Time = %g s (TimeStep %d) [Redundant]",
xi[degree + 1], Solution_S->TimeStep - 1 - j);
}
j++;
}
if(degree < degreewanted) {
Message::Info("ExtrapolatingPolynomial: " Message::Info("ExtrapolatingPolynomial: "
">> Polynomial of degree %d " "Impossible to build polynomial of degree %d "
"based on %d previous solution%s (out of %d existing)" "(%d usable previous solution%s found "
,degree, degree+1 "while %d needed for degree %d)",
, ((degree+1)==1)?"":"s" degreewanted, degree + 1, ((degree + 1) == 1) ? "" : "s",
, List_Nbr(DofData_P->Solutions)) ; degreewanted + 1, degreewanted);
}
// Build Lagrange Interpolating Polynomial Coefficients
// http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html Message::Info("ExtrapolatingPolynomial: "
LinAlg_ZeroVector(&Solution_S->x) ; ">> Polynomial of degree %d "
for (int j=0; j<=degree; j++) { "based on %d previous solution%s (out of %d existing)",
Pj=1.; degree, degree + 1, ((degree + 1) == 1) ? "" : "s",
for (int k=0;k<=degree; k++) { List_Nbr(DofData_P->Solutions));
if (j !=k){
Pj=Pj*(x-xi[k])/(xi[j]-xi[k]); // Build Lagrange Interpolating Polynomial Coefficients
} // http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html
} LinAlg_ZeroVector(&Solution_S->x);
LinAlg_AddVectorProdVectorDouble(&Solution_S->x, for(int j = 0; j <= degree; j++) {
&((struct Solution *) Pj = 1.;
List_Pointer(DofData_P->Solutions, for(int k = 0; k <= degree; k++) {
List_Nbr(DofData_P->Solutions)-1-jvalid[j]))->x, if(j != k) { Pj = Pj * (x - xi[k]) / (xi[j] - xi[k]); }
}
LinAlg_AddVectorProdVectorDouble(
&Solution_S->x,
&((struct Solution *)List_Pointer(
DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1 - jvalid[j]))
->x,
Pj, &Solution_S->x); Pj, &Solution_S->x);
} }
delete [] xi; delete[] xi;
delete [] jvalid; delete[] jvalid;
} }
void Generate_System(struct DefineSystem * DefineSystem_P, void Generate_System(struct DefineSystem *DefineSystem_P,
struct DofData * DofData_P, struct DofData *DofData_P, struct DofData *DofData_P0,
struct DofData * DofData_P0, int Flag_Jac, int Flag_Separate, int Flag_Cumulative = 0)
int Flag_Jac, int Flag_Separate,
int Flag_Cumulative = 0)
{ {
int Nbr_Formulation, Index_Formulation, i_TimeStep, iMat ; int Nbr_Formulation, Index_Formulation, i_TimeStep, iMat;
struct Solution * Solution_P, Solution_S ; struct Solution *Solution_P, Solution_S;
struct Formulation * Formulation_P ; struct Formulation *Formulation_P;
/* (Re)creation des liens entre FunctionSpace et DofData: /* (Re)creation des liens entre FunctionSpace et DofData: seuls les FS
seuls les FS n'intervenant pas dans le DD courant peuvent n'intervenant pas dans le DD courant peuvent pointer vers un autre DD */
pointer vers un autre DD */ Init_DofDataInFunctionSpace(1, DofData_P);
Init_DofDataInFunctionSpace(1, DofData_P) ;
if(!DofData_P->Solutions)
if (!DofData_P->Solutions) DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution));
DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution)) ;
i_TimeStep = (int)Current.TimeStep;
i_TimeStep = (int)Current.TimeStep ;
if(!(Solution_P = (struct Solution *)List_PQuery(DofData_P->Solutions,
if (!(Solution_P = (struct Solution*) &i_TimeStep, fcmp_int))) {
List_PQuery(DofData_P->Solutions, &i_TimeStep, fcmp_int))) { Solution_S.TimeStep = (int)Current.TimeStep;
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.Time = Current.Time;
Solution_S.Time = Current.Time ; Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeImag = Current.TimeImag ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_S.SolutionExist = 1;
Solution_S.SolutionExist = 1 ; LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof) ; if(List_Nbr(DofData_P->Solutions)) {
if (List_Nbr(DofData_P->Solutions)) { if(Flag_ExtrapolationOrder <= 0) {
if(Flag_ExtrapolationOrder <= 0){ LinAlg_CopyVector(
LinAlg_CopyVector(&((struct Solution *) &((struct Solution *)List_Pointer(DofData_P->Solutions,
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1))
List_Nbr(DofData_P->Solutions)-1))->x, ->x,
&Solution_S.x) ; &Solution_S.x);
} }
else{ else {
ExtrapolatingPolynomial(Flag_ExtrapolationOrder, DofData_P, &Solution_S) ExtrapolatingPolynomial(Flag_ExtrapolationOrder, DofData_P,
; &Solution_S);
} }
} }
else { else {
LinAlg_ZeroVector(&Solution_S.x) ; LinAlg_ZeroVector(&Solution_S.x);
} }
List_Add(DofData_P->Solutions, &Solution_S) ; List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ; DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
} }
else if (Solution_P != DofData_P->CurrentSolution && !Flag_Separate) { else if(Solution_P != DofData_P->CurrentSolution && !Flag_Separate) {
/* the test on Flag_Separate is necessary for high order time /* the test on Flag_Separate is necessary for high order time schemes, where
schemes, where InitSolution[] gets called multiple times, InitSolution[] gets called multiple times, resulting in multiple stored
resulting in multiple stored solutions with the same TimeStep solutions with the same TimeStep number. Since GenerateSeparate[] is
number. Since GenerateSeparate[] is called outside the time called outside the time loop (i.e. before TimeStep+=1), the List_PQuery
loop (i.e. before TimeStep+=1), the List_PQuery may return (in may return (in an unpredictable way) any of the initial solutions. */
an unpredictable way) any of the initial solutions. */ Message::Error("Incompatible time");
Message::Error("Incompatible time") ;
} }
else{ else {
// fix time values if we recompute the same step (with different time) // fix time values if we recompute the same step (with different time)
Solution_P->Time = Current.Time ; Solution_P->Time = Current.Time;
Solution_P->TimeImag = Current.TimeImag ; Solution_P->TimeImag = Current.TimeImag;
Free(Solution_P->TimeFunctionValues) ; Free(Solution_P->TimeFunctionValues);
Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
} }
if(Flag_Separate){ if(Flag_Separate) {
for(int i = 0; i < List_Nbr(DofData_P->TimeFunctionIndex); i++) for(int i = 0; i < List_Nbr(DofData_P->TimeFunctionIndex); i++)
if(*(int*)List_Pointer(DofData_P->TimeFunctionIndex, i) > 0) if(*(int *)List_Pointer(DofData_P->TimeFunctionIndex, i) > 0)
Message::Warning("Ignored TimeFunction in Constraint for GenerateSeparate Message::Warning(
") ; "Ignored TimeFunction in Constraint for GenerateSeparate");
for(int i = 0; i < List_Nbr(Problem_S.Expression); i++){ for(int i = 0; i < List_Nbr(Problem_S.Expression); i++) {
DofData_P->CurrentSolution->TimeFunctionValues[i] = 1. ; DofData_P->CurrentSolution->TimeFunctionValues[i] = 1.;
} }
if(Current.DofData->Flag_Init[1] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[1] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M1, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M1, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m1); LinAlg_ZeroVector(&Current.DofData->m1);
for(int i = 0; i < List_Nbr(DofData_P->m1s); i++) for(int i = 0; i < List_Nbr(DofData_P->m1s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m1s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m1s, i));
} }
if(Current.DofData->Flag_Init[2] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[2] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M2, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M2, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m2); LinAlg_ZeroVector(&Current.DofData->m2);
for(int i = 0; i < List_Nbr(DofData_P->m2s); i++) for(int i = 0; i < List_Nbr(DofData_P->m2s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m2s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m2s, i));
} }
if(Current.DofData->Flag_Init[3] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[3] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M3, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M3, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m3); LinAlg_ZeroVector(&Current.DofData->m3);
for(int i = 0; i < List_Nbr(DofData_P->m3s); i++) for(int i = 0; i < List_Nbr(DofData_P->m3s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m3s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m3s, i));
} }
if(Current.DofData->Flag_Init[4] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[4] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M4, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M4, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m4); LinAlg_ZeroVector(&Current.DofData->m4);
for(int i = 0; i < List_Nbr(DofData_P->m4s); i++) for(int i = 0; i < List_Nbr(DofData_P->m4s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m4s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m4s, i));
} }
if(Current.DofData->Flag_Init[5] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[5] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M5, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M5, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m5); LinAlg_ZeroVector(&Current.DofData->m5);
for(int i = 0; i < List_Nbr(DofData_P->m5s); i++) for(int i = 0; i < List_Nbr(DofData_P->m5s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m5s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m5s, i));
} }
if(Current.DofData->Flag_Init[6] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[6] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M6, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M6, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m6); LinAlg_ZeroVector(&Current.DofData->m6);
for(int i = 0; i < List_Nbr(DofData_P->m6s); i++) for(int i = 0; i < List_Nbr(DofData_P->m6s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m6s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m6s, i));
} }
if(Current.DofData->Flag_Init[7] && !Flag_Cumulative){ if(Current.DofData->Flag_Init[7] && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->M7, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->M7, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->m7); LinAlg_ZeroVector(&Current.DofData->m7);
for(int i = 0; i < List_Nbr(DofData_P->m7s); i++) for(int i = 0; i < List_Nbr(DofData_P->m7s); i++)
LinAlg_ZeroVector((gVector*)List_Pointer(DofData_P->m7s, i)); LinAlg_ZeroVector((gVector *)List_Pointer(DofData_P->m7s, i));
} }
} }
else{ else {
if(!Current.DofData->Flag_RHS && !Flag_Cumulative){ if(!Current.DofData->Flag_RHS && !Flag_Cumulative) {
ZeroMatrix(&Current.DofData->A, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->A, &Current.DofData->Solver,
Current.DofData->NbrDof); Current.DofData->NbrDof);
} }
if(!Flag_Cumulative){ if(!Flag_Cumulative) {
LinAlg_ZeroVector(&Current.DofData->b) ; LinAlg_ZeroVector(&Current.DofData->b);
if(Current.DofData->Flag_ListOfRHS){ if(Current.DofData->Flag_ListOfRHS) {
for(int i = 0; i < Current.DofData->TotalNumberOfRHS; i++){ for(int i = 0; i < Current.DofData->TotalNumberOfRHS; i++) {
gVector m; gVector m;
LinAlg_CreateVector(&m, &Current.DofData->Solver, LinAlg_CreateVector(&m, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&m); LinAlg_ZeroVector(&m);
Current.DofData->ListOfRHS.push_back(m); Current.DofData->ListOfRHS.push_back(m);
} }
} }
} }
if(DofData_P->Flag_Only){ if(DofData_P->Flag_Only) {
for(int i = 0 ; i < List_Nbr( DofData_P->OnlyTheseMatrices ); i++){ for(int i = 0; i < List_Nbr(DofData_P->OnlyTheseMatrices); i++) {
List_Read(DofData_P->OnlyTheseMatrices, i, &iMat); List_Read(DofData_P->OnlyTheseMatrices, i, &iMat);
if(iMat && !Flag_Cumulative){ if(iMat && !Flag_Cumulative) {
switch(iMat){ switch(iMat) {
case 1 : case 1:
ZeroMatrix(&Current.DofData->A1, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->A1, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->b1) ; LinAlg_ZeroVector(&Current.DofData->b1);
break; break;
case 2 : case 2:
ZeroMatrix(&Current.DofData->A2, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->A2, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->b2) ; LinAlg_ZeroVector(&Current.DofData->b2);
break; break;
case 3 : case 3:
ZeroMatrix(&Current.DofData->A3, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->A3, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->b3) ; LinAlg_ZeroVector(&Current.DofData->b3);
break; break;
} }
} }
} }
} }
} }
if(Flag_Jac && !Flag_Cumulative) if(Flag_Jac && !Flag_Cumulative)
ZeroMatrix(&Current.DofData->Jac, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->Jac, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ; Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex);
for (int i = 0 ; i < Nbr_Formulation ; i++) { for(int i = 0; i < Nbr_Formulation; i++) {
List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ; List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation);
Formulation_P = (struct Formulation*) Formulation_P = (struct Formulation *)List_Pointer(Problem_S.Formulation,
List_Pointer(Problem_S.Formulation, Index_Formulation) ; Index_Formulation);
Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P); Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P);
Treatment_Formulation(Formulation_P) ; Treatment_Formulation(Formulation_P);
} }
if(Flag_Separate){ if(Flag_Separate) {
DofData_P->CurrentSolution->TimeFunctionValues = Get_TimeFunctionValues(DofD DofData_P->CurrentSolution->TimeFunctionValues =
ata_P) ; Get_TimeFunctionValues(DofData_P);
if(DofData_P->Flag_Init[1]){ if(DofData_P->Flag_Init[1]) {
LinAlg_AssembleMatrix(&DofData_P->M1) ; LinAlg_AssembleMatrix(&DofData_P->M1);
LinAlg_AssembleVector(&DofData_P->m1) ; LinAlg_AssembleVector(&DofData_P->m1);
for(int i = 0; i < List_Nbr(DofData_P->m1s); i++) for(int i = 0; i < List_Nbr(DofData_P->m1s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m1s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m1s, i));
} }
if(DofData_P->Flag_Init[2]){ if(DofData_P->Flag_Init[2]) {
LinAlg_AssembleMatrix(&DofData_P->M2) ; LinAlg_AssembleMatrix(&DofData_P->M2);
LinAlg_AssembleVector(&DofData_P->m2) ; LinAlg_AssembleVector(&DofData_P->m2);
for(int i = 0; i < List_Nbr(DofData_P->m2s); i++) for(int i = 0; i < List_Nbr(DofData_P->m2s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m2s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m2s, i));
} }
if(DofData_P->Flag_Init[3]){ if(DofData_P->Flag_Init[3]) {
LinAlg_AssembleMatrix(&DofData_P->M3) ; LinAlg_AssembleMatrix(&DofData_P->M3);
LinAlg_AssembleVector(&DofData_P->m3) ; LinAlg_AssembleVector(&DofData_P->m3);
for(int i = 0; i < List_Nbr(DofData_P->m3s); i++) for(int i = 0; i < List_Nbr(DofData_P->m3s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m3s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m3s, i));
} }
if(DofData_P->Flag_Init[4]){ if(DofData_P->Flag_Init[4]) {
LinAlg_AssembleMatrix(&DofData_P->M4) ; LinAlg_AssembleMatrix(&DofData_P->M4);
LinAlg_AssembleVector(&DofData_P->m4) ; LinAlg_AssembleVector(&DofData_P->m4);
for(int i = 0; i < List_Nbr(DofData_P->m4s); i++) for(int i = 0; i < List_Nbr(DofData_P->m4s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m4s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m4s, i));
} }
if(DofData_P->Flag_Init[5]){ if(DofData_P->Flag_Init[5]) {
LinAlg_AssembleMatrix(&DofData_P->M5) ; LinAlg_AssembleMatrix(&DofData_P->M5);
LinAlg_AssembleVector(&DofData_P->m5) ; LinAlg_AssembleVector(&DofData_P->m5);
for(int i = 0; i < List_Nbr(DofData_P->m5s); i++) for(int i = 0; i < List_Nbr(DofData_P->m5s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m5s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m5s, i));
} }
if(DofData_P->Flag_Init[6]){ if(DofData_P->Flag_Init[6]) {
LinAlg_AssembleMatrix(&DofData_P->M6) ; LinAlg_AssembleMatrix(&DofData_P->M6);
LinAlg_AssembleVector(&DofData_P->m6) ; LinAlg_AssembleVector(&DofData_P->m6);
for(int i = 0; i < List_Nbr(DofData_P->m6s); i++) for(int i = 0; i < List_Nbr(DofData_P->m6s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m6s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m6s, i));
} }
if(DofData_P->Flag_Init[7]){ if(DofData_P->Flag_Init[7]) {
LinAlg_AssembleMatrix(&DofData_P->M7) ; LinAlg_AssembleMatrix(&DofData_P->M7);
LinAlg_AssembleVector(&DofData_P->m7) ; LinAlg_AssembleVector(&DofData_P->m7);
for(int i = 0; i < List_Nbr(DofData_P->m7s); i++) for(int i = 0; i < List_Nbr(DofData_P->m7s); i++)
LinAlg_AssembleVector((gVector*)List_Pointer(DofData_P->m7s, i)); LinAlg_AssembleVector((gVector *)List_Pointer(DofData_P->m7s, i));
} }
} }
else{ else {
LinAlg_AssembleMatrix(&DofData_P->A) ; LinAlg_AssembleMatrix(&DofData_P->A);
LinAlg_AssembleVector(&DofData_P->b) ; LinAlg_AssembleVector(&DofData_P->b);
if(DofData_P->Flag_ListOfRHS){ if(DofData_P->Flag_ListOfRHS) {
LinAlg_CopyVector(&DofData_P->b , &DofData_P->ListOfRHS[DofData_P->Counter LinAlg_CopyVector(&DofData_P->b,
OfRHS]); &DofData_P->ListOfRHS[DofData_P->CounterOfRHS]);
Message::Info("There are now %d RHS terms", DofData_P->CounterOfRHS + 1); Message::Info("There are now %d RHS terms", DofData_P->CounterOfRHS + 1);
} }
int i; int i;
LinAlg_GetVectorSize(&DofData_P->b, &i) ; LinAlg_GetVectorSize(&DofData_P->b, &i);
if(!i) Message::Info("Generated system is of dimension zero"); if(!i) Message::Info("Generated system is of dimension zero");
if(DofData_P->Flag_Only){ if(DofData_P->Flag_Only) {
for(int i = 0 ; i < List_Nbr( DofData_P->OnlyTheseMatrices ); i++){ for(int i = 0; i < List_Nbr(DofData_P->OnlyTheseMatrices); i++) {
List_Read(DofData_P->OnlyTheseMatrices, i, &iMat); List_Read(DofData_P->OnlyTheseMatrices, i, &iMat);
switch(iMat){ switch(iMat) {
case 1 : case 1:
LinAlg_AssembleMatrix(&Current.DofData->A1) ; LinAlg_AssembleMatrix(&Current.DofData->A1);
LinAlg_AssembleVector(&Current.DofData->b1) ; LinAlg_AssembleVector(&Current.DofData->b1);
break; break;
case 2 : case 2:
LinAlg_AssembleMatrix(&Current.DofData->A2) ; LinAlg_AssembleMatrix(&Current.DofData->A2);
LinAlg_AssembleVector(&Current.DofData->b2) ; LinAlg_AssembleVector(&Current.DofData->b2);
break; break;
case 3: case 3:
LinAlg_AssembleMatrix(&Current.DofData->A3) ; LinAlg_AssembleMatrix(&Current.DofData->A3);
LinAlg_AssembleVector(&Current.DofData->b3) ; LinAlg_AssembleVector(&Current.DofData->b3);
break; break;
} }
} }
} }
} }
if(Flag_Jac){ /* This should in fact only be done if a JacNL term if(Flag_Jac) {
exists in the formulation... */ // This should in fact only be done if a JacNL term exists in the
LinAlg_AssembleMatrix(&DofData_P->Jac) ; // formulation...
LinAlg_AssembleMatrix(&DofData_P->Jac);
} }
Free_UnusedSolutions(DofData_P); Free_UnusedSolutions(DofData_P);
} }
void ReGenerate_System(struct DefineSystem * DefineSystem_P, void ReGenerate_System(struct DefineSystem *DefineSystem_P,
struct DofData * DofData_P, struct DofData *DofData_P, struct DofData *DofData_P0,
struct DofData * DofData_P0, int Flag_Jac=0) int Flag_Jac = 0)
{ {
int Nbr_Formulation, Index_Formulation ; int Nbr_Formulation, Index_Formulation;
struct Formulation * Formulation_P ; struct Formulation *Formulation_P;
ZeroMatrix(&Current.DofData->A, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->A, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
LinAlg_ZeroVector(&Current.DofData->b) ; LinAlg_ZeroVector(&Current.DofData->b);
if(Flag_Jac) if(Flag_Jac)
ZeroMatrix(&Current.DofData->Jac, &Current.DofData->Solver, ZeroMatrix(&Current.DofData->Jac, &Current.DofData->Solver,
Current.DofData->NbrDof) ; Current.DofData->NbrDof);
Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ; Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex);
for (int i = 0 ; i < Nbr_Formulation ; i++) { for(int i = 0; i < Nbr_Formulation; i++) {
List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ; List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation);
Formulation_P = (struct Formulation*) Formulation_P = (struct Formulation *)List_Pointer(Problem_S.Formulation,
List_Pointer(Problem_S.Formulation, Index_Formulation) ; Index_Formulation);
Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P); Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P);
Treatment_Formulation(Formulation_P) ; Treatment_Formulation(Formulation_P);
} }
LinAlg_AssembleMatrix(&DofData_P->A) ; LinAlg_AssembleMatrix(&DofData_P->A);
LinAlg_AssembleVector(&DofData_P->b) ; LinAlg_AssembleVector(&DofData_P->b);
int i; int i;
LinAlg_GetVectorSize(&DofData_P->b, &i) ; LinAlg_GetVectorSize(&DofData_P->b, &i);
if(!i) Message::Info("ReGenerated system is of dimension zero"); if(!i) Message::Info("ReGenerated system is of dimension zero");
if(Flag_Jac){ /* This should in fact only be done if a JacNL term if(Flag_Jac) {
exists in the formulation... */ // This should in fact only be done if a JacNL term exists in the
LinAlg_AssembleMatrix(&DofData_P->Jac) ; // formulation...
LinAlg_AssembleMatrix(&DofData_P->Jac);
} }
} }
void Generate_Residual(gVector *x, gVector *f) void Generate_Residual(gVector *x, gVector *f)
{ {
struct DefineSystem * DefineSystem_P ; struct DefineSystem *DefineSystem_P;
struct DofData * DofData_P ; struct DofData *DofData_P;
struct DofData * DofData_P0 ; struct DofData *DofData_P0;
if(Message::GetVerbosity() == 10) if(Message::GetVerbosity() == 10)
Message::Info("Generating Residual = b(xn)-A(xn)*xn"); Message::Info("Generating Residual = b(xn)-A(xn)*xn");
DofData_P = Current.DofData ; DofData_P = Current.DofData;
DofData_P0 = Current.DofData_P0; DofData_P0 = Current.DofData_P0;
DefineSystem_P = Current.DefineSystem_P ; DefineSystem_P = Current.DefineSystem_P;
if(!DofData_P->CurrentSolution){ if(!DofData_P->CurrentSolution) {
Message::Error("No current solution available"); Message::Error("No current solution available");
return; return;
} }
// new trial solution // new trial solution
LinAlg_CopyVector(x, &DofData_P->dx); LinAlg_CopyVector(x, &DofData_P->dx);
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &DofData_P->d LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x,
x, &DofData_P->dx, -1.,
-1., &DofData_P->CurrentSolution->x); &DofData_P->CurrentSolution->x);
// calculate residual with new solution // calculate residual with new solution
ReGenerate_System(DefineSystem_P, DofData_P, DofData_P0, 1) ; ReGenerate_System(DefineSystem_P, DofData_P, DofData_P0, 1);
// calculate residual with new solution // calculate residual with new solution
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, &DofDat LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
a_P->res) ; &DofData_P->res);
// res = b(xn)-A(xn)*xn // res = b(xn)-A(xn)*xn
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res) ; LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res);
if(Message::GetVerbosity() == 10){ if(Message::GetVerbosity() == 10) {
Message::Info("dx"); LinAlg_PrintVector(stdout, &DofData_P->dx) ; Message::Info("dx");
Message::Info("A"); LinAlg_PrintMatrix(stdout, &DofData_P->A) ; LinAlg_PrintVector(stdout, &DofData_P->dx);
Message::Info("A");
LinAlg_PrintMatrix(stdout, &DofData_P->A);
} }
*f = DofData_P->res ; *f = DofData_P->res;
LinAlg_AssembleVector(f) ; LinAlg_AssembleVector(f);
} }
void Generate_FullJacobian(gVector *x, gMatrix *Jac) void Generate_FullJacobian(gVector *x, gMatrix *Jac)
{ {
struct DofData * DofData_P ; struct DofData *DofData_P;
Message::Debug("Generating Full Jacobian = A(x) + DofData_P->Jac"); Message::Debug("Generating Full Jacobian = A(x) + DofData_P->Jac");
DofData_P = Current.DofData ; DofData_P = Current.DofData;
if(!DofData_P->CurrentSolution){ if(!DofData_P->CurrentSolution) {
Message::Error("No current solution available"); Message::Error("No current solution available");
return; return;
} }
//LinAlg_CopyVector(x, &DofData_P->dx); // LinAlg_CopyVector(x, &DofData_P->dx);
LinAlg_AddVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, LinAlg_AddVectorVector(
&DofData_P->CurrentSolution->x); // updating solution s &DofData_P->CurrentSolution->x, &DofData_P->dx,
olution &DofData_P->CurrentSolution->x); // updating solution solution
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac) ; LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac);
*Jac = DofData_P->Jac ; *Jac = DofData_P->Jac;
LinAlg_AssembleMatrix(Jac) ; LinAlg_AssembleMatrix(Jac);
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* U p d a t e _ C o n s t r a i n t S y s t e m */ /* U p d a t e _ C o n s t r a i n t S y s t e m */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void UpdateConstraint_System(struct DefineSystem * DefineSystem_P, void UpdateConstraint_System(struct DefineSystem *DefineSystem_P,
struct DofData * DofData_P, struct DofData *DofData_P,
struct DofData * DofData_P0, struct DofData *DofData_P0, int GroupIndex,
int GroupIndex, int Type_Constraint, int Flag_Jac) int Type_Constraint, int Flag_Jac)
{ {
// Update constraints, i.e. new preprocessing of STATUS_CST type // Update constraints, i.e. new preprocessing of STATUS_CST type
int Nbr_Formulation, Index_Formulation, Save_TreatmentStatus ; int Nbr_Formulation, Index_Formulation, Save_TreatmentStatus;
struct Formulation * Formulation_P ; struct Formulation *Formulation_P;
// Incrementing Current.SubTimeStep, so that Generate_Link is re-triggered // Incrementing Current.SubTimeStep, so that Generate_Link is re-triggered
Current.SubTimeStep++; Current.SubTimeStep++;
Save_TreatmentStatus = TreatmentStatus ; Save_TreatmentStatus = TreatmentStatus;
TreatmentStatus = STATUS_CST ; TreatmentStatus = STATUS_CST;
Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ; Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex);
for (int k = 0; k < Nbr_Formulation; k++) { for(int k = 0; k < Nbr_Formulation; k++) {
List_Read(DefineSystem_P->FormulationIndex, k, &Index_Formulation) ; List_Read(DefineSystem_P->FormulationIndex, k, &Index_Formulation);
Formulation_P = (struct Formulation*) Formulation_P = (struct Formulation *)List_Pointer(Problem_S.Formulation,
List_Pointer(Problem_S.Formulation, Index_Formulation) ; Index_Formulation);
Message::Info("UpdateConstraint: Treatment Formulation '%s'", Formulation_P- Message::Info("UpdateConstraint: Treatment Formulation '%s'",
>Name) ; Formulation_P->Name);
Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P) ; Init_DofDataInDefineQuantity(DefineSystem_P, DofData_P0, Formulation_P);
Treatment_Formulation(Formulation_P) ; Treatment_Formulation(Formulation_P);
} }
Dof_InitDofType(DofData_P) ; /* Attention: Init for only one DofData */ Dof_InitDofType(DofData_P); /* Attention: Init for only one DofData */
TreatmentStatus = Save_TreatmentStatus ; TreatmentStatus = Save_TreatmentStatus;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* I n i t _ O p e r a t i o n O n S y s t e m */ /* I n i t _ O p e r a t i o n O n S y s t e m */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Init_OperationOnSystem(const char * Name, void Init_OperationOnSystem(const char *Name, struct Resolution *Resolution_P,
struct Resolution * Resolution_P, struct Operation *Operation_P,
struct Operation * Operation_P , struct DofData *DofData_P0,
struct DofData * DofData_P0, struct GeoData *GeoData_P0,
struct GeoData * GeoData_P0, struct DefineSystem **DefineSystem_P,
struct DefineSystem ** DefineSystem_P, struct DofData **DofData_P,
struct DofData ** DofData_P, struct Resolution *Resolution2_P)
struct Resolution * Resolution2_P)
{ {
*DefineSystem_P = (struct DefineSystem*) *DefineSystem_P = (struct DefineSystem *)List_Pointer(
List_Pointer(Resolution_P->DefineSystem,Operation_P->DefineSystemIndex) ; Resolution_P->DefineSystem, Operation_P->DefineSystemIndex);
Current.DefineSystem_P = *DefineSystem_P ; Current.DefineSystem_P = *DefineSystem_P;
*DofData_P = DofData_P0 + Operation_P->DefineSystemIndex ; *DofData_P = DofData_P0 + Operation_P->DefineSystemIndex;
Dof_SetCurrentDofData(Current.DofData = *DofData_P) ; Dof_SetCurrentDofData(Current.DofData = *DofData_P);
Current.NbrHar = Current.DofData->NbrHar ; Current.NbrHar = Current.DofData->NbrHar;
Geo_SetCurrentGeoData(Current.GeoData = GeoData_P0 + (*DofData_P)->GeoDataInde Geo_SetCurrentGeoData(Current.GeoData =
x) ; GeoData_P0 + (*DofData_P)->GeoDataIndex);
if((*DefineSystem_P)->DestinationSystemName && if((*DefineSystem_P)->DestinationSystemName &&
(*DefineSystem_P)->DestinationSystemIndex == -1){ (*DefineSystem_P)->DestinationSystemIndex == -1) {
int i;
int i ; if(Resolution2_P) { /* pre-resolution */
if(Resolution2_P){ /* pre-resolution */ if((i = List_ISearchSeq(Resolution2_P->DefineSystem,
if ((i = List_ISearchSeq(Resolution2_P->DefineSystem, (*DefineSystem_P)->DestinationSystemName,
(*DefineSystem_P)->DestinationSystemName, fcmp_DefineSystem_Name)) < 0) {
fcmp_DefineSystem_Name)) < 0){ Message::Error("Unknown DestinationSystem (%s) in System (%s)",
Message::Error("Unknown DestinationSystem (%s) in System (%s)", (*DefineSystem_P)->DestinationSystemName,
(*DefineSystem_P)->DestinationSystemName, (*DefineSystem_ (*DefineSystem_P)->Name);
P)->Name) ;
return; return;
} }
(*DefineSystem_P)->DestinationSystemIndex = i ; (*DefineSystem_P)->DestinationSystemIndex = i;
Dof_DefineUnknownDofFromSolveOrInitDof(DofData_P) ; Dof_DefineUnknownDofFromSolveOrInitDof(DofData_P);
} }
else { /* a changer !!! */ else { /* a changer !!! */
if ((i = List_ISearchSeq(Resolution_P->DefineSystem, if((i = List_ISearchSeq(Resolution_P->DefineSystem,
(*DefineSystem_P)->DestinationSystemName, (*DefineSystem_P)->DestinationSystemName,
fcmp_DefineSystem_Name)) < 0){ fcmp_DefineSystem_Name)) < 0) {
Message::Error("Unknown DestinationSystem (%s) in System (%s)", Message::Error("Unknown DestinationSystem (%s) in System (%s)",
(*DefineSystem_P)->DestinationSystemName, (*DefineSystem_ (*DefineSystem_P)->DestinationSystemName,
P)->Name) ; (*DefineSystem_P)->Name);
return; return;
} }
(*DefineSystem_P)->DestinationSystemIndex = i ; (*DefineSystem_P)->DestinationSystemIndex = i;
} }
} }
const char *str = Name ? Name : Get_StringForDefine(Operation_Type, Operation_ const char *str =
P->Type); Name ? Name : Get_StringForDefine(Operation_Type, Operation_P->Type);
Message::Info("%s[%s]", str, (*DefineSystem_P)->Name) ; Message::Info("%s[%s]", str, (*DefineSystem_P)->Name);
Message::ProgressMeter(0, 0, "Processing (%s)", str); Message::ProgressMeter(0, 0, "Processing (%s)", str);
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ O p e r a t i o n */ /* T r e a t m e n t _ O p e r a t i o n */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Treatment_Operation(struct Resolution * Resolution_P, void Treatment_Operation(struct Resolution *Resolution_P, List_T *Operation_L,
List_T * Operation_L, struct DofData *DofData_P0, struct GeoData *GeoData_P0,
struct DofData * DofData_P0, struct Resolution *Resolution2_P,
struct GeoData * GeoData_P0, struct DofData *DofData2_P0)
struct Resolution * Resolution2_P,
struct DofData * DofData2_P0)
{ {
double d, d1, d2, *Scales ; double d, d1, d2, *Scales;
int Nbr_Operation, Nbr_Sol, i_Operation, Num_Iteration ; int Nbr_Operation, Nbr_Sol, i_Operation, Num_Iteration;
int Flag_Jac, Flag_Binary = 0 ; int Flag_Jac, Flag_Binary = 0;
int Save_TypeTime ; int Save_TypeTime;
double Save_Time, Save_DTime ; double Save_Time, Save_DTime;
double Save_Iteration ; double Save_Iteration;
double MeanError, RelFactor_Modified ; double MeanError, RelFactor_Modified;
char ResName[256], ResNum[256] ; char ResName[256], ResNum[256];
char FileName[256]; char FileName[256];
char FileName_exMH[256]; char FileName_exMH[256];
gScalar tmp ; gScalar tmp;
FILE *fp = stdout; FILE *fp = stdout;
struct Operation * Operation_P ; struct Operation *Operation_P;
struct DefineSystem * DefineSystem_P ; struct DefineSystem *DefineSystem_P;
struct DofData * DofData_P, * DofData2_P ; struct DofData *DofData_P, *DofData2_P;
struct Solution * Solution_P, Solution_S ; struct Solution *Solution_P, Solution_S;
struct Dof Dof, * Dof_P ; struct Dof Dof, *Dof_P;
struct Value Value ; struct Value Value;
int N ; int N;
static int RES0 = -1 ; static int RES0 = -1;
/* adaptive relaxation */ /* adaptive relaxation */
gVector x_Save; gVector x_Save;
int NbrSteps_relax; int NbrSteps_relax;
double Norm; double Norm;
double Frelax, Frelax_Opt, Error_Prev, Frelax_Prev = 0., Fratio; double Frelax, Frelax_Opt, Error_Prev, Frelax_Prev = 0., Fratio;
int istep; int istep;
int Nbr_Formulation, Index_Formulation ; int Nbr_Formulation, Index_Formulation;
struct Formulation * Formulation_P ; struct Formulation *Formulation_P;
int iTime ; int iTime;
double *Val_Pulsation ; double *Val_Pulsation;
double hop[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC] ; double hop[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC];
double DCfactor ; double DCfactor;
int NbrHar1, NbrHar2, NbrDof1, NbrDof2 ; int NbrHar1, NbrHar2, NbrDof1, NbrDof2;
double dd ; double dd;
int NumDof, iMat ; int NumDof, iMat;
int row_old, row_new, col_old, col_new; int row_old, row_new, col_old, col_new;
double aii, ajj; double aii, ajj;
int nnz__; int nnz__;
List_T *DofList_MH_moving; List_T *DofList_MH_moving;
static int NbrDof_MH_moving; static int NbrDof_MH_moving;
static int *NumDof_MH_moving; static int *NumDof_MH_moving;
static struct Dof ** Dof_MH_moving; static struct Dof **Dof_MH_moving;
gMatrix A_MH_moving_tmp ; gMatrix A_MH_moving_tmp;
//gVector b_MH_moving_tmp ; // gVector b_MH_moving_tmp ;
Nbr_Operation = List_Nbr(Operation_L) ; Nbr_Operation = List_Nbr(Operation_L);
for (i_Operation = 0 ; i_Operation < Nbr_Operation ; i_Operation++) { for(i_Operation = 0; i_Operation < Nbr_Operation; i_Operation++) {
Operation_P = (struct Operation*)List_Pointer(Operation_L, i_Operation) ; Operation_P = (struct Operation *)List_Pointer(Operation_L, i_Operation);
Flag_Jac = 0 ; Flag_Jac = 0;
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) break; if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) break;
switch (Operation_P->Type) { switch(Operation_P->Type) {
/* --> S y s t e m C o m m a n d */ /* --> S y s t e m C o m m a n d */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SYSTEMCOMMAND : case OPERATION_SYSTEMCOMMAND:
BlockingSystemCall(Operation_P->Case.SystemCommand.String); BlockingSystemCall(Operation_P->Case.SystemCommand.String);
break ; break;
/* --> E r r o r */ /* --> E r r o r */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_ERROR : case OPERATION_ERROR:
Message::Error(Operation_P->Case.Error.String); Message::Error(Operation_P->Case.Error.String);
break ; break;
/* --> G e n e r a t e L i s t O f R H S */ /* --> G e n e r a t e L i s t O f R H S */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GENERATELISTOFRHS : case OPERATION_GENERATELISTOFRHS: {
{ Init_OperationOnSystem(
Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P-> Get_StringForDefine(Operation_Type, Operation_P->Type), Resolution_P,
Type), Operation_P, DofData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 Resolution2_P);
, Current.TypeAssembly = ASSEMBLY_AGGREGATE;
&DefineSystem_P, &DofData_P, Resolution2_P) ; DofData_P->TotalNumberOfRHS = Operation_P->Case.Generate.NumListOfRHS;
Current.TypeAssembly = ASSEMBLY_AGGREGATE; Init_SystemData(DofData_P, 0);
DofData_P->TotalNumberOfRHS = Operation_P->Case.Generate.NumListOfRHS; DofData_P->Flag_ListOfRHS = 1;
Init_SystemData(DofData_P, 0) ; // if(Operation_P->Case.Generate.GroupIndex >= 0){
DofData_P->Flag_ListOfRHS = 1; // Generate_Group = (struct Group *)
// if(Operation_P->Case.Generate.GroupIndex >= 0){ // List_Pointer(Problem_S.Group,
// Generate_Group = (struct Group *) // Operation_P->Case.Generate.GroupIndex) ;
// List_Pointer(Problem_S.Group, // }
// Operation_P->Case.Generate.GroupIndex) ; Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0, 0);
// } DofData_P->CounterOfRHS += 1;
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0, 0) ; } break;
DofData_P->CounterOfRHS += 1;
}
break ;
/* --> G e n e r a t e */ /* --> G e n e r a t e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GENERATEJAC : Flag_Jac = 1 ; case OPERATION_GENERATEJAC: Flag_Jac = 1;
case OPERATION_GENERATEJAC_CUMULATIVE : Flag_Jac = 1 ; case OPERATION_GENERATEJAC_CUMULATIVE: Flag_Jac = 1;
case OPERATION_GENERATERHS : case OPERATION_GENERATERHS:
case OPERATION_GENERATERHS_CUMULATIVE : case OPERATION_GENERATERHS_CUMULATIVE:
case OPERATION_GENERATE : case OPERATION_GENERATE:
case OPERATION_GENERATE_CUMULATIVE : case OPERATION_GENERATE_CUMULATIVE: {
{ int cumulative = (Operation_P->Type == OPERATION_GENERATEJAC_CUMULATIVE ||
#ifdef TIMER Operation_P->Type == OPERATION_GENERATERHS_CUMULATIVE ||
double tstart = MPI_Wtime(); Operation_P->Type == OPERATION_GENERATE_CUMULATIVE);
#endif
int cumulative = (Operation_P->Type == OPERATION_GENERATEJAC_CUMULATIVE
||
Operation_P->Type == OPERATION_GENERATERHS_CUMULATIVE
||
Operation_P->Type == OPERATION_GENERATE_CUMULATIVE);
Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->
Type),
Resolution_P, Operation_P, DofData_P0, GeoData_P0
,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
if(Operation_P->Type == OPERATION_GENERATERHS) DofData_P->Flag_RHS = 1; Init_OperationOnSystem(
Get_StringForDefine(Operation_Type, Operation_P->Type), Resolution_P,
Operation_P, DofData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
Resolution2_P);
Current.TypeAssembly = ASSEMBLY_AGGREGATE ; if(Operation_P->Type == OPERATION_GENERATERHS) DofData_P->Flag_RHS = 1;
Init_SystemData(DofData_P, Flag_Jac) ; Current.TypeAssembly = ASSEMBLY_AGGREGATE;
if(Operation_P->Case.Generate.GroupIndex >= 0){
Generate_Group = (struct Group *)
List_Pointer(Problem_S.Group,
Operation_P->Case.Generate.GroupIndex) ;
}
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0, Init_SystemData(DofData_P, Flag_Jac);
cumulative) ; if(Operation_P->Case.Generate.GroupIndex >= 0) {
Generate_Group = (struct Group *)List_Pointer(
Problem_S.Group, Operation_P->Case.Generate.GroupIndex);
}
if(Flag_Jac && !DofData_P->Flag_Only){ Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0,
if(Flag_AddMHMoving){ cumulative);
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving, &DofD
ata_P->A) ;
}
// compute full Jacobian J = A + JacNL, and store it in Jac
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac
) ;
// res = b(xn)-A(xn)*xn if(Flag_Jac && !DofData_P->Flag_Only) {
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, if(Flag_AddMHMoving) {
&DofData_P->res) ; LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving,
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res &DofData_P->A);
) ;
LinAlg_DummyVector(&DofData_P->res) ;
} }
// compute full Jacobian J = A + JacNL, and store it in Jac
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac);
if(Operation_P->Case.Generate.GroupIndex >= 0) // res = b(xn)-A(xn)*xn
Generate_Group = NULL ; LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
&DofData_P->res);
DofData_P->Flag_RHS = 0; LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res);
#ifdef TIMER LinAlg_DummyVector(&DofData_P->res);
double timer = MPI_Wtime() - tstart;
if(Operation_P->Type == OPERATION_GENERATERHS)
printf("Proc %d, time spent in Generate_RHS %.16g\n", Message::GetComm
Rank(), timer);
else
printf("Proc %d, time spent in Generate %.16g\n", Message::GetCommRank
(), timer);
#endif
} }
break ;
if(Operation_P->Case.Generate.GroupIndex >= 0) Generate_Group = NULL;
DofData_P->Flag_RHS = 0;
} break;
/* --> G e n e r a t e S e p a r a t e */ /* --> G e n e r a t e S e p a r a t e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GENERATESEPARATE : case OPERATION_GENERATESEPARATE:
Init_OperationOnSystem("GenerateSeparate", Init_OperationOnSystem("GenerateSeparate", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if (Operation_P->Case.Generate.GroupIndex >= 0) if(Operation_P->Case.Generate.GroupIndex >= 0)
Generate_Group = (struct Group *) List_Pointer(Problem_S.Group, Generate_Group = (struct Group *)List_Pointer(
Operation_P->Case.Generate Problem_S.Group, Operation_P->Case.Generate.GroupIndex);
.GroupIndex) ; Current.TypeAssembly = ASSEMBLY_SEPARATE;
Current.TypeAssembly = ASSEMBLY_SEPARATE ; Init_Update = 0; /* modif... ! */
Init_Update = 0 ; /* modif... ! */
Init_SystemData(DofData_P, Flag_Jac) ; Init_SystemData(DofData_P, Flag_Jac);
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 1) ; Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 1);
if (Operation_P->Case.Generate.GroupIndex >= 0) Generate_Group = NULL ; if(Operation_P->Case.Generate.GroupIndex >= 0) Generate_Group = NULL;
break ; break;
/* --> G e n e r a t e O n l y */ /* --> G e n e r a t e O n l y */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GENERATEONLYJAC : Flag_Jac = 1 ; case OPERATION_GENERATEONLYJAC: Flag_Jac = 1;
case OPERATION_GENERATEONLY: case OPERATION_GENERATEONLY:
Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->Ty Init_OperationOnSystem(
pe), Get_StringForDefine(Operation_Type, Operation_P->Type), Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution2_P);
if(DofData_P->Flag_Only < 2) DofData_P->Flag_Only += 1 ; if(DofData_P->Flag_Only < 2) DofData_P->Flag_Only += 1;
DofData_P->OnlyTheseMatrices = Operation_P->Case.GenerateOnly.MatrixIndex_ DofData_P->OnlyTheseMatrices =
L ; Operation_P->Case.GenerateOnly.MatrixIndex_L;
if (DofData_P->Flag_Only <= 2) if(DofData_P->Flag_Only <= 2)
for (int i = 0; i < List_Nbr(DofData_P->OnlyTheseMatrices); i++){ for(int i = 0; i < List_Nbr(DofData_P->OnlyTheseMatrices); i++) {
List_Read( DofData_P->OnlyTheseMatrices, i, &iMat); List_Read(DofData_P->OnlyTheseMatrices, i, &iMat);
switch(iMat){ switch(iMat) {
case 1: DofData_P->Flag_InitOnly[0] = 1 ; break ; case 1: DofData_P->Flag_InitOnly[0] = 1; break;
case 2: DofData_P->Flag_InitOnly[1] = 1 ; break ; case 2: DofData_P->Flag_InitOnly[1] = 1; break;
case 3: DofData_P->Flag_InitOnly[2] = 1 ; break ; case 3: DofData_P->Flag_InitOnly[2] = 1; break;
} }
} }
Current.TypeAssembly = ASSEMBLY_AGGREGATE ; Current.TypeAssembly = ASSEMBLY_AGGREGATE;
Init_SystemData(DofData_P, Flag_Jac) ; Init_SystemData(DofData_P, Flag_Jac);
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0) ; Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0);
break; break;
/* --> U p d a t e */ /* --> U p d a t e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_UPDATE : case OPERATION_UPDATE:
Init_OperationOnSystem("Update", Init_OperationOnSystem("Update", Resolution_P, Operation_P, DofData_P0,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution2_P);
Operation_Update(DefineSystem_P, DofData_P, DofData_P0, Operation_Update(DefineSystem_P, DofData_P, DofData_P0,
Operation_P->Case.Update.ExpressionIndex) ; Operation_P->Case.Update.ExpressionIndex);
break ; break;
/* --> U p d a t e C o n s t r a i n t */ /* --> U p d a t e C o n s t r a i n t */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_UPDATECONSTRAINT : case OPERATION_UPDATECONSTRAINT:
Init_OperationOnSystem("UpdateConstraint", Init_OperationOnSystem("UpdateConstraint", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
UpdateConstraint_System(DefineSystem_P, DofData_P, DofData_P0, UpdateConstraint_System(DefineSystem_P, DofData_P, DofData_P0,
Operation_P->Case.UpdateConstraint.GroupIndex, Operation_P->Case.UpdateConstraint.GroupIndex,
Operation_P->Case.UpdateConstraint.Type, Flag_Jac) Operation_P->Case.UpdateConstraint.Type,
; Flag_Jac);
break ; break;
/* --> S e l e c t C o r r e c t i o n */ /* --> S e l e c t C o r r e c t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SELECTCORRECTION : case OPERATION_SELECTCORRECTION:
Init_OperationOnSystem("SelectCorrection", Init_OperationOnSystem("SelectCorrection", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if (!Operation_P->Case.SelectCorrection.Iteration) { if(!Operation_P->Case.SelectCorrection.Iteration) {
/* Full solution to be considered again */ /* Full solution to be considered again */
Message::Info(" Full solution to be considered again"); Message::Info(" Full solution to be considered again");
if (DofData_P->CorrectionSolutions.Flag) { if(DofData_P->CorrectionSolutions.Flag) {
DofData_P->CorrectionSolutions.Flag = 0; DofData_P->CorrectionSolutions.Flag = 0;
DofData_P->Solutions = DofData_P->CorrectionSolutions.Save_FullSolution DofData_P->Solutions =
s ; DofData_P->CorrectionSolutions.Save_FullSolutions;
DofData_P->CurrentSolution DofData_P->CurrentSolution =
= DofData_P->CorrectionSolutions.Save_CurrentFullSolution ; DofData_P->CorrectionSolutions.Save_CurrentFullSolution;
} }
else { else {
Message::Error("SelectCorrection: DofData #%d already selected as a ful Message::Error(
l solution", "SelectCorrection: DofData #%d already selected as a full solution",
DofData_P->Num); DofData_P->Num);
} }
} }
else { else {
/* Last correction to be considered */ /* Last correction to be considered */
if (!DofData_P->CorrectionSolutions.Flag) { if(!DofData_P->CorrectionSolutions.Flag) {
DofData_P->CorrectionSolutions.Flag = 1; DofData_P->CorrectionSolutions.Flag = 1;
DofData_P->CorrectionSolutions.Save_FullSolutions = DofData_P->Solution DofData_P->CorrectionSolutions.Save_FullSolutions =
s ; DofData_P->Solutions;
DofData_P->CorrectionSolutions.Save_CurrentFullSolution DofData_P->CorrectionSolutions.Save_CurrentFullSolution =
= DofData_P->CurrentSolution ; DofData_P->CurrentSolution;
/* last correction solutions */ /* last correction solutions */
int i; int i;
if ((i = List_Nbr(DofData_P->CorrectionSolutions.AllSolutions)-1) >= 0) if((i = List_Nbr(DofData_P->CorrectionSolutions.AllSolutions) - 1) >=
{ 0) {
List_Read(DofData_P->CorrectionSolutions.AllSolutions, i, List_Read(DofData_P->CorrectionSolutions.AllSolutions, i,
&DofData_P->Solutions); &DofData_P->Solutions);
} }
else { else {
DofData_P->CorrectionSolutions.AllSolutions = DofData_P->CorrectionSolutions.AllSolutions =
List_Create(10, 10, sizeof(List_T*)); List_Create(10, 10, sizeof(List_T *));
DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution)) ; DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution));
List_Add(DofData_P->CorrectionSolutions.AllSolutions, &DofData_P->Sol List_Add(DofData_P->CorrectionSolutions.AllSolutions,
utions); &DofData_P->Solutions);
} }
/* last time step correction */ /* last time step correction */
if ((i = List_Nbr(DofData_P->Solutions)-1) >= 0) { if((i = List_Nbr(DofData_P->Solutions) - 1) >= 0) {
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution =
List_Pointer(DofData_P->Solutions, i) ; (struct Solution *)List_Pointer(DofData_P->Solutions, i);
} }
else { else {
DofData_P->CurrentSolution = NULL; DofData_P->CurrentSolution = NULL;
/* CurrentSolution will be defined later */ /* CurrentSolution will be defined later */
} }
} }
else { else {
Message::Error("SelectCorrection: DofData #%d already selected as a cor Message::Error(
rection", "SelectCorrection: DofData #%d already selected as a correction",
DofData_P->Num); DofData_P->Num);
} }
} }
break ; break;
/* --> A d d C o r r e c t i o n */ /* --> A d d C o r r e c t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_ADDCORRECTION : case OPERATION_ADDCORRECTION:
Init_OperationOnSystem("AddCorrection", Init_OperationOnSystem("AddCorrection", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if (DofData_P->CorrectionSolutions.Flag) { if(DofData_P->CorrectionSolutions.Flag) {
if(DofData_P->CorrectionSolutions.Save_CurrentFullSolution->TimeStep !=
if (DofData_P->CorrectionSolutions.Save_CurrentFullSolution->TimeStep DofData_P->CurrentSolution->TimeStep) {
!= Solution_S.TimeStep = (int)Current.TimeStep;
DofData_P->CurrentSolution->TimeStep) { Solution_S.Time = Current.Time;
Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
Solution_S.Time = Current.Time ; Solution_S.SolutionExist = 1;
Solution_S.TimeImag = Current.TimeImag ; LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver,
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; DofData_P->NbrDof);
Solution_S.SolutionExist = 1 ; LinAlg_ZeroVector(&Solution_S.x);
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDo
f) ; List_Add(DofData_P->CorrectionSolutions.Save_FullSolutions,
LinAlg_ZeroVector(&Solution_S.x) ; &Solution_S);
DofData_P->CorrectionSolutions.Save_CurrentFullSolution =
List_Add(DofData_P->CorrectionSolutions.Save_FullSolutions, &Solution_S (struct Solution *)List_Pointer(
) ; DofData_P->CorrectionSolutions.Save_FullSolutions,
DofData_P->CorrectionSolutions.Save_CurrentFullSolution = List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutions) - 1);
(struct Solution*) }
List_Pointer(DofData_P->CorrectionSolutions.Save_FullSolutions,
List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutio Cal_SolutionError(
ns)-1) ; &DofData_P->CurrentSolution->x,
} &DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x, 0,
&MeanError);
Cal_SolutionError // LinAlg_VectorNorm2(&DofData_P->CurrentSolution->x, &MeanError);
(&DofData_P->CurrentSolution->x, Message::Info("Mean error: %.3e (after %d iteration%s)", MeanError,
&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x, (int)Current.Iteration,
0, &MeanError) ; ((int)Current.Iteration == 1) ? "" : "s");
//LinAlg_VectorNorm2(&DofData_P->CurrentSolution->x, &MeanError); if(Message::GetProgressMeterStep() > 0 &&
Message::Info("Mean error: %.3e (after %d iteration%s)", Message::GetProgressMeterStep() < 100)
MeanError, (int)Current.Iteration, ((int)Current.Iteration Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
==1)?"":"s") ; "/Residual",
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep(
) < 100)
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Resi
dual",
std::vector<double>(1, MeanError)); std::vector<double>(1, MeanError));
Current.RelativeDifference += Current.RelativeDifference +=
MeanError * Operation_P->Case.AddCorrection.Alpha ; MeanError * Operation_P->Case.AddCorrection.Alpha;
LinAlg_AddVectorVector
(&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x,
&DofData_P->CurrentSolution->x,
&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x) ;
LinAlg_AddVectorVector(
&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x,
&DofData_P->CurrentSolution->x,
&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x);
} }
else { else {
Message::Error("AddCorrection: DofData #%d is not selected as a correcti Message::Error(
on", "AddCorrection: DofData #%d is not selected as a correction",
DofData_P->Num); DofData_P->Num);
} }
break ; break;
/* --> I n i t C o r r e c t i o n */ /* --> I n i t C o r r e c t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_INITCORRECTION : case OPERATION_INITCORRECTION:
Init_OperationOnSystem("InitCorrection", Init_OperationOnSystem("InitCorrection", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if (DofData_P->CorrectionSolutions.Flag) { if(DofData_P->CorrectionSolutions.Flag) {
Solution_S.TimeStep = (int)Current.TimeStep;
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.Time = Current.Time;
Solution_S.Time = Current.Time ; Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeImag = Current.TimeImag ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_S.SolutionExist = 1;
Solution_S.SolutionExist = 1 ; LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver,
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof) DofData_P->NbrDof);
;
/* The last full solution, if any, initializes the current correction */
/* The last full solution, if any, initializes the current correction */ if(List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutions)) {
if (List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutions)) { LinAlg_CopyVector(
LinAlg_CopyVector &((struct Solution *)List_Pointer(
(&((struct Solution *) DofData_P->CorrectionSolutions.Save_FullSolutions,
List_Pointer List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutions) -
(DofData_P->CorrectionSolutions.Save_FullSolutions, 1))
List_Nbr(DofData_P->CorrectionSolutions.Save_FullSolutions)-1))-> ->x,
x, &Solution_S.x);
&Solution_S.x) ; }
} else {
else { LinAlg_ZeroVector(&Solution_S.x);
LinAlg_ZeroVector(&Solution_S.x) ; }
}
List_Add(DofData_P->Solutions, &Solution_S);
List_Add(DofData_P->Solutions, &Solution_S) ; DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
DofData_P->CurrentSolution = DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
(struct Solution*)
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ;
} }
else { else {
Message::Error("InitCorrection: DofData #%d is not selected as a correct Message::Error(
ion", "InitCorrection: DofData #%d is not selected as a correction",
DofData_P->Num); DofData_P->Num);
} }
break ; break;
/* --> M u l t i p l y S o l u t i o n */ /* --> M u l t i p l y S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_MULTIPLYSOLUTION : case OPERATION_MULTIPLYSOLUTION:
Init_OperationOnSystem("MultiplySolution", Init_OperationOnSystem("MultiplySolution", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x, LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x,
Operation_P->Case.MultiplySolution.Alpha, Operation_P->Case.MultiplySolution.Alpha,
&DofData_P->CurrentSolution->x) ; &DofData_P->CurrentSolution->x);
break ; break;
/* --> A d d O p p o s i t e F u l l S o l u t i o n */ /* --> A d d O p p o s i t e F u l l S o l u t i o n */
/* -------------------------------------------------- */ /* -------------------------------------------------- */
case OPERATION_ADDOPPOSITEFULLSOLUTION : case OPERATION_ADDOPPOSITEFULLSOLUTION:
Init_OperationOnSystem("AddOppositeFullSolution", Init_OperationOnSystem("AddOppositeFullSolution", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DefineSystem_P, &DofData_P, Resolution2_P);
LinAlg_AddVectorProdVectorDouble LinAlg_AddVectorProdVectorDouble(
(&DofData_P->CurrentSolution->x, &DofData_P->CurrentSolution->x,
&DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x, -1., &DofData_P->CorrectionSolutions.Save_CurrentFullSolution->x, -1.,
&DofData_P->CurrentSolution->x) ; &DofData_P->CurrentSolution->x);
break ; break;
/* --> S e t R H S A s S o l u t i o n */ /* --> S e t R H S A s S o l u t i o n */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_SETRHSASSOLUTION : case OPERATION_SETRHSASSOLUTION: {
{ /* Compute : x <- b */
/* Compute : x <- b */ Init_OperationOnSystem("SetRHSAsSolution", Resolution_P, Operation_P,
Init_OperationOnSystem("SetRHSAsSolution", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_CopyVector(&DofData_P->b, &DofData_P->CurrentSolution->x);
if(DofData_P->CurrentSolution) else
LinAlg_CopyVector(&DofData_P->b, &DofData_P->CurrentSolution->x); Message::Error("No current solution available");
else } break;
Message::Error("No current solution available");
}
break ;
/* --> S e t S o l u t i o n A s R H S */ /* --> S e t S o l u t i o n A s R H S */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_SETSOLUTIONASRHS : case OPERATION_SETSOLUTIONASRHS: {
{ /* Compute : b <- x */
/* Compute : b <- x */ Init_OperationOnSystem("SetSolutionAsRHS", Resolution_P, Operation_P,
Init_OperationOnSystem("SetSolutionAsRHS", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->b);
if(DofData_P->CurrentSolution) else
LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->b); Message::Error("No current solution available");
else } break;
Message::Error("No current solution available");
}
break ;
/* --> S e t I n c r e m e n t A s S o l u t i o n */ /* --> S e t I n c r e m e n t A s S o l u t i o n */
/* ---------------------------------------------------- */ /* ---------------------------------------------------- */
case OPERATION_SETINCREMENTASSOLUTION : case OPERATION_SETINCREMENTASSOLUTION: {
{ /* Compute : x <- dx */
/* Compute : x <- dx */ Init_OperationOnSystem("SetIncrementAsSolution", Resolution_P,
Init_OperationOnSystem("SetIncrementAsSolution", Operation_P, DofData_P0, GeoData_P0,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DefineSystem_P, &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_CopyVector(&DofData_P->dx, &DofData_P->CurrentSolution->x);
if(DofData_P->CurrentSolution) else
LinAlg_CopyVector(&DofData_P->dx, &DofData_P->CurrentSolution->x); Message::Error("No current solution available");
else } break;
Message::Error("No current solution available");
}
break ;
/* --> S w a p S o l u t i o n */ /* --> S w a p S o l u t i o n */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_SWAPSOLUTIONANDRHS : case OPERATION_SWAPSOLUTIONANDRHS: {
{ Init_OperationOnSystem("SwapSolutionAndRHS", Resolution_P, Operation_P,
Init_OperationOnSystem("SwapSolutionAndRHS", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_SwapVector(&DofData_P->CurrentSolution->x, &DofData_P->b);
if(DofData_P->CurrentSolution) else
LinAlg_SwapVector(&DofData_P->CurrentSolution->x, &DofData_P->b); Message::Error("No current solution available");
else } break;
Message::Error("No current solution available");
}
break ;
case OPERATION_SWAPSOLUTIONANDRESIDUAL : case OPERATION_SWAPSOLUTIONANDRESIDUAL: {
{ Init_OperationOnSystem("SwapSolutionAndResidual", Resolution_P,
Init_OperationOnSystem("SwapSolutionAndResidual", Operation_P, DofData_P0, GeoData_P0,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DefineSystem_P, &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_SwapVector(&DofData_P->CurrentSolution->x, &DofData_P->res);
if(DofData_P->CurrentSolution) else
LinAlg_SwapVector(&DofData_P->CurrentSolution->x, &DofData_P->res); Message::Error("No current solution available");
else } break;
Message::Error("No current solution available");
}
break ;
/* --> C o p y V e c t o r */ /* --> C o p y V e c t o r */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_COPYSOLUTION : case OPERATION_COPYSOLUTION:
case OPERATION_COPYRHS : case OPERATION_COPYRHS:
case OPERATION_COPYRESIDUAL : case OPERATION_COPYRESIDUAL:
case OPERATION_COPYINCREMENT : case OPERATION_COPYINCREMENT:
Init_OperationOnSystem Init_OperationOnSystem(
((Operation_P->Type == OPERATION_COPYSOLUTION) ? "CopySolution" : (Operation_P->Type == OPERATION_COPYSOLUTION) ? "CopySolution" :
(Operation_P->Type == OPERATION_COPYRHS) ? "CopyRightHandSide" : (Operation_P->Type == OPERATION_COPYRHS) ? "CopyRightHandSide" :
(Operation_P->Type == OPERATION_COPYRESIDUAL) ? "CopyResidual" : (Operation_P->Type == OPERATION_COPYRESIDUAL) ? "CopyResidual" :
"CopyIncrement", Resolution_P, Operation_P, DofData_P0, GeoData_P0, "CopyIncrement",
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution_P, Operation_P, DofData_P0, GeoData_P0, &DefineSystem_P,
&DofData_P, Resolution2_P);
Operation_CopyVector(Operation_P, DofData_P); Operation_CopyVector(Operation_P, DofData_P);
break ; break;
/* --> A d d V e c t o r */ /* --> A d d V e c t o r */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_ADDVECTOR : case OPERATION_ADDVECTOR:
Init_OperationOnSystem Init_OperationOnSystem("AddVector", Resolution_P, Operation_P, DofData_P0,
("AddVector", Resolution_P, Operation_P, DofData_P0, GeoData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution2_P);
Operation_AddVector(Operation_P, DofData_P); Operation_AddVector(Operation_P, DofData_P);
break ; break;
/* --> C o p y D o f s */ /* --> C o p y D o f s */
/* ---------------------------------------- */ /* ---------------------------------------- */
case OPERATION_COPYDOFS : case OPERATION_COPYDOFS:
Init_OperationOnSystem("CopyDegreesOfFreedom", Resolution_P, Init_OperationOnSystem("CopyDegreesOfFreedom", Resolution_P, Operation_P,
Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if(!Operation_P->Case.Copy.useList || !Operation_P->Case.Copy.to){ if(!Operation_P->Case.Copy.useList || !Operation_P->Case.Copy.to) {
Message::Error("Degrees of freedom can only be copied to a list") ; Message::Error("Degrees of freedom can only be copied to a list");
} }
else{ else {
std::vector<double> v; std::vector<double> v;
for(int i = 0; i < List_Nbr(DofData_P->DofList); i++){ for(int i = 0; i < List_Nbr(DofData_P->DofList); i++) {
Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i) ; Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i);
v.push_back(Dof_P->NumType); v.push_back(Dof_P->Entity); v.push_back(Dof_P->NumType);
v.push_back(Dof_P->Harmonic); v.push_back(Dof_P->Type); v.push_back(Dof_P->Entity);
switch (Dof_P->Type) { v.push_back(Dof_P->Harmonic);
case DOF_UNKNOWN : v.push_back(Dof_P->Type);
v.push_back(Dof_P->Case.Unknown.NumDof); switch(Dof_P->Type) {
break; case DOF_UNKNOWN: v.push_back(Dof_P->Case.Unknown.NumDof); break;
case DOF_FIXEDWITHASSOCIATE : case DOF_FIXEDWITHASSOCIATE:
v.push_back(Dof_P->Case.FixedAssociate.NumDof); v.push_back(Dof_P->Case.FixedAssociate.NumDof);
break ; break;
case DOF_FIXED : case DOF_FIXED:
case DOF_FIXED_SOLVE : case DOF_FIXED_SOLVE: v.push_back(0); break;
v.push_back(0); case DOF_UNKNOWN_INIT: v.push_back(Dof_P->Case.Unknown.NumDof); break;
break ; case DOF_LINK:
case DOF_UNKNOWN_INIT : case DOF_LINKCPLX: v.push_back(Dof_P->Case.Link.EntityRef); break;
v.push_back(Dof_P->Case.Unknown.NumDof);
break ;
case DOF_LINK :
case DOF_LINKCPLX :
v.push_back(Dof_P->Case.Link.EntityRef) ;
break ;
} }
} }
GetDPNumbers[Operation_P->Case.Copy.to] = v; GetDPNumbers[Operation_P->Case.Copy.to] = v;
} }
break; break;
/* --> G e t N o r m */ /* --> G e t N o r m */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GETNORMSOLUTION : case OPERATION_GETNORMSOLUTION:
case OPERATION_GETNORMRHS : case OPERATION_GETNORMRHS:
case OPERATION_GETNORMRESIDUAL : case OPERATION_GETNORMRESIDUAL:
case OPERATION_GETNORMINCREMENT : case OPERATION_GETNORMINCREMENT: {
{ Init_OperationOnSystem(
Init_OperationOnSystem (Operation_P->Type == OPERATION_GETNORMSOLUTION) ? "GetNormSolution" :
((Operation_P->Type == OPERATION_GETNORMSOLUTION) ? "GetNormSolution" (Operation_P->Type == OPERATION_GETNORMRHS) ? "GetNormRightHandSide" :
: (Operation_P->Type == OPERATION_GETNORMRESIDUAL) ? "GetNormResidual" :
(Operation_P->Type == OPERATION_GETNORMRHS) ? "GetNormRightHandSide" "GetNormIncrement",
: Resolution_P, Operation_P, DofData_P0, GeoData_P0, &DefineSystem_P,
(Operation_P->Type == OPERATION_GETNORMRESIDUAL) ? "GetNormResidual" &DofData_P, Resolution2_P);
: double norm = 0.;
"GetNormIncrement", Resolution_P, Operation_P, DofData_P0, GeoData_P0 if(DofData_P->CurrentSolution &&
, Operation_P->Type == OPERATION_GETNORMSOLUTION)
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_VectorNorm2(&DofData_P->CurrentSolution->x, &norm);
double norm = 0.; else if(Operation_P->Type == OPERATION_GETNORMRESIDUAL)
if(DofData_P->CurrentSolution && Operation_P->Type == OPERATION_GETNORMS LinAlg_VectorNorm2(&DofData_P->res, &norm);
OLUTION) else if(Operation_P->Type == OPERATION_GETNORMRHS)
LinAlg_VectorNorm2(&DofData_P->CurrentSolution->x, &norm); LinAlg_VectorNorm2(&DofData_P->b, &norm);
else if(Operation_P->Type == OPERATION_GETNORMRESIDUAL) else if(Operation_P->Type == OPERATION_GETNORMINCREMENT)
LinAlg_VectorNorm2(&DofData_P->res, &norm); LinAlg_VectorNorm2(&DofData_P->dx, &norm);
else if(Operation_P->Type == OPERATION_GETNORMRHS) Cal_ZeroValue(&Value);
LinAlg_VectorNorm2(&DofData_P->b, &norm); Value.Type = SCALAR;
else if(Operation_P->Type == OPERATION_GETNORMINCREMENT) Value.Val[0] = norm;
LinAlg_VectorNorm2(&DofData_P->dx, &norm); Cal_StoreInVariable(&Value, Operation_P->Case.GetNorm.VariableName);
Cal_ZeroValue(&Value); } break;
Value.Type = SCALAR;
Value.Val[0] = norm;
Cal_StoreInVariable(&Value, Operation_P->Case.GetNorm.VariableName);
}
break ;
/* --> G e t R e s i d u a l */ /* --> G e t R e s i d u a l */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GETRESIDUAL : case OPERATION_GETRESIDUAL: {
{ /* Compute : res = b - A x and return ||res||_2 */
/* Compute : res = b - A x and return ||res||_2 */ Init_OperationOnSystem("GetResidual", Resolution_P, Operation_P,
Init_OperationOnSystem("GetResidual", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, if(DofData_P->CurrentSolution) {
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
if(DofData_P->CurrentSolution){ &DofData_P->res);
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res);
&DofData_P->res); double residual;
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res LinAlg_VectorNorm2(&DofData_P->res, &residual);
); Cal_ZeroValue(&Value);
double residual; Value.Type = SCALAR;
LinAlg_VectorNorm2(&DofData_P->res, &residual); Value.Val[0] = residual;
Cal_ZeroValue(&Value); Cal_StoreInVariable(&Value, Operation_P->Case.GetNorm.VariableName);
Value.Type = SCALAR; if(Message::GetProgressMeterStep() > 0 &&
Value.Val[0] = residual; Message::GetProgressMeterStep() < 100)
Cal_StoreInVariable(&Value, Operation_P->Case.GetNorm.VariableName); Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterSte "/Residual",
p() < 100) std::vector<double>(1, residual));
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Re
sidual",
std::vector<double>(1, residual));
}
else
Message::Error("No current solution available");
} }
break ; else
Message::Error("No current solution available");
} break;
/* --> A p p l y */ /* --> A p p l y */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_APPLY : case OPERATION_APPLY: {
{ /* Compute : x <- A x */
/* Compute : x <- A x */ Init_OperationOnSystem("Apply", Resolution_P, Operation_P, DofData_P0,
Init_OperationOnSystem("Apply", GeoData_P0, &DefineSystem_P, &DofData_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 Resolution2_P);
, if(DofData_P->CurrentSolution) {
&DefineSystem_P, &DofData_P, Resolution2_P) ; LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
if(DofData_P->CurrentSolution){ &DofData_P->res);
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, LinAlg_CopyVector(&DofData_P->res, &DofData_P->CurrentSolution->x);
&DofData_P->res);
LinAlg_CopyVector(&DofData_P->res, &DofData_P->CurrentSolution->x);
}
else
Message::Error("No current solution available");
} }
break ; else
Message::Error("No current solution available");
} break;
/* --> S e t S o l v e r O p t i o n s */ /* --> S e t S o l v e r O p t i o n s */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETGLOBALSOLVEROPTIONS : case OPERATION_SETGLOBALSOLVEROPTIONS: {
{ Message::Info("SetGlobalSolverOptions[\"%s\"]",
Message::Info("SetGlobalSolverOptions[\"%s\"]", Operation_P->Case.SetGlobalSolverOptions.String);
Operation_P->Case.SetGlobalSolverOptions.String); LinAlg_SetGlobalSolverOptions(
LinAlg_SetGlobalSolverOptions(Operation_P->Case.SetGlobalSolverOptions.S Operation_P->Case.SetGlobalSolverOptions.String);
tring); } break;
}
break ;
/* --> S o l v e */ /* --> S o l v e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SOLVEAGAINWITHOTHER : case OPERATION_SOLVEAGAINWITHOTHER:
case OPERATION_SOLVEAGAIN : case OPERATION_SOLVEAGAIN:
case OPERATION_SOLVE : case OPERATION_SOLVE: {
{ int again = (Operation_P->Type == OPERATION_SOLVEAGAINWITHOTHER) ? 2 :
int again = (Operation_P->Type == OPERATION_SOLVEAGAINWITHOTHER) ? 2 : (Operation_P->Type == OPERATION_SOLVEAGAIN) ? 1 :
(Operation_P->Type == OPERATION_SOLVEAGAIN) ? 1 : 0; 0;
/* Solve : A x = b */
Init_OperationOnSystem((again == 2) ? "SolveAgainWithOther" :
(again == 1) ? "SolveAgain" :
"Solve",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P);
#ifdef TIMER if(!DofData_P->CurrentSolution) {
double tstart = MPI_Wtime(); Message::Error("No current solution available");
#endif break;
/* Solve : A x = b */ }
Init_OperationOnSystem((again == 2) ? "SolveAgainWithOther" :
(again == 1) ? "SolveAgain" : "Solve",
Resolution_P, Operation_P, DofData_P0, GeoData_P0
,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
if(!DofData_P->CurrentSolution){ if(DofData_P->Flag_Only) {
Message::Error("No current solution available"); // FIXME: this should move to a separate operation, so that solve
break; // does just solve...
if(DofData_P->Flag_InitOnly[0]) {
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A);
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b);
} }
if (DofData_P->Flag_Only){ if(DofData_P->Flag_InitOnly[1]) {
// FIXME: this should move to a separate operation, so that solve LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A);
// does just solve... LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b);
if(DofData_P->Flag_InitOnly[0]){ }
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A) if(DofData_P->Flag_InitOnly[2]) {
; LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A);
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b) LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b);
;
}
if(DofData_P->Flag_InitOnly[1]){
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A)
;
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b)
;
}
if(DofData_P->Flag_InitOnly[2]){
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A)
;
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b)
;
}
LinAlg_AssembleMatrix(&DofData_P->A) ;
LinAlg_AssembleVector(&DofData_P->b) ;
} }
LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->dx); //kj+ LinAlg_AssembleMatrix(&DofData_P->A);
++ In prevision to build 'dx' in the following (needed for "IterativeLoopPro") LinAlg_AssembleVector(&DofData_P->b);
}
if(!again){ // In prevision to build 'dx' in the following (needed for
LinAlg_Solve(&DofData_P->A, &DofData_P->b, &DofData_P->Solver, // "IterativeLoopPro")
&DofData_P->CurrentSolution->x, LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->dx);
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag) ;
}
else{
DofData *d = (again == 1) ? DofData_P :
DofData_P0 + Operation_P->Case.SolveAgainWithOther.DefineSystemIndex
;
LinAlg_SolveAgain(&d->A, &DofData_P->b, &d->Solver,
&DofData_P->CurrentSolution->x,
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag) ;
}
LinAlg_SubVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, & if(!again) {
DofData_P->dx) ; //kj+++ In order to build 'dx' (needed for "IterativeLoopPro") LinAlg_Solve(&DofData_P->A, &DofData_P->b, &DofData_P->Solver,
#ifdef TIMER &DofData_P->CurrentSolution->x,
double timer = MPI_Wtime() - tstart; (Operation_P->Flag < 0) ? 0 : Operation_P->Flag);
printf("Proc %d, time spent in %s %.16g\n", again ? "SolveAgain" : "Solv
e",
Message::GetCommRank(), timer);
#endif
} }
break ; else {
DofData *d =
(again == 1) ?
DofData_P :
DofData_P0 +
Operation_P->Case.SolveAgainWithOther.DefineSystemIndex;
LinAlg_SolveAgain(&d->A, &DofData_P->b, &d->Solver,
&DofData_P->CurrentSolution->x,
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag);
}
// In order to build 'dx' (needed for "IterativeLoopPro")
LinAlg_SubVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx,
&DofData_P->dx);
} break;
/* --> S o l v e N L */ /* --> S o l v e N L */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SOLVENL : case OPERATION_SOLVENL:
Init_OperationOnSystem("Using PETSc SNES: SolveNL", Init_OperationOnSystem("Using PETSc SNES: SolveNL", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DefineSystem_P, &DofData_P, Resolution2_P);
Init_SystemData(DofData_P, 1); Init_SystemData(DofData_P, 1);
LinAlg_SolveNL(&DofData_P->A, &DofData_P->b, &DofData_P->Jac, &DofData_P-> LinAlg_SolveNL(&DofData_P->A, &DofData_P->b, &DofData_P->Jac,
res, &DofData_P->res, &DofData_P->Solver, &DofData_P->dx,
&DofData_P->Solver, &DofData_P->dx, (Operation_P->Flag < 0) ? 0 : Operation_P->Flag);
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag) ;
break ;
/*
case OPERATION_SOLVENLTS :
Init_OperationOnSystem("Using PETSc SNES and TS: SolveNL",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
LinAlg_SolveNLTS(&DofData_P->A, &DofData_P->b, &DofData_P->Jac, &DofData_P
->res,
&DofData_P->Solver, &DofData_P->dx,
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag) ;
break; break;
*/
/* --> S o l v e J a c */ /* --> S o l v e J a c */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SOLVEJACAGAIN : case OPERATION_SOLVEJACAGAIN:
case OPERATION_SOLVEJAC : case OPERATION_SOLVEJAC: {
{ /* SolveJac : J(xn) dx = b(xn) - A(xn) xn ; x = xn + dx */
/* SolveJac : J(xn) dx = b(xn) - A(xn) xn ; x = xn + dx */
int again = (Operation_P->Type == OPERATION_SOLVEJACAGAIN) ? 1 : 0; int again = (Operation_P->Type == OPERATION_SOLVEJACAGAIN) ? 1 : 0;
Flag_Jac = 1 ; Flag_Jac = 1;
Init_OperationOnSystem(again ? "SolveJacAgain" : "SolveJac", Init_OperationOnSystem(again ? "SolveJacAgain" : "SolveJac", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 Operation_P, DofData_P0, GeoData_P0,
, &DefineSystem_P, &DofData_P, Resolution2_P);
&DefineSystem_P, &DofData_P, Resolution2_P) ;
if(DofData_P->Flag_Init[0] < 2){ if(DofData_P->Flag_Init[0] < 2) {
Message::Error("Jacobian system not initialized (missing GenerateJac?) Message::Error(
"); "Jacobian system not initialized (missing GenerateJac?)");
break; break;
} }
if(!DofData_P->CurrentSolution){ if(!DofData_P->CurrentSolution) {
Message::Error("No current solution available"); Message::Error("No current solution available");
break; break;
} }
if (DofData_P->Flag_Only){ if(DofData_P->Flag_Only) {
// FIXME: this should move to a separate operation, so that solve // FIXME: this should move to a separate operation, so that solve
// does just solve... // does just solve...
if(DofData_P->Flag_InitOnly[0]){ if(DofData_P->Flag_InitOnly[0]) {
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A) LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A);
; LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b);
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b) }
;
} if(DofData_P->Flag_InitOnly[1]) {
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A);
if(DofData_P->Flag_InitOnly[1]){ LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b);
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A) }
; if(DofData_P->Flag_InitOnly[2]) {
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b) LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A);
; LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b);
} }
if(DofData_P->Flag_InitOnly[2]){ LinAlg_AssembleMatrix(&DofData_P->A);
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A) LinAlg_AssembleVector(&DofData_P->b);
;
LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b) // for normal (without Flag_Only) assemblies, the full Jacobian is
; // computed at the end of GenerateJac, as it should be.
} LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac);
LinAlg_AssembleMatrix(&DofData_P->A) ; LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
LinAlg_AssembleVector(&DofData_P->b) ; &DofData_P->res);
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res);
// for normal (without Flag_Only) assemblies, the full Jacobian is LinAlg_DummyVector(&DofData_P->res);
// computed at the end of GenerateJac, as it should be. }
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->Jac, &DofData_P->Jac
) ; if(!again)
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver,
&DofData_P->res) ; &DofData_P->dx);
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res else
) ; LinAlg_SolveAgain(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver,
LinAlg_DummyVector(&DofData_P->res) ; &DofData_P->dx);
}
if(!again) if(!Flag_IterativeLoopN) {
LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver, &Do Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0,
fData_P->dx) ; &MeanError);
else Current.Residual = MeanError;
LinAlg_SolveAgain(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver
, &DofData_P->dx) ;
if(!Flag_IterativeLoopN){ if(MeanError != MeanError) {
//kj+++ (the following lines were outside the condition (!Flag_Iterati Message::Warning("No valid solution found (NaN or Inf)!");
veLoopN) but
// Current.Residual is only needed when Flag_IterativeLoopN==0)
//....................................................................
.....................
Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0, &
MeanError) ;
// NormType: 1=LINFNORM, 2=L1NORM, 3=MEANL1NORM, 4=L2NORM, 5=MEANL2NOR
M
// Closest behaviour to old function Cal_SolutionError with MEANL2NORM
// Cal_SolutionErrorRatio(&DofData_P->dx, &DofData_P->CurrentSolution-
>x,
// reltol, abstol, MEANL2NORM, &MeanError) ;
//LinAlg_VectorNorm2(&DofData_P->dx, &MeanError);
Current.Residual = MeanError; // NB: Residual computed for classical I
terativeLoop using SolveJac
//....................................................................
.....................
if (MeanError != MeanError){
Message::Warning("No valid solution found (NaN or Inf)!");
}
else{
Message::Info("%3ld Nonlinear Residual norm %14.12e",
(int)Current.Iteration, MeanError);
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterS
tep() < 100)
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/
Residual",
std::vector<double>(1, MeanError));
}
} }
else {
Current.RelativeDifference += MeanError ; Message::Info("%3ld Nonlinear Residual norm %14.12e",
//NB: Current.RelativeDifference is what is used for classical Iterative (int)Current.Iteration, MeanError);
Loop stop criterion if(Message::GetProgressMeterStep() > 0 &&
//NB: Current.RelativeDifference is reset to 0 at the begin of every ite Message::GetProgressMeterStep() < 100)
r in iterloop Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
//NB: if only one SolveJac is done: Current.RelativeDifference = MeanErr "/Residual",
or = Current.Residual std::vector<double>(1, MeanError));
if (!Flag_IterativeLoop) {
LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor, &Dof
Data_P->dx) ;
}
else { // Attention: phase test ... Technique bricolee ... provisoire
if (Current.Iteration == 1. || MeanError < Current.RelativeDifferenceO
ld){
LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor, &D
ofData_P->dx) ;
}
else {
RelFactor_Modified = Current.RelaxationFactor /
(MeanError / Current.RelativeDifferenceOld) ;
Message::Info("RelFactor modified = %g", RelFactor_Modified) ;
LinAlg_ProdVectorDouble(&DofData_P->dx, RelFactor_Modified, &DofData
_P->dx) ;
Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0,
&MeanError) ;
//LinAlg_VectorNorm2(&DofData_P->dx, &MeanError);
Message::Info("Mean error: %.3e", MeanError) ;
}
} }
}
LinAlg_AddVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, Current.RelativeDifference += MeanError;
&DofData_P->CurrentSolution->x) ; // NB: Current.RelativeDifference is what is used for classical
// IterativeLoop stopping criterion, and is reset to 0 at the beginning of
// every iteration; if only one SolveJac is done:
// Current.RelativeDifference = MeanError = Current.Residual
if(!Flag_IterativeLoop) {
LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor,
&DofData_P->dx);
}
else { // Attention: phase test ... Technique bricolee ... provisoire
if(Current.Iteration == 1. ||
MeanError < Current.RelativeDifferenceOld) {
LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor,
&DofData_P->dx);
}
else {
RelFactor_Modified = Current.RelaxationFactor /
(MeanError / Current.RelativeDifferenceOld);
Message::Info("RelFactor modified = %g", RelFactor_Modified);
LinAlg_ProdVectorDouble(&DofData_P->dx, RelFactor_Modified,
&DofData_P->dx);
Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0,
&MeanError);
// LinAlg_VectorNorm2(&DofData_P->dx, &MeanError);
Message::Info("Mean error: %.3e", MeanError);
}
} }
break ;
LinAlg_AddVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx,
&DofData_P->CurrentSolution->x);
} break;
/* --> S o l v e J a c _ A d a p t R e l a x */ /* --> S o l v e J a c _ A d a p t R e l a x */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SOLVEJACADAPTRELAX : case OPERATION_SOLVEJACADAPTRELAX:
/* get increment dx by solving : J(xn) dx = b(xn) - A(xn) xn */ /* get increment dx by solving : J(xn) dx = b(xn) - A(xn) xn */
Flag_Jac = 1 ; Flag_Jac = 1;
Fratio = GOLDENRATIO; Fratio = GOLDENRATIO;
//Fratio = 2; // Fratio = 2;
Init_OperationOnSystem("SolveJacAdaptRelax", Init_OperationOnSystem("SolveJacAdaptRelax", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if(DofData_P->Flag_Init[0] < 2){ if(DofData_P->Flag_Init[0] < 2) {
Message::Error("Jacobian system not initialized (missing GenerateJa Message::Error(
c?)"); "Jacobian system not initialized (missing GenerateJac?)");
break; break;
} }
if(!DofData_P->CurrentSolution){ if(!DofData_P->CurrentSolution) {
Message::Error("No current solution available"); Message::Error("No current solution available");
break; break;
} }
LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver, &DofDat LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver,
a_P->dx) ; &DofData_P->dx);
/* save CurrentSolution */ /* save CurrentSolution */
LinAlg_CreateVector(&x_Save, &DofData_P->Solver, DofData_P->NbrDof) ; LinAlg_CreateVector(&x_Save, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &x_Save); LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &x_Save);
Flag_RHS = 1; Flag_RHS = 1;
/* MHBilinear-terms do not contribute to the RHS and residual, and are thu /* MHBilinear-terms do not contribute to the RHS and residual, and are
s disregarded */ * thus disregarded */
/* init dummy values */ /* init dummy values */
Error_Prev = 1e99 ; Frelax_Opt = 1. ; Error_Prev = 1e99;
//if(Current.Iteration==1) Current.Residual_Iter1=1.0; //kj+++ to compute Frelax_Opt = 1.;
a relative residual (relative to residual at iter 1) in SolveJacAdapt (not nece
ssary)
if (!(NbrSteps_relax = List_Nbr(Operation_P->Case.SolveJac_AdaptRelax.Fact if(!(NbrSteps_relax =
or_L))){ List_Nbr(Operation_P->Case.SolveJac_AdaptRelax.Factor_L))) {
Message::Error("No factors provided for Adaptive Relaxation"); Message::Error("No factors provided for Adaptive Relaxation");
break; break;
} }
/*======================================================================== /* CheckAll Meaning:
======= 0 : try first relaxation factors (from the list) and stops when the
CheckAll Meaning: residual goes up
0 : try first relaxation factors (from the list) and stops when the re 1 : try every relaxation factors (from the list) and keep the optimal
sidual goes up one 2 : find the maximum relaxation factor that decreases the residual:
1 : try every relaxation factors (from the list) and keep the optimal - start with the relaxation factor from the previous time step or
one from the previous iteration
2 : find the maximum relaxation factor that decreases the residual: - the relaxation factor is multiplied by a ratio as long as the
- start with the relaxation factor from the previous time step or residual decreases
from the previous iteration - the relaxation factor is decreased by a ratio until a decreasing
- the relaxation factor is multiplied by a ratio as long as the re residual is found */
sidual decreases
- the relaxation factor is decreased by a ratio until a decreasing if(Operation_P->Case.SolveJac_AdaptRelax.CheckAll == 2) {
residual is found Frelax = 1;
[3 : Build a parabola based on three relaxation factors and deduce a m if(Current.Iteration > 1) { Error_Prev = Current.Residual; }
inimizing relaxation factor (NOT WORKING!)] Frelax_Opt = Frelax;
========================================================================== Frelax_Prev = Frelax;
=======*/ }
/* init Frelax for CheckAll==2: */
if (Operation_P->Case.SolveJac_AdaptRelax.CheckAll==2) {
if (Current.Residual>0) {
// if the residual is not 0, a previous TimeStep or iteration has been
done
// ==> Take the last RelaxFactor used from the previous TimeStep or it
eration
Frelax=1;//Current.RelaxFac; //NEW
/*Message::Warning("Init SolveJacAdapt: Start with RelaxFactor "
"from the last iteration (or TimeStep) = %g",Frelax);
//*/
}
else {
// At first TimeStep, at first iteration ...
// ==> Init a RelaxFac
Frelax=1;
/*Message::Warning("Init SolveJacAdapt: Start with "
"a fixed RelaxFactor = %g",Frelax);//*/
}
if(Current.Iteration>1) {
// if a first iteration has been done ...
// ==> Take the Residual from the previous iteration
Error_Prev=Current.Residual;
}
Frelax_Opt = Frelax; Frelax_Prev = Frelax;
}
/* >>>>>>>>>>>>>for printing Error(Frelax) curve : <<<<<<<<<<<<<<<*/
/*
double *FrelaxAll; FrelaxAll = new double[NbrSteps_relax];
double *ErrorAll; ErrorAll = new double[NbrSteps_relax];
for( istep = 0 ; istep < NbrSteps_relax ; istep++ ){
FrelaxAll[istep]=0.;
ErrorAll[istep]=0.;
}
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//*/
/* >>>>>>>>>>>>>test for CheckAll==3(NOT WORKING!): <<<<<<<<<<<<<<<*/
/*
Frelax=1.;
//double Frelax3[3]={Frelax/Fratio, Frelax, Frelax*Fratio};
double Frelax3[3]={0, 1e-3, Frelax*Fratio};
//double Frelax3[3]={0., 1e-1, Fratio};
//List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, NbrSteps_relax
-1, &Frelax3[1]);
//List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, 0, &Frelax3[2]
);
double Error3[3]={0,0,0};
double para=1., parb=1.;
double xopt;
//if (Operation_P->Case.SolveJac_AdaptRelax.CheckAll==3)
// NbrSteps_relax=4;
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//*/
for( istep = 0 ; istep < NbrSteps_relax ; istep++ ){
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) bre
ak;
/* >>>>>>>>>>>>>test for CheckAll==3(NOT WORKING!): <<<<<<<<<<<<<<<*/
/*
if (Operation_P->Case.SolveJac_AdaptRelax.CheckAll==3) {
if (istep<3)
Frelax=Frelax3[istep];
else if (istep>3)
List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, istep, &Fr
elax);
}
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//*/
/* set Frelax : */ for(istep = 0; istep < NbrSteps_relax; istep++) {
if (Operation_P->Case.SolveJac_AdaptRelax.CheckAll<2 ) if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, istep, &Frel break;
ax);
/* new trial solution = x + Frelax * dx */ /* set Frelax : */
LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x); if(Operation_P->Case.SolveJac_AdaptRelax.CheckAll < 2)
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &DofData List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, istep,
_P->dx, &Frelax);
Frelax, &DofData_P->CurrentSolution->x);
/* LinAlg_PrintVector(stdout, &DofData_P->CurrentSolution->x); */ /* new trial solution = x + Frelax * dx */
LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x);
/* calculate residual with trial solution */ LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x,
&DofData_P->dx, Frelax,
ReGenerate_System(DefineSystem_P, DofData_P, DofData_P0) ; &DofData_P->CurrentSolution->x);
if(Flag_AddMHMoving){// Contribution of the moving band (precalcul
ated) /* calculate residual with trial solution */
// Jac does not change (Flag_Jac = 0, default argument of ReGene ReGenerate_System(DefineSystem_P, DofData_P, DofData_P0);
rate_System) if(Flag_AddMHMoving) { // Contribution of the moving band
LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving, & // (precalculated)
DofData_P->A) ; // Jac does not change (Flag_Jac = 0, default argument of
} // ReGenerate_System)
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving,
&DofData_P->res) ; &DofData_P->A);
LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res) ; }
LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x,
/* check whether norm of residual is smaller than previous ones */ &DofData_P->res);
LinAlg_VectorNorm2(&DofData_P->res, &Norm); LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res);
LinAlg_GetVectorSize(&DofData_P->res, &N);
Norm /= (double)N; /* check whether norm of residual is smaller than previous ones */
//Norm /= Current.Residual_Iter1; //kj+++ to compute a relative residual LinAlg_VectorNorm2(&DofData_P->res, &Norm);
(relative to residual at iter 1) in SolveJacAdapt (not necessary) LinAlg_GetVectorSize(&DofData_P->res, &N);
Norm /= (double)N;
Current.Residual = Norm; Current.Residual = Norm;
Current.NbrTestedFac=istep+1; //+++ Current.NbrTestedFac = istep + 1;
/* if(Norm < Error_Prev) {
if (Norm < Error_Prev ) // if the residual has decreased save the current relaxation factor as
Message::Warning("SolveJacAdapt: relaxation factor = %8f" // optimal
" Residual norm = %10.4e" Error_Prev = Norm;
" (BETTER than %10.4e)",
Frelax, Norm, Error_Prev) ;
else
Message::Warning("SolveJacAdapt: relaxation factor = %8f"
" Residual norm = %10.4e"
" (WORST than %10.4e)",
Frelax, Norm, Error_Prev) ;
//*/
/* >>>>>>>>>>>>>for printing Error(Frelax) curve : <<<<<<<<<<<<<<<*/
/*
ErrorAll[istep]=Norm;
FrelaxAll[istep]=Frelax;
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//*/
if (Norm < Error_Prev ) {
// if the residual has decreased ...............................
// save the current relaxation factor as optimal
Error_Prev = Norm;
Frelax_Opt = Frelax; Frelax_Opt = Frelax;
if ( Operation_P->Case.SolveJac_AdaptRelax.CheckAll==2 ){ if(Operation_P->Case.SolveJac_AdaptRelax.CheckAll == 2) {
/* for CheckAll==2: */ if(Frelax < Frelax_Prev && istep > 0) {
if ( Frelax<Frelax_Prev && istep > 0) // if the factor has been decreased ... and a decreasing residual
// if the factor has been decreased ... // has been found => break
// and a decreasing residual has been found => break
break; break;
// if the factor has been increased ...
// => increase the factor (as long as the residual decreases)
Frelax_Prev=Frelax;
Frelax=Frelax*Fratio;
}
}
else if ( Operation_P->Case.SolveJac_AdaptRelax.CheckAll==2 ) {
/* for CheckAll==2: */
// if the residual has increased ...............................
if ( Frelax>Frelax_Prev && istep > 0)
// if the factor has been increased ...
// but the residual has increased => break
break;
// if the factor has been decreased ...
// => decrease the factor (until a decreasing residual is found)
Frelax_Prev=Frelax;
Frelax=Frelax/Fratio;
}
else if (Operation_P->Case.SolveJac_AdaptRelax.CheckAll==0 && istep > 0
)
/* for CheckAll==0: */
// if the residual has increased ...............................
// => break
break ;
/* >>>>>>>>>>>>>test for CheckAll==3(NOT WORKING!): <<<<<<<<<<<<<<<*/
/*
if ( Operation_P->Case.SolveJac_AdaptRelax.CheckAll==3 && istep<3){
Error3[istep]=Norm;
if (istep==2) {
// find Frelax "Opt" and compute corresponding Error
para=(Error3[2]-Error3[0])/((Frelax3[2]-Frelax3[0])*(Frelax3[2]-Frel
ax3[1]))
-(Error3[1]-Error3[0])/((Frelax3[1]-Frelax3[0])*(Frelax3[2]-Fr
elax3[1]));
parb=(Error3[1]-Error3[0])/(Frelax3[1]-Frelax3[0])
-para*(Frelax3[1]+Frelax3[0]);
Frelax=-parb/(2*para);
xopt=Frelax;
if (para<0 || Frelax<0) {
//Something has to be done...
//Frelax=Frelax_Opt;
//Frelax=0.1;
} }
// if the factor has been increased ... => increase the factor (as
// long as the residual decreases)
Frelax_Prev = Frelax;
Frelax = Frelax * Fratio;
} }
} }
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< else if(Operation_P->Case.SolveJac_AdaptRelax.CheckAll == 2) {
//*/ if(Frelax > Frelax_Prev && istep > 0) {
// if the factor has been increased ... but the residual has
if (istep==NbrSteps_relax-1 && Operation_P->Case.SolveJac_AdaptRelax.Che // increased => break
ckAll!=1) { break;
Message::Warning("SolveJacAdapt: LineSearch failed at TimeStep %g iter }
%g",Current.TimeStep,Current.Iteration) ; // if the factor has been decreased ... => decrease the factor (until
Current.SolveJacAdaptFailed=1; // a decreasing residual is found)
//getchar(); Frelax_Prev = Frelax;
} Frelax = Frelax / Fratio;
} }
//Message::Info(" => optimal relaxation factor = %f", Frelax_Opt) ; else if(Operation_P->Case.SolveJac_AdaptRelax.CheckAll == 0 &&
istep > 0) {
/* >>>>>>>>>>>>>for printing Error(Frelax) curve : <<<<<<<<<<<<<<<*/ // if the residual has increased ... => break
/* break;
if ( Operation_P->Case.SolveJac_AdaptRelax.CheckAll==3){ }
printf("xi=[%.g,%g,%g];\n",Frelax3[0],Frelax3[1],Frelax3[2]); if(istep == NbrSteps_relax - 1 &&
printf("yi=[%.18g,%.18g,%.18g];\n",Error3[0],Error3[1],Error3[2]); Operation_P->Case.SolveJac_AdaptRelax.CheckAll != 1) {
printf("xoptpara=[%g];\n",xopt); Message::Warning(
} "SolveJacAdapt: LineSearch failed at TimeStep %g iter %g",
printf("xtaken=[%g];\n",Frelax_Opt); Current.TimeStep, Current.Iteration);
printf("xx=[%g",FrelaxAll[0]); Current.SolveJacAdaptFailed = 1;
for( istep = 1 ; istep < NbrSteps_relax ; istep++ ){ }
printf(",%g",FrelaxAll[istep]); }
if(istep%10==0) // Message::Info(" => optimal relaxation factor = %f", Frelax_Opt) ;
printf("...\n");
}
printf("];\n");
printf("yy=[%g",ErrorAll[0]);
for( istep = 1 ; istep < NbrSteps_relax ; istep++ ){
printf(",%g",ErrorAll[istep]);
if(istep%10==0)
printf("...\n");
}
printf("];\n");
//if (xopt<0 || xopt>2)
// getchar();
delete [] ErrorAll;
delete [] FrelaxAll;
//>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//*/
/* solution = x + Frelax_Opt * dx */ /* solution = x + Frelax_Opt * dx */
LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x); LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x);
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &DofData_ LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x,
P->dx, &DofData_P->dx, Frelax_Opt,
Frelax_Opt, &DofData_P->CurrentSolution->x &DofData_P->CurrentSolution->x);
);
MeanError = Error_Prev;
MeanError = Error_Prev ;
Current.RelaxFac = Frelax_Opt;
Current.RelaxFac = Frelax_Opt; //kj+++
// Residual computed here with SolveJacAdapt (useful to test stop
Current.Residual = MeanError; //kj+++ Residual computed here with SolveJac // criterion
Adapt (usefull to test stop criterion in classical IterativeLoop then) // in classical IterativeLoop then)
Message::Info("%3ld Nonlinear Residual norm %14.12e (optimal relaxation fa Current.Residual = MeanError;
ctor = %f)", Message::Info(
(int)Current.Iteration, MeanError, Frelax_Opt); "%3ld Nonlinear Residual norm %14.12e (optimal relaxation factor = %f)",
//if(Current.Iteration==1) Current.Residual_Iter1=MeanError; //kj+++ to co (int)Current.Iteration, MeanError, Frelax_Opt);
mpute a relative residual (relative to residual at iter 1) in SolveJacAdapt (not if(Message::GetProgressMeterStep() > 0 &&
necessary) Message::GetProgressMeterStep() < 100)
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep() Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
< 100) "/Residual",
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Residu
al",
std::vector<double>(1, MeanError)); std::vector<double>(1, MeanError));
// NB: Current.RelativeDifference is what is used for classical IterativeL // NB: Current.RelativeDifference is what is used for classical
oop stop criterion // IterativeLoop stop criterion here: Current.RelativeDifference =
// here: Current.RelativeDifference = MeanError = Current.Residual; // MeanError = Current.Residual;
Current.RelativeDifference = MeanError; Current.RelativeDifference = MeanError;
Flag_RHS = 0 ; Flag_RHS = 0;
LinAlg_DestroyVector(&x_Save); LinAlg_DestroyVector(&x_Save);
break ; break;
/* --> EigenSolve */ /* --> EigenSolve */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_EIGENSOLVE : case OPERATION_EIGENSOLVE:
Init_OperationOnSystem("EigenSolve", Init_OperationOnSystem("EigenSolve", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
// TODO : DefineOtherSystemIndex DofData2_P0 =
DofData2_P0 = DofData_P0 + Operation_P->Case.EigenSolve.DefineOtherSystemI DofData_P0 + Operation_P->Case.EigenSolve.DefineOtherSystemIndex;
ndex ;
EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues, EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues,
Operation_P->Case.EigenSolve.Shift_r, Operation_P->Case.EigenSolve.Shift_r,
Operation_P->Case.EigenSolve.Shift_i, Operation_P->Case.EigenSolve.Shift_i,
Operation_P->Case.EigenSolve.FilterExpressionIndex, Operation_P->Case.EigenSolve.FilterExpressionIndex,
Operation_P->Case.EigenSolve.RationalCoefsNum, Operation_P->Case.EigenSolve.RationalCoefsNum,
Operation_P->Case.EigenSolve.RationalCoefsDen, Operation_P->Case.EigenSolve.RationalCoefsDen,
Operation_P->Case.EigenSolve.ApplyResolventRealFreqs, Operation_P->Case.EigenSolve.ApplyResolventRealFreqs,
DofData2_P0); DofData2_P0);
break ; break;
/* --> EigenSolveJac */ /* --> EigenSolveJac */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_EIGENSOLVEJAC : case OPERATION_EIGENSOLVEJAC:
Init_OperationOnSystem("EigenSolveJac", Init_OperationOnSystem("EigenSolveJac", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues, EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues,
Operation_P->Case.EigenSolve.Shift_r, Operation_P->Case.EigenSolve.Shift_r,
Operation_P->Case.EigenSolve.Shift_i, Operation_P->Case.EigenSolve.Shift_i,
Operation_P->Case.EigenSolve.FilterExpressionIndex, Operation_P->Case.EigenSolve.FilterExpressionIndex, NULL, NULL,
NULL, NULL, NULL, NULL); NULL, NULL);
/* Insert intelligent convergence test here :-) */ /* Insert intelligent convergence test here :-) */
Current.RelativeDifference = 1.0 ; Current.RelativeDifference = 1.0;
break ; break;
/* --> H P D D M S o l v e */
/* ------------------------------------------ */
case OPERATION_HPDDMSOLVE :
Init_OperationOnSystem("HPDDMSolve",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
Operation_HPDDMSolve(Operation_P, DofData_P);
break;
/* --> Perturbation */
/* ------------------------------------------ */
case OPERATION_PERTURBATION :
Init_OperationOnSystem("Perturbation",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
/*
Perturbation(DofData_P,
DofData_P0+Operation_P->Case.Perturbation.DefineSystemIndex2,
DofData_P0+Operation_P->Case.Perturbation.DefineSystemIndex3,
Operation_P->Case.Perturbation.Size,
Operation_P->Case.Perturbation.Save,
Operation_P->Case.Perturbation.Shift,
Operation_P->Case.Perturbation.PertFreq) ;
*/
break ;
/* --> S e t C u r r e n t S y s t e m */ /* --> S e t C u r r e n t S y s t e m */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETCURRENTSYSTEM : case OPERATION_SETCURRENTSYSTEM:
Init_OperationOnSystem("SetCurrentSystem", Init_OperationOnSystem("SetCurrentSystem", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
break ; break;
/* --> C r e a t e S o l u t i o n */ /* --> C r e a t e S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_CREATESOLUTION : case OPERATION_CREATESOLUTION:
Init_OperationOnSystem("CreateSolution", Init_OperationOnSystem("CreateSolution", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if (!DofData_P->Solutions) if(!DofData_P->Solutions)
DofData_P->Solutions = List_Create( 20, 20, sizeof(struct Solution)) ; DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution));
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.TimeStep = (int)Current.TimeStep;
Solution_S.Time = Current.Time ; Solution_S.Time = Current.Time;
Solution_S.TimeImag = Current.TimeImag ; Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof) LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof);
; LinAlg_ZeroVector(&Solution_S.x);
LinAlg_ZeroVector(&Solution_S.x) ;
{ {
int ts = Operation_P->Case.CreateSolution.CopyFromTimeStep; int ts = Operation_P->Case.CreateSolution.CopyFromTimeStep;
if(ts >= 0){ // FIXME Inno: maybe better to search for the actual if(ts >= 0) { // FIXME Inno: maybe better to search for the actual
// timestep instead of assuming we provide an index // timestep instead of assuming we provide an index
if(ts < List_Nbr(DofData_P->Solutions)){ if(ts < List_Nbr(DofData_P->Solutions)) {
LinAlg_CopyVector(&((struct Solution *) LinAlg_CopyVector(
List_Pointer(DofData_P->Solutions, ts))->x, &((struct Solution *)List_Pointer(DofData_P->Solutions, ts))->x,
&Solution_S.x) ; &Solution_S.x);
} }
else{ else {
Message::Error("Solution at step %d does not exist", ts); Message::Error("Solution at step %d does not exist", ts);
} }
} }
} }
LinAlg_AssembleVector(&Solution_S.x) ; LinAlg_AssembleVector(&Solution_S.x);
List_Add(DofData_P->Solutions, &Solution_S) ; List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ; DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
break ; break;
/* --> I n i t S o l u t i o n */ /* --> I n i t S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_INITSOLUTION : case OPERATION_INITSOLUTION:
case OPERATION_INITSOLUTION1 : case OPERATION_INITSOLUTION1:
Init_OperationOnSystem("InitSolution", Init_OperationOnSystem("InitSolution", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if(Flag_RESTART){ if(Flag_RESTART) {
if (!DofData_P->Solutions){ if(!DofData_P->Solutions) {
Message::Error("No solution to restart the computation"); Message::Error("No solution to restart the computation");
break; break;
} }
for(int i = 0; i < DofData_P->NbrAnyDof; i++){ for(int i = 0; i < DofData_P->NbrAnyDof; i++) {
Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i) ; Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i);
if(Dof_P->Type == DOF_UNKNOWN_INIT) Dof_P->Type = DOF_UNKNOWN ; if(Dof_P->Type == DOF_UNKNOWN_INIT) Dof_P->Type = DOF_UNKNOWN;
} }
for(int i = 0; i < List_Nbr(DofData_P->Solutions); i++){ for(int i = 0; i < List_Nbr(DofData_P->Solutions); i++) {
Solution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, i); Solution_P = (struct Solution *)List_Pointer(DofData_P->Solutions, i);
Free(Solution_P->TimeFunctionValues); Free(Solution_P->TimeFunctionValues);
Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
/* The last solution is the current one */ /* The last solution is the current one */
if(i == List_Nbr(DofData_P->Solutions) - 1) if(i == List_Nbr(DofData_P->Solutions) - 1)
DofData_P->CurrentSolution = Solution_P; DofData_P->CurrentSolution = Solution_P;
} }
RES0 = (int)Current.TimeStep ; RES0 = (int)Current.TimeStep;
} }
else{ else {
if (!DofData_P->Solutions) if(!DofData_P->Solutions)
DofData_P->Solutions = List_Create( 20, 20, sizeof(struct Solution)) ; DofData_P->Solutions = List_Create(20, 20, sizeof(struct Solution));
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.TimeStep = (int)Current.TimeStep;
Solution_S.Time = Current.Time ; Solution_S.Time = Current.Time;
Solution_S.TimeImag = Current.TimeImag ; Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P);
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver,
) ; DofData_P->NbrDof);
/* The last solution, if any, initializes the current one. /* The last solution, if any, initializes the current one. Otherwise a
Otherwise a null solution is used. null solution is used. a revoir qd les conditions initiales
a revoir qd les conditions initiales multiples seront mieux traitees * multiples seront mieux traitees
/ */
if (List_Nbr(DofData_P->Solutions)) { if(List_Nbr(DofData_P->Solutions)) {
LinAlg_CopyVector(&((struct Solution *) LinAlg_CopyVector(
List_Pointer(DofData_P->Solutions, &((struct Solution *)List_Pointer(
List_Nbr(DofData_P->Solutions)-1))->x, DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1))
&Solution_S.x) ; ->x,
} &Solution_S.x);
else { }
LinAlg_ZeroVector(&Solution_S.x) ; else {
} LinAlg_ZeroVector(&Solution_S.x);
}
for(int i = 0; i < DofData_P->NbrAnyDof; i++){
Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i) ; for(int i = 0; i < DofData_P->NbrAnyDof; i++) {
if(Dof_P->Type == DOF_UNKNOWN_INIT){ /* Init values loaded */ Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i);
if(Operation_P->Type == OPERATION_INITSOLUTION){ if(Dof_P->Type == DOF_UNKNOWN_INIT) { /* Init values loaded */
Dof_P->Type = DOF_UNKNOWN ; if(Operation_P->Type == OPERATION_INITSOLUTION) {
LinAlg_SetScalarInVector Dof_P->Type = DOF_UNKNOWN;
(&Dof_P->Val, &Solution_S.x, Dof_P->Case.Unknown.NumDof-1) ; LinAlg_SetScalarInVector(&Dof_P->Val, &Solution_S.x,
Dof_P->Case.Unknown.NumDof - 1);
} }
else{ else {
LinAlg_SetScalarInVector LinAlg_SetScalarInVector(&Dof_P->Val2, &Solution_S.x,
(&Dof_P->Val2, &Solution_S.x, Dof_P->Case.Unknown.NumDof-1) ; Dof_P->Case.Unknown.NumDof - 1);
} }
} }
} }
LinAlg_AssembleVector(&Solution_S.x) ; LinAlg_AssembleVector(&Solution_S.x);
List_Add(DofData_P->Solutions, &Solution_S) ; List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ; DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
} }
break ; break;
/* --> S a v e S o l u t i o n */ /* --> S a v e S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SAVESOLUTION : case OPERATION_SAVESOLUTION:
Init_OperationOnSystem("SaveSolution", Init_OperationOnSystem("SaveSolution", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
strcpy(ResName, Name_Generic) ; strcpy(ResName, Name_Generic);
if(!Flag_SPLIT){ if(!Flag_SPLIT) {
strcat(ResName, ".res") ; strcat(ResName, ".res");
if(RES0 < 0){ if(RES0 < 0) {
Dof_WriteFileRES0(ResName, Flag_BIN) ; Dof_WriteFileRES0(ResName, Flag_BIN);
RES0 = 1 ; RES0 = 1;
} }
} }
else{ else {
strcat(ResName, "-") ; strcat(ResName, "-");
sprintf(ResNum, "%d.res", (int)Current.TimeStep) ; sprintf(ResNum, "%d.res", (int)Current.TimeStep);
for(int i = 0; i < 5+4-(int)strlen(ResNum); i++) strcat(ResName, "0") ; for(int i = 0; i < 5 + 4 - (int)strlen(ResNum); i++)
strcat(ResName, ResNum) ; strcat(ResName, "0");
if(RES0 != (int)Current.TimeStep){ strcat(ResName, ResNum);
Dof_WriteFileRES0(ResName, Flag_BIN) ; if(RES0 != (int)Current.TimeStep) {
RES0 = (int)Current.TimeStep ; Dof_WriteFileRES0(ResName, Flag_BIN);
} RES0 = (int)Current.TimeStep;
} }
Dof_WriteFileRES(ResName, DofData_P, Flag_BIN, Current.Time, Current.TimeI }
mag, Dof_WriteFileRES(ResName, DofData_P, Flag_BIN, Current.Time,
(int)Current.TimeStep) ; Current.TimeImag, (int)Current.TimeStep);
break ; break;
/* --> S a v e S o l u t i o n W i t h E n t i t y N u m */ /* --> S a v e S o l u t i o n W i t h E n t i t y N u m */
/* ------------------------------------------------ */ /* ------------------------------------------------ */
case OPERATION_SAVESOLUTION_WITH_ENTITY_NUM : case OPERATION_SAVESOLUTION_WITH_ENTITY_NUM:
Init_OperationOnSystem("SaveSolutionWithEntityNum", Init_OperationOnSystem("SaveSolutionWithEntityNum", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DefineSystem_P, &DofData_P, Resolution2_P);
strcpy(ResName, Name_Generic) ; strcpy(ResName, Name_Generic);
//strcat(ResName, ".txt") ; // strcat(ResName, ".txt") ;
{ {
int num = Operation_P->Case.SaveSolutionWithEntityNum.GroupIndex; int num = Operation_P->Case.SaveSolutionWithEntityNum.GroupIndex;
Group *g = 0; Group *g = 0;
if (num >= 0) g = (Group*)List_Pointer(Problem_S.Group, num); if(num >= 0) g = (Group *)List_Pointer(Problem_S.Group, num);
bool saveFixed = Operation_P->Case.SaveSolutionWithEntityNum.SaveFixed; bool saveFixed = Operation_P->Case.SaveSolutionWithEntityNum.SaveFixed;
Dof_WriteFileRES_WithEntityNum(ResName, DofData_P, GeoData_P0, g, saveFi Dof_WriteFileRES_WithEntityNum(ResName, DofData_P, GeoData_P0, g,
xed) ; saveFixed);
} }
break ; break;
/* --> S a v e S o l u t i o n s */ /* --> S a v e S o l u t i o n s */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SAVESOLUTIONS : case OPERATION_SAVESOLUTIONS:
Init_OperationOnSystem("SaveSolutions", Init_OperationOnSystem("SaveSolutions", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
strcpy(ResName, Name_Generic) ; strcpy(ResName, Name_Generic);
strcat(ResName, ".res") ; strcat(ResName, ".res");
if(RES0 < 0){ if(RES0 < 0) {
Dof_WriteFileRES0(ResName, Flag_BIN) ; Dof_WriteFileRES0(ResName, Flag_BIN);
RES0 = 1 ; RES0 = 1;
} }
for(int i = 0; i < List_Nbr(DofData_P->Solutions); i++){ for(int i = 0; i < List_Nbr(DofData_P->Solutions); i++) {
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution =
List_Pointer(DofData_P->Solutions, i) ; (struct Solution *)List_Pointer(DofData_P->Solutions, i);
if (!DofData_P->CurrentSolution->SolutionExist) if(!DofData_P->CurrentSolution->SolutionExist)
Message::Warning("Solution #%d doesn't exist anymore: skipping", i) ; Message::Warning("Solution #%d doesn't exist anymore: skipping", i);
else else
Dof_WriteFileRES(ResName, DofData_P, Flag_BIN, Dof_WriteFileRES(ResName, DofData_P, Flag_BIN,
DofData_P->CurrentSolution->Time, DofData_P->CurrentSolution->Time,
DofData_P->CurrentSolution->TimeImag, i) ; DofData_P->CurrentSolution->TimeImag, i);
} }
break ; break;
/* --> M o v i n g B a n d */ /* --> M o v i n g B a n d */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_INIT_MOVINGBAND2D : case OPERATION_INIT_MOVINGBAND2D:
Message::Info("InitMovingBand2D") ; Message::Info("InitMovingBand2D");
Init_MovingBand2D( (struct Group *) Init_MovingBand2D((struct Group *)List_Pointer(
List_Pointer(Problem_S.Group, Problem_S.Group, Operation_P->Case.Init_MovingBand2D.GroupIndex));
Operation_P->Case.Init_MovingBand2D.GroupIn break;
dex)) ;
break ;
case OPERATION_MESH_MOVINGBAND2D : case OPERATION_MESH_MOVINGBAND2D:
if(Message::GetVerbosity() == 10) // +++ if(Message::GetVerbosity() == 10) // +++
Message::Info("MeshMovingBand2D") ; Message::Info("MeshMovingBand2D");
Mesh_MovingBand2D( (struct Group *) Mesh_MovingBand2D((struct Group *)List_Pointer(
List_Pointer(Problem_S.Group, Problem_S.Group, Operation_P->Case.Mesh_MovingBand2D.GroupIndex));
Operation_P->Case.Mesh_MovingBand2D.GroupIn break;
dex)) ;
break ;
case OPERATION_GENERATE_MH_MOVING :
Init_OperationOnSystem("GenerateMHMoving",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
if(gSCALAR_SIZE == 2){ case OPERATION_GENERATE_MH_MOVING:
Message::Error("FIXME: GenerateMHMoving will not work in complex arithmet Init_OperationOnSystem("GenerateMHMoving", Resolution_P, Operation_P,
ic"); DofData_P0, GeoData_P0, &DefineSystem_P,
&DofData_P, Resolution2_P);
if(gSCALAR_SIZE == 2) {
Message::Error(
"FIXME: GenerateMHMoving will not work in complex arithmetic");
break; break;
} }
if (!(Val_Pulsation = Current.DofData->Val_Pulsation)){ if(!(Val_Pulsation = Current.DofData->Val_Pulsation)) {
Message::Error("GenerateMHMoving can only be used for harmonic problems") Message::Error(
; "GenerateMHMoving can only be used for harmonic problems");
break; break;
} }
Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ; Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex);
Generate_Group = (struct Group *) Generate_Group = (struct Group *)List_Pointer(
List_Pointer(Problem_S.Group, Operation_P->Case.Generate_MH_Moving.GroupI Problem_S.Group, Operation_P->Case.Generate_MH_Moving.GroupIndex);
ndex) ;
MH_Moving_Matrix = (double **) Malloc(Current.NbrHar*sizeof(double *)) ; MH_Moving_Matrix = (double **)Malloc(Current.NbrHar * sizeof(double *));
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
MH_Moving_Matrix[k] = (double *) Malloc(Current.NbrHar*sizeof(double)) ; MH_Moving_Matrix[k] = (double *)Malloc(Current.NbrHar * sizeof(double));
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int l = 0; l < Current.NbrHar; l++) for(int l = 0; l < Current.NbrHar; l++) hop[k][l] = 0.;
hop[k][l] = 0.;
Save_Time = Current.Time; Save_Time = Current.Time;
Save_DTime = Current.DTime; Save_DTime = Current.DTime;
MHMoving_assemblyType = 1; // Assembly done in current system: A, b MHMoving_assemblyType = 1; // Assembly done in current system: A, b
for (iTime = 0 ; iTime < Operation_P->Case.Generate_MH_Moving.NbrStep ; iT for(iTime = 0; iTime < Operation_P->Case.Generate_MH_Moving.NbrStep;
ime++) { iTime++) {
Current.Time = (double)iTime /
Current.Time = (double)iTime/(double)Operation_P->Case.Generate_MH_Moving (double)Operation_P->Case.Generate_MH_Moving.NbrStep *
.NbrStep * Operation_P->Case.Generate_MH_Moving.Period;
Operation_P->Case.Generate_MH_Moving.Period ; Current.DTime = 1. /
Current.DTime = 1./(double)Operation_P->Case.Generate_MH_Moving.NbrStep * (double)Operation_P->Case.Generate_MH_Moving.NbrStep *
Operation_P->Case.Generate_MH_Moving.Period ; Operation_P->Case.Generate_MH_Moving.Period;
Current.TimeStep = iTime; Current.TimeStep = iTime;
if(Message::GetVerbosity() == 10) if(Message::GetVerbosity() == 10)
Message::Info("GenerateMHMoving: Step %d/%d (Time = %e DTime %e)", Message::Info("GenerateMHMoving: Step %d/%d (Time = %e DTime %e)",
(int)(Current.TimeStep+1), Operation_P->Case.Generate_MH (int)(Current.TimeStep + 1),
_Moving.NbrStep, Operation_P->Case.Generate_MH_Moving.NbrStep,
Current.Time, Current.DTime) ; Current.Time, Current.DTime);
Treatment_Operation(Resolution_P, Operation_P->Case.Generate_MH_Moving.Op Treatment_Operation(Resolution_P,
eration, Operation_P->Case.Generate_MH_Moving.Operation,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int l = 0; l < Current.NbrHar; l++) { for(int l = 0; l < Current.NbrHar; l++) {
if (Val_Pulsation[k/2]) DCfactor = 2. ; else DCfactor = 1. ; if(Val_Pulsation[k / 2])
MH_Moving_Matrix[k][l] = DCfactor / DCfactor = 2.;
(double)Operation_P->Case.Generate_MH_Moving.NbrStep * else
( fmod(k,2) ? -sin(Val_Pulsation[k/2]*Current.Time) : DCfactor = 1.;
cos(Val_Pulsation[k/2]*Current.Time) ) * MH_Moving_Matrix[k][l] =
( fmod(l,2) ? -sin(Val_Pulsation[l/2]*Current.Time) : DCfactor / (double)Operation_P->Case.Generate_MH_Moving.NbrStep *
cos(Val_Pulsation[l/2]*Current.Time) ) ; (fmod(k, 2) ? -sin(Val_Pulsation[k / 2] * Current.Time) :
hop[k][l] += MH_Moving_Matrix[k][l] ; cos(Val_Pulsation[k / 2] * Current.Time)) *
} (fmod(l, 2) ? -sin(Val_Pulsation[l / 2] * Current.Time) :
cos(Val_Pulsation[l / 2] * Current.Time));
for (int k = 0; k < Current.NbrHar/2; k++) hop[k][l] += MH_Moving_Matrix[k][l];
if (!Val_Pulsation[k]) MH_Moving_Matrix[2*k+1][2*k+1] = 1. ; }
for (int i = 0; i < Nbr_Formulation; i++) { for(int k = 0; k < Current.NbrHar / 2; k++)
List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ; if(!Val_Pulsation[k]) MH_Moving_Matrix[2 * k + 1][2 * k + 1] = 1.;
Formulation_P = (struct Formulation*)
List_Pointer(Problem_S.Formulation, Index_Formulation) ; for(int i = 0; i < Nbr_Formulation; i++) {
Treatment_Formulation(Formulation_P) ; List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation);
} Formulation_P = (struct Formulation *)List_Pointer(
Problem_S.Formulation, Index_Formulation);
Treatment_Formulation(Formulation_P);
}
} }
Current.Time = Save_Time; Current.Time = Save_Time;
Current.DTime = Save_DTime; Current.DTime = Save_DTime;
for (int k = 0; k < Current.NbrHar; k++) Free(MH_Moving_Matrix[k]) ; for(int k = 0; k < Current.NbrHar; k++) Free(MH_Moving_Matrix[k]);
Free(MH_Moving_Matrix) ; Free(MH_Moving_Matrix);
MH_Moving_Matrix = NULL ; MH_Moving_Matrix = NULL;
Generate_Group = NULL; Generate_Group = NULL;
LinAlg_AssembleMatrix(&DofData_P->A) ; LinAlg_AssembleMatrix(&DofData_P->A);
LinAlg_AssembleVector(&DofData_P->b) ; LinAlg_AssembleVector(&DofData_P->b);
LinAlg_AssembleMatrix(&DofData_P->Jac) ; LinAlg_AssembleMatrix(&DofData_P->Jac);
MHMoving_assemblyType = 0; MHMoving_assemblyType = 0;
break ; break;
case OPERATION_GENERATE_MH_MOVING_S : case OPERATION_GENERATE_MH_MOVING_S:
Init_OperationOnSystem("GenerateMHMovingSeparate", Init_OperationOnSystem("GenerateMHMovingSeparate", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DefineSystem_P, &DofData_P, Resolution2_P);
if(gSCALAR_SIZE == 2){ if(gSCALAR_SIZE == 2) {
Message::Error("FIXME: GenerateMHMovingSeparate will not work in complex Message::Error("FIXME: GenerateMHMovingSeparate will not work in "
arithmetic"); "complex arithmetic");
break; break;
} }
if (!(Val_Pulsation = Current.DofData->Val_Pulsation)){ if(!(Val_Pulsation = Current.DofData->Val_Pulsation)) {
Message::Error("GenerateMHMovingSeparate can only be used for harmonic pr Message::Error(
oblems"); "GenerateMHMovingSeparate can only be used for harmonic problems");
break; break;
} }
Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ; Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex);
Generate_Group = (struct Group *) Generate_Group = (struct Group *)List_Pointer(
List_Pointer(Problem_S.Group, Operation_P->Case.Generate_MH_Moving_S.Grou Problem_S.Group, Operation_P->Case.Generate_MH_Moving_S.GroupIndex);
pIndex) ;
MH_Moving_Matrix = (double **) Malloc(Current.NbrHar*sizeof(double *)) ; MH_Moving_Matrix = (double **)Malloc(Current.NbrHar * sizeof(double *));
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
MH_Moving_Matrix[k] = (double *) Malloc(Current.NbrHar*sizeof(double)) ; MH_Moving_Matrix[k] = (double *)Malloc(Current.NbrHar * sizeof(double));
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int l = 0; l < Current.NbrHar; l++) for(int l = 0; l < Current.NbrHar; l++) hop[k][l] = 0.;
hop[k][l] = 0.;
DummyDof = DofData_P->DummyDof ; DummyDof = DofData_P->DummyDof;
DofData_P->DummyDof = NULL ; DofData_P->DummyDof = NULL;
Save_Time = Current.Time; Save_Time = Current.Time;
Save_DTime = Current.DTime; Save_DTime = Current.DTime;
for (iTime = 0 ; iTime < Operation_P->Case.Generate_MH_Moving_S.NbrStep ; for(iTime = 0; iTime < Operation_P->Case.Generate_MH_Moving_S.NbrStep;
iTime++) { iTime++) {
Current.Time = (double)iTime/(double)Operation_P->Case.Generate_MH_Moving Current.Time = (double)iTime /
_S.NbrStep * (double)Operation_P->Case.Generate_MH_Moving_S.NbrStep *
Operation_P->Case.Generate_MH_Moving_S.Period ; Operation_P->Case.Generate_MH_Moving_S.Period;
Current.DTime = 1./(double)Operation_P->Case.Generate_MH_Moving_S.NbrStep Current.DTime = 1. /
* (double)Operation_P->Case.Generate_MH_Moving_S.NbrStep *
Operation_P->Case.Generate_MH_Moving_S.Period ; Operation_P->Case.Generate_MH_Moving_S.Period;
Current.TimeStep = iTime; Current.TimeStep = iTime;
if (!iTime) { if(!iTime) {
//Message::Info("GenerateMHMovingSeparate: probing for any degrees of f // Message::Info("GenerateMHMovingSeparate: probing for any degrees of
reedom"); // freedom");
DofTree_MH_moving = Tree_Create(sizeof(struct Dof), fcmp_Dof) ; DofTree_MH_moving = Tree_Create(sizeof(struct Dof), fcmp_Dof);
// probing assembly // probing assembly
MHMoving_assemblyType = 3; // Constraints - Dofs: Unknown or Link MHMoving_assemblyType = 3; // Constraints - Dofs: Unknown or Link
for (int i = 0; i < Nbr_Formulation; i++) { for(int i = 0; i < Nbr_Formulation; i++) {
List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ; List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation);
Formulation_P = (struct Formulation*) Formulation_P = (struct Formulation *)List_Pointer(
List_Pointer(Problem_S.Formulation, Index_Formulation) ; Problem_S.Formulation, Index_Formulation);
Treatment_Formulation(Formulation_P) ; Treatment_Formulation(Formulation_P);
} }
DofList_MH_moving = Tree2List(DofTree_MH_moving) ; DofList_MH_moving = Tree2List(DofTree_MH_moving);
Tree_Delete(DofTree_MH_moving) ; Tree_Delete(DofTree_MH_moving);
NbrDof_MH_moving = List_Nbr(DofList_MH_moving) ; NbrDof_MH_moving = List_Nbr(DofList_MH_moving);
Message::Info("GenerateMHMovingSeparate: NbrDof_MHMoving = %d", NbrDof_ Message::Info("GenerateMHMovingSeparate: NbrDof_MHMoving = %d",
MH_moving); NbrDof_MH_moving);
Dof_MH_moving = (struct Dof **)Malloc(NbrDof_MH_moving * sizeof(struct Dof_MH_moving =
Dof *)) ; (struct Dof **)Malloc(NbrDof_MH_moving * sizeof(struct Dof *));
NumDof_MH_moving = (int *)Malloc(NbrDof_MH_moving * sizeof(int)) ; NumDof_MH_moving = (int *)Malloc(NbrDof_MH_moving * sizeof(int));
for (int i = 0; i < NbrDof_MH_moving; i++) { for(int i = 0; i < NbrDof_MH_moving; i++) {
Dof_P = (struct Dof*)List_Pointer(DofList_MH_moving,i) ; Dof_P = (struct Dof *)List_Pointer(DofList_MH_moving, i);
if (Dof_P->Type != DOF_UNKNOWN){ if(Dof_P->Type != DOF_UNKNOWN) {
Message::Error("Dof_MH_moving not of type unknown !?"); Message::Error("Dof_MH_moving not of type unknown !?");
break; break;
} }
NumDof_MH_moving[i] = Dof_P->Case.Unknown.NumDof; NumDof_MH_moving[i] = Dof_P->Case.Unknown.NumDof;
if(!(Dof_MH_moving[i] = (struct Dof *)List_PQuery(Current.DofData->Do if(!(Dof_MH_moving[i] = (struct Dof *)List_PQuery(
fList, Current.DofData->DofList, Dof_P, fcmp_Dof))) {
Dof_P, fcmp_Dof))){ Message::Error("GenerateMHMovingSeparate: Dof_MH_moving[%d]=%d "
Message::Error("GenerateMHMovingSeparate: Dof_MH_moving[%d]=%d not "not in Current.DofData->DofList!!!",
in Current.DofData->DofList!!!", i, Dof_MH_moving[i]);
i, Dof_MH_moving[i]) ;
break; break;
} }
for (int k = 0; k < Current.NbrHar; k++) { for(int k = 0; k < Current.NbrHar; k++) {
(Dof_MH_moving[i]+k)->Case.Unknown.NumDof = i*Current.NbrHar+k+1 ; (Dof_MH_moving[i] + k)->Case.Unknown.NumDof =
} i * Current.NbrHar + k + 1;
} /* if (!iTime) */ }
} /* if (!iTime) */
LinAlg_CreateMatrix(&DofData_P->A_MH_moving, &DofData_P->Solver,
NbrDof_MH_moving*Current.NbrHar, LinAlg_CreateMatrix(&DofData_P->A_MH_moving, &DofData_P->Solver,
NbrDof_MH_moving*Current.NbrHar) ; NbrDof_MH_moving * Current.NbrHar,
LinAlg_ZeroMatrix(&DofData_P->A_MH_moving) ; NbrDof_MH_moving * Current.NbrHar);
LinAlg_ZeroMatrix(&DofData_P->A_MH_moving);
/* /*
LinAlg_CreateVector(&DofData_P->b_MH_moving, &DofData_P->Solver, LinAlg_CreateVector(&DofData_P->b_MH_moving, &DofData_P->Solver,
NbrDof_MH_moving*Current.NbrHar) ; NbrDof_MH_moving*Current.NbrHar) ;
LinAlg_ZeroVector(&DofData_P->b_MH_moving) ; LinAlg_ZeroVector(&DofData_P->b_MH_moving) ;
*/ */
} }
if(Message::GetVerbosity() == 10) if(Message::GetVerbosity() == 10)
Message::Info("GenerateMHMovingSeparate : Step %d/%d (Time = %e DTime Message::Info(
%e)", "GenerateMHMovingSeparate : Step %d/%d (Time = %e DTime %e)",
(int)(Current.TimeStep+1), (int)(Current.TimeStep + 1),
Operation_P->Case.Generate_MH_Moving_S.NbrStep, Operation_P->Case.Generate_MH_Moving_S.NbrStep, Current.Time,
Current.Time, Current.DTime) ; Current.DTime);
Treatment_Operation(Resolution_P, Operation_P->Case.Generate_MH_Moving.Op Treatment_Operation(Resolution_P,
eration, Operation_P->Case.Generate_MH_Moving.Operation,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int l = 0; l < Current.NbrHar; l++) { for(int l = 0; l < Current.NbrHar; l++) {
if (Val_Pulsation[k/2]) DCfactor = 2. ; else DCfactor = 1. ; if(Val_Pulsation[k / 2])
MH_Moving_Matrix[k][l] = DCfactor / DCfactor = 2.;
(double)Operation_P->Case.Generate_MH_Moving.NbrStep * else
( fmod(k,2) ? -sin(Val_Pulsation[k/2]*Current.Time) : DCfactor = 1.;
cos(Val_Pulsation[k/2]*Current.Time) ) * MH_Moving_Matrix[k][l] =
( fmod(l,2) ? -sin(Val_Pulsation[l/2]*Current.Time) : DCfactor / (double)Operation_P->Case.Generate_MH_Moving.NbrStep *
cos(Val_Pulsation[l/2]*Current.Time) ) ; (fmod(k, 2) ? -sin(Val_Pulsation[k / 2] * Current.Time) :
hop[k][l] += MH_Moving_Matrix[k][l] ; cos(Val_Pulsation[k / 2] * Current.Time)) *
} (fmod(l, 2) ? -sin(Val_Pulsation[l / 2] * Current.Time) :
cos(Val_Pulsation[l / 2] * Current.Time));
hop[k][l] += MH_Moving_Matrix[k][l];
}
for (int k = 0; k < Current.NbrHar/2; k++) for(int k = 0; k < Current.NbrHar / 2; k++)
if (!Val_Pulsation[k]) MH_Moving_Matrix[2*k+1][2*k+1] = 1. ; if(!Val_Pulsation[k]) MH_Moving_Matrix[2 * k + 1][2 * k + 1] = 1.;
// Assembly in dedicated system: A_MH_Moving, b_MH_moving // Assembly in dedicated system: A_MH_Moving, b_MH_moving
MHMoving_assemblyType = 2; MHMoving_assemblyType = 2;
for (int i = 0; i < Nbr_Formulation; i++) { for(int i = 0; i < Nbr_Formulation; i++) {
List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ; List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation);
Formulation_P = (struct Formulation*) Formulation_P = (struct Formulation *)List_Pointer(
List_Pointer(Problem_S.Formulation, Index_Formulation) ; Problem_S.Formulation, Index_Formulation);
Treatment_Formulation(Formulation_P) ; Treatment_Formulation(Formulation_P);
} }
} /* for iTime */ } /* for iTime */
LinAlg_AssembleMatrix(&DofData_P->A_MH_moving) ; LinAlg_AssembleMatrix(&DofData_P->A_MH_moving);
// LinAlg_AssembleVector(&DofData_P->b_MH_moving) ; // LinAlg_AssembleVector(&DofData_P->b_MH_moving) ;
for (int k = 0; k < Current.NbrHar; k++) Free(MH_Moving_Matrix[k]) ; for(int k = 0; k < Current.NbrHar; k++) Free(MH_Moving_Matrix[k]);
Free(MH_Moving_Matrix) ; Free(MH_Moving_Matrix);
MH_Moving_Matrix = NULL ; MH_Moving_Matrix = NULL;
Generate_Group = NULL; Generate_Group = NULL;
for (int i = 0; i < NbrDof_MH_moving; i++) { for(int i = 0; i < NbrDof_MH_moving; i++) {
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
(Dof_MH_moving[i]+k)->Case.Unknown.NumDof = NumDof_MH_moving[i] + k ; (Dof_MH_moving[i] + k)->Case.Unknown.NumDof = NumDof_MH_moving[i] + k;
} }
LinAlg_CreateMatrix(&A_MH_moving_tmp, &DofData_P->Solver, LinAlg_CreateMatrix(&A_MH_moving_tmp, &DofData_P->Solver,
DofData_P->NbrDof, DofData_P->NbrDof) ; DofData_P->NbrDof, DofData_P->NbrDof);
LinAlg_ZeroMatrix(&A_MH_moving_tmp) ; LinAlg_ZeroMatrix(&A_MH_moving_tmp);
// LinAlg_CreateVector(&b_MH_moving_tmp, &DofData_P->Solver, // LinAlg_CreateVector(&b_MH_moving_tmp, &DofData_P->Solver,
// Current.DofData->NbrDof) ; // Current.DofData->NbrDof) ;
// LinAlg_ZeroVector(&b_MH_moving_tmp) ; // LinAlg_ZeroVector(&b_MH_moving_tmp) ;
nnz__=0; nnz__ = 0;
for (int i = 0; i < NbrDof_MH_moving; i++) { for(int i = 0; i < NbrDof_MH_moving; i++) {
for (int k = 0; k < Current.NbrHar; k++) { for(int k = 0; k < Current.NbrHar; k++) {
row_old = Current.NbrHar*i+k ; row_old = Current.NbrHar * i + k;
row_new = NumDof_MH_moving[i]+k-1 ; row_new = NumDof_MH_moving[i] + k - 1;
// LinAlg_GetDoubleInVector(&d, &DofData_P->b_MH_moving, row_old) ; // LinAlg_GetDoubleInVector(&d, &DofData_P->b_MH_moving, row_old) ;
// LinAlg_SetDoubleInVector( d, &b_MH_moving_tmp, row_new) ; // LinAlg_SetDoubleInVector( d, &b_MH_moving_tmp, row_new) ;
for (int j = 0; j < NbrDof_MH_moving; j++) { for(int j = 0; j < NbrDof_MH_moving; j++) {
for (int l = 0; l < Current.NbrHar; l++) { for(int l = 0; l < Current.NbrHar; l++) {
col_old = Current.NbrHar*j+l ; col_old = Current.NbrHar * j + l;
col_new = NumDof_MH_moving[j]+l-1 ; col_new = NumDof_MH_moving[j] + l - 1;
LinAlg_GetDoubleInMatrix(&d, &DofData_P->A_MH_moving, col_old, r LinAlg_GetDoubleInMatrix(&d, &DofData_P->A_MH_moving, col_old,
ow_old) ; row_old);
LinAlg_GetDoubleInMatrix(&aii, &DofData_P->A_MH_moving, row_old, r LinAlg_GetDoubleInMatrix(&aii, &DofData_P->A_MH_moving, row_old,
ow_old) ; row_old);
LinAlg_GetDoubleInMatrix(&ajj, &DofData_P->A_MH_moving, col_old, c LinAlg_GetDoubleInMatrix(&ajj, &DofData_P->A_MH_moving, col_old,
ol_old) ; col_old);
if(DummyDof==NULL){ if(DummyDof == NULL) {
if(d*d > 1e-12*aii*ajj){ if(d * d > 1e-12 * aii * ajj) {
LinAlg_AddDoubleInMatrix(d, &A_MH_moving_tmp, col_new, row_new LinAlg_AddDoubleInMatrix(d, &A_MH_moving_tmp, col_new,
) ; row_new);
nnz__++; nnz__++;
} }
} }
else{ else {
if(d*d > 1e-12*aii*ajj && if(d * d > 1e-12 * aii * ajj &&
( (DummyDof[row_new]==0 && DummyDof[col_new] == 0) || (row_ne ((DummyDof[row_new] == 0 && DummyDof[col_new] == 0) ||
w == col_new) ) ){ (row_new == col_new))) {
LinAlg_AddDoubleInMatrix(d, &A_MH_moving_tmp, col_new, row_new LinAlg_AddDoubleInMatrix(d, &A_MH_moving_tmp, col_new,
) ; row_new);
nnz__++; nnz__++;
} }
} }
} }
} }
} }
} }
LinAlg_DestroyMatrix(&DofData_P->A_MH_moving); LinAlg_DestroyMatrix(&DofData_P->A_MH_moving);
// LinAlg_DestroyVector(&DofData_P->b_MH_moving); // LinAlg_DestroyVector(&DofData_P->b_MH_moving);
DofData_P->A_MH_moving = A_MH_moving_tmp; DofData_P->A_MH_moving = A_MH_moving_tmp;
// DofData_P->b_MH_moving = b_MH_moving_tmp; // DofData_P->b_MH_moving = b_MH_moving_tmp;
LinAlg_AssembleMatrix(&DofData_P->A_MH_moving); LinAlg_AssembleMatrix(&DofData_P->A_MH_moving);
// LinAlg_AssembleVector(&DofData_P->b_MH_moving); // LinAlg_AssembleVector(&DofData_P->b_MH_moving);
// LinAlg_PrintVector(stdout, &DofData_P->b_MH_moving); // LinAlg_PrintVector(stdout, &DofData_P->b_MH_moving);
Current.Time = Save_Time; Current.Time = Save_Time;
Current.DTime = Save_DTime; Current.DTime = Save_DTime;
Current.TimeStep = 0; // Inner time iteration for integral, no solution in Current.TimeStep =
time 0; // Inner time iteration for integral, no solution in time
DofData_P->DummyDof = DummyDof ; DofData_P->DummyDof = DummyDof;
MHMoving_assemblyType = 0; MHMoving_assemblyType = 0;
Flag_AddMHMoving = 1; Flag_AddMHMoving = 1;
Message::Info("GenerateMHMovingSeparate, contrib. precalculated & assemble Message::Info("GenerateMHMovingSeparate, contrib. precalculated & "
d: Flag_AddMHMoving = %d", Flag_AddMHMoving); "assembled: Flag_AddMHMoving = %d",
Flag_AddMHMoving);
break; break;
case OPERATION_DOFSFREQUENCYSPECTRUM : case OPERATION_DOFSFREQUENCYSPECTRUM:
Dof_GetDummies(DefineSystem_P, DofData_P); Dof_GetDummies(DefineSystem_P, DofData_P);
Message::Info("DofsFrequencySpectrum... DummyDofs"); Message::Info("DofsFrequencySpectrum... DummyDofs");
// FIXME: Name is misleading // FIXME: Name is misleading
// what is taken care of by this function is the Dofs linked to the harmon // what is taken care of by this function is the Dofs linked to the
ics that are not considered // harmonics that are not considered in a particular region (e.g. when
// in a particular region (e.g. when rotor and stator have different spect // rotor and stator have different spectrum) dummydofs ==
rum) // DOFS_NOT_IN_FREQUENCYSPECTRUM_OF_QUANTITY
// dummydofs == DOFS_NOT_IN_FREQUENCYSPECTRUM_OF_QUANTITY break;
break ;
case OPERATION_ADDMHMOVING : Flag_AddMHMoving = 1; case OPERATION_ADDMHMOVING:
Flag_AddMHMoving = 1;
// I think this operation could be merged with GenerateMHMovingSeparate // I think this operation could be merged with GenerateMHMovingSeparate
// LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving, &DofData // LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving,
_P->A) ; // &DofData_P->A) ;
Message::Info("AddMHMoving: contribution of moving band precalculated"); Message::Info("AddMHMoving: contribution of moving band precalculated");
break ; break;
/* --> S a v e S o l u t i o n E x t e n d e d M H */ /* --> S a v e S o l u t i o n E x t e n d e d M H */
/* ----------------------------------------------------------- */ /* ----------------------------------------------------------- */
case OPERATION_SAVESOLUTIONEXTENDEDMH : case OPERATION_SAVESOLUTIONEXTENDEDMH:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
Message::Warning("ExtendSolutionMH can only be used with multi-harmonics" Message::Warning(
) ; "ExtendSolutionMH can only be used with multi-harmonics");
break ; break;
} }
else if (!List_Nbr(DofData_P->Solutions)) { else if(!List_Nbr(DofData_P->Solutions)) {
Message::Warning("No solution available for ExtendSolutionMH"); Message::Warning("No solution available for ExtendSolutionMH");
break ; break;
} }
else if (List_Nbr(DofData_P->Solutions) > 1) { else if(List_Nbr(DofData_P->Solutions) > 1) {
Message::Warning("Only last solution will be extended multi-harmonically Message::Warning(
and saved"); "Only last solution will be extended multi-harmonically and saved");
} }
Init_OperationOnSystem("SaveSolutionExtendedMH", Init_OperationOnSystem("SaveSolutionExtendedMH", Resolution_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DefineSystem_P, &DofData_P, Resolution2_P);
strcpy(FileName_exMH, Name_Generic) ; strcpy(FileName_exMH, Name_Generic);
strcat(FileName_exMH, Operation_P->Case.SaveSolutionExtendedMH.ResFile) ; strcat(FileName_exMH, Operation_P->Case.SaveSolutionExtendedMH.ResFile);
strcat(FileName_exMH, ".res") ; strcat(FileName_exMH, ".res");
Dof_WriteFileRES0(FileName_exMH, Flag_BIN) ; Dof_WriteFileRES0(FileName_exMH, Flag_BIN);
Dof_WriteFileRES_ExtendMH(FileName_exMH, DofData_P, Flag_BIN, Dof_WriteFileRES_ExtendMH(
Current.NbrHar + FileName_exMH, DofData_P, Flag_BIN,
2*Operation_P->Case.SaveSolutionExtendedMH.NbrFre Current.NbrHar + 2 * Operation_P->Case.SaveSolutionExtendedMH.NbrFreq);
q);
Message::Direct(" > '%s' (%d to %d frequencies)", FileName_exMH, Message::Direct(" > '%s' (%d to %d frequencies)", FileName_exMH,
Current.NbrHar/2, Current.NbrHar/2 + Current.NbrHar / 2,
Operation_P->Case.SaveSolutionExtendedMH.NbrFreq) ; Current.NbrHar / 2 +
Operation_P->Case.SaveSolutionExtendedMH.NbrFreq);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
break ; break;
/* --> S a v e S o l u t i o n M H T o T i m e */ /* --> S a v e S o l u t i o n M H T o T i m e */
/* ----------------------------------------------------------- */ /* ----------------------------------------------------------- */
case OPERATION_SAVESOLUTIONMHTOTIME : case OPERATION_SAVESOLUTIONMHTOTIME:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
Message::Warning("SaveSolutionMHtoTime can only to be used with multi-har Message::Warning(
monics") ; "SaveSolutionMHtoTime can only to be used with multi-harmonics");
break ; break;
} }
else if (!List_Nbr(DofData_P->Solutions)) { else if(!List_Nbr(DofData_P->Solutions)) {
Message::Warning("No solution available for SaveSolutionMHtoTime"); Message::Warning("No solution available for SaveSolutionMHtoTime");
break ; break;
} }
else if (List_Nbr(DofData_P->Solutions) > 1) { else if(List_Nbr(DofData_P->Solutions) > 1) {
Message::Warning("Only last mult-harmonic solution will be saved for time Message::Warning(
X"); "Only last mult-harmonic solution will be saved for time X");
} }
Init_OperationOnSystem("SaveSolutionMHtoTime", Init_OperationOnSystem("SaveSolutionMHtoTime", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
strcpy(FileName_exMH, Name_Generic) ; strcpy(FileName_exMH, Name_Generic);
strcat(FileName_exMH, Operation_P->Case.SaveSolutionMHtoTime.ResFile) ; strcat(FileName_exMH, Operation_P->Case.SaveSolutionMHtoTime.ResFile);
strcat(FileName_exMH, ".res") ; strcat(FileName_exMH, ".res");
Dof_WriteFileRES0(FileName_exMH, Flag_BIN) ; Dof_WriteFileRES0(FileName_exMH, Flag_BIN);
Dof_WriteFileRES_MHtoTime(FileName_exMH, DofData_P, Flag_BIN, Dof_WriteFileRES_MHtoTime(FileName_exMH, DofData_P, Flag_BIN,
Operation_P->Case.SaveSolutionMHtoTime.Time) ; Operation_P->Case.SaveSolutionMHtoTime.Time);
Message::Direct(" > '%s' (time = %e)", FileName_exMH, Message::Direct(" > '%s' (time = %e)", FileName_exMH,
Operation_P->Case.SaveSolutionMHtoTime.Time) ; Operation_P->Case.SaveSolutionMHtoTime.Time);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
break ; break;
/* --> R e a d S o l u t i o n */ /* --> R e a d S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_READSOLUTION : case OPERATION_READSOLUTION: {
{ Init_OperationOnSystem("ReadSolution", Resolution_P, Operation_P,
Init_OperationOnSystem("ReadSolution", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, int i = 0;
&DefineSystem_P, &DofData_P, Resolution2_P) ; while(Name_ResFile[i]) {
int i = 0 ; Message::Info("Loading Processing data '%s'", Name_ResFile[i]);
while(Name_ResFile[i]){ Dof_OpenFile(DOF_TMP, Name_ResFile[i], "rb");
Message::Info("Loading Processing data '%s'", Name_ResFile[i]) ; Dof_ReadFileRES(NULL, DofData_P, DofData_P->Num, &Current.Time,
Dof_OpenFile(DOF_TMP, Name_ResFile[i], "rb"); &Current.TimeImag, &Current.TimeStep);
Dof_ReadFileRES(NULL, DofData_P, DofData_P->Num, &Current.Time, &Curre Dof_CloseFile(DOF_TMP);
nt.TimeImag, i++;
&Current.TimeStep) ; }
Dof_CloseFile(DOF_TMP); if(!List_Nbr(DofData_P->Solutions)) {
i++ ; Message::Error("No valid data found for ReadSolution[%s]",
} DefineSystem_P->Name);
if(!List_Nbr(DofData_P->Solutions)){ break;
Message::Error("No valid data found for ReadSolution[%s]", DefineSyste
m_P->Name);
break;
}
DofData_P->CurrentSolution = (struct Solution*)
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ;
Free(DofData_P->CurrentSolution->TimeFunctionValues);
DofData_P->CurrentSolution->TimeFunctionValues = Get_TimeFunctionValues(
DofData_P) ;
} }
break ;
DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
Free(DofData_P->CurrentSolution->TimeFunctionValues);
DofData_P->CurrentSolution->TimeFunctionValues =
Get_TimeFunctionValues(DofData_P);
} break;
/* --> G m s h R e a d */ /* --> G m s h R e a d */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_GMSHREAD : case OPERATION_GMSHREAD:
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
if(Operation_P->Case.GmshRead.RunTimeVar){ if(Operation_P->Case.GmshRead.RunTimeVar) {
// FIXME: well, this is reaaally ugly and unsafe - we should sanitize // FIXME: well, this is reaaally ugly and unsafe - we should sanitize
// the string and verify that it contains a valid format specification : // the string and verify that it contains a valid format specification
-) // :-)
struct Value val; struct Value val;
Cal_GetValueSaved(&val, Operation_P->Case.GmshRead.RunTimeVar); Cal_GetValueSaved(&val, Operation_P->Case.GmshRead.RunTimeVar);
char tmp[256]; char tmp[256];
sprintf(tmp, Operation_P->Case.GmshRead.FileName, val.Val[0]); sprintf(tmp, Operation_P->Case.GmshRead.FileName, val.Val[0]);
Message::Info("GmshRead[%s]", tmp); Message::Info("GmshRead[%s]", tmp);
GmshMergePostProcessingFile(tmp); GmshMergePostProcessingFile(tmp);
} }
else{ else {
if(Operation_P->Case.GmshRead.ViewTag >= 0){ if(Operation_P->Case.GmshRead.ViewTag >= 0) {
PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag); PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag);
Message::Info("GmshRead[%s] -> View[%d]", Operation_P->Case.GmshRead.F Message::Info("GmshRead[%s] -> View[%d]",
ileName, Operation_P->Case.GmshRead.FileName,
Operation_P->Case.GmshRead.ViewTag); Operation_P->Case.GmshRead.ViewTag);
} }
else { else {
Message::Info("GmshRead[%s]", Operation_P->Case.GmshRead.FileName); Message::Info("GmshRead[%s]", Operation_P->Case.GmshRead.FileName);
} }
GmshMergePostProcessingFile(Operation_P->Case.GmshRead.FileName); GmshMergePostProcessingFile(Operation_P->Case.GmshRead.FileName);
} }
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'GmshRe Message::Error(
ad'"); "You need to compile GetDP with Gmsh support to use 'GmshRead'");
#endif #endif
break ; break;
case OPERATION_GMSHMERGE : case OPERATION_GMSHMERGE:
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
if(Operation_P->Case.GmshRead.ViewTag >= 0){ if(Operation_P->Case.GmshRead.ViewTag >= 0) {
PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag); PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag);
Message::Info("GmshMerge[%s] -> View[%d]", Operation_P->Case.GmshRead.Fi Message::Info("GmshMerge[%s] -> View[%d]",
leName, Operation_P->Case.GmshRead.FileName,
Operation_P->Case.GmshRead.ViewTag); Operation_P->Case.GmshRead.ViewTag);
} }
else{ else {
Message::Info("GmshMerge[%s]", Operation_P->Case.GmshRead.FileName); Message::Info("GmshMerge[%s]", Operation_P->Case.GmshRead.FileName);
} }
GmshMergeFile(Operation_P->Case.GmshRead.FileName); GmshMergeFile(Operation_P->Case.GmshRead.FileName);
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'GmshMe Message::Error(
rge'"); "You need to compile GetDP with Gmsh support to use 'GmshMerge'");
#endif #endif
break ; break;
case OPERATION_GMSHOPEN : case OPERATION_GMSHOPEN:
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
if(Operation_P->Case.GmshRead.ViewTag >= 0){ if(Operation_P->Case.GmshRead.ViewTag >= 0) {
PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag); PView::setGlobalTag(Operation_P->Case.GmshRead.ViewTag);
Message::Info("GmshOpen[%s] -> View[%d]", Operation_P->Case.GmshRead.Fil Message::Info("GmshOpen[%s] -> View[%d]",
eName, Operation_P->Case.GmshRead.FileName,
Operation_P->Case.GmshRead.ViewTag); Operation_P->Case.GmshRead.ViewTag);
} }
else{ else {
Message::Info("GmshOpen[%s]", Operation_P->Case.GmshRead.FileName); Message::Info("GmshOpen[%s]", Operation_P->Case.GmshRead.FileName);
} }
GmshOpenProject(Operation_P->Case.GmshRead.FileName); GmshOpenProject(Operation_P->Case.GmshRead.FileName);
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'GmshOp Message::Error(
en'"); "You need to compile GetDP with Gmsh support to use 'GmshOpen'");
#endif #endif
break ; break;
case OPERATION_GMSHCLEARALL : case OPERATION_GMSHCLEARALL:
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
Message::Info("GmshClearAll[]"); Message::Info("GmshClearAll[]");
while(PView::list.size()) delete PView::list[0]; while(PView::list.size()) delete PView::list[0];
PView::setGlobalTag(0); PView::setGlobalTag(0);
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'GmshCl Message::Error(
earAll'"); "You need to compile GetDP with Gmsh support to use 'GmshClearAll'");
#endif #endif
break ; break;
case OPERATION_GMSHWRITE : case OPERATION_GMSHWRITE:
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
{ {
Message::Info("GmshWrite[%s]", Operation_P->Case.GmshRead.FileName); Message::Info("GmshWrite[%s]", Operation_P->Case.GmshRead.FileName);
PView *view = PView::getViewByTag(Operation_P->Case.GmshRead.ViewTag); PView *view = PView::getViewByTag(Operation_P->Case.GmshRead.ViewTag);
if(view) if(view)
view->write(Operation_P->Case.GmshRead.FileName, 10); view->write(Operation_P->Case.GmshRead.FileName, 10);
else else
Message::Error("View %d does not exist"); Message::Error("View %d does not exist");
} }
#else #else
Message::Error("You need to compile GetDP with Gmsh support to use 'GmshWr Message::Error(
ite'"); "You need to compile GetDP with Gmsh support to use 'GmshWrite'");
#endif #endif
break ; break;
/* --> S a v e M e s h */ /* --> S a v e M e s h */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SAVEMESH : case OPERATION_SAVEMESH:
Init_OperationOnSystem("SaveMesh", Init_OperationOnSystem("SaveMesh", Resolution_P, Operation_P, DofData_P0,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution2_P);
// FIXME: wrong on Windows -- see Fix_RelativePath // FIXME: wrong on Windows -- see Fix_RelativePath
if(Operation_P->Case.SaveMesh.FileName[0] == '/' || if(Operation_P->Case.SaveMesh.FileName[0] == '/' ||
Operation_P->Case.SaveMesh.FileName[0] == '\\'){ Operation_P->Case.SaveMesh.FileName[0] == '\\') {
strcpy(FileName, Operation_P->Case.SaveMesh.FileName); strcpy(FileName, Operation_P->Case.SaveMesh.FileName);
} }
else { else {
strcpy(FileName, Name_Path); strcpy(FileName, Name_Path);
strcat(FileName, Operation_P->Case.SaveMesh.FileName); strcat(FileName, Operation_P->Case.SaveMesh.FileName);
} }
if (Operation_P->Case.SaveMesh.ExprIndex >= 0) { if(Operation_P->Case.SaveMesh.ExprIndex >= 0) {
Get_ValueOfExpressionByIndex(Operation_P->Case.SaveMesh.ExprIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.SaveMesh.ExprIndex, NULL,
NULL, 0., 0., 0., &Value) ; 0., 0., 0., &Value);
char fmt[256]; strcpy(fmt, FileName); char fmt[256];
sprintf(FileName, fmt, Value.Val[0]); strcpy(fmt, FileName);
sprintf(FileName, fmt, Value.Val[0]);
} }
Geo_SaveMesh(Current.GeoData, Geo_SaveMesh(Current.GeoData,
((struct Group*) ((struct Group *)List_Pointer(
List_Pointer(Problem_S.Group, Problem_S.Group, Operation_P->Case.SaveMesh.GroupIndex))
Operation_P->Case.SaveMesh.GroupIndex))->Initial ->InitialList,
List, FileName);
FileName) ; break;
break ;
/* --> T r a n s f e r S o l u t i o n */ /* --> T r a n s f e r S o l u t i o n */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_TRANSFERSOLUTION : case OPERATION_TRANSFERSOLUTION:
Init_OperationOnSystem("TransferSolution", Init_OperationOnSystem("TransferSolution", Resolution_P, Operation_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P0, GeoData_P0, &DefineSystem_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; &DofData_P, Resolution2_P);
if(Resolution2_P){ /* pre-resolution */ if(Resolution2_P) { /* pre-resolution */
DofData2_P = DofData2_P0 + DefineSystem_P->DestinationSystemIndex ; DofData2_P = DofData2_P0 + DefineSystem_P->DestinationSystemIndex;
Dof_TransferDof(DofData_P, &DofData2_P); Dof_TransferDof(DofData_P, &DofData2_P);
} }
else{ else {
/* a changer!!! Il faut se mettre d'accord sur ce que doit faire /* a changer!!! Il faut se mettre d'accord sur ce que doit faire
Dof_TransferDof. Ceci sert a transferer la derniere solution d'un Dof_TransferDof. Ceci sert a transferer la derniere solution d'un
DofData dans un autre (ds la meme resolution), base sur le meme DofData dans un autre (ds la meme resolution), base sur le meme
espace fonctionnel. */ espace fonctionnel. */
DofData2_P = DofData_P0 + DefineSystem_P->DestinationSystemIndex ; DofData2_P = DofData_P0 + DefineSystem_P->DestinationSystemIndex;
if(DofData_P->NbrAnyDof != DofData2_P->NbrAnyDof){ if(DofData_P->NbrAnyDof != DofData2_P->NbrAnyDof) {
Message::Error("Dimensions do not match for TransferSolution"); Message::Error("Dimensions do not match for TransferSolution");
break; break;
} }
Solution_S.TimeStep = (int)Current.TimeStep ; Solution_S.TimeStep = (int)Current.TimeStep;
Solution_S.Time = Current.Time ; Solution_S.Time = Current.Time;
Solution_S.TimeImag = Current.TimeImag ; Solution_S.TimeImag = Current.TimeImag;
Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData2_P) ; Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData2_P);
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver, DofData2_P->NbrD LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver,
of) ; DofData2_P->NbrDof);
LinAlg_ZeroVector(&Solution_S.x) ; LinAlg_ZeroVector(&Solution_S.x);
if (List_Nbr(DofData_P->Solutions)) { if(List_Nbr(DofData_P->Solutions)) {
Solution_P = (struct Solution *)List_Pointer(
Solution_P = (struct Solution *)List_Pointer(DofData_P->Solutions, DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
List_Nbr(DofData_P->Soluti for(int i = 0; i < DofData_P->NbrAnyDof; i++) {
ons)-1) ; Dof = *(struct Dof *)List_Pointer(DofData_P->DofList, i);
for(int i = 0; i < DofData_P->NbrAnyDof; i++){ if(Dof.Type == DOF_UNKNOWN) {
Dof = *(struct Dof *)List_Pointer(DofData_P->DofList, i) ; LinAlg_GetScalarInVector(&tmp, &Solution_P->x,
if(Dof.Type == DOF_UNKNOWN){ Dof.Case.Unknown.NumDof - 1);
LinAlg_GetScalarInVector(&tmp, &Solution_P->x, Dof.Case.Unknown.Num
Dof-1) ; if((Dof_P = (struct Dof *)List_PQuery(DofData2_P->DofList, &Dof,
fcmp_Dof))) {
if((Dof_P = (struct Dof*)List_PQuery(DofData2_P->DofList, &Dof, fcm LinAlg_SetScalarInVector(&tmp, &Solution_S.x,
p_Dof))){ Dof_P->Case.Unknown.NumDof - 1);
LinAlg_SetScalarInVector(&tmp, &Solution_S.x, Dof_P->Case.Unknown Dof_P->Type = DOF_UNKNOWN;
.NumDof-1) ; }
Dof_P->Type = DOF_UNKNOWN ; else {
} Message::Warning("Unknown Dof in TransferSolution");
else{ }
Message::Warning("Unknown Dof in TransferSolution") ; }
} else {
} // Message::Warning("Trying to transfer a non symmetrical Dof
else{ // (type %d)", Dof.Type);
// Message::Warning("Trying to transfer a non symmetrical Dof (ty }
pe %d)", }
// Dof.Type); LinAlg_AssembleVector(&Solution_S.x);
}
} if(!DofData2_P->Solutions)
LinAlg_AssembleVector(&Solution_S.x) ; DofData2_P->Solutions =
List_Create(20, 20, sizeof(struct Solution));
if (!DofData2_P->Solutions)
DofData2_P->Solutions = List_Create(20, 20, sizeof(struct Solution)) List_Add(DofData2_P->Solutions, &Solution_S);
; DofData2_P->CurrentSolution = (struct Solution *)List_Pointer(
DofData2_P->Solutions, List_Nbr(DofData2_P->Solutions) - 1);
List_Add(DofData2_P->Solutions, &Solution_S) ; }
DofData2_P->CurrentSolution = (struct Solution*) }
List_Pointer(DofData2_P->Solutions, List_Nbr(DofData2_P->Solutions)-1 break;
) ;
} case OPERATION_REMOVELASTSOLUTION:
} Init_OperationOnSystem("RemoveLastSolution", Resolution_P, Operation_P,
break ; DofData_P0, GeoData_P0, &DefineSystem_P,
&DofData_P, Resolution2_P);
case OPERATION_REMOVELASTSOLUTION : if(List_Nbr(DofData_P->Solutions)) {
Init_OperationOnSystem("RemoveLastSolution", Solution_P = (struct Solution *)List_Pointer(
Resolution_P, Operation_P, DofData_P0, GeoData_P0, DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
&DefineSystem_P, &DofData_P, Resolution2_P) ; if(Solution_P->SolutionExist) {
if(List_Nbr(DofData_P->Solutions)){
Solution_P = (struct Solution *)List_Pointer(DofData_P->Solutions,
List_Nbr(DofData_P->Solutio
ns) - 1);
if(Solution_P->SolutionExist){
Message::Info("Freeing Solution %d", Solution_P->TimeStep); Message::Info("Freeing Solution %d", Solution_P->TimeStep);
LinAlg_DestroyVector(&Solution_P->x); LinAlg_DestroyVector(&Solution_P->x);
Free(Solution_P->TimeFunctionValues) ; Free(Solution_P->TimeFunctionValues);
Solution_P->SolutionExist = 0; Solution_P->SolutionExist = 0;
} }
List_Pop(DofData_P->Solutions); List_Pop(DofData_P->Solutions);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1) DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
;
} }
else{ else {
DofData_P->CurrentSolution = 0; DofData_P->CurrentSolution = 0;
} }
break; break;
/* --> E v a l u a t e */ /* --> E v a l u a t e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_EVALUATE : case OPERATION_EVALUATE:
for(int i = 0 ; i < List_Nbr(Operation_P->Case.Evaluate.Expressions); i++) for(int i = 0; i < List_Nbr(Operation_P->Case.Evaluate.Expressions);
{ i++) {
int j; int j;
List_Read(Operation_P->Case.Evaluate.Expressions, i, &j) ; List_Read(Operation_P->Case.Evaluate.Expressions, i, &j);
Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value) ; Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value);
} }
break ; break;
/* --> S e t T i m e */ /* --> S e t T i m e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETTIME : case OPERATION_SETTIME:
Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Current.Time = Value.Val[0] ; Current.Time = Value.Val[0];
break ; break;
/* --> S e t D T i m e */ /* --> S e t D T i m e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETDTIME : case OPERATION_SETDTIME:
Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Current.DTime = Value.Val[0] ; Current.DTime = Value.Val[0];
break ; break;
/* --> S e t T i m e S t e p */ /* --> S e t T i m e S t e p */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETTIMESTEP : case OPERATION_SETTIMESTEP:
Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.SetTime.ExpressionIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Current.TimeStep = Value.Val[0] ; Current.TimeStep = Value.Val[0];
break ; break;
/* --> S e t F r e q u e n c y */ /* --> S e t F r e q u e n c y */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETFREQUENCY : case OPERATION_SETFREQUENCY:
DefineSystem_P = (struct DefineSystem*) DefineSystem_P = (struct DefineSystem *)List_Pointer(
List_Pointer(Resolution_P->DefineSystem, Resolution_P->DefineSystem, Operation_P->DefineSystemIndex);
Operation_P->DefineSystemIndex) ; DofData_P = DofData_P0 + Operation_P->DefineSystemIndex;
DofData_P = DofData_P0 + Operation_P->DefineSystemIndex ;
if(DefineSystem_P->Type == VAL_COMPLEX) {
if (DefineSystem_P->Type == VAL_COMPLEX){ if(DefineSystem_P->FrequencyValue)
if(DefineSystem_P->FrequencyValue) List_Reset(DefineSystem_P->FrequencyValue);
List_Reset(DefineSystem_P->FrequencyValue); else
else DefineSystem_P->FrequencyValue = List_Create(1, 1, sizeof(double));
DefineSystem_P->FrequencyValue = List_Create(1, 1, sizeof(double)) ; /* Provisoire: une seule frequence */
/* Provisoire: une seule frequence */ Get_ValueOfExpressionByIndex(
Get_ValueOfExpressionByIndex(Operation_P->Case.SetFrequency.ExpressionInd Operation_P->Case.SetFrequency.ExpressionIndex, NULL, 0., 0., 0.,
ex, &Value);
NULL, 0., 0., 0., &Value) ; List_Add(DefineSystem_P->FrequencyValue, &Value.Val[0]);
List_Add(DefineSystem_P->FrequencyValue, &Value.Val[0]); if(DofData_P->Pulsation == NULL)
if (DofData_P->Pulsation == NULL) DofData_P->Pulsation = List_Create(1, 2, sizeof(double));
DofData_P->Pulsation = List_Create(1, 2, sizeof(double)) ; List_Reset(DofData_P->Pulsation);
List_Reset(DofData_P->Pulsation); Init_HarInDofData(DefineSystem_P, DofData_P);
Init_HarInDofData(DefineSystem_P, DofData_P) ;
} }
else else
Message::Error("Invalid SetFrequency for real system '%s'", DefineSystem_ Message::Error("Invalid SetFrequency for real system '%s'",
P->Name) ; DefineSystem_P->Name);
break; break;
/* --> T i m e L o o p T h e t a */ /* --> T i m e L o o p T h e t a */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_TIMELOOPTHETA : case OPERATION_TIMELOOPTHETA:
if(!List_Nbr(Current.DofData->Solutions)){ if(!List_Nbr(Current.DofData->Solutions)) {
Message::Error("Not enough initial solutions for TimeLoopTheta"); Message::Error("Not enough initial solutions for TimeLoopTheta");
break; break;
} }
Message::Info("TimeLoopTheta ...") ; Message::Info("TimeLoopTheta ...");
Save_TypeTime = Current.TypeTime ; Save_TypeTime = Current.TypeTime;
Save_DTime = Current.DTime ; Save_DTime = Current.DTime;
Flag_NextThetaFixed = 0 ; /* Attention: Test */ Flag_NextThetaFixed = 0; /* Attention: Test */
Current.TypeTime = TIME_THETA ; Current.TypeTime = TIME_THETA;
if(Flag_RESTART) { if(Flag_RESTART) {
if (Current.Time < Operation_P->Case.TimeLoopTheta.TimeMax * 0.999999) if(Current.Time < Operation_P->Case.TimeLoopTheta.TimeMax * 0.999999)
Flag_RESTART = 0 ; Flag_RESTART = 0;
} }
else else
Current.Time = Operation_P->Case.TimeLoopTheta.Time0 ; Current.Time = Operation_P->Case.TimeLoopTheta.Time0;
Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopTheta.ThetaIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopTheta.ThetaIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Current.Theta = Value.Val[0] ; Current.Theta = Value.Val[0];
Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopTheta.DTimeIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopTheta.DTimeIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Current.DTime = Value.Val[0] ; Current.DTime = Value.Val[0];
while (Current.Time < Operation_P->Case.TimeLoopTheta.TimeMax * 0.999999)
{
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) bre while(Current.Time < Operation_P->Case.TimeLoopTheta.TimeMax * 0.999999) {
ak; if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
break;
if (!Flag_NextThetaFixed) { // Attention: Test if(!Flag_NextThetaFixed) { // Attention: Test
Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopTheta.ThetaInde Get_ValueOfExpressionByIndex(
x, Operation_P->Case.TimeLoopTheta.ThetaIndex, NULL, 0., 0., 0.,
NULL, 0., 0., 0., &Value) ; &Value);
Current.Theta = Value.Val[0] ; Current.Theta = Value.Val[0];
} }
Expression *DTimeExpr = (struct Expression *)List_Pointer Expression *DTimeExpr = (struct Expression *)List_Pointer(
(Problem_S.Expression, Operation_P->Case.TimeLoopTheta.DTimeIndex); Problem_S.Expression, Operation_P->Case.TimeLoopTheta.DTimeIndex);
// don't reevaluate if constant (it could have been changed via SetDTime // don't reevaluate if constant (it could have been changed via SetDTime
// in a previous operation) // in a previous operation)
if(!Is_ExpressionConstant(DTimeExpr) && Flag_NextThetaFixed != 2) { // A if(!Is_ExpressionConstant(DTimeExpr) &&
ttention: Test Flag_NextThetaFixed != 2) { // Attention: Test
Get_ValueOfExpression(DTimeExpr, NULL, 0., 0., 0., &Value) ; Get_ValueOfExpression(DTimeExpr, NULL, 0., 0., 0., &Value);
Current.DTime = Value.Val[0] ; Current.DTime = Value.Val[0];
} }
Flag_NextThetaFixed = 0 ; Flag_NextThetaFixed = 0;
Current.Time += Current.DTime ;
Current.TimeStep += 1. ;
Current.SolveJacAdaptFailed = 0;
Message::Info(3, "Theta Time = %.8g s (TimeStep %d, DTime %g)", Current.
Time,
(int)Current.TimeStep, Current.DTime) ;
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep( Current.Time += Current.DTime;
) < 100) Current.TimeStep += 1.;
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Time Current.SolveJacAdaptFailed = 0;
",
Message::Info(3, "Theta Time = %.8g s (TimeStep %d, DTime %g)",
Current.Time, (int)Current.TimeStep, Current.DTime);
if(Message::GetProgressMeterStep() > 0 &&
Message::GetProgressMeterStep() < 100)
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
"/Time",
std::vector<double>(1, Current.Time)); std::vector<double>(1, Current.Time));
// Save_Time = Current.Time ; // removed: prevents using SetTime in the // Save_Time = Current.Time ; // removed: prevents using SetTime in the
operation // operation
Treatment_Operation(Resolution_P, Operation_P->Case.TimeLoopTheta.Operati Treatment_Operation(Resolution_P,
on, Operation_P->Case.TimeLoopTheta.Operation,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
// Current.Time = Save_Time ; // removed: prevents using SetTime in the o // Current.Time = Save_Time ; // removed: prevents using SetTime in the
peration // operation
if(Flag_Break){ Flag_Break = 0; break; } if(Flag_Break) {
Flag_Break = 0;
break;
}
} }
Current.TypeTime = Save_TypeTime ; Current.TypeTime = Save_TypeTime;
Current.DTime = Save_DTime ; Current.DTime = Save_DTime;
break ; break;
/* --> T i m e L o o p N e w m a r k */ /* --> T i m e L o o p N e w m a r k */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_TIMELOOPNEWMARK : case OPERATION_TIMELOOPNEWMARK:
if(List_Nbr(Current.DofData->Solutions) < 2){ if(List_Nbr(Current.DofData->Solutions) < 2) {
Message::Error("Not enough initial solutions for TimeLoopNewmark"); Message::Error("Not enough initial solutions for TimeLoopNewmark");
break; break;
} }
Message::Info("TimeLoopNewmark ...") ; Message::Info("TimeLoopNewmark ...");
Save_TypeTime = Current.TypeTime ; Save_TypeTime = Current.TypeTime;
Save_DTime = Current.DTime ; Save_DTime = Current.DTime;
Current.Beta = Operation_P->Case.TimeLoopNewmark.Beta ; Current.Beta = Operation_P->Case.TimeLoopNewmark.Beta;
Current.Gamma = Operation_P->Case.TimeLoopNewmark.Gamma ; Current.Gamma = Operation_P->Case.TimeLoopNewmark.Gamma;
Current.TypeTime = TIME_NEWMARK ; Current.TypeTime = TIME_NEWMARK;
if(Flag_RESTART){ if(Flag_RESTART) {
if (Current.Time < Operation_P->Case.TimeLoopNewmark.TimeMax * 0.999999) if(Current.Time < Operation_P->Case.TimeLoopNewmark.TimeMax * 0.999999)
Flag_RESTART = 0 ; Flag_RESTART = 0;
} }
else else
Current.Time = Operation_P->Case.TimeLoopNewmark.Time0 ; Current.Time = Operation_P->Case.TimeLoopNewmark.Time0;
while (Current.Time < Operation_P->Case.TimeLoopNewmark.TimeMax * 0.999999
) {
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) bre while(Current.Time <
ak; Operation_P->Case.TimeLoopNewmark.TimeMax * 0.999999) {
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
break;
Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopNewmark.DTimeIndex Get_ValueOfExpressionByIndex(
, Operation_P->Case.TimeLoopNewmark.DTimeIndex, NULL, 0., 0., 0.,
NULL, 0., 0., 0., &Value) ; &Value);
Current.DTime = Value.Val[0] ; Current.DTime = Value.Val[0];
Current.Time += Current.DTime ; Current.Time += Current.DTime;
Current.TimeStep += 1. ; Current.TimeStep += 1.;
Current.SolveJacAdaptFailed = 0 ; Current.SolveJacAdaptFailed = 0;
Message::Info(3, "Newmark Time = %.8g s (TimeStep %d)", Current.Time, Message::Info(3, "Newmark Time = %.8g s (TimeStep %d)", Current.Time,
(int)Current.TimeStep) ; (int)Current.TimeStep);
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep( if(Message::GetProgressMeterStep() > 0 &&
) < 100) Message::GetProgressMeterStep() < 100)
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Time Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
", "/Time",
std::vector<double>(1, Current.Time)); std::vector<double>(1, Current.Time));
Treatment_Operation(Resolution_P, Operation_P->Case.TimeLoopNewmark.Opera Treatment_Operation(Resolution_P,
tion, Operation_P->Case.TimeLoopNewmark.Operation,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
if(Flag_Break){ Flag_Break = 0; break; } if(Flag_Break) {
Flag_Break = 0;
break;
}
} }
Current.TypeTime = Save_TypeTime ; Current.TypeTime = Save_TypeTime;
Current.DTime = Save_DTime ; Current.DTime = Save_DTime;
break ; break;
/* --> I t e r a t i v e L o o p */ /* --> I t e r a t i v e L o o p */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_ITERATIVELOOP : case OPERATION_ITERATIVELOOP:
Message::Info("IterativeLoop ...") ; Message::Info("IterativeLoop ...");
Save_Iteration = Current.Iteration ; Save_Iteration = Current.Iteration;
//Current.Residual_Iter1=1.0; //kj+++ to compute a relative residual (rela for(Num_Iteration = 1;
tive to residual at iter 1) (not necessary) Num_Iteration <= Operation_P->Case.IterativeLoop.NbrMaxIteration;
// abstol = Operation_P->Case.IterativeLoop.Criterion ; Num_Iteration++) {
// reltol = Operation_P->Case.IterativeLoop.Criterion/100 ; if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
break;
for (Num_Iteration = 1 ;
Num_Iteration <= Operation_P->Case.IterativeLoop.NbrMaxIteration ;
Num_Iteration++) {
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) bre Current.Iteration = (double)Num_Iteration;
ak; Current.RelativeDifference = 0.;
Current.Iteration = (double)Num_Iteration ; Get_ValueOfExpressionByIndex(
Current.RelativeDifference = 0. ; Operation_P->Case.IterativeLoop.RelaxationFactorIndex, NULL, 0., 0.,
0., &Value);
if(Current.RelaxationFactor != Value.Val[0] || Num_Iteration == 1) {
Current.RelaxationFactor = Value.Val[0];
Message::Info("Nonlinear Iteration Relaxation %g",
Current.RelaxationFactor);
}
Flag_IterativeLoop =
Operation_P->Case.IterativeLoop.Flag; /* Attention: Test */
// NB: SolveJac OR SolveJacAdapt are called here
// Resolution2_P and DofData2_P0 added as arguments for allowing
// TransferSolution of a nonlinear resolution
Treatment_Operation(Resolution_P,
Operation_P->Case.IterativeLoop.Operation,
DofData_P0, GeoData_P0, Resolution2_P, DofData2_P0);
if(Current.RelaxFac == 0) {
// SolveJacAdapt has not been called
// ==> Copy the default RelaxationFactor in RelaxFac
Current.RelaxFac = Current.RelaxationFactor;
}
// NB: Current.RelativeDifference is what is used for classical
// IterativeLoop stop criterion NB: In SolveJac:
// (Current.RelativeDifference+=Current.Residual) NB: In SolveJacAdapt:
// (Current.RelativeDifference=Current.Residual)
if((Current.RelativeDifference <=
Operation_P->Case.IterativeLoop.Criterion) ||
(Current.RelativeDifference !=
Current.RelativeDifference)) // NaN or Inf
break;
Get_ValueOfExpressionByIndex if(Flag_Break) {
(Operation_P->Case.IterativeLoop.RelaxationFactorIndex, Flag_Break = 0;
NULL, 0., 0., 0., &Value) ; break;
if(Current.RelaxationFactor != Value.Val[0] || Num_Iteration == 1){
Current.RelaxationFactor = Value.Val[0] ;
Message::Info("Nonlinear Iteration Relaxation %g", Current.Relaxation
Factor) ;
} }
Flag_IterativeLoop = Operation_P->Case.IterativeLoop.Flag ; /* Attention: Current.RelativeDifferenceOld =
Test */ Current.RelativeDifference; /* Attention: pt */
//NB: SolveJac OR SolveJacAdapt are called here
// Resolution2_P and DofData2_P0 added as arguments for allowing
// TransferSolution of a nonlinear resolution
Treatment_Operation(Resolution_P, Operation_P->Case.IterativeLoop.Operati
on,
DofData_P0, GeoData_P0, Resolution2_P, DofData2_P0) ;
if (Current.RelaxFac==0)
// SolveJacAdapt has not been called
// ==> Copy the default RelaxationFactor in RelaxFac
Current.RelaxFac = Current.RelaxationFactor ; // +++
//if (Current.Iteration==1) Current.Residual_Iter1=Current.RelativeDifference;
//kj+++ to compute a relative residual (relative to residual at iter 1) (not ne
cessary)
//NB: Current.RelativeDifference is what is used for classical IterativeLoop s
top criterion
//NB: In SolveJac: (Current.RelativeDifference+=Current.Residual)
//NB: In SolveJacAdapt: (Current.RelativeDifference=Current.Residual)
if ( (Current.RelativeDifference <= Operation_P->Case.IterativeLoop.Cri
terion) ||
//(Current.RelativeDifference/Current.Residual_Iter1 <= Operation_P->Case
.IterativeLoop.Criterion) || //kj+++ to compute a relative residual (relative to
residual at iter 1) (not necessary)
(Current.RelativeDifference != Current.RelativeDifference) ) // NaN or
Inf
break ;
if(Flag_Break){ Flag_Break = 0; break; }
Current.RelativeDifferenceOld = Current.RelativeDifference ; /* Attention
: pt */
} }
if ( (Num_Iteration > Operation_P->Case.IterativeLoop.NbrMaxIteration) if((Num_Iteration > Operation_P->Case.IterativeLoop.NbrMaxIteration) ||
|| (Current.RelativeDifference != Current.RelativeDifference) ) // Na (Current.RelativeDifference !=
N or Inf Current.RelativeDifference)) // NaN or Inf
{ {
// Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration ; // Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration ;
Flag_IterativeLoopConverged = 0; Flag_IterativeLoopConverged = 0;
//Message::Info(3, "IterativeLoop did NOT converge (%d iterations, residual %g // Message::Info(3, "IterativeLoop did NOT converge (%d iterations,
)", // residual %g)",
Message::Warning("IterativeLoop did NOT converge (%d iterations, residua Message::Warning(
l %g)", "IterativeLoop did NOT converge (%d iterations, residual %g)",
Num_Iteration, Current.RelativeDifference); Num_Iteration, Current.RelativeDifference);
// Either it has reached the max num of iterations or a NaN at a given // Either it has reached the max num of iterations or a NaN at a given
iteration // iteration
Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration ; Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration;
} }
else{ else {
Message::Info(3, "IterativeLoop converged (%d iteration%s, residual %g)" Message::Info(3,
, "IterativeLoop converged (%d iteration%s, residual %g)",
Num_Iteration, Num_Iteration > 1 ? "s" : "", Num_Iteration, Num_Iteration > 1 ? "s" : "",
Current.RelativeDifference); Current.RelativeDifference);
} }
/* kj+++ to compute a relative residual (relative to residual at iter 1) ( Current.Iteration = Save_Iteration;
not necessary) break;
Message::Info(3, "IterativeLoop did NOT converge (%d iterations, residual
%g)\n rel %g",
Num_Iteration, Current.RelativeDifference, Current.Relat
iveDifference/Current.Residual_Iter1);
// Either it has reached the max num of iterations or a NaN at a given
iteration
Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration ;
}
else{
Message::Info(3, "IterativeLoop converged (%d iteration%s, residual %g)\n
rel %g",
Num_Iteration, Num_Iteration > 1 ? "s" : "",
Current.RelativeDifference, Current.RelativeDifference/C
urrent.Residual_Iter1);
}
*/
Current.Iteration = Save_Iteration ;
break ;
case OPERATION_ITERATIVELINEARSOLVER : case OPERATION_ITERATIVELINEARSOLVER:
Message::Info("IterativeLinearSolver ...") ; Message::Info("IterativeLinearSolver ...");
Operation_IterativeLinearSolver Operation_IterativeLinearSolver(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_BROADCASTFIELDS : case OPERATION_BROADCASTFIELDS:
Message::Info("BroadcastFields ...") ; Message::Info("BroadcastFields ...");
Operation_BroadcastFields Operation_BroadcastFields(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_BROADCASTVARIABLES : case OPERATION_BROADCASTVARIABLES:
Message::Info("BroadcastVariables ...") ; Message::Info("BroadcastVariables ...");
Operation_BroadcastVariables Operation_BroadcastVariables(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_GATHERVARIABLES : case OPERATION_GATHERVARIABLES:
Message::Info("GatherVariables ...") ; Message::Info("GatherVariables ...");
Operation_GatherVariables Operation_GatherVariables(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_SCATTERVARIABLES : case OPERATION_SCATTERVARIABLES:
Message::Info("ScatterVariables ...") ; Message::Info("ScatterVariables ...");
Operation_ScatterVariables Operation_ScatterVariables(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_CLEARVARIABLES : case OPERATION_CLEARVARIABLES: {
{ Message::Info("ClearVariables ...");
Message::Info("ClearVariables ...") ;
if (List_Nbr(Operation_P->Case.ClearVariables.Names)==0) if(List_Nbr(Operation_P->Case.ClearVariables.Names) == 0) {
{
Message::Info("ClearVariables: Clear All Run-time Variables"); Message::Info("ClearVariables: Clear All Run-time Variables");
Get_AllValueSaved().clear(); Get_AllValueSaved().clear();
} }
else else {
{
std::map<std::string, struct Value> &values = Get_AllValueSaved(); std::map<std::string, struct Value> &values = Get_AllValueSaved();
for(int i = 0; i < List_Nbr(Operation_P->Case.ClearVariables.Names); i++ for(int i = 0; i < List_Nbr(Operation_P->Case.ClearVariables.Names);
){ i++) {
char *s; char *s;
List_Read(Operation_P->Case.ClearVariables.Names, i, &s); List_Read(Operation_P->Case.ClearVariables.Names, i, &s);
if(values.find(s) != values.end()) if(values.find(s) != values.end()) {
{
Message::Info("ClearVariables: Clear Run-time Variable %s", s); Message::Info("ClearVariables: Clear Run-time Variable %s", s);
values.erase(s); values.erase(s);
} }
else else
Message::Info("ClearVariables: Unknown Run-time Variable %s", s); Message::Info("ClearVariables: Unknown Run-time Variable %s", s);
} }
} }
} } break;
break;
case OPERATION_CHECKVARIABLES : case OPERATION_CHECKVARIABLES:
Message::Info("CheckVariables ...") ; Message::Info("CheckVariables ...");
Operation_CheckVariables Operation_CheckVariables(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break; break;
case OPERATION_CLEARVECTORS : case OPERATION_CLEARVECTORS: {
{ Message::Info("ClearVectors ...");
Message::Info("ClearVectors ...") ; Operation_ClearVectors(Operation_P, DofData_P);
Operation_ClearVectors } break;
( Operation_P, DofData_P);
}
break;
/* --> I t e r a t i v e T i m e R e d u c t i o n */ /* --> I t e r a t i v e T i m e R e d u c t i o n */
/* ------------------------------------------------ */ /* ------------------------------------------------ */
case OPERATION_ITERATIVETIMEREDUCTION : case OPERATION_ITERATIVETIMEREDUCTION:
Message::Info("IterativeTimeReduction ...") ; Message::Info("IterativeTimeReduction ...");
Operation_IterativeTimeReduction Operation_IterativeTimeReduction(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break ; break;
/* --> T e s t */ /* --> T e s t */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_TEST : case OPERATION_TEST:
Message::Info("Test") ; Message::Info("Test");
Get_ValueOfExpressionByIndex(Operation_P->Case.Test.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.Test.ExpressionIndex, NULL,
NULL, 0., 0., 0., &Value) ; 0., 0., 0., &Value);
if(Value.Val[0]){ if(Value.Val[0]) {
Treatment_Operation(Resolution_P, Operation_P->Case.Test.Operation_True, Treatment_Operation(Resolution_P, Operation_P->Case.Test.Operation_True,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
} }
else{ else {
if(Operation_P->Case.Test.Operation_False) if(Operation_P->Case.Test.Operation_False)
Treatment_Operation(Resolution_P, Operation_P->Case.Test.Operation_Fals Treatment_Operation(Resolution_P,
e, Operation_P->Case.Test.Operation_False,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
} }
break ; break;
/* --> W h i l e */ /* --> W h i l e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_WHILE : case OPERATION_WHILE:
Message::Info("While...") ; Message::Info("While...");
while(1){ while(1) {
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) bre if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
ak; break;
Get_ValueOfExpressionByIndex(Operation_P->Case.While.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.While.ExpressionIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
if(!Value.Val[0]) break; if(!Value.Val[0]) break;
Treatment_Operation(Resolution_P, Operation_P->Case.While.Operation, Treatment_Operation(Resolution_P, Operation_P->Case.While.Operation,
DofData_P0, GeoData_P0, NULL, NULL) ; DofData_P0, GeoData_P0, NULL, NULL);
if(Flag_Break){ Flag_Break = 0; break; } if(Flag_Break) {
Flag_Break = 0;
break;
}
} }
break ; break;
/* --> F o u r i e r T r a n s f o r m */ /* --> F o u r i e r T r a n s f o r m */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_FOURIERTRANSFORM2 : case OPERATION_FOURIERTRANSFORM2:
Message::Info("FourierTransform") ; Message::Info("FourierTransform");
if(gSCALAR_SIZE == 2){ if(gSCALAR_SIZE == 2) {
Message::Error("FIXME: FourierTransform2 will not work in complex arithme Message::Error(
tic"); "FIXME: FourierTransform2 will not work in complex arithmetic");
break; break;
} }
DofData_P = DofData_P0 + Operation_P->Case.FourierTransform2.DefineSystem DofData_P =
Index[0] ; DofData_P0 + Operation_P->Case.FourierTransform2.DefineSystemIndex[0];
DofData2_P = DofData_P0 + Operation_P->Case.FourierTransform2.DefineSystem DofData2_P =
Index[1] ; DofData_P0 + Operation_P->Case.FourierTransform2.DefineSystemIndex[1];
NbrHar1 = DofData_P->NbrHar;
NbrDof1 = List_Nbr(DofData_P->DofList);
NbrHar2 = DofData2_P->NbrHar;
NbrDof2 = List_Nbr(DofData2_P->DofList);
NbrHar1 = DofData_P->NbrHar ; if(NbrHar1 != 1 || NbrHar2 < 2 || NbrDof2 != (NbrDof1 * NbrHar2)) {
NbrDof1 = List_Nbr(DofData_P->DofList) ; Message::Error("Uncompatible System definitions for FourierTransform"
NbrHar2 = DofData2_P->NbrHar ;
NbrDof2 = List_Nbr(DofData2_P->DofList) ;
if (NbrHar1 != 1 || NbrHar2 < 2 || NbrDof2 != (NbrDof1*NbrHar2)){
Message::Error("Uncompatible System definitions for FourierTransform"
" (NbrHar = %d|%d NbrDof = %d|%d)", " (NbrHar = %d|%d NbrDof = %d|%d)",
NbrHar1, NbrHar2, NbrDof1, NbrDof2) ; NbrHar1, NbrHar2, NbrDof1, NbrDof2);
break; break;
} }
if(!DofData2_P->Solutions){ if(!DofData2_P->Solutions) {
DofData2_P->Solutions = List_Create(1, 1, sizeof(struct Solution)) ; DofData2_P->Solutions = List_Create(1, 1, sizeof(struct Solution));
Operation_P->Case.FourierTransform2.Scales = (double *)Malloc(NbrHar2*siz Operation_P->Case.FourierTransform2.Scales =
eof(double)) ; (double *)Malloc(NbrHar2 * sizeof(double));
} }
Nbr_Sol = List_Nbr(DofData2_P->Solutions) ; Nbr_Sol = List_Nbr(DofData2_P->Solutions);
Scales = Operation_P->Case.FourierTransform2.Scales ; Scales = Operation_P->Case.FourierTransform2.Scales;
if ( (Operation_P->Case.FourierTransform2.Period_sofar + Current.DTime > if((Operation_P->Case.FourierTransform2.Period_sofar + Current.DTime >
Operation_P->Case.FourierTransform2.Period) && Nbr_Sol ) { Operation_P->Case.FourierTransform2.Period) &&
Message::Info("Normalizing and finalizing Fourier Analysis" Nbr_Sol) {
Message::Info("Normalizing and finalizing Fourier Analysis"
" (solution %d) (Period: %e out of %e)", " (solution %d) (Period: %e out of %e)",
Nbr_Sol, Operation_P->Case.FourierTransform2.Period_sofar, Nbr_Sol, Operation_P->Case.FourierTransform2.Period_sofar,
Operation_P->Case.FourierTransform2.Period); Operation_P->Case.FourierTransform2.Period);
for (int i = 0; i < NbrHar2; i++) for(int i = 0; i < NbrHar2; i++)
Message::Info("Har %d : Scales %e ", i, Scales[i]) ; Message::Info("Har %d : Scales %e ", i, Scales[i]);
Solution_P = (struct Solution*)List_Pointer(DofData2_P->Solutions, Nbr_So Solution_P =
l-1); (struct Solution *)List_Pointer(DofData2_P->Solutions, Nbr_Sol - 1);
for(int j = 0; j<DofData2_P->NbrDof; j += NbrHar2){ for(int j = 0; j < DofData2_P->NbrDof; j += NbrHar2) {
NumDof = ((struct Dof *)List_Pointer(DofData2_P->DofList,j))->Case.Unkn NumDof = ((struct Dof *)List_Pointer(DofData2_P->DofList, j))
own.NumDof - 1 ; ->Case.Unknown.NumDof -
for(int k = 0; k<NbrHar2; k++){ 1;
LinAlg_GetDoubleInVector(&d1, &Solution_P->x, NumDof+k) ; for(int k = 0; k < NbrHar2; k++) {
if (Scales[k]) d1 /= Scales[k] ; LinAlg_GetDoubleInVector(&d1, &Solution_P->x, NumDof + k);
LinAlg_SetDoubleInVector(d1, &Solution_P->x, NumDof+k) ; if(Scales[k]) d1 /= Scales[k];
} LinAlg_SetDoubleInVector(d1, &Solution_P->x, NumDof + k);
} }
Operation_P->Case.FourierTransform2.Period_sofar = 0 ; }
break; Operation_P->Case.FourierTransform2.Period_sofar = 0;
} break;
}
if (Operation_P->Case.FourierTransform2.Period_sofar == 0) {
Message::Info("Starting new Fourier Analysis : solution %d ", Nbr_Sol); if(Operation_P->Case.FourierTransform2.Period_sofar == 0) {
Solution_S.TimeStep = Nbr_Sol; Message::Info("Starting new Fourier Analysis : solution %d ", Nbr_Sol);
Solution_S.Time = Nbr_Sol; Solution_S.TimeStep = Nbr_Sol;
Solution_S.Time = Nbr_Sol;
Solution_S.TimeFunctionValues = NULL; Solution_S.TimeFunctionValues = NULL;
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver, DofData2_P->NbrDo LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver,
f) ; DofData2_P->NbrDof);
LinAlg_ZeroVector(&Solution_S.x) ; LinAlg_ZeroVector(&Solution_S.x);
List_Add(DofData2_P->Solutions, &Solution_S) ; List_Add(DofData2_P->Solutions, &Solution_S);
Nbr_Sol++ ; Nbr_Sol++;
for (int k = 0; k<NbrHar2; k++) Scales[k] = 0 ; for(int k = 0; k < NbrHar2; k++) Scales[k] = 0;
} }
DofData2_P->CurrentSolution = Solution_P = DofData2_P->CurrentSolution = Solution_P =
(struct Solution*)List_Pointer(DofData2_P->Solutions, Nbr_Sol-1) ; (struct Solution *)List_Pointer(DofData2_P->Solutions, Nbr_Sol - 1);
for (int k = 0; k<NbrHar2; k+=2) { for(int k = 0; k < NbrHar2; k += 2) {
d = DofData2_P->Val_Pulsation[k/2] * Current.Time ; d = DofData2_P->Val_Pulsation[k / 2] * Current.Time;
Scales[k ] += cos(d) * cos(d) * Current.DTime ; Scales[k] += cos(d) * cos(d) * Current.DTime;
Scales[k+1] += sin(d) * sin(d) * Current.DTime ; Scales[k + 1] += sin(d) * sin(d) * Current.DTime;
} }
for(int j = 0; j < NbrDof1; j++){ for(int j = 0; j < NbrDof1; j++) {
Dof_GetRealDofValue(DofData_P, (struct Dof *)List_Pointer(DofData_P->DofL Dof_GetRealDofValue(
ist,j), &dd) ; DofData_P, (struct Dof *)List_Pointer(DofData_P->DofList, j), &dd);
NumDof = ((struct Dof *)List_Pointer(DofData2_P->DofList, NumDof = ((struct Dof *)List_Pointer(DofData2_P->DofList, j * NbrHar2))
j*NbrHar2))->Case.Unknown.NumDof - 1 ->Case.Unknown.NumDof -
; 1;
if (((struct Dof *)List_Pointer(DofData2_P->DofList,j*NbrHar2))->Type != if(((struct Dof *)List_Pointer(DofData2_P->DofList, j * NbrHar2))
DOF_UNKNOWN) ->Type != DOF_UNKNOWN)
Message::Info("Dof not unknown %d", j) ; Message::Info("Dof not unknown %d", j);
for (int k = 0; k < NbrHar2; k+=2) { for(int k = 0; k < NbrHar2; k += 2) {
d = DofData2_P->Val_Pulsation[k/2] * Current.Time ; d = DofData2_P->Val_Pulsation[k / 2] * Current.Time;
LinAlg_AddDoubleInVector( dd*cos(d)*Current.DTime, &Solution_P->x, NumD LinAlg_AddDoubleInVector(dd * cos(d) * Current.DTime, &Solution_P->x,
of+k ) ; NumDof + k);
LinAlg_AddDoubleInVector(-dd*sin(d)*Current.DTime, &Solution_P->x, NumD LinAlg_AddDoubleInVector(-dd * sin(d) * Current.DTime, &Solution_P->x,
of+k+1) ; NumDof + k + 1);
} }
} }
Operation_P->Case.FourierTransform2.Period_sofar += Current.DTime ; Operation_P->Case.FourierTransform2.Period_sofar += Current.DTime;
break; break;
case OPERATION_FOURIERTRANSFORM : case OPERATION_FOURIERTRANSFORM:
Message::Info("FourierTransform") ; Message::Info("FourierTransform");
DofData_P = DofData_P0 + Operation_P->Case.FourierTransform.DefineSystemIn DofData_P =
dex[0] ; DofData_P0 + Operation_P->Case.FourierTransform.DefineSystemIndex[0];
DofData2_P = DofData_P0 + Operation_P->Case.FourierTransform.DefineSystemI DofData2_P =
ndex[1] ; DofData_P0 + Operation_P->Case.FourierTransform.DefineSystemIndex[1];
if(!DofData2_P->Solutions){ if(!DofData2_P->Solutions) {
int k = List_Nbr(Operation_P->Case.FourierTransform.Frequency) ; int k = List_Nbr(Operation_P->Case.FourierTransform.Frequency);
if(DofData2_P->NbrDof != gCOMPLEX_INCREMENT * DofData_P->NbrDof){ if(DofData2_P->NbrDof != gCOMPLEX_INCREMENT * DofData_P->NbrDof) {
Message::Error("Uncompatible System definitions for FourierTransform") Message::Error(
; "Uncompatible System definitions for FourierTransform");
break; break;
} }
DofData2_P->Solutions = List_Create(k, 1, sizeof(struct Solution)) ; DofData2_P->Solutions = List_Create(k, 1, sizeof(struct Solution));
for(int i = 0; i < k; i++){ for(int i = 0; i < k; i++) {
List_Read(Operation_P->Case.FourierTransform.Frequency, i, &d) ; List_Read(Operation_P->Case.FourierTransform.Frequency, i, &d);
Solution_S.TimeStep = i ; Solution_S.TimeStep = i;
Solution_S.Time = TWO_PI * d; Solution_S.Time = TWO_PI * d;
Solution_S.TimeImag = 0.; Solution_S.TimeImag = 0.;
Solution_S.TimeFunctionValues = NULL; Solution_S.TimeFunctionValues = NULL;
Solution_S.SolutionExist = 1 ; Solution_S.SolutionExist = 1;
LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver, DofData2_P->Nbr LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver,
Dof) ; DofData2_P->NbrDof);
LinAlg_ZeroVector(&Solution_S.x) ; LinAlg_ZeroVector(&Solution_S.x);
List_Add(DofData2_P->Solutions, &Solution_S) ; List_Add(DofData2_P->Solutions, &Solution_S);
} }
DofData2_P->CurrentSolution = DofData2_P->CurrentSolution =
(struct Solution*)List_Pointer(DofData2_P->Solutions, k/2) ; (struct Solution *)List_Pointer(DofData2_P->Solutions, k / 2);
} }
for(int i = 0; i < List_Nbr(DofData2_P->Solutions); i++){ for(int i = 0; i < List_Nbr(DofData2_P->Solutions); i++) {
Solution_P = (struct Solution*)List_Pointer(DofData2_P->Solutions, i); Solution_P = (struct Solution *)List_Pointer(DofData2_P->Solutions, i);
d = Solution_P->Time * Current.Time ; d = Solution_P->Time * Current.Time;
for(int j=0,k=0 ; j<DofData_P->NbrDof ; j++,k+=gCOMPLEX_INCREMENT){ for(int j = 0, k = 0; j < DofData_P->NbrDof;
LinAlg_GetDoubleInVector(&d2, &DofData_P->CurrentSolution->x, j); j++, k += gCOMPLEX_INCREMENT) {
LinAlg_AddComplexInVector( d2 * cos(d) * Current.DTime, LinAlg_GetDoubleInVector(&d2, &DofData_P->CurrentSolution->x, j);
-d2 * sin(d) * Current.DTime, LinAlg_AddComplexInVector(d2 * cos(d) * Current.DTime,
&Solution_P->x, k, k+1) ; -d2 * sin(d) * Current.DTime,
} &Solution_P->x, k, k + 1);
}
} }
break; break;
/* --> P r i n t / W r i t e */ /* --> P r i n t / W r i t e */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_WRITE : Flag_Binary = 1 ; case OPERATION_WRITE: Flag_Binary = 1;
case OPERATION_PRINT : case OPERATION_PRINT:
if(Operation_P->Case.Print.FileOut){ if(Operation_P->Case.Print.FileOut) {
if(Operation_P->Case.Print.FileOut[0] == '/' || if(Operation_P->Case.Print.FileOut[0] == '/' ||
Operation_P->Case.Print.FileOut[0] == '\\'){ Operation_P->Case.Print.FileOut[0] == '\\') {
strcpy(FileName, Operation_P->Case.Print.FileOut); strcpy(FileName, Operation_P->Case.Print.FileOut);
} }
else{ else {
strcpy(FileName, Name_Path); strcpy(FileName, Name_Path);
strcat(FileName, Operation_P->Case.Print.FileOut); strcat(FileName, Operation_P->Case.Print.FileOut);
} }
if(!(fp = FOpen(FileName, "ab"))){ if(!(fp = FOpen(FileName, "ab"))) {
Message::Error("Unable to open file '%s'", FileName) ; Message::Error("Unable to open file '%s'", FileName);
break; break;
} }
Message::Info("Print -> '%s'", FileName) ; Message::Info("Print -> '%s'", FileName);
} }
else{ else {
fp = stdout ; fp = stdout;
Message::Info("Print") ; Message::Info("Print");
} }
if(Operation_P->Case.Print.Expressions){ if(Operation_P->Case.Print.Expressions) {
List_T *list = 0; List_T *list = 0;
if(Operation_P->Case.Print.FormatString) if(Operation_P->Case.Print.FormatString)
list = List_Create(10, 10, sizeof(double)); list = List_Create(10, 10, sizeof(double));
for(int i = 0; i < List_Nbr(Operation_P->Case.Print.Expressions); i++){ for(int i = 0; i < List_Nbr(Operation_P->Case.Print.Expressions); i++) {
int j; int j;
List_Read(Operation_P->Case.Print.Expressions, i, &j) ; List_Read(Operation_P->Case.Print.Expressions, i, &j);
Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value) ; Get_ValueOfExpressionByIndex(j, NULL, 0., 0., 0., &Value);
if(list) if(list)
List_Add(list, &Value.Val[0]); List_Add(list, &Value.Val[0]);
else else
Print_Value(&Value, fp) ; Print_Value(&Value, fp);
} }
if(list){ if(list) {
char buffer[1024]; char buffer[1024];
Print_ListOfDouble(Operation_P->Case.Print.FormatString, list, buffer) Print_ListOfDouble(Operation_P->Case.Print.FormatString, list,
; buffer);
Message::Direct(3, buffer); Message::Direct(3, buffer);
if(fp != stdout) fprintf(fp, "%s\n", buffer); if(fp != stdout) fprintf(fp, "%s\n", buffer);
List_Delete(list); List_Delete(list);
} }
} }
else if (Operation_P->Case.Print.DofNumber){ else if(Operation_P->Case.Print.DofNumber) {
DofData_P = DofData_P0 + Operation_P->DefineSystemIndex ; DofData_P = DofData_P0 + Operation_P->DefineSystemIndex;
for(int i = 0; i < List_Nbr(Operation_P->Case.Print.DofNumber); i++){ for(int i = 0; i < List_Nbr(Operation_P->Case.Print.DofNumber); i++) {
int j = *(int*)List_Pointer(Operation_P->Case.Print.DofNumber, i) ; int j = *(int *)List_Pointer(Operation_P->Case.Print.DofNumber, i);
if(j >= 0 && j < DofData_P->NbrDof){ if(j >= 0 && j < DofData_P->NbrDof) {
if(Operation_P->Case.Print.TimeStep) if(Operation_P->Case.Print.TimeStep)
for(int k = 0 ; k < List_Nbr(Operation_P->Case.Print.TimeStep); k++ for(int k = 0; k < List_Nbr(Operation_P->Case.Print.TimeStep);
){ k++) {
int l = *(int*)List_Pointer(Operation_P->Case.Print.TimeStep, k) int l =
; *(int *)List_Pointer(Operation_P->Case.Print.TimeStep, k);
if(l >= 0 && l < List_Nbr(DofData_P->Solutions)){ if(l >= 0 && l < List_Nbr(DofData_P->Solutions)) {
Solution_P = (struct Solution*)List_Pointer(DofData_P->Solution Solution_P =
s, l) ; (struct Solution *)List_Pointer(DofData_P->Solutions, l);
LinAlg_GetScalarInVector(&tmp, &Solution_P->x, j) ; LinAlg_GetScalarInVector(&tmp, &Solution_P->x, j);
if(Flag_Binary){ if(Flag_Binary) { LinAlg_WriteScalar(fp, &tmp); }
LinAlg_WriteScalar(fp, &tmp) ; else {
} LinAlg_PrintScalar(fp, &tmp);
else{ fprintf(fp, " ");
LinAlg_PrintScalar(fp, &tmp) ; }
fprintf(fp, " ") ; }
} else
} Message::Warning("Print of Dof out of TimeStep range [0,%d]",
else Message::Warning("Print of Dof out of TimeStep range [0,%d]" List_Nbr(DofData_P->Solutions) - 1);
, }
List_Nbr(DofData_P->Solutions)-1); else {
} LinAlg_GetScalarInVector(&tmp, &DofData_P->CurrentSolution->x, j);
else{ if(Flag_Binary) { LinAlg_WriteScalar(fp, &tmp); }
LinAlg_GetScalarInVector(&tmp, &DofData_P->CurrentSolution->x, j) ; else {
if(Flag_Binary){ LinAlg_PrintScalar(fp, &tmp);
LinAlg_WriteScalar(fp, &tmp) ; fprintf(fp, " ");
} }
else{ }
LinAlg_PrintScalar(fp, &tmp) ; }
fprintf(fp, " ") ; else
} Message::Warning(
} "Wrong number of Dof to Print (%d is out of [0,%d])", j,
} DofData_P->NbrDof - 1);
else Message::Warning("Wrong number of Dof to Print (%d is out of [0,%d }
])", fprintf(fp, "\n");
j, DofData_P->NbrDof-1); }
} else {
fprintf(fp, "\n") ; DofData_P = DofData_P0 + Operation_P->DefineSystemIndex;
} if(Flag_Binary) {
else{ LinAlg_WriteMatrix(fp, &DofData_P->A);
DofData_P = DofData_P0 + Operation_P->DefineSystemIndex ; LinAlg_WriteVector(fp, &DofData_P->b);
if(Flag_Binary){ }
LinAlg_WriteMatrix(fp, &DofData_P->A) ; else {
LinAlg_WriteVector(fp, &DofData_P->b) ;
}
else{
// use matlab format if available // use matlab format if available
DefineSystem_P = (struct DefineSystem*) DefineSystem_P = (struct DefineSystem *)List_Pointer(
List_Pointer(Resolution_P->DefineSystem, Operation_P->DefineSystemIn Resolution_P->DefineSystem, Operation_P->DefineSystemIndex);
dex) ;
std::string path(Name_Path), file("file_"); std::string path(Name_Path), file("file_");
std::string mat("mat_"), vec("vec_"), sol("sol_"); std::string mat("mat_"), vec("vec_"), sol("sol_");
std::string jac("jac_"), res("res_"), dx("dx_"); std::string jac("jac_"), res("res_"), dx("dx_");
std::string name(Operation_P->Case.Print.FileOut ? Operation_P->Case.P std::string name(Operation_P->Case.Print.FileOut ?
rint.FileOut : Operation_P->Case.Print.FileOut :
DefineSystem_P->Name); DefineSystem_P->Name);
if(DofData_P->Flag_Init[1] || DofData_P->Flag_Init[2] || if(DofData_P->Flag_Init[1] || DofData_P->Flag_Init[2] ||
DofData_P->Flag_Init[3] || DofData_P->Flag_Init[4] || DofData_P->Flag_Init[3] || DofData_P->Flag_Init[4] ||
DofData_P->Flag_Init[5] || DofData_P->Flag_Init[6] || DofData_P->Flag_Init[5] || DofData_P->Flag_Init[6] ||
DofData_P->Flag_Init[7]){ DofData_P->Flag_Init[7]) {
if(DofData_P->Flag_Init[1]){ if(DofData_P->Flag_Init[1]) {
std::string name1 = name + "1"; std::string name1 = name + "1";
LinAlg_PrintMatrix(fp, &DofData_P->M1, true, LinAlg_PrintMatrix(fp, &DofData_P->M1, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name).c_str()) ; (mat + name).c_str());
LinAlg_PrintVector(fp, &DofData_P->m1, true, LinAlg_PrintVector(fp, &DofData_P->m1, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[2]){ if(DofData_P->Flag_Init[2]) {
std::string name1 = name + "2"; std::string name1 = name + "2";
LinAlg_PrintMatrix(fp, &DofData_P->M2, true, LinAlg_PrintMatrix(fp, &DofData_P->M2, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m2, true, LinAlg_PrintVector(fp, &DofData_P->m2, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[3]){ if(DofData_P->Flag_Init[3]) {
std::string name1 = name + "3"; std::string name1 = name + "3";
LinAlg_PrintMatrix(fp, &DofData_P->M3, true, LinAlg_PrintMatrix(fp, &DofData_P->M3, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m3, true, LinAlg_PrintVector(fp, &DofData_P->m3, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[4]){ if(DofData_P->Flag_Init[4]) {
std::string name1 = name + "4"; std::string name1 = name + "4";
LinAlg_PrintMatrix(fp, &DofData_P->M4, true, LinAlg_PrintMatrix(fp, &DofData_P->M4, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m4, true, LinAlg_PrintVector(fp, &DofData_P->m4, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[5]){ if(DofData_P->Flag_Init[5]) {
std::string name1 = name + "5"; std::string name1 = name + "5";
LinAlg_PrintMatrix(fp, &DofData_P->M5, true, LinAlg_PrintMatrix(fp, &DofData_P->M5, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m5, true, LinAlg_PrintVector(fp, &DofData_P->m5, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[6]){ if(DofData_P->Flag_Init[6]) {
std::string name1 = name + "6"; std::string name1 = name + "6";
LinAlg_PrintMatrix(fp, &DofData_P->M6, true, LinAlg_PrintMatrix(fp, &DofData_P->M6, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m6, true, LinAlg_PrintVector(fp, &DofData_P->m6, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
if(DofData_P->Flag_Init[7]){ if(DofData_P->Flag_Init[7]) {
std::string name1 = name + "7"; std::string name1 = name + "7";
LinAlg_PrintMatrix(fp, &DofData_P->M7, true, LinAlg_PrintMatrix(fp, &DofData_P->M7, true,
(path + file + mat + name1 + ".m").c_str(), (path + file + mat + name1 + ".m").c_str(),
(mat + name1).c_str()) ; (mat + name1).c_str());
LinAlg_PrintVector(fp, &DofData_P->m7, true, LinAlg_PrintVector(fp, &DofData_P->m7, true,
(path + file + vec + name1 + ".m").c_str(), (path + file + vec + name1 + ".m").c_str(),
(vec + name1).c_str()) ; (vec + name1).c_str());
} }
} }
else{ else {
if(DofData_P->Flag_Init[0]){ if(DofData_P->Flag_Init[0]) {
LinAlg_PrintMatrix(fp, &DofData_P->A, true, LinAlg_PrintMatrix(fp, &DofData_P->A, true,
(path + file + mat + name + ".m").c_str(), (path + file + mat + name + ".m").c_str(),
(mat + name).c_str()) ; (mat + name).c_str());
LinAlg_PrintVector(fp, &DofData_P->b, true, LinAlg_PrintVector(fp, &DofData_P->b, true,
(path + file + vec + name + ".m").c_str(), (path + file + vec + name + ".m").c_str(),
(vec + name).c_str()) ; (vec + name).c_str());
} }
if(DofData_P->Flag_Init[0] == 2){ if(DofData_P->Flag_Init[0] == 2) {
LinAlg_PrintMatrix(fp, &DofData_P->Jac, true, LinAlg_PrintMatrix(fp, &DofData_P->Jac, true,
(path + file + jac + name + ".m").c_str(), (path + file + jac + name + ".m").c_str(),
(jac + name).c_str()) ; (jac + name).c_str());
LinAlg_PrintVector(fp, &DofData_P->res, true, LinAlg_PrintVector(fp, &DofData_P->res, true,
(path + file + res + name + ".m").c_str(), (path + file + res + name + ".m").c_str(),
(res + name).c_str()) ; (res + name).c_str());
LinAlg_PrintVector(fp, &DofData_P->dx, true, LinAlg_PrintVector(fp, &DofData_P->dx, true,
(path + file + dx + name + ".m").c_str(), (path + file + dx + name + ".m").c_str(),
(dx + name).c_str()) ; (dx + name).c_str());
} }
} }
if(DofData_P->CurrentSolution) if(DofData_P->CurrentSolution)
LinAlg_PrintVector(fp, &DofData_P->CurrentSolution->x, true, LinAlg_PrintVector(fp, &DofData_P->CurrentSolution->x, true,
(path + file + sol + name + ".m").c_str(), (path + file + sol + name + ".m").c_str(),
(sol + name).c_str()) ; (sol + name).c_str());
} }
} }
fflush(fp); fflush(fp);
if(Operation_P->Case.Print.FileOut){ if(Operation_P->Case.Print.FileOut) {
fclose(fp); fclose(fp);
fp = stdout ; fp = stdout;
} }
Flag_Binary = 0; Flag_Binary = 0;
break; break;
case OPERATION_DEBUG : case OPERATION_DEBUG:
Init_OperationOnSystem("Debug", Init_OperationOnSystem("Debug", Resolution_P, Operation_P, DofData_P0,
Resolution_P, Operation_P, DofData_P0, GeoData_P0, GeoData_P0, &DefineSystem_P, &DofData_P,
&DefineSystem_P, &DofData_P, Resolution2_P) ; Resolution2_P);
Operation_Debug(Operation_P, DofData_P); Operation_Debug(Operation_P, DofData_P);
break; break;
/* --> C h a n g e O f C o o r d i n a t e s */ /* --> C h a n g e O f C o o r d i n a t e s */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_CHANGEOFCOORDINATES : case OPERATION_CHANGEOFCOORDINATES:
if(Message::GetVerbosity() == 10) // +++ if(Message::GetVerbosity() == 10) // +++
Message::Info("ChangeOfCoordinates") ; Message::Info("ChangeOfCoordinates");
/* Geo_SetCurrentGeoData(Current.GeoData = GeoData_P0) ; */ /* Geo_SetCurrentGeoData(Current.GeoData = GeoData_P0) ; */
Operation_ChangeOfCoordinates Operation_ChangeOfCoordinates(Resolution_P, Operation_P, DofData_P0,
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ; GeoData_P0);
break ; break;
/* --> D e f o r m e M e s h */ /* --> D e f o r m e M e s h */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_DEFORMMESH : case OPERATION_DEFORMMESH: {
{ if(Operation_P->Case.DeformMesh.Name_MshFile == NULL)
if (Operation_P->Case.DeformMesh.Name_MshFile == NULL) Operation_P->Case.DeformMesh.Name_MshFile = Name_MshFile;
Operation_P->Case.DeformMesh.Name_MshFile = Name_MshFile ; Message::Info(
Message::Info("DeformMesh[%s, %s, '%s']", "DeformMesh[%s, %s, '%s']",
((struct DefineSystem *) ((struct DefineSystem *)List_Pointer(Resolution_P->DefineSystem,
List_Pointer(Resolution_P->DefineSystem, Operation_P->DefineSystemIndex))
Operation_P->DefineSystemIndex))->Name, ->Name,
Operation_P->Case.DeformMesh.Quantity, Operation_P->Case.DeformMesh.Quantity,
Operation_P->Case.DeformMesh.Name_MshFile) ; Operation_P->Case.DeformMesh.Name_MshFile);
int i; int i;
if ((i = List_ISearchSeq(GeoData_L, Operation_P->Case.DeformMesh.Name_Ms if((i = List_ISearchSeq(GeoData_L,
hFile, Operation_P->Case.DeformMesh.Name_MshFile,
fcmp_GeoData_Name)) < 0){ fcmp_GeoData_Name)) < 0) {
Message::Error("DeformMesh: Wrong NameOfMeshFile %s", Message::Error("DeformMesh: Wrong NameOfMeshFile %s",
Operation_P->Case.DeformMesh.Name_MshFile); Operation_P->Case.DeformMesh.Name_MshFile);
break; break;
}
Operation_P->Case.DeformMesh.GeoDataIndex = i ;
Operation_DeformMesh
(Resolution_P, Operation_P, DofData_P0, GeoData_P0) ;
} }
break; Operation_P->Case.DeformMesh.GeoDataIndex = i;
Operation_DeformMesh(Resolution_P, Operation_P, DofData_P0, GeoData_P0);
} break;
/* --> P o s t O p e r a t i o n */ /* --> P o s t O p e r a t i o n */
/* ------------------------------- */ /* ------------------------------- */
case OPERATION_POSTOPERATION : case OPERATION_POSTOPERATION:
#ifdef TIMER Message::Info("PostOperation");
{double tstart = MPI_Wtime();
#endif
Message::Info("PostOperation") ;
Operation_PostOperation(Resolution_P, DofData_P0, GeoData_P0, Operation_PostOperation(Resolution_P, DofData_P0, GeoData_P0,
Operation_P->Case.PostOperation.PostOperations); Operation_P->Case.PostOperation.PostOperations);
#ifdef TIMER break;
double timer = MPI_Wtime() - tstart;
printf("Proc %d, time spent in PostOperation %.16g\n", Message::GetCommRan
k(), timer);
}
#endif
break ;
/* --> D e l e t e F i l e */ /* --> D e l e t e F i l e */
/* ------------------------- */ /* ------------------------- */
case OPERATION_DELETEFILE : case OPERATION_DELETEFILE:
Message::Info("DeleteFile[%s]", Operation_P->Case.DeleteFile.FileName) ; Message::Info("DeleteFile[%s]", Operation_P->Case.DeleteFile.FileName);
RemoveFile(Operation_P->Case.DeleteFile.FileName); RemoveFile(Operation_P->Case.DeleteFile.FileName);
break ; break;
/* --> R e n a m e F i l e */ /* --> R e n a m e F i l e */
/* ------------------------- */ /* ------------------------- */
case OPERATION_RENAMEFILE : case OPERATION_RENAMEFILE:
Message::Info("RenameFile[%s, %s]", Operation_P->Case.RenameFile.OldFileNa Message::Info("RenameFile[%s, %s]",
me, Operation_P->Case.RenameFile.OldFileName,
Operation_P->Case.RenameFile.NewFileName) ; Operation_P->Case.RenameFile.NewFileName);
RenameFile(Operation_P->Case.RenameFile.OldFileName, RenameFile(Operation_P->Case.RenameFile.OldFileName,
Operation_P->Case.RenameFile.NewFileName); Operation_P->Case.RenameFile.NewFileName);
break ; break;
/* --> C r e a t e D i r */ /* --> C r e a t e D i r */
/* ------------------------ */ /* ------------------------ */
case OPERATION_CREATEDIR : case OPERATION_CREATEDIR:
Message::Info("CreateDir[%s]", Operation_P->Case.CreateDir.DirName) ; Message::Info("CreateDir[%s]", Operation_P->Case.CreateDir.DirName);
CreateDirs(Operation_P->Case.CreateDir.DirName); CreateDirs(Operation_P->Case.CreateDir.DirName);
break ; break;
/* --> T i m e L o o p A d a p t i v e */ /* --> T i m e L o o p A d a p t i v e */
/* ------------------------------------- */ /* ------------------------------------- */
case OPERATION_TIMELOOPADAPTIVE : case OPERATION_TIMELOOPADAPTIVE:
Message::Info("TimeLoopAdaptve ...") ; Message::Info("TimeLoopAdaptve ...");
Save_TypeTime = Current.TypeTime ; Save_TypeTime = Current.TypeTime;
Save_DTime = Current.DTime ; Save_DTime = Current.DTime;
Operation_TimeLoopAdaptive(Resolution_P, Operation_P, DofData_P0, GeoData_ Operation_TimeLoopAdaptive(Resolution_P, Operation_P, DofData_P0,
P0, GeoData_P0, &Flag_Break);
&Flag_Break) ; Current.TypeTime = Save_TypeTime;
Current.TypeTime = Save_TypeTime ; Current.DTime = Save_DTime;
Current.DTime = Save_DTime ;
break; break;
/* --> I t e r a t i v e L o o p N */ /* --> I t e r a t i v e L o o p N */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_ITERATIVELOOPN : case OPERATION_ITERATIVELOOPN:
Message::Info("IterativeLoopN ...") ; Message::Info("IterativeLoopN ...");
Save_Iteration = Current.Iteration ; Save_Iteration = Current.Iteration;
Operation_IterativeLoopN(Resolution_P, Operation_P, DofData_P0, GeoData_P0 Operation_IterativeLoopN(Resolution_P, Operation_P, DofData_P0,
, GeoData_P0, Resolution2_P, DofData2_P0,
Resolution2_P, DofData2_P0, &Flag_Break) ; &Flag_Break);
Current.Iteration = Save_Iteration ; Current.Iteration = Save_Iteration;
break; break;
/* --> S e t E x t r a p o l a t i o n O r d e r */ /* --> S e t E x t r a p o l a t i o n O r d e r */
/* --------------------------------------------- */ /* --------------------------------------------- */
case OPERATION_SETEXTRAPOLATIONORDER : case OPERATION_SETEXTRAPOLATIONORDER:
Message::Info("SetExtrapolationOrder [%d]...", Message::Info("SetExtrapolationOrder [%d]...",
Operation_P->Case.SetExtrapolationOrder.order) ; Operation_P->Case.SetExtrapolationOrder.order);
Flag_ExtrapolationOrder = Operation_P->Case.SetExtrapolationOrder.order; Flag_ExtrapolationOrder = Operation_P->Case.SetExtrapolationOrder.order;
if(Flag_ExtrapolationOrder < 0) Flag_ExtrapolationOrder = 0; if(Flag_ExtrapolationOrder < 0) Flag_ExtrapolationOrder = 0;
break; break;
/* --> T i m e L o o p R u n g e K u t t a */ /* --> T i m e L o o p R u n g e K u t t a */
/* ----------------------------------------- */ /* ----------------------------------------- */
case OPERATION_TIMELOOPRUNGEKUTTA : case OPERATION_TIMELOOPRUNGEKUTTA: {
{ Init_OperationOnSystem("TimeLoopRungeKutta", Resolution_P, Operation_P,
Init_OperationOnSystem("TimeLoopRungeKutta", DofData_P0, GeoData_P0, &DefineSystem_P,
Resolution_P, Operation_P, DofData_P0, GeoData_P0 &DofData_P, Resolution2_P);
, int numStepRK = List_Nbr(Operation_P->Case.TimeLoopRungeKutta.ButcherC);
&DefineSystem_P, &DofData_P, Resolution2_P) ; if(numStepRK != List_Nbr(Operation_P->Case.TimeLoopRungeKutta.ButcherB) ||
int numStepRK = List_Nbr(Operation_P->Case.TimeLoopRungeKutta.ButcherC); numStepRK * numStepRK !=
if(numStepRK != List_Nbr(Operation_P->Case.TimeLoopRungeKutta.ButcherB) List_Nbr(Operation_P->Case.TimeLoopRungeKutta.ButcherA)) {
|| Message::Error("Incompatible sizes of Butcher Tableaux");
numStepRK * numStepRK != List_Nbr(Operation_P->Case.TimeLoopRungeKutt break;
a.ButcherA)){ }
Message::Error("Incompatible sizes of Butcher Tableaux"); Current.Time = Operation_P->Case.TimeLoopRungeKutta.Time0;
gVector xn, rhs;
LinAlg_CreateVector(&xn, &DofData_P->Solver, Current.DofData->NbrDof);
LinAlg_CreateVector(&rhs, &DofData_P->Solver, Current.DofData->NbrDof);
std::vector<gVector> ki(numStepRK);
for(int i = 0; i < numStepRK; i++)
LinAlg_CreateVector(&ki[i], &DofData_P->Solver,
Current.DofData->NbrDof);
while(Current.Time <
Operation_P->Case.TimeLoopRungeKutta.TimeMax * 0.9999999) {
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount())
break; break;
}
Current.Time = Operation_P->Case.TimeLoopRungeKutta.Time0 ; double tn = Current.Time;
gVector xn, rhs; LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &xn);
LinAlg_CreateVector(&xn, &DofData_P->Solver, Current.DofData->NbrDof); Get_ValueOfExpressionByIndex(
LinAlg_CreateVector(&rhs, &DofData_P->Solver, Current.DofData->NbrDof); Operation_P->Case.TimeLoopRungeKutta.DTimeIndex, NULL, 0., 0., 0.,
std::vector<gVector> ki(numStepRK); &Value);
for(int i = 0; i < numStepRK; i++) Current.DTime = Value.Val[0];
LinAlg_CreateVector(&ki[i], &DofData_P->Solver, Current.DofData->NbrDo Current.TimeStep += 1.;
f); for(int i = 0; i < numStepRK; i++) {
double ci;
while (Current.Time < Operation_P->Case.TimeLoopRungeKutta.TimeMax * 0.9 List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherC, i, &ci);
999999) { Current.Time = tn + ci * Current.DTime;
if(Message::GetOnelabAction() == "stop" || Message::GetErrorCount()) b
reak;
double tn = Current.Time;
LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &xn);
Get_ValueOfExpressionByIndex(Operation_P->Case.TimeLoopRungeKutta.DTim
eIndex,
NULL, 0., 0., 0., &Value) ;
Current.DTime = Value.Val[0];
Current.TimeStep += 1.;
for(int i = 0; i < numStepRK; i++){
double ci;
List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherC, i, &ci);
Current.Time = tn + ci * Current.DTime;
LinAlg_CopyVector(&xn, &DofData_P->CurrentSolution->x);
// FIXME: warning, this assumes an explicit RK scheme!
for(int j = 0; j < i; j++){
double aij;
List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherA, i * numSt
epRK + j, &aij);
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &
ki[j], aij,
&DofData_P->CurrentSolution->x);
}
Current.TypeAssembly = ASSEMBLY_SEPARATE ;
Init_SystemData(DofData_P, Flag_Jac) ;
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 1);
LinAlg_ProdMatrixVector(&DofData_P->M1, &DofData_P->CurrentSolution-
>x, &rhs);
LinAlg_ProdVectorDouble(&rhs, -1., &rhs);
LinAlg_AddVectorProdVectorDouble(&rhs, &DofData_P->b, 1., &rhs)
;
LinAlg_ProdVectorDouble(&rhs, Current.DTime, &rhs);
LinAlg_Solve(&DofData_P->M2, &rhs, &DofData_P->Solver, &ki[i]) ;
}
// restore previous time step
LinAlg_CopyVector(&xn, &((struct Solution*)
List_Pointer(DofData_P->Solutions,
List_Nbr(DofData_P->Solutions)-2
))->x) ;
LinAlg_CopyVector(&xn, &DofData_P->CurrentSolution->x); LinAlg_CopyVector(&xn, &DofData_P->CurrentSolution->x);
for(int i = 0; i < numStepRK; i++){ // FIXME: warning, this assumes an explicit RK scheme!
double bi; for(int j = 0; j < i; j++) {
List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherB, i, &bi); double aij;
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &ki List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherA,
[i], bi, i * numStepRK + j, &aij);
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x,
&ki[j], aij,
&DofData_P->CurrentSolution->x); &DofData_P->CurrentSolution->x);
} }
Current.TypeAssembly = ASSEMBLY_SEPARATE;
Init_SystemData(DofData_P, Flag_Jac);
Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 1);
LinAlg_ProdMatrixVector(&DofData_P->M1,
&DofData_P->CurrentSolution->x, &rhs);
LinAlg_ProdVectorDouble(&rhs, -1., &rhs);
LinAlg_AddVectorProdVectorDouble(&rhs, &DofData_P->b, 1., &rhs);
LinAlg_ProdVectorDouble(&rhs, Current.DTime, &rhs);
LinAlg_Solve(&DofData_P->M2, &rhs, &DofData_P->Solver, &ki[i]);
}
// restore previous time step
LinAlg_CopyVector(
&xn, &((struct Solution *)List_Pointer(
DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 2))
->x);
LinAlg_CopyVector(&xn, &DofData_P->CurrentSolution->x);
for(int i = 0; i < numStepRK; i++) {
double bi;
List_Read(Operation_P->Case.TimeLoopRungeKutta.ButcherB, i, &bi);
LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x,
&ki[i], bi,
&DofData_P->CurrentSolution->x);
}
Current.Time = tn + Current.DTime; Current.Time = tn + Current.DTime;
if(Flag_Break){ Flag_Break = 0; break; } if(Flag_Break) {
Flag_Break = 0;
break;
} }
} }
break ; } break;
/* case OPERATION_BREAK: Flag_Break = 1; break;
case OPERATION_GETSIZE:
{
Init_OperationOnSystem("GetSize",
Resolution_P, Operation_P, DofData_P0, GeoData_P0
,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
Cal_ZeroValue(&Value);
Value.Type = SCALAR;
Value.Val[0] = DofData_P->NbrDof;
Cal_StoreInVariable(&Value, Operation_P->Case.GetSize.VariableName);
}
break;
*/
case OPERATION_BREAK :
Flag_Break = 1;
break ;
case OPERATION_EXIT : case OPERATION_EXIT: Message::Exit(0); break;
Message::Exit(0);
break ;
case OPERATION_SLEEP : case OPERATION_SLEEP:
Get_ValueOfExpressionByIndex(Operation_P->Case.Sleep.ExpressionIndex, Get_ValueOfExpressionByIndex(Operation_P->Case.Sleep.ExpressionIndex,
NULL, 0., 0., 0., &Value) ; NULL, 0., 0., 0., &Value);
Message::Info("Sleeping for %g seconds", Value.Val[0]); Message::Info("Sleeping for %g seconds", Value.Val[0]);
SleepSeconds(Value.Val[0]); SleepSeconds(Value.Val[0]);
break ; break;
/* --> P a r a l l e l C o m p u t i n g */ /* --> P a r a l l e l C o m p u t i n g */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_SETCOMMSELF : case OPERATION_SETCOMMSELF: LinAlg_SetCommSelf(); break;
LinAlg_SetCommSelf();
break ;
case OPERATION_SETCOMMWORLD :
LinAlg_SetCommWorld();
break ;
case OPERATION_BARRIER : case OPERATION_SETCOMMWORLD: LinAlg_SetCommWorld(); break;
case OPERATION_BARRIER:
#if defined(HAVE_PETSC) #if defined(HAVE_PETSC)
Message::Info("Barrier: waiting"); Message::Info("Barrier: waiting");
MPI_Barrier(PETSC_COMM_WORLD); MPI_Barrier(PETSC_COMM_WORLD);
Message::Info("Barrier: let's continue"); Message::Info("Barrier: let's continue");
#endif #endif
break ; break;
/* --> O p t i m i z e r */ /* --> O p t i m i z e r */
/* ------------------------------------------ */ /* ------------------------------------------ */
case OPERATION_OPTIMIZER_INITIALIZE : case OPERATION_OPTIMIZER_INITIALIZE:
Operation_OptimizerInitialize(Operation_P); Operation_OptimizerInitialize(Operation_P);
break ; break;
case OPERATION_OPTIMIZER_UPDATE : case OPERATION_OPTIMIZER_UPDATE:
Operation_OptimizerUpdate(Operation_P); Operation_OptimizerUpdate(Operation_P);
break ; break;
case OPERATION_OPTIMIZER_FINALIZE : case OPERATION_OPTIMIZER_FINALIZE:
Operation_OptimizerFinalize(Operation_P); Operation_OptimizerFinalize(Operation_P);
break ; break;
/* --> O t h e r */ /* --> O t h e r */
/* ------------------------------------------ */ /* ------------------------------------------ */
default : default: Message::Warning("Operation: ? ? ?"); break;
Message::Warning("Operation: ? ? ?") ;
break ;
} }
} }
Message::Barrier(); Message::Barrier();
} }
 End of changes. 555 change blocks. 
2863 lines changed or deleted 2459 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)