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 <cmath>
11 #include <iostream>
12 #include <fstream>
13 #include <string>
14 #include <vector>
15 #include <map>
16 #include <algorithm>
17 #include <cctype>
18
19 enum var_t { W, T, R, E, X, Y };
20
21
main(int argc,char ** argv)22 int main(int argc, char** argv)
23 {
24
25 // This test problem is an OUU example from Applied Research Associates
26 // (42nd AIAA SDM conference, April 2001).
27
28 std::ifstream fin(argv[1]);
29 if (!fin) {
30 std::cerr << "\nError: failure opening " << argv[1] << std::endl;
31 exit(-1);
32 }
33 size_t i, j, num_vars, num_fns, num_deriv_vars;
34 std::string vars_text, fns_text, dvv_text;
35
36 // define the std::string to enumeration map
37 std::map<std::string, var_t> var_t_map;
38 var_t_map["w"] = W; var_t_map["t"] = T;
39 var_t_map["r"] = R; var_t_map["e"] = E;
40 var_t_map["x"] = X; var_t_map["y"] = Y;
41
42 // Get the parameter std::vector and ignore the labels
43 fin >> num_vars >> vars_text;
44 std::map<var_t, double> vars;
45 std::vector<var_t> labels(num_vars);
46 double var_i; std::string label_i; var_t v_i;
47 std::map<std::string, var_t>::iterator v_iter;
48 for (i=0; i<num_vars; i++) {
49 fin >> var_i >> label_i;
50 transform(label_i.begin(), label_i.end(), label_i.begin(),
51 (int(*)(int))tolower);
52 v_iter = var_t_map.find(label_i);
53 if (v_iter == var_t_map.end()) {
54 std::cerr << "Error: label \"" << label_i << "\" not supported in analysis "
55 << "driver." << std::endl;
56 exit(-1);
57 }
58 else
59 v_i = v_iter->second;
60 vars[v_i] = var_i;
61 labels[i] = v_i;
62 }
63
64 // Get the ASV std::vector and ignore the labels
65 fin >> num_fns >> fns_text;
66 std::vector<short> ASV(num_fns);
67 for (i=0; i<num_fns; i++) {
68 fin >> ASV[i];
69 fin.ignore(256, '\n');
70 }
71
72 // Get the DVV std::vector and ignore the labels
73 fin >> num_deriv_vars >> dvv_text;
74 std::vector<var_t> DVV(num_deriv_vars);
75 unsigned int dvv_i;
76 for (i=0; i<num_deriv_vars; i++) {
77 fin >> dvv_i;
78 fin.ignore(256, '\n');
79 DVV[i] = labels[dvv_i-1];
80 }
81
82 if (num_vars != 4 && num_vars != 6) {
83 std::cerr << "Error: Wrong number of variables in cantilever test fn." << std::endl;
84 exit(-1);
85 }
86 if (num_fns < 2 || num_fns > 3) {
87 std::cerr << "Error: wrong number of response functions in cantilever test fn."
88 << std::endl;
89 exit(-1);
90 }
91
92 // Compute the cross-sectional area, stress, and displacement of the
93 // cantilever beam. This simulator is unusual in that it supports both
94 // the case of design variable insertion and the case of design variable
95 // augmentation. It does not support mixed insertion/augmentation. In
96 // the 6 variable case, w,t,R,E,X,Y are all passed in; in the 4 variable
97 // case, w,t assume local values.
98 std::map<var_t, double>::iterator m_iter = vars.find(W);
99 double w = (m_iter == vars.end()) ? 2.5 : m_iter->second; // beam width
100 m_iter = vars.find(T);
101 double t = (m_iter == vars.end()) ? 2.5 : m_iter->second; // beam thickness
102 double r = vars[R]; // yield strength
103 double e = vars[E]; // Young's modulus
104 double x = vars[X]; // horizontal load
105 double y = vars[Y]; // vertical load
106
107 // allow f,c1,c2 (optimization) or just c1,c2 (calibration)
108 bool objective; size_t c1i, c2i;
109 if (num_fns == 2) { objective = false; c1i = 0; c2i = 1; }
110 else { objective = true; c1i = 1; c2i = 2; }
111
112 // optimization inequality constraint: <= 0 and scaled O(1)
113 //Real g_stress = stress/R - 1.0;
114 //Real g_disp = disp/D0 - 1.0;
115
116 // Compute the results and output them directly to argv[2] (the NO_FILTER
117 // option is used). Response tags are optional; output them for ease of
118 // results readability.
119 double D0 = 2.2535, L = 100., area = w*t, w_sq = w*w, t_sq = t*t,
120 r_sq = r*r, x_sq = x*x, y_sq = y*y;
121 double stress = 600.*y/w/t_sq + 600.*x/w_sq/t;
122 double D1 = 4.*pow(L,3)/e/area, D2 = pow(y/t_sq, 2)+pow(x/w_sq, 2);
123 double D3 = D1/sqrt(D2)/D0, D4 = D1*sqrt(D2)/D0;
124
125 std::ofstream fout(argv[2]); // do not instantiate until ready to write results
126 if (!fout) {
127 std::cerr << "\nError: failure creating " << argv[2] << std::endl;
128 exit(-1);
129 }
130 fout.precision(15); // 16 total digits
131 fout.setf(std::ios::scientific);
132 fout.setf(std::ios::right);
133
134 // **** f:
135 if (objective && (ASV[0] & 1))
136 fout << " " << area << '\n';
137
138 // **** c1:
139 if (ASV[c1i] & 1)
140 fout << " " << stress/r - 1. << '\n';
141
142 // **** c2:
143 if (ASV[c2i] & 1)
144 fout << " " << D4 - 1. << '\n';
145
146 // **** df/dx:
147 if (objective && (ASV[0] & 2)) {
148 fout << "[ ";
149 for (i=0; i<num_deriv_vars; i++)
150 switch (DVV[i]) {
151 case W: fout << t << ' '; break; // design var derivative
152 case T: fout << w << ' '; break; // design var derivative
153 default: fout << "0. "; break; // uncertain var derivative
154 }
155 fout << "]\n";
156 }
157
158 // **** dc1/dx:
159 if (ASV[c1i] & 2) {
160 fout << "[ ";
161 for (i=0; i<num_deriv_vars; i++)
162 switch (DVV[i]) {
163 case W: fout << -600.*(y/t + 2.*x/w)/w_sq/t/r << ' '; break; // des var
164 case T: fout << -600.*(2.*y/t + x/w)/w/t_sq/r << ' '; break; // des var
165 case R: fout << -stress/r_sq << ' '; break; // uncertain var deriv
166 case E: fout << "0. "; break; // uncertain var deriv
167 case X: fout << 600./w_sq/t/r << ' '; break; // uncertain var deriv
168 case Y: fout << 600./w/t_sq/r << ' '; break; // uncertain var deriv
169 }
170 fout << "]\n";
171 }
172
173 // **** dc2/dx:
174 if (ASV[c2i] & 2) {
175 fout << "[ ";
176 for (i=0; i<num_deriv_vars; i++)
177 switch (DVV[i]) {
178 case W: fout << -D3*2.*x_sq/w_sq/w_sq/w - D4/w << ' '; break; // des var
179 case T: fout << -D3*2.*y_sq/t_sq/t_sq/t - D4/t << ' '; break; // des var
180 case R: fout << "0. "; break; // unc var deriv
181 case E: fout << -D4/e << ' '; break; // unc var deriv
182 case X: fout << D3*x/w_sq/w_sq << ' '; break; // unc var deriv
183 case Y: fout << D3*y/t_sq/t_sq << ' '; break; // unc var deriv
184 }
185 fout << "]\n";
186 }
187
188 // **** d^2f/dx^2:
189 if (objective && (ASV[0] & 4)) {
190 fout << "[[ ";
191 for (i=0; i<num_deriv_vars; i++)
192 for (j=0; j<num_deriv_vars; j++)
193 if ( (DVV[i] == W && DVV[j] == T) || (DVV[i] == T && DVV[j] == W) )
194 fout << "1. ";
195 else
196 fout << "0. ";
197 fout << "]]\n";
198 }
199
200 // **** d^2c1/dx^2:
201 if (ASV[c1i] & 4) {
202 fout << "[[ ";
203 for (i=0; i<num_deriv_vars; i++)
204 for (j=0; j<num_deriv_vars; j++)
205 if (DVV[i] == W && DVV[j] == W) // d^2g/dw^2
206 fout << 1200.*(y/t + 3.*x/w)/w_sq/area/r << ' ';
207 else if (DVV[i] == T && DVV[j] == T) // d^2g/dt^2
208 fout << 1200.*(3.*y/t + x/w)/t_sq/area/r << ' ';
209 else if (DVV[i] == R && DVV[j] == R) // d^2g/dr^2
210 fout << 2.*stress/pow(r, 3) << ' ';
211 else if ( (DVV[i] == W && DVV[j] == T) ||
212 (DVV[i] == T && DVV[j] == W) ) // d^2g/dwdt
213 fout << 1200.*(y/t + x/w)/w_sq/t_sq/r << ' ';
214 else if ( (DVV[i] == W && DVV[j] == R) ||
215 (DVV[i] == R && DVV[j] == W) ) // d^2g/dwdr
216 fout << 600.*(y/t + 2.*x/w)/w_sq/t/r_sq << ' ';
217 else if ( (DVV[i] == W && DVV[j] == X) ||
218 (DVV[i] == X && DVV[j] == W) ) // d^2g/dwdx
219 fout << -1200./w_sq/w/t/r << ' ';
220 else if ( (DVV[i] == W && DVV[j] == Y) ||
221 (DVV[i] == Y && DVV[j] == W) ) // d^2g/dwdy
222 fout << -600./w_sq/t_sq/r << ' ';
223 else if ( (DVV[i] == T && DVV[j] == R) ||
224 (DVV[i] == R && DVV[j] == T) ) // d^2g/dtdr
225 fout << 600.*(2.*y/t + x/w)/w/t_sq/r_sq << ' ';
226 else if ( (DVV[i] == T && DVV[j] == X) ||
227 (DVV[i] == X && DVV[j] == T) ) // d^2g/dtdx
228 fout << -600./w_sq/t_sq/r << ' ';
229 else if ( (DVV[i] == T && DVV[j] == Y) ||
230 (DVV[i] == Y && DVV[j] == T) ) // d^2g/dtdy
231 fout << -1200./w/t_sq/t/r << ' ';
232 else if ( (DVV[i] == R && DVV[j] == X) ||
233 (DVV[i] == X && DVV[j] == R) ) // d^2g/drdx
234 fout << -600./w_sq/t/r_sq << ' ';
235 else if ( (DVV[i] == R && DVV[j] == Y) ||
236 (DVV[i] == Y && DVV[j] == R) ) // d^2g/drdy
237 fout << -600./w/t_sq/r_sq << ' ';
238 else
239 fout << "0. ";
240 fout << "]]\n";
241 }
242
243 // **** d^2c2/dx^2:
244 if (ASV[c2i] & 4) {
245 double D5 = 1./sqrt(D2)/D0, D6 = -D1/2./D0/pow(D2,1.5);
246 double D7 = sqrt(D2)/D0, D8 = D1/2./D0/sqrt(D2);
247 double dD2_dx = 2.*x/w_sq/w_sq, dD3_dx = D6*dD2_dx, dD4_dx = D8*dD2_dx;
248 double dD2_dy = 2.*y/t_sq/t_sq, dD3_dy = D6*dD2_dy, dD4_dy = D8*dD2_dy;
249 double dD1_dw = -D1/w, dD2_dw = -4.*x_sq/w_sq/w_sq/w,
250 dD3_dw = D5*dD1_dw + D6*dD2_dw, dD4_dw = D7*dD1_dw + D8*dD2_dw;
251 double dD1_dt = -D1/t, dD2_dt = -4.*y_sq/t_sq/t_sq/t,
252 dD3_dt = D5*dD1_dt + D6*dD2_dt, dD4_dt = D7*dD1_dt + D8*dD2_dt;
253 fout << "[[ ";
254 for (i=0; i<num_deriv_vars; i++)
255 for (j=0; j<num_deriv_vars; j++)
256 if (DVV[i] == W && DVV[j] == W) // d^2g/dw^2
257 fout << D3*10.*x_sq/pow(w_sq,3)
258 - 2.*x_sq/w_sq/w_sq/w*dD3_dw + D4/w_sq - dD4_dw/w << ' ';
259 else if (DVV[i] == T && DVV[j] == T) // d^2g/dt^2
260 fout << D3*10.*y_sq/pow(t_sq,3)
261 - 2.*y_sq/t_sq/t_sq/t*dD3_dt + D4/t_sq - dD4_dt/t << ' ';
262 else if (DVV[i] == E && DVV[j] == E) { // d^2g/de^2
263 double dD1_de = -D1/e, dD4_de = D7*dD1_de;
264 fout << D4/e/e - dD4_de/e << ' ';
265 }
266 else if (DVV[i] == X && DVV[j] == X) // d^2g/dx^2
267 fout << D3/w_sq/w_sq + x/w_sq/w_sq*dD3_dx << ' ';
268 else if (DVV[i] == Y && DVV[j] == Y) // d^2g/dy^2
269 fout << D3/t_sq/t_sq + y/t_sq/t_sq*dD3_dy << ' ';
270 else if ( (DVV[i] == W && DVV[j] == T) ||
271 (DVV[i] == T && DVV[j] == W) ) // d^2g/dwdt
272 fout << -2.*x_sq/w_sq/w_sq/w*dD3_dt - dD4_dt/w << ' ';
273 else if ( (DVV[i] == W && DVV[j] == E) ||
274 (DVV[i] == E && DVV[j] == W) ) // d^2g/dwde
275 fout << -dD4_dw/e << ' ';
276 else if ( (DVV[i] == W && DVV[j] == X) ||
277 (DVV[i] == X && DVV[j] == W) ) // d^2g/dwdx
278 fout << -4.*x*D3/w_sq/w_sq/w + x/w_sq/w_sq*dD3_dw << ' ';
279 else if ( (DVV[i] == W && DVV[j] == Y) ||
280 (DVV[i] == Y && DVV[j] == W) ) // d^2g/dwdy
281 fout << y/t_sq/t_sq*dD3_dw << ' ';
282 else if ( (DVV[i] == T && DVV[j] == E) ||
283 (DVV[i] == E && DVV[j] == T) ) // d^2g/dtde
284 fout << -dD4_dt/e << ' ';
285 else if ( (DVV[i] == T && DVV[j] == X) ||
286 (DVV[i] == X && DVV[j] == T) ) // d^2g/dtdx
287 fout << x/w_sq/w_sq*dD3_dt << ' ';
288 else if ( (DVV[i] == T && DVV[j] == Y) ||
289 (DVV[i] == Y && DVV[j] == T) ) // d^2g/dtdy
290 fout << -4.*y*D3/t_sq/t_sq/t + y/t_sq/t_sq*dD3_dt << ' ';
291 else if ( (DVV[i] == E && DVV[j] == X) ||
292 (DVV[i] == X && DVV[j] == E) ) // d^2g/dedx
293 fout << -dD4_dx/e << ' ';
294 else if ( (DVV[i] == E && DVV[j] == Y) ||
295 (DVV[i] == Y && DVV[j] == E) ) // d^2g/dedy
296 fout << -dD4_dy/e << ' ';
297 else if ( (DVV[i] == X && DVV[j] == Y) ||
298 (DVV[i] == Y && DVV[j] == X) ) // d^2g/dxdy
299 fout << x/w_sq/w_sq*dD3_dy << ' ';
300 else
301 fout << "0. ";
302 fout << "]]\n";
303 }
304
305 fout.flush();
306 fout.close();
307 return 0;
308 }
309