1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19
20 #ifndef LIBMESH_SENSITIVITY_DATA_H
21 #define LIBMESH_SENSITIVITY_DATA_H
22
23
24 // Local Includes
25 #include "libmesh/libmesh_common.h"
26 #include "libmesh/parameter_vector.h"
27 #include "libmesh/system.h"
28
29 // C++ Includes
30 #include <vector>
31
32 namespace libMesh
33 {
34
35 // Forward declarations
36 class QoISet;
37
38 /**
39 * Data structure for holding completed parameter sensitivity
40 * calculations.
41 *
42 * \author Roy Stogner
43 * \date 2009
44 * \brief Holds completed parameter sensitivity calculations.
45 */
46 class SensitivityData
47 {
48 public:
49 class Row
50 {
51 public:
Row(SensitivityData & sd,unsigned int qoi)52 Row(SensitivityData & sd, unsigned int qoi) : _sd(sd), _qoi(qoi) {}
53
54 Number & operator[] (unsigned int parameter) { return _sd.derivative(_qoi, parameter); }
55 private:
56 SensitivityData & _sd;
57 unsigned int _qoi;
58 };
59
60 class ConstRow
61 {
62 public:
ConstRow(const SensitivityData & sd,unsigned int qoi)63 ConstRow(const SensitivityData & sd, unsigned int qoi) : _sd(sd), _qoi(qoi) {}
64
65 const Number & operator[] (unsigned int parameter) { return _sd.derivative(_qoi, parameter); }
66 private:
67 const SensitivityData & _sd;
68 unsigned int _qoi;
69 };
70
71 /**
72 * Default constructor: empty data set
73 */
SensitivityData()74 SensitivityData() {}
75
76 /**
77 * Constructor from QoISet and ParameterVector: allocates space
78 * for all required sensitivities
79 */
80 SensitivityData(const QoISet & qoi_indices,
81 const System & sys,
82 const ParameterVector & parameter_vector);
83
84 /**
85 * Clears and deallocates all data
86 */
clear()87 void clear() { _grad_data.clear(); }
88
89 /**
90 * Given QoISet and ParameterVector, allocates space
91 * for all required first derivative data
92 */
93 void allocate_data(const QoISet & qoi_indices,
94 const System & sys,
95 const ParameterVector & parameter_vector);
96
97 /**
98 * Given QoISet and ParameterVector, allocates space
99 * for all required second derivative data
100 */
101 void allocate_hessian_data(const QoISet & qoi_indices,
102 const System & sys,
103 const ParameterVector & parameter_vector);
104
105 /**
106 * \returns The parameter sensitivity derivative for the specified
107 * quantity of interest for the specified parameter
108 */
109 const Number & derivative (unsigned int qoi_index,
110 unsigned int parameter_index) const;
111
112 /**
113 * \returns The parameter sensitivity derivative for the specified
114 * quantity of interest for the specified pair of parameters
115 */
116 const Number & second_derivative (unsigned int qoi_index,
117 unsigned int parameter_index1,
118 unsigned int parameter_index2) const;
119
120 /**
121 * Gets/sets the parameter sensitivity derivative for the specified
122 * quantity of interest for the specified parameter
123 */
124 Number & derivative (unsigned int qoi_index,
125 unsigned int parameter_index);
126
127 /**
128 * Gets/sets the parameter sensitivity second derivative for the
129 * specified quantity of interest for the specified pair of
130 * parameters
131 */
132 Number & second_derivative (unsigned int qoi_index,
133 unsigned int parameter_index1,
134 unsigned int parameter_index2);
135
136 /**
137 * Vector address type operator: sd[q][p] is an alias for
138 * sd.derivative(q,p)
139 */
140 ConstRow operator[] (unsigned int qoi) const { return ConstRow(*this, qoi); }
141
142 Row operator[] (unsigned int qoi) { return Row(*this, qoi); }
143
144 private:
145 /**
146 * Data storage; currently pretty trivial
147 */
148 std::vector<std::vector<Number>> _grad_data;
149 std::vector<std::vector<std::vector<Number>>> _hess_data;
150 };
151
152
153
154 // ------------------------------------------------------------
155 // SensitivityData inline methods
156
157
158
159 inline
SensitivityData(const QoISet & qoi_indices,const System & sys,const ParameterVector & parameter_vector)160 SensitivityData::SensitivityData(const QoISet & qoi_indices,
161 const System & sys,
162 const ParameterVector & parameter_vector)
163 {
164 this->allocate_data(qoi_indices, sys, parameter_vector);
165 }
166
167
168
169 inline
allocate_data(const QoISet & qoi_indices,const System & sys,const ParameterVector & parameter_vector)170 void SensitivityData::allocate_data(const QoISet & qoi_indices,
171 const System & sys,
172 const ParameterVector & parameter_vector)
173 {
174 const std::size_t Np = parameter_vector.size();
175 const unsigned int Nq = sys.n_qois();
176
177 if (_grad_data.size() < Nq)
178 _grad_data.resize(Nq);
179
180 for (unsigned int i=0; i != Nq; ++i)
181 if (qoi_indices.has_index(i))
182 {
183 _grad_data[i].clear();
184 _grad_data[i].resize(Np);
185 }
186 }
187
188
189
190 inline
allocate_hessian_data(const QoISet & qoi_indices,const System & sys,const ParameterVector & parameter_vector)191 void SensitivityData::allocate_hessian_data(const QoISet & qoi_indices,
192 const System & sys,
193 const ParameterVector & parameter_vector)
194 {
195 const std::size_t Np = parameter_vector.size();
196 const unsigned int Nq = sys.n_qois();
197
198 if (_hess_data.size() < Nq)
199 _hess_data.resize(Nq);
200
201 for (unsigned int i=0; i != Nq; ++i)
202 if (qoi_indices.has_index(i))
203 {
204 _hess_data[i].clear();
205 _hess_data[i].resize(Np);
206 for (std::size_t j=0; j != Np; ++j)
207 _hess_data[i][j].resize(Np);
208 }
209 }
210
211
212
213 inline
derivative(unsigned int qoi_index,unsigned int parameter_index)214 const Number & SensitivityData::derivative(unsigned int qoi_index,
215 unsigned int parameter_index) const
216 {
217 libmesh_assert_less (qoi_index, _grad_data.size());
218 libmesh_assert_less (parameter_index, _grad_data[qoi_index].size());
219
220 return _grad_data[qoi_index][parameter_index];
221 }
222
223
224
225 inline
derivative(unsigned int qoi_index,unsigned int parameter_index)226 Number & SensitivityData::derivative(unsigned int qoi_index,
227 unsigned int parameter_index)
228 {
229 libmesh_assert_less (qoi_index, _grad_data.size());
230 libmesh_assert_less (parameter_index, _grad_data[qoi_index].size());
231
232 return _grad_data[qoi_index][parameter_index];
233 }
234
235
236
237 inline
second_derivative(unsigned int qoi_index,unsigned int parameter_index1,unsigned int parameter_index2)238 const Number & SensitivityData::second_derivative(unsigned int qoi_index,
239 unsigned int parameter_index1,
240 unsigned int parameter_index2) const
241 {
242 libmesh_assert_less (qoi_index, _hess_data.size());
243 libmesh_assert_less (parameter_index1, _hess_data[qoi_index].size());
244 libmesh_assert_less (parameter_index2, _hess_data[qoi_index][parameter_index1].size());
245
246 return _hess_data[qoi_index][parameter_index1][parameter_index2];
247 }
248
249
250
251 inline
second_derivative(unsigned int qoi_index,unsigned int parameter_index1,unsigned int parameter_index2)252 Number & SensitivityData::second_derivative(unsigned int qoi_index,
253 unsigned int parameter_index1,
254 unsigned int parameter_index2)
255 {
256 libmesh_assert_less (qoi_index, _hess_data.size());
257 libmesh_assert_less (parameter_index1, _hess_data[qoi_index].size());
258 libmesh_assert_less (parameter_index2, _hess_data[qoi_index][parameter_index1].size());
259
260 return _hess_data[qoi_index][parameter_index1][parameter_index2];
261 }
262
263 } // namespace libMesh
264
265 #endif // LIBMESH_SENSITIVITY_DATA_H
266