Commit ddd79a61 authored by Christophe Geuzaine's avatar Christophe Geuzaine

- new ValueFromTable function to access runtime maps created by…

- new ValueFromTable function to access runtime maps created by NodeTable/ElementTable post-operations
- more work on OptimizerInitialize/Update
parent 3be75593
Pipeline #1066 passed with stage
in 10 minutes and 5 seconds
......@@ -1156,12 +1156,13 @@ struct Operation {
} Copy;
struct {
char *algorithm;
int numConstraints;
List_T *lowerBounds, *upperBounds;
List_T *currentPointLowerBounds, *currentPointUpperBounds;
} OptimizerInitialize;
struct {
char *performanceVariable, *sensitivityVariable, *residualVariable;
int sensitivityField;
char *currentPoint; // input and ouput
char *objective, *constraints; // input
char *objectiveSensitivity, *constraintsSensitivity; // input
char *residual;
} OptimizerUpdate;
} Case;
......
......@@ -1046,6 +1046,7 @@ struct StringXFunction2Nbr F_Function[] = { /* #Par #Arg */
{"SetVariable" , (CAST)F_SetVariable , -1, -1 },
{"SetCumulativeVariable" , (CAST)F_SetCumulativeVariable , -1, -1 },
{"GetVariable" , (CAST)F_GetVariable , -1, -1 },
{"ValueFromTable" , (CAST)F_ValueFromTable , -1, -1 },
{"VirtualWork" , (CAST)F_VirtualWork , 0, 1 },
{"Felec" , (CAST)F_Felec , 0, 1 },
......
......@@ -465,6 +465,7 @@ extern std::map<std::string, std::vector<double> > CommandLineNumbers;
extern std::map<std::string, std::vector<std::string> > CommandLineStrings;
extern std::map<std::string, std::vector<double> > GetDPNumbers;
extern std::map<std::string, std::vector<std::string> > GetDPStrings;
extern std::map<std::string, std::map<int, std::vector<double> > > GetDPNumbersMap;
int getdp_yyparse();
void getdp_yyrestart(FILE*);
......
This diff is collapsed.
......@@ -814,7 +814,7 @@
#if ! defined YYSTYPE && ! defined YYSTYPE_IS_DECLARED
typedef union YYSTYPE
#line 186 "ProParser.y"
#line 187 "ProParser.y"
{
char *c;
int i;
......
......@@ -44,6 +44,7 @@ std::map<std::string, std::vector<double> > CommandLineNumbers;
std::map<std::string, std::vector<std::string> > CommandLineStrings;
std::map<std::string, std::vector<double> > GetDPNumbers;
std::map<std::string, std::vector<std::string> > GetDPStrings;
std::map<std::string, std::map<int, std::vector<double> > > GetDPNumbersMap;
// Static parser variables (accessible only in this file)
......@@ -1782,7 +1783,6 @@ ParametersForFunction :
| '{' '$' String__Index '}'
{ $$ = NULL; StringForParameter = $3; }
;
/* ------------------------------------------------------------------------ */
......@@ -5833,40 +5833,28 @@ OperationTerm :
Operation_P->Case.Copy.from = $3 ;
}
| tOptimizerInitialize '[' CharExpr ',' FExpr ','
ListOfFExpr ',' ListOfFExpr ']' tEND
| tOptimizerInitialize '[' CharExpr ',' ListOfFExpr ',' ListOfFExpr ']' tEND
{
Operation_P = (struct Operation*)
List_Pointer(Operation_L, List_Nbr(Operation_L)-1) ;
Operation_P->Type = OPERATION_OPTIMIZER_INITIALIZE;
Operation_P->Case.OptimizerInitialize.algorithm = $3;
Operation_P->Case.OptimizerInitialize.numConstraints = (int)$5;
Operation_P->Case.OptimizerInitialize.lowerBounds = $7;
Operation_P->Case.OptimizerInitialize.upperBounds = $9;
}
| tOptimizerUpdate '[' '$' String__Index ',' '$' String__Index ','
'$' String__Index ']' tEND
{
Operation_P = (struct Operation*)
List_Pointer(Operation_L, List_Nbr(Operation_L)-1) ;
Operation_P->Type = OPERATION_OPTIMIZER_UPDATE;
Operation_P->Case.OptimizerUpdate.performanceVariable = $4;
Operation_P->Case.OptimizerUpdate.sensitivityVariable = $7;
Operation_P->Case.OptimizerUpdate.residualVariable = $10;
Operation_P->Case.OptimizerUpdate.sensitivityField = -1;
Operation_P->Case.OptimizerInitialize.currentPointLowerBounds = $5;
Operation_P->Case.OptimizerInitialize.currentPointUpperBounds = $7;
}
| tOptimizerUpdate '[' '$' String__Index ',' FExpr ','
'$' String__Index ']' tEND
| tOptimizerUpdate '[' CharExpr ',' CharExpr ',' CharExpr ','
CharExpr ',' CharExpr ',' '$' String__Index ']' tEND
{
Operation_P = (struct Operation*)
List_Pointer(Operation_L, List_Nbr(Operation_L)-1) ;
Operation_P->Type = OPERATION_OPTIMIZER_UPDATE;
Operation_P->Case.OptimizerUpdate.performanceVariable = $4;
Operation_P->Case.OptimizerUpdate.sensitivityVariable = 0;
Operation_P->Case.OptimizerUpdate.residualVariable = $9;
Operation_P->Case.OptimizerUpdate.sensitivityField = $6;
Operation_P->Case.OptimizerUpdate.currentPoint = $3;
Operation_P->Case.OptimizerUpdate.objective = $5;
Operation_P->Case.OptimizerUpdate.constraints = $7;
Operation_P->Case.OptimizerUpdate.objectiveSensitivity = $9;
Operation_P->Case.OptimizerUpdate.constraintsSensitivity = $11;
Operation_P->Case.OptimizerUpdate.residual = $14;
}
| Loop
......
......@@ -210,6 +210,7 @@ void F_GetNumberRunTime(F_ARG) ;
void F_SetVariable (F_ARG) ;
void F_SetCumulativeVariable (F_ARG) ;
void F_GetVariable (F_ARG) ;
void F_ValueFromTable (F_ARG) ;
void F_VirtualWork (F_ARG) ;
void F_Felec (F_ARG) ;
......
......@@ -352,7 +352,7 @@ bool Fi_dInterpolationBilinear(double *x, double *y, double *M, int NL, int NC,
bool Fi_InterpolationTrilinear (double *x, double *y, double *z, double *M, int NX, int NY, int NZ,
double xp, double yp, double zp, double *vp)
{
/*
/*
*
* a122 **************************** a222
* * | * *
......@@ -375,62 +375,62 @@ bool Fi_InterpolationTrilinear (double *x, double *y, double *z, double *M, int
* ****************************
* a111 a211
*/
double a111, a121, a211, a221, a112, a122, a212, a222;
int i, j, k;
// Interpolate point (xp,yp,zp) in a regular grid
// x[i] <= xp < x[i+1]
// y[j] <= yp < y[j+1]
// z[k] <= zp < z[k+1]
*vp = 0.0 ;
// When (xp,yp,zp) lays outside the boundaries of the table:
// the nearest border is taken
if (xp < x[0]) xp = x[0];
else if (xp > x[NX-1]) xp = x[NX-1];
for (i=0 ; i<NX-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
i = (i >= NX) ? NX-1 : i;
if (yp < y[0]) yp = y[0];
else if (yp > y[NY-1]) yp = y[NY-1];
for (j=0 ; j<NY-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
j = (j >= NY) ? NY-1 : j;
if (zp < z[0]) zp = z[0];
else if (zp > z[NZ-1]) zp = z[NZ-1];
for (k=0 ; k<NZ-1 ; ++k) if (z[k+1] >= zp && zp >= z[k]) break;
k = (k >= NZ) ? NZ-1 : k;
a111 = M[ i + NX * j + NX * NY * k ];
a211 = M[(1+i) + NX * j + NX * NY * k ];
a121 = M[ i + NX * (1+j) + NX * NY * k ];
a221 = M[(1+i) + NX * (1+j) + NX * NY * k ];
a112 = M[ i + NX * j + NX * NY * (k+1)];
a212 = M[(1+i) + NX * j + NX * NY * (k+1)];
a122 = M[ i + NX * (1+j) + NX * NY * (k+1)];
a222 = M[(1+i) + NX * (1+j) + NX * NY * (k+1)];
double xd, yd, zd;
xd = (xp-x[i])/(x[i+1]-x[i]);
yd = (yp-y[j])/(y[j+1]-y[j]);
zd = (zp-z[k])/(z[k+1]-z[k]);
double a11, a12, a21, a22;
a11 = a111*(1-xd) + a211*xd;
a12 = a112*(1-xd) + a212*xd;
a21 = a121*(1-xd) + a221*xd;
a22 = a122*(1-xd) + a222*xd;
double a1, a2;
a1 = a11*(1-yd) + a21*yd;
a2 = a12*(1-yd) + a22*yd;
*vp = a1*(1-zd) + a2*zd;
return true ;
}
......@@ -460,57 +460,57 @@ bool Fi_dInterpolationTrilinear (double *x, double *y, double *z, double *M, in
* ****************************
* a111 a211
*/
double a111, a121, a211, a221, a112, a122, a212, a222;
int i, j, k;
// Interpolate point (xp,yp,zp) in a regular grid
// x[i] <= xp < x[i+1]
// y[j] <= yp < y[j+1]
// z[k] <= zp < z[k+1]
*dvp_dx = 0.0 ;
*dvp_dy = 0.0 ;
*dvp_dz = 0.0 ;
// When (xp,yp,zp) lays outside the boundaries of the table:
// the nearest border is taken
if (xp < x[0]) xp = x[0];
else if (xp > x[NX-1]) xp = x[NX-1];
for (i=0 ; i<NX-1 ; ++i) if (x[i+1] >= xp && xp >= x[i]) break;
i = (i >= NX) ? NX-1 : i;
if (yp < y[0]) yp = y[0];
else if (yp > y[NY-1]) yp = y[NY-1];
for (j=0 ; j<NY-1 ; ++j) if (y[j+1] >= yp && yp >= y[j]) break;
j = (j >= NY) ? NY-1 : j;
if (zp < z[0]) zp = z[0];
else if (zp > z[NZ-1]) zp = z[NZ-1];
for (k=0 ; k<NZ-1 ; ++k) if (z[k+1] >= zp && zp >= z[j]) break;
k = (k >= NZ) ? NZ-1 : k;
a111 = M[ i + NX * j + NX * NY * k ];
a211 = M[(1+i) + NX * j + NX * NY * k ];
a121 = M[ i + NX * (1+j) + NX * NY * k ];
a221 = M[(1+i) + NX * (1+j) + NX * NY * k ];
a112 = M[ i + NX * j + NX * NY * (k+1)];
a212 = M[(1+i) + NX * j + NX * NY * (k+1)];
a122 = M[ i + NX * (1+j) + NX * NY * (k+1)];
a222 = M[(1+i) + NX * (1+j) + NX * NY * (k+1)];
double xd, yd, zd;
xd = (xp-x[i])/(x[i+1]-x[i]);
yd = (yp-y[j])/(y[j+1]-y[j]);
zd = (zp-z[k])/(z[k+1]-z[k]);
double dxd, dyd, dzd;
dxd = 1./(x[i+1]-x[i]);
dyd = 1./(y[j+1]-y[j]);
dzd = 1./(z[k+1]-z[k]);
double a11, a12, a21, a22, dxa11, dxa12, dxa21, dxa22;
a11 = a111*(1-xd) + a211*xd;
a12 = a112*(1-xd) + a212*xd;
......@@ -520,7 +520,7 @@ bool Fi_dInterpolationTrilinear (double *x, double *y, double *z, double *M, in
dxa12 = -a112*dxd + a212*dxd;
dxa21 = -a121*dxd + a221*dxd;
dxa22 = -a122*dxd + a222*dxd;
double a1, a2, dya1, dya2, dxa1, dxa2;
a1 = a11*(1-yd) + a21*yd;
a2 = a12*(1-yd) + a22*yd;
......@@ -528,11 +528,11 @@ bool Fi_dInterpolationTrilinear (double *x, double *y, double *z, double *M, in
dya2 = -a12*dyd + a22*dyd;
dxa1 = dxa11*(1-yd) + dxa21*yd;
dxa2 = dxa12*(1-yd) + dxa22*yd;
*dvp_dx = dxa1*(1-zd) + dxa2*zd;
*dvp_dy = dya1*(1-zd) + dya2*zd;
*dvp_dz = -a1*dzd + a2*dzd;
return true ;
}
......@@ -634,7 +634,7 @@ void F_InterpolationTrilinear(F_ARG)
/*
It performs a trilinear interpolation at point (xp,yp,zp) based
on a three-dimensional table (sorted grid).
Input parameters:
NX Number of lines X
NY Number of lines Y
......@@ -643,40 +643,40 @@ void F_InterpolationTrilinear(F_ARG)
y values (ascending order) linked to the NY lines of the table
z values (ascending order) linked to the NZ lines of the table
M Matrix M(x,y,z) = M[x+NX*y+NX*NY*z]
xp x coordinate of interpolation point
yp y coordinate of interpolation point
zp z coordinate of interpolation point
A. Royer
*/
int NX, NY, NZ;
double xp, yp, zp, vp=0., *x, *y, *z, *M;
struct FunctionActive * D;
if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
Message::Error("Three Scalar arguments required!");
if (!Fct->Active) Fi_InitListMatrix3D (Fct, A, V) ;
D = Fct->Active ;
NX = D->Case.ListMatrix3D.NbrLinesX ;
NY = D->Case.ListMatrix3D.NbrLinesY ;
NZ = D->Case.ListMatrix3D.NbrLinesZ ;
x = D->Case.ListMatrix3D.x ;
y = D->Case.ListMatrix3D.y ;
z = D->Case.ListMatrix3D.z ;
M = D->Case.ListMatrix3D.data ;
xp = (A+0)->Val[0] ;
yp = (A+1)->Val[0] ;
zp = (A+2)->Val[0] ;
bool IsInGrid = Fi_InterpolationTrilinear (x, y, z, M, NX, NY, NZ, xp, yp, zp, &vp);
if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g, zp=%g)", xp, yp, zp) ;
V->Type = SCALAR ;
V->Val[0] = vp ;
}
......@@ -686,7 +686,7 @@ void F_dInterpolationTrilinear(F_ARG)
/*
It delivers the derivative of the bilinear interpolation at point (xp, yp, zp)
based on a three-dimensional table (sorted grid).
Input parameters:
NX Number of lines X
NY Number of lines Y
......@@ -695,40 +695,40 @@ void F_dInterpolationTrilinear(F_ARG)
y values (ascending order) linked to the NY lines of the table
z values (ascending order) linked to the NZ lines of the table
M Matrix M(x,y,z) = M[x+NX*y+NX*NY*z]
xp x coordinate of interpolation point
yp y coordinate of interpolation point
zp z coordinate of interpolation point
A. Royer
*/
int NX, NY, NZ;
double xp, yp, zp, dvp_dx =0., dvp_dy =0., dvp_dz =0., *x, *y, *z, *M;
struct FunctionActive * D;
if( (A+0)->Type != SCALAR || (A+1)->Type != SCALAR || (A+2)->Type != SCALAR)
Message::Error("Three Scalar arguments required!");
if (!Fct->Active) Fi_InitListMatrix3D (Fct, A, V) ;
D = Fct->Active ;
NX = D->Case.ListMatrix3D.NbrLinesX ;
NY = D->Case.ListMatrix3D.NbrLinesY ;
NZ = D->Case.ListMatrix3D.NbrLinesZ ;
x = D->Case.ListMatrix3D.x ;
y = D->Case.ListMatrix3D.y ;
z = D->Case.ListMatrix3D.z ;
M = D->Case.ListMatrix3D.data ;
xp = (A+0)->Val[0] ;
yp = (A+1)->Val[0] ;
zp = (A+2)->Val[0] ;
bool IsInGrid = Fi_dInterpolationTrilinear (x, y, z, M, NX, NY, NZ, xp, yp, zp, &dvp_dx, &dvp_dy, &dvp_dz);
if (!IsInGrid) Message::Error("Extrapolation not allowed (xp=%g ; yp=%g, zp=%g)", xp, yp, zp) ;
V->Type = VECTOR ;
V->Val[0] = dvp_dx ;
V->Val[1] = dvp_dy ;
......@@ -786,28 +786,28 @@ void Fi_InitListMatrix3D(F_ARG)
{
int i=0, k, NX, NY, NZ, sz ;
struct FunctionActive * D ;
/*
The original table structure data(i,j,k) is furnished with the following format:
[ NX, NY, NZ, x(1..NX), y(1..NY), y(1..NZ), data(1..NX*NY*NZ) ]
A. Royer
*/
D = Fct->Active =
(struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
NX = Fct->Para[i++];
NY = Fct->Para[i++];
NZ = Fct->Para[i++];
sz = 3 + NX + NY + NZ + NX*NY*NZ ; // expected size of list matrix
if (Fct->NbrParameters != sz)
Message::Error("Bad size of input data (expected = %d ; found = %d). "
"List with format: x(NbrLines=%d), y(NbrLines=%d), z(NbrLines=%d), "
"matrix(NbrLinesX*NbrLinesY*NbrLinesZ=%d)",
sz, Fct->NbrParameters, NX, NY, NZ, NX*NY*NZ);
// Initialize structure and allocate memory
D->Case.ListMatrix3D.NbrLinesX = NX;
D->Case.ListMatrix3D.NbrLinesY = NY;
......@@ -816,7 +816,7 @@ void Fi_InitListMatrix3D(F_ARG)
D->Case.ListMatrix3D.y = (double *) malloc (sizeof(double)*NY);
D->Case.ListMatrix3D.z = (double *) malloc (sizeof(double)*NZ);
D->Case.ListMatrix3D.data = (double *) malloc (sizeof(double)*NX*NY*NZ);
// Assign values
for (k=0 ; k<NX ; ++k) D->Case.ListMatrix3D.x[k] = Fct->Para[i++];
for (k=0 ; k<NY ; ++k) D->Case.ListMatrix3D.y[k] = Fct->Para[i++];
......
......@@ -9,6 +9,7 @@
#include "GetDPConfig.h"
#include "ProData.h"
#include "ProDefine.h"
#include "ProParser.h"
#include "F.h"
#include "Message.h"
#include "Cal_Value.h"
......@@ -209,6 +210,42 @@ void F_GetVariable (F_ARG)
Cal_GetValueSaved(V, tmp);
}
void F_ValueFromTable (F_ARG)
{
if(!Fct->String){
Message::Error("Missing ElementTable or NodeTable name: use ValueFromTable[...]{name}");
return;
}
std::map<int, std::vector<double> > &table(GetDPNumbersMap[Fct->String]);
std::vector<double> &val(table[Current.NumEntity]);
if(val.size() == 1){
V->Val[0] = val[0];
V->Type = SCALAR ;
return;
}
Message::Debug("Unknown entity index %d in ValueFromTable",
Current.NumEntity);
for(int i = 0; i < Fct->NbrArguments; i++){
if(i != 0){
Message::Warning("Ignoring %d-th argument in ValueFromTable");
continue;
}
if((A+i)->Type != SCALAR){
Message::Error("Non-scalar default argument in ValueFromTable");
return;
}
Cal_CopyValue(A + i, V);
return;
}
Message::Warning("Missing table data or default value in ValueFromTable");
V->Val[0] = 0. ;
V->Type = SCALAR ;
}
void F_VirtualWork (F_ARG)
{
MATRIX3x3 Jac ;
......
......@@ -9,7 +9,8 @@
#include <string>
#include <vector>
#include "ProData.h"
#include "DofData.h"
#include "ProParser.h"
#include "Cal_Quantity.h"
#include "Message.h"
#include "GetDPConfig.h"
......@@ -17,12 +18,72 @@
void Operation_OptimizerInitialize(struct Operation *Operation_P)
{
printf("opti init\n");
printf("Opti init:\n");
printf(" - algorithm: %s\n", Operation_P->Case.OptimizerInitialize.algorithm);
List_T *lb = Operation_P->Case.OptimizerInitialize.currentPointLowerBounds;
printf(" - lower bounds: ");
for(int i = 0; i < List_Nbr(lb); i++){
double d;
List_Read(lb, i, &d);
printf("%g ", d);
}
printf("\n");
List_T *ub = Operation_P->Case.OptimizerInitialize.currentPointUpperBounds;
printf(" - upper bounds: ");
for(int i = 0; i < List_Nbr(ub); i++){
double d;
List_Read(ub, i, &d);
printf("%g ", d);
}
printf("\n");
}
static void debugInput(const std::string &type, const std::string &name)
{
std::map<std::string, std::map<int, std::vector<double> > >::iterator itv =
GetDPNumbersMap.find(name);
if(itv != GetDPNumbersMap.end()){
printf(" - table %s:\n", type.c_str());
for(std::map<int, std::vector<double> >::iterator it = itv->second.begin();
it != itv->second.end(); it++){
printf(" ele %d: ", it->first);
for(unsigned int i = 0; i < it->second.size(); i++)
printf("%g ", it->second[i]);
printf("\n");
}
}
else{
std::map<std::string, std::vector<double> >::iterator its =
GetDPNumbers.find(name);
if(its != GetDPNumbers.end()){
printf(" - scalar variable %s: ", type.c_str());
for(unsigned int i = 0; i < its->second.size(); i++)
printf("%g ", its->second[i]);
printf("\n");
}
else{
Message::Warning("Unknown %s: %s", type.c_str(), name.c_str());
}
}
}
void Operation_OptimizerUpdate(struct Operation *Operation_P)
{
printf("opti update\n");
printf("Opti update:\n");
debugInput("currentPoint", Operation_P->Case.OptimizerUpdate.currentPoint);
debugInput("objective", Operation_P->Case.OptimizerUpdate.objective);
debugInput("constraints", Operation_P->Case.OptimizerUpdate.constraints);
debugInput("objectiveSensitivity", Operation_P->Case.OptimizerUpdate.objectiveSensitivity);
debugInput("constraintsSensitivity", Operation_P->Case.OptimizerUpdate.constraintsSensitivity);
Value v;
v.Type = SCALAR;
v.Val[0] = 1e-12;
Cal_StoreInVariable(&v, Operation_P->Case.OptimizerUpdate.residual);
}
#else
......
......@@ -1421,6 +1421,7 @@ void Format_PostFooter(struct PostSubOperation *PSO_P, int Store)
exp.push_back(it->second[i]);
}
GetDPNumbers[CurrentName] = exp;
GetDPNumbersMap[CurrentName] = NodeTable;
if(PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No"))
Message::AddOnelabNumberChoice(PSO_P->SendToServer, exp, PSO_P->Color,
PSO_P->Units, PSO_P->Label, PSO_P->Visible,
......@@ -1448,6 +1449,7 @@ void Format_PostFooter(struct PostSubOperation *PSO_P, int Store)
exp.push_back(it->second[i]);
}
GetDPNumbers[CurrentName] = exp;
GetDPNumbersMap[CurrentName] = ElementTable;
if(PSO_P->SendToServer && strcmp(PSO_P->SendToServer, "No"))
Message::AddOnelabNumberChoice(PSO_P->SendToServer, exp, PSO_P->Color,
PSO_P->Units, PSO_P->Label, PSO_P->Visible,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment