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 { B, H, P, M, Y };
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["b"] = B; var_t_map["h"] = H;
35 var_t_map["p"] = P; var_t_map["m"] = M; var_t_map["y"] = Y;
36
37 // Get the parameter std::vector and ignore the labels
38 fin >> num_vars >> vars_text;
39 std::map<var_t, double> vars;
40 std::vector<var_t> labels(num_vars);
41 double var_i; std::string label_i; var_t v_i;
42 std::map<std::string, var_t>::iterator v_iter;
43 for (i=0; i<num_vars; i++) {
44 fin >> var_i >> label_i;
45 transform(label_i.begin(), label_i.end(), label_i.begin(),
46 (int(*)(int))tolower);
47 v_iter = var_t_map.find(label_i);
48 if (v_iter == var_t_map.end()) {
49 std::cerr << "Error: label \"" << label_i << "\" not supported in analysis "
50 << "driver." << std::endl;
51 exit(-1);
52 }
53 else
54 v_i = v_iter->second;
55 vars[v_i] = var_i;
56 labels[i] = v_i;
57 }
58
59 // Get the ASV std::vector and ignore the labels
60 fin >> num_fns >> fns_text;
61 std::vector<short> ASV(num_fns);
62 for (i=0; i<num_fns; i++) {
63 fin >> ASV[i];
64 fin.ignore(256, '\n');
65 }
66
67 // Get the DVV std::vector and ignore the labels
68 fin >> num_deriv_vars >> dvv_text;
69 std::vector<var_t> DVV(num_deriv_vars);
70 unsigned int dvv_i;
71 for (i=0; i<num_deriv_vars; i++) {
72 fin >> dvv_i;
73 fin.ignore(256, '\n');
74 DVV[i] = labels[dvv_i-1];
75 }
76
77 if (num_vars != 5 || num_fns != 2) {
78 std::cerr << "Error: wrong number of inputs/outputs in short_column." << std::endl;
79 exit(-1);
80 }
81
82 // Compute the results and output them directly to argv[2] (the NO_FILTER
83 // option is used). Response tags are optional; output them for ease
84 // of results readability.
85 std::ofstream fout(argv[2]);
86 if (!fout) {
87 std::cerr << "\nError: failure creating " << argv[2] << std::endl;
88 exit(-1);
89 }
90 fout.precision(15); // 16 total digits
91 fout.setf(std::ios::scientific);
92 fout.setf(std::ios::right);
93
94 // b = vars[B] = column base (design var.)
95 // h = vars[H] = column height (design var.)
96 // p = vars[P] (normal uncertain var.)
97 // m = vars[M] (normal uncertain var.)
98 // y = vars[Y] (lognormal uncertain var.)
99 double b = vars[B], h = vars[H], p = vars[P], m = vars[M], y = vars[Y],
100 b_sq = b*b, h_sq = h*h, p_sq = p*p, y_sq = y*y;
101
102 // **** f (objective = bh = cross sectional area):
103 if (ASV[0] & 1)
104 fout << " " << b*h << " f\n";
105
106 // **** g (limit state = short column response):
107 if (ASV[1] & 1)
108 fout << " "
109 << 1. - 4.*m/(b*h_sq*y) - p_sq/(b_sq*h_sq*y_sq) << " g\n";
110
111 // **** df/dx (w.r.t. active/uncertain variables):
112 if (ASV[0] & 2) {
113 fout << "[ ";
114 for (i=0; i<num_deriv_vars; i++)
115 switch (DVV[i]) {
116 case B: // design variable derivative
117 fout << h << ' ';
118 break;
119 case H: // design variable derivative
120 fout << b << ' ';
121 break;
122 default: // uncertain variable derivative
123 fout << "0. ";
124 break;
125 }
126 fout << "]\n";
127 }
128
129 // **** dg/dx (w.r.t. active/uncertain variables):
130 if (ASV[1] & 2) {
131 fout << "[ ";
132 for (i=0; i<num_deriv_vars; i++)
133 switch (DVV[i]) {
134 case B: // design variable derivative
135 fout << 4.*m/(b_sq*h_sq*y) + 2.*p_sq/(b_sq*b*h_sq*y_sq) << ' ';
136 break;
137 case H: // design variable derivative
138 fout << 8.*m/(b*h_sq*h*y) + 2.*p_sq/(b_sq*h_sq*h*y_sq) << ' ';
139 break;
140 case P: // uncertain variable derivative
141 fout << -2.*p/(b_sq*h_sq*y_sq) << ' ';
142 break;
143 case M: // uncertain variable derivative
144 fout << -4./(b*h_sq*y) << ' ';
145 break;
146 case Y: // uncertain variable derivative
147 fout << 4.*m/(b*h_sq*y_sq) + 2.*p_sq/(b_sq*h_sq*y_sq*y) << ' ';
148 break;
149 }
150 fout << "]\n";
151 }
152
153 // **** d^2f/dx^2: (SORM)
154 if (ASV[0] & 4) {
155 fout << "[[ ";
156 for (i=0; i<num_deriv_vars; i++)
157 for (j=0; j<num_deriv_vars; j++)
158 if ( (DVV[i] == B && DVV[j] == H) || (DVV[i] == H && DVV[j] == B) )
159 fout << "1. ";
160 else
161 fout << "0. ";
162 fout << "]]\n";
163 }
164
165 // **** d^2g/dx^2: (SORM)
166 if (ASV[1] & 4) {
167 fout << "[[ ";
168 for (i=0; i<num_deriv_vars; i++)
169 for (j=0; j<num_deriv_vars; j++)
170 if (DVV[i] == B && DVV[j] == B) // d^2g/db^2
171 fout << -8.*m/(b_sq*b*h_sq*y) - 6.*p_sq/(b_sq*b_sq*h_sq*y_sq) << ' ';
172 else if ( (DVV[i] == B && DVV[j] == H) ||
173 (DVV[i] == H && DVV[j] == B) ) // d^2g/dbdh
174 fout << -8.*m/(b_sq*h_sq*h*y) - 4.*p_sq/(b_sq*b*h_sq*h*y_sq) << ' ';
175 else if (DVV[i] == H && DVV[j] == H) // d^2g/dh^2
176 fout << -24.*m/(b*h_sq*h_sq*y) - 6.*p_sq/(b_sq*h_sq*h_sq*y_sq) << ' ';
177 else if (DVV[i] == P && DVV[j] == P) // d^2g/dp^2
178 fout << -2./(b_sq*h_sq*y_sq) << ' ';
179 else if ( (DVV[i] == P && DVV[j] == M) ||
180 (DVV[i] == M && DVV[j] == P) ) // d^2g/dpdm
181 fout << "0. ";
182 else if ( (DVV[i] == P && DVV[j] == Y) ||
183 (DVV[i] == Y && DVV[j] == P) ) // d^2g/dpdy
184 fout << 4.*p/(b_sq*h_sq*y_sq*y) << ' ';
185 else if (DVV[i] == M && DVV[j] == M) // d^2g/dm^2
186 fout << "0. ";
187 else if ( (DVV[i] == M && DVV[j] == Y) ||
188 (DVV[i] == Y && DVV[j] == M) ) // d^2g/dmdy
189 fout << 4./(b*h_sq*y_sq) << ' ';
190 else if (DVV[i] == Y && DVV[j] == Y) // d^2g/dy^2
191 fout << -8.*m/(b*h_sq*y_sq*y) - 6.*p_sq/(b_sq*h_sq*y_sq*y_sq) << ' ';
192 else { // unsupported cross-derivative
193 std::cerr << "Error: unsupported Hessian cross term in short_column."
194 << std::endl;
195 exit(-1);
196 }
197 fout << "]]\n";
198 }
199
200 fout.flush();
201 fout.close();
202 return 0;
203 }
204