1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 #include <cstdlib>
10 #include <iostream>
11 #include <fstream>
12 #include <string>
13 #include <vector>
14 #include <map>
15 #include <algorithm>
16 #include <cctype>
17
18 enum var_t { FS, P1, P2, P3, B, D, H, F0, E }; // order in Kuschel & Rackwitz
19
20
main(int argc,char ** argv)21 int main(int argc, char** argv)
22 {
23
24 std::ifstream fin(argv[1]);
25 if (!fin) {
26 std::cerr << "\nError: failure opening " << argv[1] << std::endl;
27 exit(-1);
28 }
29 size_t i, j, num_vars, num_fns, num_deriv_vars;
30 std::string vars_text, fns_text, dvv_text;
31
32 // define the std::string to enumeration map
33 std::map<std::string, var_t> var_t_map;
34 var_t_map["fs"] = FS; var_t_map["p1"] = P1; var_t_map["p2"] = P2;
35 var_t_map["p3"] = P3; var_t_map["b"] = B; var_t_map["d"] = D;
36 var_t_map["h"] = H; var_t_map["f0"] = F0; var_t_map["e"] = E;
37
38 // Get the parameter std::vector and ignore the labels
39 fin >> num_vars >> vars_text;
40 std::map<var_t, double> vars;
41 std::vector<var_t> labels(num_vars);
42 double var_i; std::string label_i; var_t v_i;
43 std::map<std::string, var_t>::iterator v_iter;
44 for (i=0; i<num_vars; i++) {
45 fin >> var_i >> label_i;
46 transform(label_i.begin(), label_i.end(), label_i.begin(),
47 (int(*)(int))tolower);
48 v_iter = var_t_map.find(label_i);
49 if (v_iter == var_t_map.end()) {
50 std::cerr << "Error: label \"" << label_i << "\" not supported in analysis "
51 << "driver." << std::endl;
52 exit(-1);
53 }
54 else
55 v_i = v_iter->second;
56 vars[v_i] = var_i;
57 labels[i] = v_i;
58 }
59
60 // Get the ASV std::vector and ignore the labels
61 fin >> num_fns >> fns_text;
62 std::vector<short> ASV(num_fns);
63 for (i=0; i<num_fns; i++) {
64 fin >> ASV[i];
65 fin.ignore(256, '\n');
66 }
67
68 // Get the DVV std::vector and ignore the labels
69 fin >> num_deriv_vars >> dvv_text;
70 std::vector<var_t> DVV(num_deriv_vars);
71 unsigned int dvv_i;
72 for (i=0; i<num_deriv_vars; i++) {
73 fin >> dvv_i;
74 fin.ignore(256, '\n');
75 DVV[i] = labels[dvv_i-1];
76 }
77
78 if (num_vars != 9 || num_fns != 1) {
79 std::cerr << "Error: wrong number of inputs/outputs in steel_column_perf."<<std::endl;
80 exit(-1);
81 }
82
83 // Compute the results and output them directly to argv[2] (the NO_FILTER
84 // option is used). Response tags are optional; output them for ease
85 // of results readability.
86 std::ofstream fout(argv[2]);
87 if (!fout) {
88 std::cerr << "\nError: failure creating " << argv[2] << std::endl;
89 exit(-1);
90 }
91 fout.precision(15); // 16 total digits
92 fout.setf(std::ios::scientific);
93 fout.setf(std::ios::right);
94
95 // In the steel column description in Kuschel & Rackwitz, Cost is _not_
96 // defined as a random variable. That is Cost is not a fn(B, D, H), but
97 // is rather defined as a fn(b, d, h). Since dCost/dX|_{X=mean} is not the
98 // same as dCost/dmean for non-normal X (jacobian_dX_dS is not 1), dCost/dX
99 // may not be used and an optional interface must be defined for Cost.
100
101 // set effective length s based on assumed boundary conditions
102 // actual length of the column is 7500 mm
103 double s = 7500;
104
105 // fs = yield stress (lognormal unc. var.)
106 // p1 = dead weight load (normal unc. var.)
107 // p2 = variable load (gumbel unc. var.)
108 // p3 = variable load (gumbel unc. var.)
109 // b = flange breadth (lognormal unc. var., mean is design var.)
110 // d = flange thickness (lognormal unc. var., mean is design var.)
111 // h = profile height (lognormal unc. var., mean is design var.)
112 // f0 = init. deflection (normal unc. var.)
113 // e = elastic modulus (weibull unc. var.)
114
115 double f0 = vars[F0], b = vars[B], d = vars[D], h = vars[H], fs = vars[FS],
116 e = vars[E], p = vars[P1]+vars[P2]+vars[P3], Pi = 3.1415926535897932385,
117 Pi2 = Pi*Pi, Pi4 = Pi2*Pi2, Pi6 = Pi2*Pi4, b2 = b*b, d2 = d*d, h2 = h*h,
118 h3 = h*h2, h5 = h2*h3, e2 = e*e, e3 = e*e2, s2 = s*s,
119 X = Pi2*e*b*d*h2 - 2.*s2*p, X2 = X*X, X3 = X*X2;
120
121 // **** g (limit state):
122 if (ASV[0] & 1)
123 fout << " "
124 << fs - p*(1./2./b/d + Pi2*f0*e*h/X) << " g\n";
125
126 // **** dg/dx (w.r.t. active/uncertain variables):
127 if (ASV[0] & 2) {
128 fout << "[ ";
129 for (i=0; i<num_deriv_vars; i++)
130 switch (DVV[i]) {
131 case F0: // df/df0
132 fout << -e*h*p*Pi2/X << ' ';
133 break;
134 case P1: case P2: case P3: // df/dp1, df/dp2, df/dp3
135 fout << -1./2./b/d - b*d*e2*f0*h3*Pi4/X2 << ' ';
136 break;
137 case B: // df/db
138 fout << p*(1./2./b2/d + d*e2*f0*h3*Pi4/X2) << ' ';
139 break;
140 case D: // df/dd
141 fout << p*(1./2./b/d2 + b*e2*f0*h3*Pi4/X2) << ' ';
142 break;
143 case H: // df/dh
144 fout << e*f0*p*Pi2*(X + 4.*p*s2)/X2 << ' ';
145 break;
146 case FS: // df/dfs
147 fout << "1. ";
148 break;
149 case E: // df/de
150 fout << 2.*f0*h*p*p*Pi2*s2/X2 << ' ';
151 break;
152 }
153 fout << "]\n";
154 }
155
156 // **** d^2g/dx^2: (SORM)
157 if (ASV[0] & 4) {
158 fout << "[[ ";
159 for (i=0; i<num_deriv_vars; i++)
160 for (j=0; j<num_deriv_vars; j++)
161 if (DVV[i] == FS || DVV[j] == FS) // d^2g/dfs^2
162 fout << "0. ";
163 else if ( (DVV[i] == P1 && DVV[j] == P1) ||
164 (DVV[i] == P1 && DVV[j] == P2) ||
165 (DVV[i] == P1 && DVV[j] == P3) ||
166 (DVV[i] == P2 && DVV[j] == P1) ||
167 (DVV[i] == P2 && DVV[j] == P2) ||
168 (DVV[i] == P2 && DVV[j] == P3) ||
169 (DVV[i] == P3 && DVV[j] == P1) ||
170 (DVV[i] == P3 && DVV[j] == P2) ||
171 (DVV[i] == P3 && DVV[j] == P3) ) // d^2g/dpdp
172 fout << -4.*b*d*e2*f0*h3*Pi4*s2/X3 << ' ';
173 else if ( (DVV[i] == P1 && DVV[j] == B) ||
174 (DVV[i] == P2 && DVV[j] == B) ||
175 (DVV[i] == P3 && DVV[j] == B) ||
176 (DVV[i] == B && DVV[j] == P1) ||
177 (DVV[i] == B && DVV[j] == P2) ||
178 (DVV[i] == B && DVV[j] == P3) ) // d^2g/dpdb
179 fout << 1./2./b2/d + d*e2*f0*h3*Pi4/X2*(2.*b*d*e*h2*Pi2/X - 1.) <<' ';
180 else if ( (DVV[i] == P1 && DVV[j] == D) ||
181 (DVV[i] == P2 && DVV[j] == D) ||
182 (DVV[i] == P3 && DVV[j] == D) ||
183 (DVV[i] == D && DVV[j] == P1) ||
184 (DVV[i] == D && DVV[j] == P2) ||
185 (DVV[i] == D && DVV[j] == P3) ) // d^2g/dpdd
186 fout << 1./2./b/d2 + b*e2*f0*h3*Pi4/X2*(2.*b*d*e*h2*Pi2/X - 1.) <<' ';
187 else if ( (DVV[i] == P1 && DVV[j] == H) ||
188 (DVV[i] == P2 && DVV[j] == H) ||
189 (DVV[i] == P3 && DVV[j] == H) ||
190 (DVV[i] == H && DVV[j] == P1) ||
191 (DVV[i] == H && DVV[j] == P2) ||
192 (DVV[i] == H && DVV[j] == P3) ) // d^2g/dpdh
193 fout << b*d*e2*f0*h2*Pi4*(X+8.*p*s2)/X3 << ' ';
194 else if ( (DVV[i] == P1 && DVV[j] == F0) ||
195 (DVV[i] == P2 && DVV[j] == F0) ||
196 (DVV[i] == P3 && DVV[j] == F0) ||
197 (DVV[i] == F0 && DVV[j] == P1) ||
198 (DVV[i] == F0 && DVV[j] == P2) ||
199 (DVV[i] == F0 && DVV[j] == P3) ) // d^2g/dpdf0
200 fout << -b*d*e2*h3*Pi4/X2 << ' ';
201 else if ( (DVV[i] == P1 && DVV[j] == E) ||
202 (DVV[i] == P2 && DVV[j] == E) ||
203 (DVV[i] == P3 && DVV[j] == E) ||
204 (DVV[i] == E && DVV[j] == P1) ||
205 (DVV[i] == E && DVV[j] == P2) ||
206 (DVV[i] == E && DVV[j] == P3) ) // d^2g/dpde
207 fout << 4.*b*d*e*f0*h3*p*Pi4*s2/X3 << ' ';
208 else if (DVV[i] == B && DVV[j] == B) // d^2g/db^2
209 fout << -p*(1./b2/b/d + 2.*d2*e3*f0*h5*Pi6/X3) << ' ';
210 else if ( (DVV[i] == B && DVV[j] == D) ||
211 (DVV[i] == D && DVV[j] == B) ) // d^2g/dbdd
212 fout << -p*(1./2./b2/d2 + e2*f0*h3*Pi4/X2*(2.*b*d*e*h2*Pi2/X - 1.))
213 << ' ';
214 else if ( (DVV[i] == B && DVV[j] == H) ||
215 (DVV[i] == H && DVV[j] == B) ) // d^2g/dbdh
216 fout << -d*e2*f0*h2*p*Pi4*(X + 8.*p*s2)/X3 << ' ';
217 else if ( (DVV[i] == F0 && DVV[j] == B) ||
218 (DVV[i] == B && DVV[j] == F0) ) // d^2g/dbdf0
219 fout << d*e2*h3*p*Pi4/X2 << ' ';
220 else if ( (DVV[i] == B && DVV[j] == E) ||
221 (DVV[i] == E && DVV[j] == B) ) // d^2g/dbde
222 fout << -4.*d*e*f0*h3*p*p*Pi4*s2/X3 << ' ';
223 else if (DVV[i] == D && DVV[j] == D) // d^2g/dd^2
224 fout << -p*(1./b/d2/d + 2.*b2*e3*f0*h5*Pi6/X3) << ' ';
225 else if ( (DVV[i] == D && DVV[j] == H) ||
226 (DVV[i] == H && DVV[j] == D) ) // d^2g/dddh
227 fout << -b*e2*f0*h2*p*Pi4*(X + 8.*p*s2)/X3 << ' ';
228 else if ( (DVV[i] == F0 && DVV[j] == D) ||
229 (DVV[i] == D && DVV[j] == F0) ) // d^2g/dddf0
230 fout << b*e2*h3*p*Pi4/X2 << ' ';
231 else if ( (DVV[i] == D && DVV[j] == E) ||
232 (DVV[i] == E && DVV[j] == D) ) // d^2g/ddde
233 fout << -4.*b*e*f0*h3*p*p*Pi4*s2/X3 << ' ';
234 else if (DVV[i] == H && DVV[j] == H) // d^2g/dh^2
235 fout << -2.*b*d*e2*f0*h*p*Pi4*(X + 8.*p*s2)/X3 << ' ';
236 else if ( (DVV[i] == F0 && DVV[j] == H) ||
237 (DVV[i] == H && DVV[j] == F0) ) // d^2g/dhdf0
238 fout << e*p*Pi2*(X + 4.*p*s2)/X2 << ' ';
239 else if ( (DVV[i] == H && DVV[j] == E) ||
240 (DVV[i] == E && DVV[j] == H) ) // d^2g/dhde
241 fout << -2.*f0*p*p*Pi2*s2*(3.*X + 8.*p*s2)/X3 << ' ';
242 else if (DVV[i] == F0 && DVV[j] == F0) // d^2g/df0^2
243 fout << "0. ";
244 else if ( (DVV[i] == F0 && DVV[j] == E) ||
245 (DVV[i] == E && DVV[j] == F0) ) // d^2g/df0de
246 fout << 2.*h*p*p*Pi2*s2/X2 << ' ';
247 else if (DVV[i] == E && DVV[j] == E) // d^2g/de^2
248 fout << -4.*b*d*f0*h3*p*p*Pi4*s2/X3 << ' ';
249 else { // unsupported derivative
250 std::cerr << "Error: unsupported Hessian cross term in steel_column."
251 << std::endl;
252 exit(-1);
253 }
254 fout << "]]\n";
255 }
256
257 fout.flush();
258 fout.close();
259 return 0;
260 }
261