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