"Fossies" - the Fresh Open Source Software Archive  

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

Cal_Quantity.cpp  (getdp-3.4.0-source.tgz):Cal_Quantity.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include <map> #include <map>
#include "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "DofData.h" #include "DofData.h"
skipping to change at line 28 skipping to change at line 28
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
#include "GeoData.h" #include "GeoData.h"
#include "Get_Geometry.h" #include "Get_Geometry.h"
#include "Pos_FemInterpolation.h" #include "Pos_FemInterpolation.h"
#include "Pos_Search.h" #include "Pos_Search.h"
#include "Get_FunctionValue.h" #include "Get_FunctionValue.h"
#include "Pos_Format.h" #include "Pos_Format.h"
#endif #endif
extern struct Problem Problem_S ; extern struct Problem Problem_S;
extern struct CurrentData Current ; extern struct CurrentData Current;
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
extern int TreatmentStatus ; extern int TreatmentStatus;
#else #else
const int TreatmentStatus = 0; // FIXME const int TreatmentStatus = 0; // FIXME
#endif #endif
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G e t _ V a l u e O f E x p r e s s i o n */ /* G e t _ V a l u e O f E x p r e s s i o n */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Get_ValueOfExpression(struct Expression * Expression_P, void Get_ValueOfExpression(struct Expression *Expression_P,
struct QuantityStorage * QuantityStorage_P0, struct QuantityStorage *QuantityStorage_P0, double u,
double u, double v, double w, double v, double w, struct Value *Value,
struct Value * Value, int NbrArguments, char *CallingExpressionName)
int NbrArguments,
char *CallingExpressionName)
{ {
static char *Flag_WarningUndefined = NULL ; static char *Flag_WarningUndefined = NULL;
switch (Expression_P->Type) { switch(Expression_P->Type) {
case CONSTANT:
case CONSTANT : if(Current.NbrHar == 1) { Value->Val[0] = Expression_P->Case.Constant; }
if (Current.NbrHar == 1) {
Value->Val[0] = Expression_P->Case.Constant ;
}
else { else {
for (int k = 0 ; k < Current.NbrHar ; k += 2) { for(int k = 0; k < Current.NbrHar; k += 2) {
Value->Val[MAX_DIM* k ] = Expression_P->Case.Constant ; Value->Val[MAX_DIM * k] = Expression_P->Case.Constant;
Value->Val[MAX_DIM*(k+1)] = 0. ; Value->Val[MAX_DIM * (k + 1)] = 0.;
} }
} }
Value->Type = SCALAR ; Value->Type = SCALAR;
break ; break;
case WHOLEQUANTITY : case WHOLEQUANTITY:
Cal_WholeQuantity(Current.Element, QuantityStorage_P0, Cal_WholeQuantity(
Expression_P->Case.WholeQuantity, Current.Element, QuantityStorage_P0, Expression_P->Case.WholeQuantity, u,
u,v,w, -1, 0, Value, NbrArguments, v, w, -1, 0, Value, NbrArguments,
CallingExpressionName ? CallingExpressionName ? CallingExpressionName : Expression_P->Name);
CallingExpressionName : Expression_P->Name) ; break;
break ;
case PIECEWISEFUNCTION :
struct Expression * NextExpression_P;
struct ExpressionPerRegion * ExpressionPerRegion_P ;
if (Current.Region == Expression_P->Case.PieceWiseFunction.NumLastRegion) { case PIECEWISEFUNCTION:
NextExpression_P = Expression_P->Case.PieceWiseFunction.ExpressionForLastR struct Expression *NextExpression_P;
egion; struct ExpressionPerRegion *ExpressionPerRegion_P;
if(Current.Region == Expression_P->Case.PieceWiseFunction.NumLastRegion) {
NextExpression_P =
Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion;
} }
else { else {
if ((ExpressionPerRegion_P = (struct ExpressionPerRegion*) if((ExpressionPerRegion_P = (struct ExpressionPerRegion *)List_PQuery(
List_PQuery(Expression_P->Case.PieceWiseFunction.ExpressionPerRegion, Expression_P->Case.PieceWiseFunction.ExpressionPerRegion,
&Current.Region, fcmp_int))) { &Current.Region, fcmp_int))) {
Expression_P->Case.PieceWiseFunction.NumLastRegion = Current.Region ; Expression_P->Case.PieceWiseFunction.NumLastRegion = Current.Region;
NextExpression_P = NextExpression_P =
Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion = Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion =
(struct Expression*)List_Pointer(Problem_S.Expression, (struct Expression *)List_Pointer(
ExpressionPerRegion_P->ExpressionInde Problem_S.Expression, ExpressionPerRegion_P->ExpressionIndex);
x) ;
} }
else if (Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >= 0 else if(Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >=
) { 0) {
NextExpression_P = NextExpression_P = (struct Expression *)List_Pointer(
(struct Expression*) Problem_S.Expression,
List_Pointer(Problem_S.Expression, Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default);
Expression_P->Case.PieceWiseFunction.ExpressionIndex_Defa
ult) ;
} }
else { else {
NextExpression_P = NULL; NextExpression_P = NULL;
if(Current.Region == NO_REGION) if(Current.Region == NO_REGION)
Message::Error("Function '%s' undefined in expressions without support Message::Error(
", "Function '%s' undefined in expressions without support",
Expression_P->Name); Expression_P->Name);
else else
Message::Error("Function '%s' undefined in Region %d", Message::Error("Function '%s' undefined in Region %d",
Expression_P->Name, Current.Region); Expression_P->Name, Current.Region);
} }
} }
Get_ValueOfExpression Get_ValueOfExpression(NextExpression_P, QuantityStorage_P0, u, v, w, Value,
(NextExpression_P, NbrArguments, Expression_P->Name);
QuantityStorage_P0, u, v, w, Value, NbrArguments, Expression_P->Name) ; break;
break ;
case PIECEWISEFUNCTION2:
case PIECEWISEFUNCTION2 : struct ExpressionPerRegion2 *ExpressionPerRegion2_P;
// struct Expression * NextExpression_P;
struct ExpressionPerRegion2 * ExpressionPerRegion2_P ;
int twoRegions[2]; int twoRegions[2];
if (Current.Region == Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] if(Current.Region ==
&& Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] &&
Current.SubRegion == Expression_P->Case.PieceWiseFunction2.NumLastRegion Current.SubRegion ==
[1]) { Expression_P->Case.PieceWiseFunction2.NumLastRegion[1]) {
NextExpression_P = Expression_P->Case.PieceWiseFunction2.ExpressionForLast NextExpression_P =
Region; Expression_P->Case.PieceWiseFunction2.ExpressionForLastRegion;
} }
else { else {
twoRegions[0] = Current.Region; twoRegions[0] = Current.Region;
twoRegions[1] = Current.SubRegion; twoRegions[1] = Current.SubRegion;
if ((ExpressionPerRegion2_P = (struct ExpressionPerRegion2*) if((ExpressionPerRegion2_P = (struct ExpressionPerRegion2 *)List_PQuery(
List_PQuery(Expression_P->Case.PieceWiseFunction2.ExpressionPerRegion Expression_P->Case.PieceWiseFunction2.ExpressionPerRegion,
, twoRegions, fcmp_Integer2))) {
twoRegions, fcmp_Integer2))) { Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] = Current.Region;
Expression_P->Case.PieceWiseFunction2.NumLastRegion[0] = Current.Region Expression_P->Case.PieceWiseFunction2.NumLastRegion[1] =
; Current.SubRegion;
Expression_P->Case.PieceWiseFunction2.NumLastRegion[1] = Current.SubRegi
on ;
NextExpression_P = NextExpression_P =
Expression_P->Case.PieceWiseFunction2.ExpressionForLastRegion = Expression_P->Case.PieceWiseFunction2.ExpressionForLastRegion =
(struct Expression*)List_Pointer(Problem_S.Expression, (struct Expression *)List_Pointer(
ExpressionPerRegion2_P->ExpressionInd Problem_S.Expression, ExpressionPerRegion2_P->ExpressionIndex);
ex) ;
} }
else if (Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >= 0 else if(Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default >=
) { 0) {
NextExpression_P = NextExpression_P = (struct Expression *)List_Pointer(
(struct Expression*) Problem_S.Expression,
List_Pointer(Problem_S.Expression, Expression_P->Case.PieceWiseFunction.ExpressionIndex_Default);
Expression_P->Case.PieceWiseFunction.ExpressionIndex_Defa
ult) ;
} }
else { else {
NextExpression_P = NULL; NextExpression_P = NULL;
if(Current.Region == NO_REGION) if(Current.Region == NO_REGION)
Message::Error("Function '%s' undefined in expressions without support Message::Error(
", "Function '%s' undefined in expressions without support",
Expression_P->Name); Expression_P->Name);
else else
Message::Error("Function '%s' undefined in [ Region %d, SubRegion %d ] Message::Error(
", "Function '%s' undefined in [ Region %d, SubRegion %d ]",
Expression_P->Name, Current.Region, Current.SubRegion); Expression_P->Name, Current.Region, Current.SubRegion);
} }
} }
Get_ValueOfExpression Get_ValueOfExpression(NextExpression_P, QuantityStorage_P0, u, v, w, Value,
(NextExpression_P, NbrArguments, Expression_P->Name);
QuantityStorage_P0, u, v, w, Value, NbrArguments, Expression_P->Name) ; break;
break ;
case UNDEFINED_EXP:
case UNDEFINED_EXP : if(!Flag_WarningUndefined ||
if(!Flag_WarningUndefined || strcmp(Flag_WarningUndefined, Expression_P->Nam strcmp(Flag_WarningUndefined, Expression_P->Name)) {
e)){ Message::Warning("Undefined expression '%s' (assuming zero)",
Message::Warning("Undefined expression '%s' (assuming zero)", Expression_P Expression_P->Name);
->Name) ;
Flag_WarningUndefined = Expression_P->Name; Flag_WarningUndefined = Expression_P->Name;
} }
Cal_ZeroValue(Value); Cal_ZeroValue(Value);
Value->Type = SCALAR ; Value->Type = SCALAR;
break; break;
default : default:
Message::Error("Unknown type (%d) of Expression (%s)", Message::Error("Unknown type (%d) of Expression (%s)", Expression_P->Type,
Expression_P->Type, Expression_P->Name) ; Expression_P->Name);
break; break;
} }
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* G e t _ V a l u e O f E x p r e s s i o n B y I n d e x */ /* G e t _ V a l u e O f E x p r e s s i o n B y I n d e x */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Get_ValueOfExpressionByIndex(int Index_Expression, void Get_ValueOfExpressionByIndex(int Index_Expression,
struct QuantityStorage * QuantityStorage_P0, struct QuantityStorage *QuantityStorage_P0,
double u, double v, double w, double u, double v, double w,
struct Value * Value) { struct Value *Value)
Get_ValueOfExpression {
((struct Expression *)List_Pointer(Problem_S.Expression, Index_Expression), Get_ValueOfExpression(
QuantityStorage_P0, u, v, w, Value) ; (struct Expression *)List_Pointer(Problem_S.Expression, Index_Expression),
QuantityStorage_P0, u, v, w, Value);
} }
bool Is_ExpressionConstant(struct Expression * Expression_P) bool Is_ExpressionConstant(struct Expression *Expression_P)
{ {
if(Expression_P->Type == CONSTANT) return true; if(Expression_P->Type == CONSTANT) return true;
if(Expression_P->Type == WHOLEQUANTITY){ if(Expression_P->Type == WHOLEQUANTITY) {
for(int i = 0; i < List_Nbr(Expression_P->Case.WholeQuantity); i++){ for(int i = 0; i < List_Nbr(Expression_P->Case.WholeQuantity); i++) {
struct WholeQuantity *WholeQuantity_P = (struct WholeQuantity*) struct WholeQuantity *WholeQuantity_P =
List_Pointer(Expression_P->Case.WholeQuantity, i); (struct WholeQuantity *)List_Pointer(Expression_P->Case.WholeQuantity,
if(WholeQuantity_P->Type != WQ_CONSTANT) i);
return false; if(WholeQuantity_P->Type != WQ_CONSTANT) return false;
} }
return true; return true;
} }
return false; return false;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C a l _ S o l i d A n g l e */ /* C a l _ S o l i d A n g l e */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Cal_SolidAngle(int Source, struct Element *Element, void Cal_SolidAngle(int Source, struct Element *Element,
struct QuantityStorage *QuantityStorage, struct QuantityStorage *QuantityStorage, int Nbr_Dof,
int Nbr_Dof, int Index, int Index, struct Value **Stack)
struct Value **Stack)
{ {
#if !defined(HAVE_KERNEL) #if !defined(HAVE_KERNEL)
Message::Error("Cal_SolidAngle requires Kernel"); Message::Error("Cal_SolidAngle requires Kernel");
#else #else
struct Element *Elt ; struct Element *Elt;
struct Geo_Element *GeoElement ; struct Geo_Element *GeoElement;
struct Geo_Node *GeoNode1, *GeoNode2, *GeoNode3 ; struct Geo_Node *GeoNode1, *GeoNode2, *GeoNode3;
double X, Y, Atan ; double X, Y, Atan;
int In, TypEnt, NumNode, NbrElements, *NumElements ; int In, TypEnt, NumNode, NbrElements, *NumElements;
int i, j ; int i, j;
if(Nbr_Dof != QuantityStorage->NbrElementaryBasisFunction){ if(Nbr_Dof != QuantityStorage->NbrElementaryBasisFunction) {
Message::Error("Uncompatible Quantity (%s) in SolidAngle computation", Message::Error("Uncompatible Quantity (%s) in SolidAngle computation",
QuantityStorage->DefineQuantity->Name); QuantityStorage->DefineQuantity->Name);
return; return;
} }
if(Source){ if(Source) {
Elt = Element->ElementSource ; Elt = Element->ElementSource;
In = Current.SourceIntegrationSupportIndex ; In = Current.SourceIntegrationSupportIndex;
} }
else{ else {
Elt = Element ; Elt = Element;
In = Current.IntegrationSupportIndex ; In = Current.IntegrationSupportIndex;
} }
if (Elt->NumLastElementForSolidAngle == Elt->Num) { if(Elt->NumLastElementForSolidAngle == Elt->Num) {
for(i=0 ; i<Nbr_Dof ; i++) { for(i = 0; i < Nbr_Dof; i++) {
Stack[i][Index].Type = SCALAR ; Stack[i][Index].Type = SCALAR;
Stack[i][Index].Val[0] = Elt->angle[i] ; Stack[i][Index].Val[0] = Elt->angle[i];
} }
return; return;
} }
for(i=0 ; i<Nbr_Dof ; i++){ for(i = 0; i < Nbr_Dof; i++) {
Stack[i][Index].Type = SCALAR;
Stack[i][Index].Type = SCALAR ; TypEnt = ((struct Group *)List_Pointer(
Problem_S.Group,
QuantityStorage->BasisFunction[i].BasisFunction->EntityIndex))
->FunctionType;
TypEnt = ((struct Group*) if(TypEnt != NODESOF) {
List_Pointer(Problem_S.Group, if(Elt->Type == LINE) { Stack[i][Index].Val[0] = M_PI; }
QuantityStorage->BasisFunction[i]. else {
BasisFunction->EntityIndex))->FunctionType ; Stack[i][Index].Val[0] = 2. * M_PI;
if(TypEnt != NODESOF){
if(Elt->Type == LINE){
Stack[i][Index].Val[0] = M_PI ;
}
else{
Stack[i][Index].Val[0] = 2.*M_PI ;
} }
} }
else{ else {
NumNode =
NumNode = Elt->GeoElement-> Elt->GeoElement
NumNodes[QuantityStorage->BasisFunction[i].NumEntityInElement] ; ->NumNodes[QuantityStorage->BasisFunction[i].NumEntityInElement];
Geo_CreateNodesXElements(NumNode, In, &NbrElements, &NumElements) ; Geo_CreateNodesXElements(NumNode, In, &NbrElements, &NumElements);
if(NbrElements != 2){ if(NbrElements != 2) {
Message::Error("SolidAngle not done for incidence != 2 (%d)", NbrElement Message::Error("SolidAngle not done for incidence != 2 (%d)",
s); NbrElements);
return; return;
} }
GeoNode2 = Geo_GetGeoNodeOfNum(NumNode) ; GeoNode2 = Geo_GetGeoNodeOfNum(NumNode);
GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[0])) ; GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[0]));
if(GeoElement->Type != LINE){ if(GeoElement->Type != LINE) {
Message::Error("SolidAngle not done for Elements other than LINE"); Message::Error("SolidAngle not done for Elements other than LINE");
return; return;
} }
if(NumElements[0]>0){ if(NumElements[0] > 0) {
GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]) ; GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]);
GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[1])) ; GeoElement = Geo_GetGeoElementOfNum(abs(NumElements[1]));
GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]) ; GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]);
} }
else{ else {
GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]) ; GeoNode3 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[1]);
GeoElement = Geo_GetGeoElementOfNum(NumElements[1]) ; GeoElement = Geo_GetGeoElementOfNum(NumElements[1]);
GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]) ; GeoNode1 = Geo_GetGeoNodeOfNum(GeoElement->NumNodes[0]);
} }
Y = (GeoNode1->y - GeoNode2->y) * (GeoNode3->x - GeoNode2->x) - Y = (GeoNode1->y - GeoNode2->y) * (GeoNode3->x - GeoNode2->x) -
(GeoNode1->x - GeoNode2->x) * (GeoNode3->y - GeoNode2->y) ; (GeoNode1->x - GeoNode2->x) * (GeoNode3->y - GeoNode2->y);
X = (GeoNode1->x - GeoNode2->x) * (GeoNode3->x - GeoNode2->x) + X = (GeoNode1->x - GeoNode2->x) * (GeoNode3->x - GeoNode2->x) +
(GeoNode1->y - GeoNode2->y) * (GeoNode3->y - GeoNode2->y) ; (GeoNode1->y - GeoNode2->y) * (GeoNode3->y - GeoNode2->y);
Atan = atan2(Y,X) ; Atan = atan2(Y, X);
Stack[i][Index].Val[0] = (Atan >= 0) ? Atan : (Atan+2.*M_PI) ; Stack[i][Index].Val[0] = (Atan >= 0) ? Atan : (Atan + 2. * M_PI);
if(Message::GetVerbosity() > 5){ if(Message::GetVerbosity() > 5) {
printf("Solid angle=%g node=%d, region=%s, elms=", printf("Solid angle=%g node=%d, region=%s, elms=",
Stack[i][Index].Val[0] * 180/M_PI, NumNode, Stack[i][Index].Val[0] * 180 / M_PI, NumNode,
((struct Group*)List_Pointer(Problem_S.Group, In))->Name); ((struct Group *)List_Pointer(Problem_S.Group, In))->Name);
for(j=0 ; j<NbrElements ; j++) printf("%d ", NumElements[j]); for(j = 0; j < NbrElements; j++) printf("%d ", NumElements[j]);
printf("\n"); printf("\n");
} }
} }
} }
if (Elt->NumLastElementForSolidAngle != Elt->Num) { if(Elt->NumLastElementForSolidAngle != Elt->Num) {
Elt->NumLastElementForSolidAngle = Elt->Num ; Elt->NumLastElementForSolidAngle = Elt->Num;
for(i=0 ; i<Nbr_Dof ; i++) Elt->angle[i] = Stack[i][Index].Val[0] ; for(i = 0; i < Nbr_Dof; i++) Elt->angle[i] = Stack[i][Index].Val[0];
} }
#endif #endif
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C a l _ W h o l e Q u a n t i t y */ /* C a l _ W h o l e Q u a n t i t y */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
#define CAST3V void(*)(struct Value*, struct Value*, struct Value*) #define CAST3V void (*)(struct Value *, struct Value *, struct Value *)
#define CAST1V void(*)(struct Value*) #define CAST1V void (*)(struct Value *)
#define CASTF2V void(*)(struct Function*, struct Value*, struct Value*) #define CASTF2V void (*)(struct Function *, struct Value *, struct Value *)
// There can be at max one "Dof{op qty}" per WholeQuantity, but as // There can be at max one "Dof{op qty}" per WholeQuantity, but as
// many {op qty} as you want. // many {op qty} as you want.
static std::map<int, struct Value> ValueSaved ; static std::map<int, struct Value> ValueSaved;
static std::map<std::string, struct Value> NamedValueSaved ; static std::map<std::string, struct Value> NamedValueSaved;
void Cal_WholeQuantity(struct Element * Element, void Cal_WholeQuantity(struct Element *Element,
struct QuantityStorage * QuantityStorage_P0, struct QuantityStorage *QuantityStorage_P0,
List_T * WholeQuantity_L, List_T *WholeQuantity_L, double u, double v, double w,
double u, double v, double w, int DofIndexInWholeQuantity, int Nbr_Dof,
int DofIndexInWholeQuantity, struct Value DofValue[], int NbrArguments,
int Nbr_Dof, struct Value DofValue[], char *ExpressionName)
int NbrArguments, char *ExpressionName) { {
static int Flag_WarningMissSolForDt = 0;
static int Flag_WarningMissSolForDt = 0 ; static int Flag_WarningMissSolForTime_ntime = 0;
static int Flag_WarningMissSolForTime_ntime = 0 ; static int Flag_InfoForTime_ntime = 0;
static int Flag_InfoForTime_ntime = 0 ;
int i_WQ, j, k, Flag_True, Index, DofIndex, Multi[MAX_STACK_SIZE];
int i_WQ, j, k, Flag_True, Index, DofIndex, Multi[MAX_STACK_SIZE] ; int Save_NbrHar, Save_Region, Type_Dimension, ntime;
int Save_NbrHar, Save_Region, Type_Dimension, ntime ; double Save_Time, Save_TimeImag, Save_TimeStep, X, Y, Z, Order;
double Save_Time, Save_TimeImag, Save_TimeStep, X, Y, Z, Order ; double Save_x, Save_y, Save_z;
double Save_x, Save_y, Save_z ;
struct WholeQuantity *WholeQuantity_P0, *WholeQuantity_P;
struct WholeQuantity *WholeQuantity_P0, *WholeQuantity_P ; struct DofData *Save_DofData;
struct DofData *Save_DofData ; struct Solution *Solution_P0, *Solution_PN;
struct Solution *Solution_P0, *Solution_PN ;
struct Element* Save_CurrentElement ; struct Element *Save_CurrentElement;
// we could make this dynamic (with std::vector) to reduce stack usage, but // we could make this dynamic (with std::vector) to reduce stack usage, but
// the performance hit is important // the performance hit is important
// ==> forcing a reduced size of stack for MH case... // ==> forcing a reduced size of stack for MH case...
// MAX_STACK_SIZE0 = 8 by default, 2 for MH // MAX_STACK_SIZE0 = 8 by default, 2 for MH
// segmentation violation and out of memory with high number of harmonics // segmentation violation and out of memory with high number of harmonics
// MAX_STACK_SIZE at least MAX_HAR_SIZE if harmonic function in formulation te // MAX_STACK_SIZE at least MAX_HAR_SIZE if harmonic function in formulation
rm // term
struct Value Stack[MAX_STACK_SIZE0][MAX_STACK_SIZE] ; struct Value Stack[MAX_STACK_SIZE0][MAX_STACK_SIZE];
WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ; WholeQuantity_P0 = (struct WholeQuantity *)List_Pointer(WholeQuantity_L, 0);
Index = 0 ; Index = 0;
DofIndex = -1 ; DofIndex = -1;
for (i_WQ = 0 ; i_WQ < List_Nbr(WholeQuantity_L) ; i_WQ++) { for(i_WQ = 0; i_WQ < List_Nbr(WholeQuantity_L); i_WQ++) {
if(Index >= MAX_STACK_SIZE) {
if(Index >= MAX_STACK_SIZE){
Message::Error("Stack size exceeded (%d>%d)", Index, MAX_STACK_SIZE); Message::Error("Stack size exceeded (%d>%d)", Index, MAX_STACK_SIZE);
return; return;
} }
WholeQuantity_P = WholeQuantity_P0 + i_WQ ; WholeQuantity_P = WholeQuantity_P0 + i_WQ;
switch (WholeQuantity_P->Type) {
case WQ_OPERATORANDQUANTITY : /* {op qty} Dof{op qty} BF{op qty} */ switch(WholeQuantity_P->Type) {
Save_Region = Current.Region ; case WQ_OPERATORANDQUANTITY: /* {op qty} Dof{op qty} BF{op qty} */
Save_CurrentElement = Current.Element ; Save_Region = Current.Region;
if (i_WQ != DofIndexInWholeQuantity){ /* Attention!!! || TreatmentStatus = Save_CurrentElement = Current.Element;
= STATUS_POS){ */ if(i_WQ != DofIndexInWholeQuantity) {
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
Pos_FemInterpolation Pos_FemInterpolation(
(Element, Element, QuantityStorage_P0,
QuantityStorage_P0, QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index, WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity, WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0, u, v,
WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1);
u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
#else #else
Message::Error("TODO Post_FemInterpolation"); Message::Error("TODO Post_FemInterpolation");
#endif #endif
Multi[Index] = 0 ; Multi[Index] = 0;
} }
else { else {
DofIndex = Index ; DofIndex = Index;
} }
Index++ ; Index++;
Current.Element = Save_CurrentElement ; Current.Element = Save_CurrentElement;
Current.Region = Save_Region ; Current.Region = Save_Region;
break ; break;
case WQ_ORDER : /* Order[{qty}] */ case WQ_ORDER: /* Order[{qty}] */
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
Order = Cal_InterpolationOrder Order = Cal_InterpolationOrder(
(Element, QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity Element,
.Index) ; QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index);
#else #else
Message::Error("TODO Cal_InterpolationOrder"); Message::Error("TODO Cal_InterpolationOrder");
#endif #endif
for (k = 0 ; k < Current.NbrHar ; k += 2) { for(k = 0; k < Current.NbrHar; k += 2) {
Stack[0][Index].Val[MAX_DIM* k ] = Order ; Stack[0][Index].Val[MAX_DIM * k] = Order;
Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ; Stack[0][Index].Val[MAX_DIM * (k + 1)] = 0.;
} }
Stack[0][Index].Type = SCALAR ; Stack[0][Index].Type = SCALAR;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_OPERATORANDQUANTITYEVAL : case WQ_OPERATORANDQUANTITYEVAL:
Save_Region = Current.Region ; Save_Region = Current.Region;
Save_CurrentElement = Current.Element ; Save_CurrentElement = Current.Element;
/* {op qty}[x,y,z], {op qty}[x,y,z,dimension] /* {op qty}[x,y,z], {op qty}[x,y,z,dimension]
or {op qty}[Vector[x,y,x],dimension] or {op qty}[Vector[x,y,x],dimension]
or {op qty}[ntime] */ or {op qty}[ntime] */
if (i_WQ != DofIndexInWholeQuantity || TreatmentStatus == STATUS_POS){ if(i_WQ != DofIndexInWholeQuantity || TreatmentStatus == STATUS_POS) {
j = WholeQuantity_P->Case.OperatorAndQuantity.NbrArguments; j = WholeQuantity_P->Case.OperatorAndQuantity.NbrArguments;
if (j == 2 || j == 3 || j == 4) { if(j == 2 || j == 3 || j == 4) {
if (j == 3 || j == 4) { if(j == 3 || j == 4) {
Index -= j ; Index -= j;
X = Stack[0][Index ].Val[0] ; X = Stack[0][Index].Val[0];
Y = Stack[0][Index+1].Val[0] ; Y = Stack[0][Index + 1].Val[0];
Z = Stack[0][Index+2].Val[0] ; Z = Stack[0][Index + 2].Val[0];
if(j == 4) if(j == 4)
Type_Dimension = (int)Stack[0][Index+3].Val[0] ; Type_Dimension = (int)Stack[0][Index + 3].Val[0];
else else
Type_Dimension = -1 ; Type_Dimension = -1;
} }
else { /* j==2 */ else { /* j==2 */
Index -= j ; Index -= j;
X = Stack[0][Index ].Val[0] ; X = Stack[0][Index].Val[0];
Y = Stack[0][Index ].Val[1] ; Y = Stack[0][Index].Val[1];
Z = Stack[0][Index ].Val[2] ; Z = Stack[0][Index].Val[2];
Type_Dimension = (int)Stack[0][Index+1].Val[0] ; Type_Dimension = (int)Stack[0][Index + 1].Val[0];
} }
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
Pos_FemInterpolation Pos_FemInterpolation(
(Element, Element, QuantityStorage_P0,
QuantityStorage_P0, QuantityStorage_P0 +
QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Inde WholeQuantity_P->Case.OperatorAndQuantity.Index,
x, WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity, WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator,
WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, Type_Dimens Type_Dimension, 1, u, v, w, X, Y, Z, Stack[0][Index].Val,
ion, 1, &Stack[0][Index].Type, 1);
u, v, w, X, Y, Z, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
#else #else
Message::Error("TODO Post_FemInterpolation"); Message::Error("TODO Post_FemInterpolation");
#endif #endif
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
else if (j == 1) { else if(j == 1) {
Index -= j ; Index -= j;
ntime = (int)Stack[0][Index].Val[0] ; ntime = (int)Stack[0][Index].Val[0];
for (k = 0 ; k < Current.NbrSystem ; k++){ for(k = 0; k < Current.NbrSystem; k++) {
if(!List_Nbr((Current.DofData_P0+k)->Solutions)) continue; if(!List_Nbr((Current.DofData_P0 + k)->Solutions)) continue;
Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)- Solution_P0 = (struct Solution *)List_Pointer(
>Solutions, 0); (Current.DofData_P0 + k)->Solutions, 0);
if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime) if(((Current.DofData_P0 + k)->CurrentSolution - Solution_P0) >=
{ ntime) {
((Current.DofData_P0+k)->CurrentSolution) -= ntime ; ((Current.DofData_P0 + k)->CurrentSolution) -= ntime;
if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Sol if(Flag_InfoForTime_ntime !=
utions)) { List_Nbr((Current.DofData_P0 + k)->Solutions)) {
Message::Debug("Accessing solution from %d time steps ago", ntim Message::Debug("Accessing solution from %d time steps ago",
e); ntime);
Message::Debug(" -> System %d/%d: TimeStep = %d, Time = %g + i Message::Debug(
* %g", " -> System %d/%d: TimeStep = %d, Time = %g + i * %g", k + 1,
k+1, Current.NbrSystem, Current.NbrSystem,
(Current.DofData_P0+k)->CurrentSolution->TimeStep (Current.DofData_P0 + k)->CurrentSolution->TimeStep,
, (Current.DofData_P0 + k)->CurrentSolution->Time,
(Current.DofData_P0+k)->CurrentSolution->Time, (Current.DofData_P0 + k)->CurrentSolution->TimeImag);
(Current.DofData_P0+k)->CurrentSolution->TimeImag Flag_InfoForTime_ntime =
); List_Nbr((Current.DofData_P0 + k)->Solutions);
Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Soluti
ons);
} }
} }
else { else {
if (!Flag_WarningMissSolForTime_ntime) { if(!Flag_WarningMissSolForTime_ntime) {
Message::Warning("Missing solution for time step -%d computation Message::Warning("Missing solution for time step -%d "
(System #%d/%d)", "computation (System #%d/%d)",
ntime, k+1, Current.NbrSystem); ntime, k + 1, Current.NbrSystem);
Flag_WarningMissSolForTime_ntime = 1 ; Flag_WarningMissSolForTime_ntime = 1;
} }
} }
} }
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
Pos_FemInterpolation Pos_FemInterpolation(
(Element, Element, QuantityStorage_P0,
QuantityStorage_P0, QuantityStorage_P0 +
QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Inde WholeQuantity_P->Case.OperatorAndQuantity.Index,
x, WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity, WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0, u, v,
WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1);
u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
#else #else
Message::Error("TODO Post_FemInterpolation"); Message::Error("TODO Post_FemInterpolation");
#endif #endif
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
for (k = 0 ; k < Current.NbrSystem ; k++){ for(k = 0; k < Current.NbrSystem; k++) {
if(!List_Nbr((Current.DofData_P0+k)->Solutions)) continue; if(!List_Nbr((Current.DofData_P0 + k)->Solutions)) continue;
Solution_PN = (struct Solution*) Solution_PN = (struct Solution *)List_Pointer(
List_Pointer((Current.DofData_P0+k)->Solutions, (Current.DofData_P0 + k)->Solutions,
List_Nbr((Current.DofData_P0+k)->Solutions)-1); List_Nbr((Current.DofData_P0 + k)->Solutions) - 1);
if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime) if((Solution_PN - (Current.DofData_P0 + k)->CurrentSolution) >=
((Current.DofData_P0+k)->CurrentSolution) += ntime ; ntime)
((Current.DofData_P0 + k)->CurrentSolution) += ntime;
} }
} }
else else
Message::Error("Explicit (x,y,z,time) evaluation not implemented"); Message::Error("Explicit (x,y,z,time) evaluation not implemented");
} }
else{ else {
Message::Error("Explicit Dof{} evaluation out of context"); Message::Error("Explicit Dof{} evaluation out of context");
} }
Current.Element = Save_CurrentElement ; Current.Element = Save_CurrentElement;
Current.Region = Save_Region ; Current.Region = Save_Region;
break ; break;
case WQ_TRACE : /* Trace[WholeQuantity, Group] */ case WQ_TRACE: /* Trace[WholeQuantity, Group] */
Save_Region = Current.Region ; Save_Region = Current.Region;
if(!Element->ElementTrace){ if(!Element->ElementTrace) {
Message::Error("Trace must act on discrete quantity (and not in post-pro Message::Error(
cessing)"); "Trace must act on discrete quantity (and not in post-processing)");
break; break;
} }
Current.Region = Element->ElementTrace->Region ; Current.Region = Element->ElementTrace->Region;
if(WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity >= 0){ if(WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity >= 0) {
Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0, Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
WholeQuantity_P->Case.Trace.WholeQuantity, WholeQuantity_P->Case.Trace.WholeQuantity, Current.ut,
Current.ut, Current.vt, Current.wt, Current.vt, Current.wt,
WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity, WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity,
Nbr_Dof, DofValue, NbrArguments, ExpressionName) ; Nbr_Dof, DofValue, NbrArguments, ExpressionName);
DofIndexInWholeQuantity = DofIndex = Index ; DofIndexInWholeQuantity = DofIndex = Index;
} }
else{ else {
Current.x = Current.y = Current.z = 0. ; Current.x = Current.y = Current.z = 0.;
for (j = 0 ; j < Element->GeoElement->NbrNodes ; j++) { for(j = 0; j < Element->GeoElement->NbrNodes; j++) {
Current.x += Element->x[j] * Element->n[j] ; Current.x += Element->x[j] * Element->n[j];
Current.y += Element->y[j] * Element->n[j] ; Current.y += Element->y[j] * Element->n[j];
Current.z += Element->z[j] * Element->n[j] ; Current.z += Element->z[j] * Element->n[j];
} }
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
xyz2uvwInAnElement(Element->ElementTrace, Current.x, Current.y, Current. xyz2uvwInAnElement(Element->ElementTrace, Current.x, Current.y,
z, Current.z, &Current.ut, &Current.vt, &Current.wt);
&Current.ut, &Current.vt, &Current.wt) ;
#else #else
Message::Error("TODO xyz2uvwInAnElement"); Message::Error("TODO xyz2uvwInAnElement");
#endif #endif
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0, Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
WholeQuantity_P->Case.Trace.WholeQuantity, WholeQuantity_P->Case.Trace.WholeQuantity, Current.ut,
Current.ut, Current.vt, Current.wt, Current.vt, Current.wt, -1, 0, &Stack[0][Index],
-1, 0, &Stack[0][Index], NbrArguments, ExpressionName) NbrArguments, ExpressionName);
;
} }
Current.Region = Save_Region ; Current.Region = Save_Region;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_SOLIDANGLE : /* SolidAngle[{qty}] */ case WQ_SOLIDANGLE: /* SolidAngle[{qty}] */
Cal_SolidAngle(0, Element, QuantityStorage_P0 + Cal_SolidAngle(0, Element,
WholeQuantity_P->Case.OperatorAndQuantity.Index, QuantityStorage_P0 +
WholeQuantity_P->Case.OperatorAndQuantity.Index,
Nbr_Dof, Index, (struct Value **)Stack); Nbr_Dof, Index, (struct Value **)Stack);
Multi[Index] = 1 ; Multi[Index] = 1;
Index++ ; Index++;
break ; break;
case WQ_BINARYOPERATOR : /* + - * x / % ^ < > <= >= == != && || */ case WQ_BINARYOPERATOR: /* + - * x / % ^ < > <= >= == != && || */
if (Index-2 != DofIndex && Index-1 != DofIndex){ if(Index - 2 != DofIndex && Index - 1 != DofIndex) {
if(!Multi[Index-1] && !Multi[Index-2]) if(!Multi[Index - 1] && !Multi[Index - 2])
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[0][Index-2], &Stack[0][Index-1], &Stack[0][Index-2]) ; &Stack[0][Index - 2], &Stack[0][Index - 1], &Stack[0][Index - 2]);
else if(Multi[Index-1] && Multi[Index-2]) else if(Multi[Index - 1] && Multi[Index - 2])
for(j=0 ; j<Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[j][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ; &Stack[j][Index - 2], &Stack[j][Index - 1], &Stack[j][Index - 2]);
else if(Multi[Index-2]) else if(Multi[Index - 2])
for(j=0 ; j<Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[j][Index-2], &Stack[0][Index-1], &Stack[j][Index-2]) ; &Stack[j][Index - 2], &Stack[0][Index - 1], &Stack[j][Index - 2]);
else { else {
for(j=0 ; j<Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[0][Index-2], &Stack[j][Index-1], &Stack[j][Index-2]) ; &Stack[0][Index - 2], &Stack[j][Index - 1], &Stack[j][Index - 2]);
Multi[Index-2] = 1 ; Multi[Index - 2] = 1;
} }
} }
else if (Index-1 == DofIndex) { else if(Index - 1 == DofIndex) {
if(Multi[Index-2]) if(Multi[Index - 2])
for (j = 0 ; j < Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[j][Index-2], &DofValue[j], &DofValue[j]) ; &Stack[j][Index - 2], &DofValue[j], &DofValue[j]);
else else
for (j = 0 ; j < Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&Stack[0][Index-2], &DofValue[j], &DofValue[j]) ; &Stack[0][Index - 2], &DofValue[j], &DofValue[j]);
DofIndex-- ; DofIndex--;
} }
else { /* Index-2 == DofIndex */ else { /* Index-2 == DofIndex */
if(Multi[Index-1]) if(Multi[Index - 1])
for (j = 0 ; j < Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&DofValue[j], &Stack[j][Index-1], &DofValue[j]) ; &DofValue[j], &Stack[j][Index - 1], &DofValue[j]);
else else
for (j = 0 ; j < Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST3V)WholeQuantity_P->Case.Operator.Function) ((CAST3V)WholeQuantity_P->Case.Operator.Function)(
(&DofValue[j], &Stack[0][Index-1], &DofValue[j]) ; &DofValue[j], &Stack[0][Index - 1], &DofValue[j]);
} }
Index-- ; Index--;
break ; break;
case WQ_UNARYOPERATOR : /* + - ! */ case WQ_UNARYOPERATOR: /* + - ! */
if (Index-1 == DofIndex) if(Index - 1 == DofIndex)
for (j = 0 ; j < Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST1V)WholeQuantity_P->Case.Operator.Function)(&DofValue[j]) ; ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&DofValue[j]);
else if(Multi[Index-1]) else if(Multi[Index - 1])
for(j=0 ; j<Nbr_Dof ; j++) for(j = 0; j < Nbr_Dof; j++)
((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[j][Index-1]) ((CAST1V)WholeQuantity_P->Case.Operator.Function)(
; &Stack[j][Index - 1]);
else else
((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[0][Index-1]) ; ((CAST1V)WholeQuantity_P->Case.Operator.Function)(&Stack[0][Index - 1]);
break ; break;
/* WARNING: all the rest assumes 0 multi status */ /* WARNING: all the rest assumes 0 multi status */
case WQ_TEST : case WQ_TEST:
Flag_True = (Stack[0][Index-1].Val[0] != 0.) ; Flag_True = (Stack[0][Index - 1].Val[0] != 0.);
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index-1]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index - 1]);
Multi[Index-1] = 0 ; Multi[Index - 1] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(
(Flag_True) ? WholeQuantity_P->Case.Test.WholeQuantity_T Element, QuantityStorage_P0,
rue : (Flag_True) ? WholeQuantity_P->Case.Test.WholeQuantity_True :
WholeQuantity_P->Case.Test.WholeQuantity_F WholeQuantity_P->Case.Test.WholeQuantity_False,
alse, u, v, w, -1, 0, &Stack[0][Index - 1], NbrArguments, ExpressionName);
u, v, w, -1, 0, &Stack[0][Index-1], NbrArguments, Expres break;
sionName) ;
break ; case WQ_EXPRESSION:
Index -= WholeQuantity_P->Case.Expression.NbrArguments;
case WQ_EXPRESSION : Get_ValueOfExpression(
Index -= WholeQuantity_P->Case.Expression.NbrArguments ; (struct Expression *)List_Pointer(
Get_ValueOfExpression Problem_S.Expression, WholeQuantity_P->Case.Expression.Index),
((struct Expression*)List_Pointer QuantityStorage_P0, u, v, w, &Stack[0][Index],
(Problem_S.Expression, WholeQuantity_P->Case.Expression.Index), WholeQuantity_P->Case.Expression.NbrArguments);
QuantityStorage_P0, u, v, w, &Stack[0][Index], Multi[Index] = 0;
WholeQuantity_P->Case.Expression.NbrArguments) ; Index++;
Multi[Index] = 0 ; break;
Index++ ;
break ; case WQ_BUILTINFUNCTION:
Index -= WholeQuantity_P->Case.Function.NbrArguments;
case WQ_BUILTINFUNCTION :
Index -= WholeQuantity_P->Case.Function.NbrArguments ; if(Index != DofIndex)
((CASTF2V)WholeQuantity_P->Case.Function.Fct)(
if (Index != DofIndex) &WholeQuantity_P->Case.Function, &Stack[0][Index], &Stack[0][Index]);
((CASTF2V)WholeQuantity_P->Case.Function.Fct)
(&WholeQuantity_P->Case.Function, &Stack[0][Index], &Stack[0][Index])
;
else /* Dof must be the only argument, only valid with linear functions */ else /* Dof must be the only argument, only valid with linear functions */
for (j = 0 ; j < Nbr_Dof ; j++) { for(j = 0; j < Nbr_Dof; j++) {
Current.NumEntity = Current.NumEntities[j]; /* temp */ Current.NumEntity = Current.NumEntities[j]; /* temp */
((CASTF2V)WholeQuantity_P->Case.Function.Fct) ((CASTF2V)WholeQuantity_P->Case.Function.Fct)(
(&WholeQuantity_P->Case.Function, &DofValue[j], &DofValue[j]) ; &WholeQuantity_P->Case.Function, &DofValue[j], &DofValue[j]);
} }
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_EXTERNBUILTINFUNCTION : case WQ_EXTERNBUILTINFUNCTION:
((CASTF2V)WholeQuantity_P->Case.Function.Fct) ((CASTF2V)WholeQuantity_P->Case.Function.Fct)(
(&WholeQuantity_P->Case.Function, DofValue, &Stack[0][Index]) ; &WholeQuantity_P->Case.Function, DofValue, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_CONSTANT : case WQ_CONSTANT:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
Stack[0][Index].Val[0] = WholeQuantity_P->Case.Constant ; Stack[0][Index].Val[0] = WholeQuantity_P->Case.Constant;
} }
else { else {
for (k = 0 ; k < Current.NbrHar ; k += 2) { for(k = 0; k < Current.NbrHar; k += 2) {
Stack[0][Index].Val[MAX_DIM* k ] = WholeQuantity_P->Case.Constant ; Stack[0][Index].Val[MAX_DIM * k] = WholeQuantity_P->Case.Constant;
Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ; Stack[0][Index].Val[MAX_DIM * (k + 1)] = 0.;
} }
} }
Stack[0][Index].Type = SCALAR ; Stack[0][Index].Type = SCALAR;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_CURRENTVALUE : case WQ_CURRENTVALUE:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
Stack[0][Index].Val[0] = *(WholeQuantity_P->Case.CurrentValue.Value) ; Stack[0][Index].Val[0] = *(WholeQuantity_P->Case.CurrentValue.Value);
} }
else { else {
for (k = 0 ; k < Current.NbrHar ; k += 2) { for(k = 0; k < Current.NbrHar; k += 2) {
Stack[0][Index].Val[MAX_DIM* k ] = *(WholeQuantity_P->Case.CurrentVa Stack[0][Index].Val[MAX_DIM * k] =
lue.Value) ; *(WholeQuantity_P->Case.CurrentValue.Value);
Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ; Stack[0][Index].Val[MAX_DIM * (k + 1)] = 0.;
} }
} }
Stack[0][Index].Type = SCALAR ; Stack[0][Index].Type = SCALAR;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_ARGUMENT : case WQ_ARGUMENT:
if (WholeQuantity_P->Case.Argument.Index > NbrArguments){ if(WholeQuantity_P->Case.Argument.Index > NbrArguments) {
Message::Error("Function %s called with too few arguments.", ExpressionN Message::Error("Function %s called with too few arguments.",
ame); ExpressionName);
} }
Cal_CopyValue(DofValue + WholeQuantity_P->Case.Argument.Index - 1, Cal_CopyValue(DofValue + WholeQuantity_P->Case.Argument.Index - 1,
&Stack[0][Index]) ; &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_TIMEDERIVATIVE : case WQ_TIMEDERIVATIVE:
if (Current.TypeTime == TIME_GEAR) { if(Current.TypeTime == TIME_GEAR) {
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.TimeDerivative.WholeQuantity, WholeQuantity_P->Case.TimeDerivative.WholeQuantity, u,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expres v, w, -1, 0, &Stack[0][Index], NbrArguments,
sionName); ExpressionName);
for (k = 0 ; k < Current.NbrSystem ; k++) for(k = 0; k < Current.NbrSystem; k++)
(Current.DofData_P0+k)->Save_CurrentSolution = (Current.DofData_P0 + k)->Save_CurrentSolution =
(Current.DofData_P0+k)->CurrentSolution; (Current.DofData_P0 + k)->CurrentSolution;
Save_TimeStep = Current.TimeStep ; Save_TimeStep = Current.TimeStep;
Save_Time = Current.Time ; Save_Time = Current.Time;
Save_TimeImag = Current.TimeImag ; Save_TimeImag = Current.TimeImag;
for (int n = 0; n < Current.CorrOrder; n++) { for(int n = 0; n < Current.CorrOrder; n++) {
for(k = 0; k < Current.NbrSystem; k++) {
for (k = 0 ; k < Current.NbrSystem ; k++){ Solution_P0 = (struct Solution *)List_Pointer(
Solution_P0 = (struct Solution*)List_Pointer (Current.DofData_P0 + k)->Solutions, 0);
((Current.DofData_P0+k)->Solutions, 0); if(((Current.DofData_P0 + k)->CurrentSolution - Solution_P0) > 0)
if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) > 0) ((Current.DofData_P0 + k)->CurrentSolution) -= 1;
((Current.DofData_P0+k)->CurrentSolution) -= 1 ; else {
else{
Message::Error("Too few solutions for Dt with Gear's method"); Message::Error("Too few solutions for Dt with Gear's method");
break; break;
} }
} }
Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ; Current.TimeStep = Current.DofData->CurrentSolution->TimeStep;
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index+1]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index + 1]);
Multi[Index+1] = 0 ; Multi[Index + 1] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.TimeDerivative.WholeQuantity, WholeQuantity_P->Case.TimeDerivative.WholeQuantity,
u, v, w, -1, 0, &Stack[0][Index+1], NbrArguments, Ex u, v, w, -1, 0, &Stack[0][Index + 1], NbrArguments,
pressionName); ExpressionName);
Cal_AddMultValue(&Stack[0][Index], &Stack[0][Index+1], Cal_AddMultValue(&Stack[0][Index], &Stack[0][Index + 1],
-Current.aCorrCoeff[n], &Stack[0][Index]); -Current.aCorrCoeff[n], &Stack[0][Index]);
} }
Cal_MultValue(&Stack[0][Index], 1./(Current.bCorrCoeff*Current.DTime), Cal_MultValue(&Stack[0][Index],
1. / (Current.bCorrCoeff * Current.DTime),
&Stack[0][Index]); &Stack[0][Index]);
for (k = 0 ; k < Current.NbrSystem ; k++) for(k = 0; k < Current.NbrSystem; k++)
(Current.DofData_P0+k)->CurrentSolution = (Current.DofData_P0 + k)->CurrentSolution =
(Current.DofData_P0+k)->Save_CurrentSolution; (Current.DofData_P0 + k)->Save_CurrentSolution;
Current.TimeStep = Save_TimeStep ; Current.TimeStep = Save_TimeStep;
Current.Time = Save_Time ; Current.Time = Save_Time;
Current.TimeImag = Save_TimeImag ; Current.TimeImag = Save_TimeImag;
} }
else if (Current.NbrHar == 1) { else if(Current.NbrHar == 1) {
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.TimeDerivative.WholeQuantity, WholeQuantity_P->Case.TimeDerivative.WholeQuantity, u,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expres v, w, -1, 0, &Stack[0][Index], NbrArguments,
sionName) ; ExpressionName);
for (k = 0 ; k < Current.NbrSystem ; k++){ for(k = 0; k < Current.NbrSystem; k++) {
(Current.DofData_P0+k)->Save_CurrentSolution = (Current.DofData_P0 + k)->Save_CurrentSolution =
(Current.DofData_P0+k)->CurrentSolution; (Current.DofData_P0 + k)->CurrentSolution;
if(List_Nbr((Current.DofData_P0+k)->Solutions) > 1){ if(List_Nbr((Current.DofData_P0 + k)->Solutions) > 1) {
Solution_P0 = (struct Solution*)List_Pointer Solution_P0 = (struct Solution *)List_Pointer(
((Current.DofData_P0+k)->Solutions, 0); (Current.DofData_P0 + k)->Solutions, 0);
if ((Current.DofData_P0+k)->CurrentSolution != Solution_P0) if((Current.DofData_P0 + k)->CurrentSolution != Solution_P0)
((Current.DofData_P0+k)->CurrentSolution) -- ; ((Current.DofData_P0 + k)->CurrentSolution)--;
} }
else { else {
if (!Flag_WarningMissSolForDt) { if(!Flag_WarningMissSolForDt) {
Message::Warning("Missing solution for time derivative computation Message::Warning(
" "Missing solution for time derivative computation "
"(Sys#%d/%d)", k+1, Current.NbrSystem); "(Sys#%d/%d)",
Flag_WarningMissSolForDt = 1 ; k + 1, Current.NbrSystem);
Flag_WarningMissSolForDt = 1;
} }
} }
} }
Save_TimeStep = Current.TimeStep ; Save_TimeStep = Current.TimeStep;
Save_Time = Current.Time ; Save_Time = Current.Time;
Save_TimeImag = Current.TimeImag ; Save_TimeImag = Current.TimeImag;
Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ; Current.TimeStep = Current.DofData->CurrentSolution->TimeStep;
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index+1]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index + 1]);
Multi[Index+1] = 0 ; Multi[Index + 1] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.TimeDerivative.WholeQuantity, WholeQuantity_P->Case.TimeDerivative.WholeQuantity, u,
u, v, w, -1, 0, &Stack[0][Index+1], NbrArguments, Expr v, w, -1, 0, &Stack[0][Index + 1], NbrArguments,
essionName) ; ExpressionName);
Cal_SubstractValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index Cal_SubstractValue(&Stack[0][Index], &Stack[0][Index + 1],
]) ; &Stack[0][Index]);
Stack[0][Index+1].Val[0] = Save_Time - Current.Time ; Stack[0][Index + 1].Val[0] = Save_Time - Current.Time;
Stack[0][Index+1].Type = SCALAR ; Stack[0][Index + 1].Type = SCALAR;
if(Stack[0][Index+1].Val[0]) if(Stack[0][Index + 1].Val[0])
Cal_DivideValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index] Cal_DivideValue(&Stack[0][Index], &Stack[0][Index + 1],
) ; &Stack[0][Index]);
else else
Cal_ZeroValue(&Stack[0][Index]); Cal_ZeroValue(&Stack[0][Index]);
for (k = 0 ; k < Current.NbrSystem ; k++) for(k = 0; k < Current.NbrSystem; k++)
(Current.DofData_P0+k)->CurrentSolution = (Current.DofData_P0 + k)->CurrentSolution =
(Current.DofData_P0+k)->Save_CurrentSolution; (Current.DofData_P0 + k)->Save_CurrentSolution;
Current.TimeStep = Save_TimeStep ; Current.TimeStep = Save_TimeStep;
Current.Time = Save_Time ; Current.Time = Save_Time;
Current.TimeImag = Save_TimeImag ; Current.TimeImag = Save_TimeImag;
} }
else { else {
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.TimeDerivative.WholeQuantity, WholeQuantity_P->Case.TimeDerivative.WholeQuantity, u,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expres v, w, -1, 0, &Stack[0][Index], NbrArguments,
sionName) ; ExpressionName);
for (k = 0 ; k < Current.NbrHar ; k += 2) { for(k = 0; k < Current.NbrHar; k += 2) {
Stack[0][Index+1].Val[MAX_DIM* k ] = 0. ; Stack[0][Index + 1].Val[MAX_DIM * k] = 0.;
Stack[0][Index+1].Val[MAX_DIM*(k+1)] = Current.DofData->Val_Pulsation[ Stack[0][Index + 1].Val[MAX_DIM * (k + 1)] =
k/2] ; Current.DofData->Val_Pulsation[k / 2];
} }
Stack[0][Index+1].Type = SCALAR ; Stack[0][Index + 1].Type = SCALAR;
Cal_ProductValue(&Stack[0][Index], &Stack[0][Index+1], &Stack[0][Index]) Cal_ProductValue(&Stack[0][Index], &Stack[0][Index + 1],
; &Stack[0][Index]);
} }
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_ATANTERIORTIMESTEP : case WQ_ATANTERIORTIMESTEP:
ntime = WholeQuantity_P->Case.AtAnteriorTimeStep.TimeStep ; ntime = WholeQuantity_P->Case.AtAnteriorTimeStep.TimeStep;
for (k = 0 ; k < Current.NbrSystem ; k++){ for(k = 0; k < Current.NbrSystem; k++) {
Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Sol Solution_P0 = (struct Solution *)List_Pointer(
utions, 0); (Current.DofData_P0 + k)->Solutions, 0);
if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime){ if(((Current.DofData_P0 + k)->CurrentSolution - Solution_P0) >= ntime) {
((Current.DofData_P0+k)->CurrentSolution) -= ntime ; ((Current.DofData_P0 + k)->CurrentSolution) -= ntime;
if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Solutio if(Flag_InfoForTime_ntime !=
ns)) { List_Nbr((Current.DofData_P0 + k)->Solutions)) {
Message::Info("Accessing solution from %d time steps ago", ntime); Message::Info("Accessing solution from %d time steps ago", ntime);
Message::Info(" -> System %d/%d: TimeStep = %d, Time = %g + i * %g" Message::Info(
, " -> System %d/%d: TimeStep = %d, Time = %g + i * %g", k + 1,
k+1, Current.NbrSystem, Current.NbrSystem,
(Current.DofData_P0+k)->CurrentSolution->TimeStep, (Current.DofData_P0 + k)->CurrentSolution->TimeStep,
(Current.DofData_P0+k)->CurrentSolution->Time, (Current.DofData_P0 + k)->CurrentSolution->Time,
(Current.DofData_P0+k)->CurrentSolution->TimeImag); (Current.DofData_P0 + k)->CurrentSolution->TimeImag);
Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Solutions) Flag_InfoForTime_ntime =
; List_Nbr((Current.DofData_P0 + k)->Solutions);
} }
} }
else { else {
if (!Flag_WarningMissSolForTime_ntime) { if(!Flag_WarningMissSolForTime_ntime) {
Message::Warning("Missing solution for time step -%d computation " Message::Warning("Missing solution for time step -%d computation "
"(System #%d/%d)", ntime, k+1, Current.NbrSystem); "(System #%d/%d)",
Flag_WarningMissSolForTime_ntime = 1 ; ntime, k + 1, Current.NbrSystem);
Flag_WarningMissSolForTime_ntime = 1;
} }
} }
} }
Save_TimeStep = Current.TimeStep ; Save_TimeStep = Current.TimeStep;
Save_Time = Current.Time ; Save_Time = Current.Time;
Save_TimeImag = Current.TimeImag ; Save_TimeImag = Current.TimeImag;
Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ; Current.TimeStep = Current.DofData->CurrentSolution->TimeStep;
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.AtAnteriorTimeStep.WholeQuantity, WholeQuantity_P->Case.AtAnteriorTimeStep.WholeQuantity,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expressi u, v, w, -1, 0, &Stack[0][Index], NbrArguments,
onName) ; ExpressionName);
Current.TimeStep = Save_TimeStep;
Current.Time = Save_Time;
Current.TimeImag = Save_TimeImag;
for(k = 0; k < Current.NbrSystem; k++) {
Solution_PN = (struct Solution *)List_Pointer(
(Current.DofData_P0 + k)->Solutions,
List_Nbr((Current.DofData_P0 + k)->Solutions) - 1);
if((Solution_PN - (Current.DofData_P0 + k)->CurrentSolution) >= ntime)
((Current.DofData_P0 + k)->CurrentSolution) += ntime;
}
Current.TimeStep = Save_TimeStep ; Multi[Index] = 0;
Current.Time = Save_Time ; Index++;
Current.TimeImag = Save_TimeImag ; break;
for (k = 0 ; k < Current.NbrSystem ; k++){
Solution_PN = (struct Solution*)
List_Pointer((Current.DofData_P0+k)->Solutions,
List_Nbr((Current.DofData_P0+k)->Solutions)-1);
if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime)
((Current.DofData_P0+k)->CurrentSolution) += ntime ;
}
Multi[Index] = 0 ;
Index++ ;
break ;
case WQ_MAXOVERTIME : case WQ_MAXOVERTIME:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
double time_init = WholeQuantity_P->Case.MaxOverTime.TimeInit; double time_init = WholeQuantity_P->Case.MaxOverTime.TimeInit;
double time_final = WholeQuantity_P->Case.MaxOverTime.TimeFinal; double time_final = WholeQuantity_P->Case.MaxOverTime.TimeFinal;
/* /*
for (k = 0 ; k < Current.NbrSystem ; k++) for (k = 0 ; k < Current.NbrSystem ; k++)
(Current.DofData_P0+k)->Save_CurrentSolution = (Current.DofData_P0+k)->Save_CurrentSolution =
(Current.DofData_P0+k)->CurrentSolution; (Current.DofData_P0+k)->CurrentSolution;
*/ */
Save_TimeStep = Current.TimeStep ; Save_TimeStep = Current.TimeStep;
Save_Time = Current.Time ; Save_Time = Current.Time;
Save_TimeImag = Current.TimeImag ; Save_TimeImag = Current.TimeImag;
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
double val_maxInTime = -1.e99; double val_maxInTime = -1.e99;
for (j=1; j<List_Nbr((Current.DofData)->Solutions); j++) { for(j = 1; j < List_Nbr((Current.DofData)->Solutions); j++) {
Current.DofData->CurrentSolution =
Current.DofData->CurrentSolution = (struct Solution *)List_Pointer((Current.DofData)->Solutions, j);
(struct Solution*)List_Pointer((Current.DofData)->Solutions, j);
//++++ Add: also for other systems! //++++ Add: also for other systems!
Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ; Current.TimeStep = Current.DofData->CurrentSolution->TimeStep;
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
//++++ test to do more accurately! //++++ test to do more accurately!
if (Current.Time >= time_init && Current.Time <= time_final) { if(Current.Time >= time_init && Current.Time <= time_final) {
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.MaxOverTime.WholeQuantity, WholeQuantity_P->Case.MaxOverTime.WholeQuantity,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Ex u, v, w, -1, 0, &Stack[0][Index], NbrArguments,
pressionName) ; ExpressionName);
if (Stack[0][Index].Type == SCALAR) { if(Stack[0][Index].Type == SCALAR) {
if (Stack[0][Index].Val[0] > val_maxInTime){ if(Stack[0][Index].Val[0] > val_maxInTime) {
val_maxInTime = Stack[0][Index].Val[0]; val_maxInTime = Stack[0][Index].Val[0];
} }
} }
else { else {
Message::Error("MaxOverTime can only be applied on scalar values") Message::Error(
; "MaxOverTime can only be applied on scalar values");
} }
} }
} }
Stack[0][Index].Val[0] = val_maxInTime; Stack[0][Index].Val[0] = val_maxInTime;
/* /*
for (k = 0 ; k < Current.NbrSystem ; k++) for (k = 0 ; k < Current.NbrSystem ; k++)
(Current.DofData_P0+k)->CurrentSolution = (Current.DofData_P0+k)->CurrentSolution =
(Current.DofData_P0+k)->Save_CurrentSolution; (Current.DofData_P0+k)->Save_CurrentSolution;
*/ */
Current.TimeStep = Save_TimeStep ; Current.TimeStep = Save_TimeStep;
Current.Time = Save_Time ; Current.Time = Save_Time;
Current.TimeImag = Save_TimeImag ; Current.TimeImag = Save_TimeImag;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
else { else {
Message::Error("MaxOverTime can only be used in time domain") ; Message::Error("MaxOverTime can only be used in time domain");
break; break;
} }
break ; break;
case WQ_FOURIERSTEINMETZ : case WQ_FOURIERSTEINMETZ:
if (Current.NbrHar == 1) { if(Current.NbrHar == 1) {
double time_init = WholeQuantity_P->Case.FourierSteinmetz.TimeInit; double time_init = WholeQuantity_P->Case.FourierSteinmetz.TimeInit;
double time_final = WholeQuantity_P->Case.FourierSteinmetz.TimeFinal; double time_final = WholeQuantity_P->Case.FourierSteinmetz.TimeFinal;
int nbrFrequencyInFormula = WholeQuantity_P->Case.FourierSteinmetz.NbrFr int nbrFrequencyInFormula =
equency; WholeQuantity_P->Case.FourierSteinmetz.NbrFrequency;
double exponent_f = WholeQuantity_P->Case.FourierSteinmetz.Exponent_f; double exponent_f = WholeQuantity_P->Case.FourierSteinmetz.Exponent_f;
double exponent_b = WholeQuantity_P->Case.FourierSteinmetz.Exponent_b; double exponent_b = WholeQuantity_P->Case.FourierSteinmetz.Exponent_b;
/* /*
for (k = 0 ; k < Current.NbrSystem ; k++) for (k = 0 ; k < Current.NbrSystem ; k++)
(Current.DofData_P0+k)->Save_CurrentSolution = (Current.DofData_P0+k)->Save_CurrentSolution =
(Current.DofData_P0+k)->CurrentSolution; (Current.DofData_P0+k)->CurrentSolution;
*/ */
Save_TimeStep = Current.TimeStep ; Save_TimeStep = Current.TimeStep;
Save_Time = Current.Time ; Save_Time = Current.Time;
Save_TimeImag = Current.TimeImag ; Save_TimeImag = Current.TimeImag;
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
int i_Solution_init = -1, i_Solution_final = -1; int i_Solution_init = -1, i_Solution_final = -1;
int NbrTimeStep, Size=-1; int NbrTimeStep, Size = -1;
for (j=0; j<List_Nbr((Current.DofData)->Solutions); j++) { for(j = 0; j < List_Nbr((Current.DofData)->Solutions); j++) {
Current.DofData->CurrentSolution = Current.DofData->CurrentSolution =
(struct Solution*)List_Pointer((Current.DofData)->Solutions, j); (struct Solution *)List_Pointer((Current.DofData)->Solutions, j);
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
if (Current.Time >= time_init && i_Solution_init < 0) if(Current.Time >= time_init && i_Solution_init < 0)
i_Solution_init = j; i_Solution_init = j;
if (Current.Time <= time_final) { if(Current.Time <= time_final) { i_Solution_final = j; }
i_Solution_final = j; if(Current.Time > time_final) { break; }
}
if (Current.Time > time_final) {
break;
}
} }
NbrTimeStep = i_Solution_final-i_Solution_init+1; NbrTimeStep = i_Solution_final - i_Solution_init + 1;
if (NbrTimeStep < 2) if(NbrTimeStep < 2)
Message::Error("Wrong time interval in Function FourierSteinmetz (%d,% Message::Error(
d)", "Wrong time interval in Function FourierSteinmetz (%d,%d)",
i_Solution_init, i_Solution_final) ; i_Solution_init, i_Solution_final);
double *Times = (double *)Malloc(NbrTimeStep*sizeof(double)); double *Times = (double *)Malloc(NbrTimeStep * sizeof(double));
struct Value *TmpValues = (struct Value *)Malloc(NbrTimeStep*sizeof(stru struct Value *TmpValues =
ct Value)); (struct Value *)Malloc(NbrTimeStep * sizeof(struct Value));
for (j=0; j<NbrTimeStep; j++) { for(j = 0; j < NbrTimeStep; j++) {
Current.DofData->CurrentSolution = Current.DofData->CurrentSolution = (struct Solution *)List_Pointer(
(struct Solution*)List_Pointer((Current.DofData)->Solutions, (Current.DofData)->Solutions, i_Solution_init + j);
i_Solution_init+j);
//++++ Add: also for other systems! //++++ Add: also for other systems!
Current.TimeStep = Current.DofData->CurrentSolution->TimeStep ; Current.TimeStep = Current.DofData->CurrentSolution->TimeStep;
Current.Time = Current.DofData->CurrentSolution->Time ; Current.Time = Current.DofData->CurrentSolution->Time;
Current.TimeImag = Current.DofData->CurrentSolution->TimeImag ; Current.TimeImag = Current.DofData->CurrentSolution->TimeImag;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.MaxOverTime.WholeQuantity, WholeQuantity_P->Case.MaxOverTime.WholeQuantity, u,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expr v, w, -1, 0, &Stack[0][Index], NbrArguments,
essionName) ; ExpressionName);
Times[j] = Current.Time ; Times[j] = Current.Time;
Cal_CopyValue(&Stack[0][Index], &TmpValues[j]) ; Cal_CopyValue(&Stack[0][Index], &TmpValues[j]);
if (Stack[0][Index].Type == SCALAR) if(Stack[0][Index].Type == SCALAR)
Size = 1; Size = 1;
else if (Stack[0][Index].Type == VECTOR) else if(Stack[0][Index].Type == VECTOR)
Size = 3; Size = 3;
else else
Message::Error("FourierSteinmetz can only be applied on scalar or ve Message::Error("FourierSteinmetz can only be applied on scalar or "
ctor values") ; "vector values");
} }
// FourierTransform // FourierTransform
int NbrFreq ; int NbrFreq;
double *Frequencies; double *Frequencies;
struct Value *FourierValues; struct Value *FourierValues;
#if defined(HAVE_KERNEL) #if defined(HAVE_KERNEL)
Pos_FourierTransform(NbrTimeStep, 1, Times, TmpValues, Size, Pos_FourierTransform(NbrTimeStep, 1, Times, TmpValues, Size, 2,
2, nbrFrequencyInFormula, nbrFrequencyInFormula, &NbrFreq, &Frequencies,
&NbrFreq, &Frequencies, &FourierValues); &FourierValues);
#else #else
Message::Error("TODO Pos_FourierTransform"); Message::Error("TODO Pos_FourierTransform");
#endif #endif
/* /*
we calculate the Sum for all frequencies of frequency_i^exponent_f * we calculate the Sum for all frequencies of frequency_i^exponent_f *
b_i^exponent_b b_i^exponent_b
*/ */
if (nbrFrequencyInFormula > NbrFreq-1) if(nbrFrequencyInFormula > NbrFreq - 1)
Message::Error("FourierSteinmetz: too many frequencies asked " Message::Error("FourierSteinmetz: too many frequencies asked "
"(%d asked and only %d available)", "(%d asked and only %d available)",
nbrFrequencyInFormula, NbrFreq-1) ; nbrFrequencyInFormula, NbrFreq - 1);
double val=0.; double val = 0.;
for (j=0; j<nbrFrequencyInFormula; j++) { for(j = 0; j < nbrFrequencyInFormula; j++) {
if (Size==1) { if(Size == 1) {
val += pow(Frequencies[j+1], exponent_f) * val += pow(Frequencies[j + 1], exponent_f) *
pow(FourierValues[j+1].Val[0], exponent_b); pow(FourierValues[j + 1].Val[0], exponent_b);
} }
else { else {
val += val +=
pow(Frequencies[j+1], exponent_f) * pow(Frequencies[j + 1], exponent_f) *
pow(sqrt(FourierValues[j+1].Val[0]*FourierValues[j+1].Val[0] pow(
+ FourierValues[j+1].Val[1]*FourierValues[j+1].Val[1] sqrt(FourierValues[j + 1].Val[0] * FourierValues[j + 1].Val[0] +
+ FourierValues[j+1].Val[2]*FourierValues[j+1].Val[2]), FourierValues[j + 1].Val[1] * FourierValues[j + 1].Val[1] +
exponent_b); FourierValues[j + 1].Val[2] * FourierValues[j + 1].Val[2]),
exponent_b);
} }
} }
Stack[0][Index].Type = SCALAR; Stack[0][Index].Type = SCALAR;
Stack[0][Index].Val[0] = val; Stack[0][Index].Val[0] = val;
Free(Frequencies); Free(Frequencies);
Free(FourierValues); Free(FourierValues);
Free(Times); Free(TmpValues) ; Free(Times);
Free(TmpValues);
/* /*
for (k = 0 ; k < Current.NbrSystem ; k++) for (k = 0 ; k < Current.NbrSystem ; k++)
(Current.DofData_P0+k)->CurrentSolution = (Current.DofData_P0+k)->CurrentSolution =
(Current.DofData_P0+k)->Save_CurrentSolution; (Current.DofData_P0+k)->Save_CurrentSolution;
*/ */
Current.TimeStep = Save_TimeStep ; Current.TimeStep = Save_TimeStep;
Current.Time = Save_Time ; Current.Time = Save_Time;
Current.TimeImag = Save_TimeImag ; Current.TimeImag = Save_TimeImag;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
else { else {
Message::Error("FourierSteinmetz can only be used in time domain") ; Message::Error("FourierSteinmetz can only be used in time domain");
break; break;
} }
break ; break;
case WQ_MHTRANSFORM : case WQ_MHTRANSFORM:
if(Current.NbrHar == 1){ if(Current.NbrHar == 1) {
Message::Error("MHTransform can only be used in complex (multi-harmonic) Message::Error(
" "MHTransform can only be used in complex (multi-harmonic)"
" calculations") ; " calculations");
break; break;
} }
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
{ {
int N = List_Nbr(WholeQuantity_P->Case.MHTransform.WholeQuantity_L); int N = List_Nbr(WholeQuantity_P->Case.MHTransform.WholeQuantity_L);
std::vector<struct Value> MH_Values(N); std::vector<struct Value> MH_Values(N);
for(int j = 0; j < N; j++){ for(int j = 0; j < N; j++) {
List_T *WQ; List_Read(WholeQuantity_P->Case.MHTransform.WholeQuantity_ List_T *WQ;
L, j, &WQ); List_Read(WholeQuantity_P->Case.MHTransform.WholeQuantity_L, j, &WQ);
Cal_WholeQuantity(Element, QuantityStorage_P0, WQ, u, v, w, -1, 0, Cal_WholeQuantity(Element, QuantityStorage_P0, WQ, u, v, w, -1, 0,
&MH_Values[j], NbrArguments, ExpressionName) ; &MH_Values[j], NbrArguments, ExpressionName);
} }
MHTransform(Element, QuantityStorage_P0, u, v, w, MH_Values, MHTransform(
(struct Expression *)List_Pointer(Problem_S.Expression, Element, QuantityStorage_P0, u, v, w, MH_Values,
WholeQuantity_P->Case.MHTr (struct Expression *)List_Pointer(
ansform.Index), Problem_S.Expression, WholeQuantity_P->Case.MHTransform.Index),
WholeQuantity_P->Case.MHTransform.NbrPoints, Stack[0][Index] WholeQuantity_P->Case.MHTransform.NbrPoints, Stack[0][Index]);
) ; }
} Multi[Index] = 0;
Multi[Index] = 0 ; Index++;
Index++ ; break;
break ;
case WQ_CAST : case WQ_CAST:
/* This should be changed... */ /* This should be changed... */
Save_NbrHar = Current.NbrHar ; Save_NbrHar = Current.NbrHar;
Save_DofData = Current.DofData ; Save_DofData = Current.DofData;
if (!WholeQuantity_P->Case.Cast.NbrHar){ if(!WholeQuantity_P->Case.Cast.NbrHar) {
Current.DofData = Current.DofData =
((struct FunctionSpace *) ((struct FunctionSpace *)List_Pointer(
List_Pointer(Problem_S.FunctionSpace, Problem_S.FunctionSpace,
WholeQuantity_P->Case.Cast.FunctionSpaceIndexForType)) WholeQuantity_P->Case.Cast.FunctionSpaceIndexForType))
->DofData ; ->DofData;
Current.NbrHar = Current.DofData->NbrHar ; Current.NbrHar = Current.DofData->NbrHar;
} }
else{ else {
Current.NbrHar = WholeQuantity_P->Case.Cast.NbrHar ; Current.NbrHar = WholeQuantity_P->Case.Cast.NbrHar;
} }
for (j=0; j<NbrArguments; j++) { for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index]) ; Cal_CopyValue(DofValue + j, &Stack[0][Index]);
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
} }
Index -= NbrArguments ; Index -= NbrArguments;
Cal_WholeQuantity(Element, QuantityStorage_P0, Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.Cast.WholeQuantity, WholeQuantity_P->Case.Cast.WholeQuantity, u, v, w, -1,
u, v, w, -1, 0, &Stack[0][Index], NbrArguments, Expressi 0, &Stack[0][Index], NbrArguments, ExpressionName);
onName) ;
if (Current.NbrHar < Save_NbrHar) /* ne plus a completer ...?? */ if(Current.NbrHar < Save_NbrHar) /* ne plus a completer ...?? */
Cal_SetZeroHarmonicValue(&Stack[0][Index], Save_NbrHar) ; Cal_SetZeroHarmonicValue(&Stack[0][Index], Save_NbrHar);
Current.NbrHar = Save_NbrHar ; Current.NbrHar = Save_NbrHar;
Current.DofData = Save_DofData ; Current.DofData = Save_DofData;
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_CHANGECURRENTPOSITION :
Save_x = Current.x ; Save_y = Current.y ; Save_z = Current.z ;
Current.x = Stack[0][Index-1].Val[0] ;
Current.y = Stack[0][Index-1].Val[1] ;
Current.z = Stack[0][Index-1].Val[2] ;
for (j=0; j<NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index-1]) ;
Multi[Index-1] = 0 ;
Index++ ;
}
Index -= NbrArguments ;
Cal_WholeQuantity(Element, QuantityStorage_P0,
WholeQuantity_P->Case.ChangeCurrentPosition.WholeQuantit
y,
u, v, w, -1, 0, &Stack[0][Index-1], NbrArguments, Expres
sionName) ;
Current.x = Save_x ; Current.y = Save_y ; Current.z = Save_z ; case WQ_CHANGECURRENTPOSITION:
break ; Save_x = Current.x;
Save_y = Current.y;
Save_z = Current.z;
Current.x = Stack[0][Index - 1].Val[0];
Current.y = Stack[0][Index - 1].Val[1];
Current.z = Stack[0][Index - 1].Val[2];
for(j = 0; j < NbrArguments; j++) {
Cal_CopyValue(DofValue + j, &Stack[0][Index - 1]);
Multi[Index - 1] = 0;
Index++;
}
Index -= NbrArguments;
Cal_WholeQuantity(
Element, QuantityStorage_P0,
WholeQuantity_P->Case.ChangeCurrentPosition.WholeQuantity, u, v, w, -1,
0, &Stack[0][Index - 1], NbrArguments, ExpressionName);
Current.x = Save_x;
Current.y = Save_y;
Current.z = Save_z;
break;
case WQ_SAVEVALUE : case WQ_SAVEVALUE:
Cal_CopyValue(&Stack[0][Index-1], Cal_CopyValue(&Stack[0][Index - 1],
&ValueSaved[WholeQuantity_P->Case.SaveValue.Index]) ; &ValueSaved[WholeQuantity_P->Case.SaveValue.Index]);
break ; break;
case WQ_VALUESAVED : case WQ_VALUESAVED:
if(ValueSaved.count(WholeQuantity_P->Case.ValueSaved.Index)) if(ValueSaved.count(WholeQuantity_P->Case.ValueSaved.Index))
Cal_CopyValue(&ValueSaved[WholeQuantity_P->Case.ValueSaved.Index], Cal_CopyValue(&ValueSaved[WholeQuantity_P->Case.ValueSaved.Index],
&Stack[0][Index]) ; &Stack[0][Index]);
else{ else {
if(TreatmentStatus != STATUS_PRE) if(TreatmentStatus != STATUS_PRE)
Message::Warning("Empty register %d: assuming zero value", Message::Warning("Empty register %d: assuming zero value",
WholeQuantity_P->Case.ValueSaved.Index + 1); WholeQuantity_P->Case.ValueSaved.Index + 1);
Cal_ZeroValue(&Stack[0][Index]); Cal_ZeroValue(&Stack[0][Index]);
Stack[0][Index].Type = SCALAR ; Stack[0][Index].Type = SCALAR;
} }
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_SAVENAMEDVALUE : case WQ_SAVENAMEDVALUE:
Cal_CopyValue(&Stack[0][Index-1], Cal_CopyValue(&Stack[0][Index - 1],
&NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name]) ; &NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name]);
break ; break;
case WQ_NAMEDVALUESAVED : case WQ_NAMEDVALUESAVED:
if(NamedValueSaved.count(WholeQuantity_P->Case.NamedValue.Name)) if(NamedValueSaved.count(WholeQuantity_P->Case.NamedValue.Name))
Cal_CopyValue(&NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name], Cal_CopyValue(&NamedValueSaved[WholeQuantity_P->Case.NamedValue.Name],
&Stack[0][Index]) ; &Stack[0][Index]);
else{ else {
if(TreatmentStatus != STATUS_PRE) if(TreatmentStatus != STATUS_PRE)
Message::Warning("Unknown current value '$%s': assuming zero value", Message::Warning("Unknown current value '$%s': assuming zero value",
WholeQuantity_P->Case.NamedValue.Name); WholeQuantity_P->Case.NamedValue.Name);
Cal_ZeroValue(&Stack[0][Index]); Cal_ZeroValue(&Stack[0][Index]);
Stack[0][Index].Type = SCALAR ; Stack[0][Index].Type = SCALAR;
} }
Multi[Index] = 0 ; Multi[Index] = 0;
Index++ ; Index++;
break ; break;
case WQ_SHOWVALUE : case WQ_SHOWVALUE:
if (Index-1 == DofIndex) { if(Index - 1 == DofIndex) {
for(j=0 ; j<Nbr_Dof ; j++){ for(j = 0; j < Nbr_Dof; j++) {
fprintf(stderr, "##%d Dof %d ", WholeQuantity_P->Case.ShowValue.Index, fprintf(stderr, "##%d Dof %d ", WholeQuantity_P->Case.ShowValue.Index,
j+1); j + 1);
Show_Value(&DofValue[j]); Show_Value(&DofValue[j]);
} }
} }
else { else {
fprintf(stderr, "##%d ", WholeQuantity_P->Case.ShowValue.Index); fprintf(stderr, "##%d ", WholeQuantity_P->Case.ShowValue.Index);
Show_Value(&Stack[0][Index-1]); Show_Value(&Stack[0][Index - 1]);
} }
break ; break;
default : default:
Message::Error("Unknown type of WholeQuantity (%d)", WholeQuantity_P->Type Message::Error("Unknown type of WholeQuantity (%d)",
); WholeQuantity_P->Type);
break; break;
} }
} }
if (DofIndexInWholeQuantity < 0) Cal_CopyValue(&Stack[0][0], &DofValue[0]) ; if(DofIndexInWholeQuantity < 0) Cal_CopyValue(&Stack[0][0], &DofValue[0]);
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C a l _ S t o r e I n R e g i s t e r */ /* C a l _ S t o r e I n R e g i s t e r */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Cal_StoreInRegister(struct Value *Value, int RegisterIndex) void Cal_StoreInRegister(struct Value *Value, int RegisterIndex)
{ {
Cal_CopyValue(Value, &ValueSaved[RegisterIndex]) ; Cal_CopyValue(Value, &ValueSaved[RegisterIndex]);
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C a l _ S t o r e I n V a r i a b l e */ /* C a l _ S t o r e I n V a r i a b l e */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void Cal_StoreInVariable(struct Value *Value, const char *name) void Cal_StoreInVariable(struct Value *Value, const char *name)
{ {
Cal_CopyValue(Value, &NamedValueSaved[name]) ; Cal_CopyValue(Value, &NamedValueSaved[name]);
Export_Value(Value, GetDPNumbers[name], 0, false); // don't append Export_Value(Value, GetDPNumbers[name], 0, false); // don't append
} }
void Cal_GetValueSaved(struct Value *Value, const char *name) void Cal_GetValueSaved(struct Value *Value, const char *name)
{ {
if(NamedValueSaved.count(name)) if(NamedValueSaved.count(name))
Cal_CopyValue(&NamedValueSaved[name], Value) ; Cal_CopyValue(&NamedValueSaved[name], Value);
else{ else {
Message::Warning("Unknown current value '$%s': assuming zero value", name); Message::Warning("Unknown current value '$%s': assuming zero value", name);
Cal_ZeroValue(Value); Cal_ZeroValue(Value);
Value->Type = SCALAR ; Value->Type = SCALAR;
} }
} }
std::map<std::string, struct Value> &Get_AllValueSaved() std::map<std::string, struct Value> &Get_AllValueSaved()
{ {
return NamedValueSaved; return NamedValueSaved;
} }
 End of changes. 213 change blocks. 
897 lines changed or deleted 857 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)