b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2taylor_expansions_buckling_solver.H
1//------------------------------------------------------------------------
2// b2taylor_expansions_buckling_solver.H --
3//
4//
5// written by Mathias Doreille
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7// Thomas Blome <thomas.blome@dlr.de>
8//
9// (c) 2012-2013 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef B2LINEARIZE_PREBUCKLING_SOLVER_H_
22#define B2LINEARIZE_PREBUCKLING_SOLVER_H_
23
24#include <iostream>
25
26#include "b2ppconfig.h"
27#include "b2static_nonlinear_utile.H"
28#include "model/b2solution.H"
29#include "model/b2solver.H"
30#include "utils/b2linear_algebra.H"
31#include "utils/b2logging.H"
32#include "utils/b2nonlinear_system_solver.H"
33#include "utils/b2numerical_differentiation.H"
34#include "utils/b2object.H"
35
36#define NEW_METHOD0
37#define NEW_METHOD1
38
39namespace b2000::solver {
40
41template <typename T, typename MATRIX_FORMAT>
42class TaylorExpansionsBucklingSolver : public Solver {
43public:
44 using type_t = ObjectTypeComplete<TaylorExpansionsBucklingSolver, Solver::type_t>;
45
46 void solve() override {
47 while (solve_iteration()) { ; }
48 }
49
50 bool solve_iteration() override {
51 if (model_ == nullptr) { Exception() << THROW; }
52 if (case_iterator_.get() == nullptr) {
53 case_iterator_ = model_->get_case_list().get_case_iterator();
54 }
55 case_ = case_iterator_->next();
56 if (case_) {
57 if (case_->get_number_of_stage() != 1) {
58 logging::get_logger("solver.linearised_prebuckling")
59 << logging::warning
60 << "The Taylor expansion buckling solver cannot hand a case with multiple "
61 "stages. The extra stages are ignored."
62 << logging::LOGFLUSH;
63 }
64 solve_on_case(*case_);
65 return true;
66 }
67 return false;
68 }
69
70 static type_t type;
71
72protected:
73 void solve_on_case(Case& case_);
74
75 void solve_nonlinear_refinement(
76 Case& case_, DefaultStaticResidueFunction<T, MATRIX_FORMAT>& residue_function,
77 b2linalg::Vector<double, b2linalg::Vdense_ref> value, double& time,
78 b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector,
79 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K);
80
81 void solve_nonlinear_refinement_sha2(
82 Case& case_, DefaultStaticResidueFunction<T, MATRIX_FORMAT>& residue_function,
83 b2linalg::Vector<double, b2linalg::Vdense_ref> value, double& time,
84 b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector,
85 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K) {
86 NonLinearSystem nls(case_, residue_function, K);
87 b2linalg::Vector<T> x;
88 nls.init(value, time, eigenvector, x);
89 Sha2NonlinearSolver solver(5);
90 solver.solve(nls, x);
91 nls.set_result(x, value, time, eigenvector);
92 }
93
94 class NonLinearSystem {
95 public:
96 typedef b2linalg::Vector<T> x_type;
97 typedef b2linalg::Vector<T> r_type;
98
99 NonLinearSystem(
100 Case& case_, DefaultStaticResidueFunction<T, MATRIX_FORMAT>& residue_function_,
101 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_)
102 : residue_function(residue_function_), K(K_) {
103 approximation_order = case_.get_int(
104 "APPROXIMATION_ORDER_REFINEMENT", case_.get_int("APPROXIMATION_ORDER", 1));
105 h = std::pow(1e-15, 1.0 / (approximation_order + 1))
106 * case_.get_double("H_SCALE_REFINEMENT", case_.get_double("H_SCALE", 1.0));
107 get_forward_coefficient(1, 0, approximation_order, diff_coeff);
108 eps = 1e-5;
109 use_constraint_with_time = false;
110
111 s = residue_function.get_number_of_dof();
112 residue.resize(s);
113 d_residue_d_time.resize(s);
114 Z.set_same_structure(K);
115 TMP.set_same_structure(K);
116 tmp.resize(s);
117 qfp.resize(s, 2);
118 phifp.resize(s, 2);
119 }
120
121 void init(
122 b2linalg::Vector<double, b2linalg::Vdense_ref> value, double& time,
123 b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector_,
124 b2linalg::Vector<double>& x) {
125 x.resize(2 * s + 1);
126 x(b2linalg::Interval(0, s)) = value;
127 x(b2linalg::Interval(s, 2 * s)) = eigenvector_;
128 x[2 * s] = time;
129
130 b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector =
131 x(b2linalg::Interval(s, 2 * s));
132 eigenvector *=
133 1.0 / (eigenvector.norm2() * (use_constraint_with_time ? std::sqrt(time) : 1));
134 }
135
136 void set_result(
137 b2linalg::Vector<double>& x, b2linalg::Vector<double, b2linalg::Vdense_ref> value,
138 double& time, b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector_) {
139 value = x(b2linalg::Interval(0, s));
140 eigenvector_ = x(b2linalg::Interval(s, 2 * s));
141 time = x[2 * s];
142 }
143
144 void operator()(const x_type& x, r_type& r, bool update_der = false) {
145 EquilibriumSolution esf(false);
146 b2linalg::Vector<double, b2linalg::Vdense_constref> value = x(b2linalg::Interval(0, s));
147 b2linalg::Vector<double, b2linalg::Vdense_constref> eigenvector =
148 x(b2linalg::Interval(s, 2 * s));
149 double time = x[2 * s];
150
151 r.resize(x.size());
152 b2linalg::Vector<double, b2linalg::Vdense_ref> residue = r(b2linalg::Interval(0, s));
153 b2linalg::Vector<double, b2linalg::Vdense_ref> Ke = r(b2linalg::Interval(s, 2 * s));
154
155 residue.set_zero();
156 K.set_zero_same_struct();
157 d_residue_d_time.set_zero();
158 if (update_der) {
159 residue_function.get_residue(
160 value, time, 0, esf, residue, K, d_residue_d_time, 0, 0, 0);
161 } else {
162 residue_function.get_residue(
163 value, time, 0, esf, residue,
164 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
165 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
166 }
167 Ke = (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)K * eigenvector;
168 const double eigenvector_eigenvector = eigenvector * eigenvector;
169 r[2 * s] = eigenvector_eigenvector * (use_constraint_with_time ? time : 1) - 1;
170 }
171
172 void operator()(const x_type& x, r_type& r, const x_type& dx, r_type& dr) {
173 EquilibriumSolution esf(false);
174 b2linalg::Vector<double, b2linalg::Vdense_constref> value = x(b2linalg::Interval(0, s));
175 b2linalg::Vector<double, b2linalg::Vdense_constref> eigenvector =
176 x(b2linalg::Interval(s, 2 * s));
177 double time = x[2 * s];
178
179 r.resize(x.size());
180 dr.resize(x.size());
181 b2linalg::Vector<double, b2linalg::Vdense_ref> residue = r(b2linalg::Interval(0, s));
182 b2linalg::Vector<double, b2linalg::Vdense_ref> Ke = r(b2linalg::Interval(s, 2 * s));
183
184 residue.set_zero();
185 K.set_zero_same_struct();
186 d_residue_d_time.set_zero();
187 residue_function.get_residue(
188 value, time, 0, esf, residue, K, d_residue_d_time, 0, 0, 0);
189 Ke = (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)K * eigenvector;
190 const double eigenvector_eigenvector = eigenvector * eigenvector;
191 r[2 * s] = eigenvector_eigenvector * (use_constraint_with_time ? time : 1) - 1;
192
193 dr(b2linalg::Interval(0, s)) = K * dx(b2linalg::Interval(0, s));
194 dr(b2linalg::Interval(0, s)) += d_residue_d_time * dx[2 * s];
195 dr(b2linalg::Interval(s, 2 * s)) = K * dx(b2linalg::Interval(s, 2 * s));
196
197 {
198 const double hh = h / eigenvector.norm2();
199 Z.set_zero_same_struct();
200 Z += (diff_coeff[0] / hh) * K;
201 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
202 tmp = value;
203 tmp += (i * hh) * eigenvector;
204 TMP.set_zero_same_struct();
205 residue_function.get_residue(
206 tmp, time, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
207 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
208 Z += (diff_coeff[i] / hh) * TMP;
209 }
210 }
211 dr(b2linalg::Interval(s, 2 * s)) += Z * dx(b2linalg::Interval(0, s));
212
213 {
214 tmp = (diff_coeff[0] / h) * Ke;
215 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
216 TMP.set_zero_same_struct();
217 residue_function.get_residue(
218 value, time + i * h, 0, esf,
219 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
220 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
221 tmp += (diff_coeff[i] / h)
222 * (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)TMP
223 * eigenvector;
224 }
225 }
226 dr(b2linalg::Interval(s, 2 * s)) += tmp * dx[2 * s];
227
228 if (use_constraint_with_time) {
229 dr[2 * s] = 2 * x[2 * s] * (eigenvector * dx(b2linalg::Interval(s, 2 * s)))
230 + (eigenvector * eigenvector) * dx[2 * s];
231 } else {
232 dr[2 * s] = 2 * (eigenvector * dx(b2linalg::Interval(s, 2 * s)));
233 }
234
235 dr *= -1;
236 }
237
238 bool nd(const x_type& x, const r_type& r, x_type& d) {
239 EquilibriumSolution esf(false);
240 b2linalg::Vector<double, b2linalg::Vdense_constref> value = x(b2linalg::Interval(0, s));
241 b2linalg::Vector<double, b2linalg::Vdense_constref> eigenvector =
242 x(b2linalg::Interval(s, 2 * s));
243 double time = x[2 * s];
244 b2linalg::Vector<double, b2linalg::Vdense_constref> residue =
245 r(b2linalg::Interval(0, s));
246 b2linalg::Vector<double, b2linalg::Vdense_constref> Ke =
247 r(b2linalg::Interval(s, 2 * s));
248 d.resize(x.size());
249
250 qfp[0] = b2linalg::inverse(K) * residue;
251 qfp[1] = b2linalg::inverse(K) * d_residue_d_time;
252 qfp *= -1;
253
254 {
255 const double hh = h / eigenvector.norm2();
256 Z.set_zero_same_struct();
257 Z += (diff_coeff[0] / hh) * K;
258 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
259 tmp = value;
260 tmp += (i * hh) * eigenvector;
261 TMP.set_zero_same_struct();
262 residue_function.get_residue(
263 tmp, time, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
264 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
265 Z += (diff_coeff[i] / hh) * TMP;
266 }
267 }
268
269 phifp = (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)Z * qfp;
270 phifp[0] += Ke;
271
272 {
273 tmp = (diff_coeff[0] / h) * Ke;
274 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
275 TMP.set_zero_same_struct();
276 residue_function.get_residue(
277 value, time + i * h, 0, esf,
278 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
279 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
280 tmp += (diff_coeff[i] / h)
281 * (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)TMP
282 * eigenvector;
283 }
284 }
285
286 phifp[1] += tmp;
287
288 phifp *= -1;
289 phifp[0] = inverse(K) * phifp[0];
290 phifp[1] = inverse(K) * phifp[1];
291
292 // computation of delta lambda with constraint := eigenvector * eigenvector - 1
293 double dc;
294 const double eigenvector_eigenvector = eigenvector * eigenvector;
295
296 if (use_constraint_with_time) {
297 dc = -(r[2 * s] + 2 * time * (eigenvector * phifp[0]))
298 / (eigenvector_eigenvector + 2 * time * (eigenvector * phifp[1]));
299 } else {
300 dc = -(r[2 * s] + 2 * (eigenvector * phifp[0])) / (2 * (eigenvector * phifp[1]));
301 }
302
303 qfp[0] += dc * qfp[1];
304 phifp[0] += dc * phifp[1];
305
306 d(b2linalg::Interval(0, s)) = qfp[0];
307 d(b2linalg::Interval(s, 2 * s)) = phifp[0];
308 d[2 * s] = dc;
309 return true;
310 }
311
312 void gd(const x_type& x, const r_type& r, x_type& d) {
313 EquilibriumSolution esf(false);
314 b2linalg::Vector<double, b2linalg::Vdense_constref> value = x(b2linalg::Interval(0, s));
315 b2linalg::Vector<double, b2linalg::Vdense_constref> eigenvector =
316 x(b2linalg::Interval(s, 2 * s));
317 double time = x[2 * s];
318 b2linalg::Vector<double, b2linalg::Vdense_constref> Ke =
319 r(b2linalg::Interval(s, 2 * s));
320 d.resize(x.size());
321
322 K.set_zero_same_struct();
323 d_residue_d_time.set_zero();
324 residue_function.get_residue(
325 value, time, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, K,
326 d_residue_d_time, 0, 0, 0);
327
328 d(b2linalg::Interval(0, s)) = K * r(b2linalg::Interval(0, s));
329 d(b2linalg::Interval(0, s)) += d_residue_d_time * r[2 * s];
330 d(b2linalg::Interval(s, 2 * s)) = K * r(b2linalg::Interval(s, 2 * s));
331
332 {
333 const double hh = h / eigenvector.norm2();
334 Z.set_zero_same_struct();
335 Z += (diff_coeff[0] / hh) * K;
336 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
337 tmp = value;
338 tmp += (i * hh) * eigenvector;
339 TMP.set_zero_same_struct();
340 residue_function.get_residue(
341 tmp, time, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
342 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
343 Z += (diff_coeff[i] / hh) * TMP;
344 }
345 }
346 d(b2linalg::Interval(s, 2 * s)) += Z * r(b2linalg::Interval(0, s));
347
348 {
349 tmp = (diff_coeff[0] / h) * Ke;
350 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
351 TMP.set_zero_same_struct();
352 residue_function.get_residue(
353 value, time + i * h, 0, esf,
354 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
355 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
356 tmp += (diff_coeff[i] / h)
357 * (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)TMP
358 * eigenvector;
359 }
360 }
361 d(b2linalg::Interval(s, 2 * s)) += tmp * r[2 * s];
362
363 if (use_constraint_with_time) {
364 d[2 * s] = 2 * x[2 * s] * (eigenvector * r(b2linalg::Interval(s, 2 * s)))
365 + (eigenvector * eigenvector) * r[2 * s];
366 } else {
367 d[2 * s] = 2 * (eigenvector * r(b2linalg::Interval(s, 2 * s)));
368 }
369
370 d *= -1;
371 }
372
373 private:
374 DefaultStaticResidueFunction<T, MATRIX_FORMAT>& residue_function;
375 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K;
376 size_t s;
377
378 b2linalg::Vector<T, b2linalg::Vdense> residue;
379 b2linalg::Vector<T, b2linalg::Vdense> d_residue_d_time;
380 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> Z;
381 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> TMP;
382 b2linalg::Vector<T, b2linalg::Vdense> tmp;
383 b2linalg::Matrix<T, b2linalg::Mrectangle> qfp;
384 b2linalg::Matrix<T, b2linalg::Mrectangle> phifp;
385
386 double approximation_order;
387 double h;
388 std::vector<double> diff_coeff;
389 double eps;
390 bool use_constraint_with_time;
391 };
392
393 SolverUtils<T, MATRIX_FORMAT> solver_utile;
394};
395
396template <typename T, typename MATRIX_FORMAT>
397typename TaylorExpansionsBucklingSolver<T, MATRIX_FORMAT>::type_t
398 TaylorExpansionsBucklingSolver<T, MATRIX_FORMAT>::type(
399 "TaylorExpansionsBucklingSolver", type_name<T, MATRIX_FORMAT>(),
400 StringList("TAYLOR_EXPANSIONS_BUCKLING", "P_BIFURCATION"), module, &Solver::type);
401
402} // namespace b2000::solver
403
409template <typename T, typename MATRIX_FORMAT>
410void b2000::solver::TaylorExpansionsBucklingSolver<T, MATRIX_FORMAT>::solve_nonlinear_refinement(
411 Case& case_, DefaultStaticResidueFunction<T, MATRIX_FORMAT>& residue_function,
412 b2linalg::Vector<double, b2linalg::Vdense_ref> value, double& time,
413 b2linalg::Vector<double, b2linalg::Vdense_ref> eigenvector,
414 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K) {
415 const double approximation_order =
416 case_.get_int("APPROXIMATION_ORDER_REFINEMENT", case_.get_int("APPROXIMATION_ORDER", 1));
417 const double h = std::pow(1e-15, 1.0 / (approximation_order + 1))
418 * case_.get_double("H_SCALE_REFINEMENT", case_.get_double("H_SCALE", 1.0));
419 std::vector<double> diff_coeff;
420 get_forward_coefficient(1, 0, approximation_order, diff_coeff);
421 const double eps = 1e-5;
422
423 const bool use_constraint_with_time = false;
424
425 const size_t s = value.size();
426 b2linalg::Vector<T, b2linalg::Vdense> residue(s);
427 b2linalg::Vector<T, b2linalg::Vdense> d_residue_d_time(s);
428 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> Z;
429 Z.set_same_structure(K);
430 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> TMP;
431 TMP.set_same_structure(K);
432 b2linalg::Vector<T, b2linalg::Vdense> tmp(s);
433 b2linalg::Vector<T, b2linalg::Vdense> Ke(s);
434 b2linalg::Matrix<T, b2linalg::Mrectangle> qfp(s, 2);
435 b2linalg::Matrix<T, b2linalg::Mrectangle> phifp(s, 2);
436 eigenvector *= 1.0 / (eigenvector.norm2() * (use_constraint_with_time ? std::sqrt(time) : 1));
437
438 for (;;) {
439 EquilibriumSolution esf(false);
440 residue.set_zero();
441 K.set_zero_same_struct();
442 d_residue_d_time.set_zero();
443 residue_function.get_residue(value, time, 0, esf, residue, K, d_residue_d_time, 0, 0, 0);
444 Ke = (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)K * eigenvector;
445 const double eigenvector_eigenvector = eigenvector * eigenvector;
446 double c = eigenvector_eigenvector * (use_constraint_with_time ? time : 1) - 1;
447 // termination test
448 const double res1 = residue.norm2();
449 const double res2 = Ke.norm2();
450 const double res3 = c;
451 std::cout << res1 << " " << res2 << " " << res3 << " " << time << std::endl;
452 if (res1 < eps && res2 < eps && res3 < eps) { break; }
453
454 qfp[0] = inverse(K) * residue;
455 qfp[1] = inverse(K) * d_residue_d_time;
456 qfp *= -1;
457
458#ifndef NEW_METHOD0
459 {
460 double hh = h / eigenvector.norm2();
461 tmp = value;
462 tmp += hh * eigenvector;
463 Z.set_zero_same_struct();
464 residue_function.get_residue(
465 tmp, time, 0, esf, b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref>::null,
466 Z, b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref>::null, 0, 0, 0);
467 Z -= K;
468 Z *= 1 / hh;
469 }
470#endif
471
472#ifdef NEW_METHOD0
473 {
474 const double hh = h / eigenvector.norm2();
475 Z.set_zero_same_struct();
476 Z += (diff_coeff[0] / hh) * K;
477 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
478 tmp = value;
479 tmp += (i * hh) * eigenvector;
480 TMP.set_zero_same_struct();
481 residue_function.get_residue(
482 tmp, time, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, TMP,
483 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
484 Z += (diff_coeff[i] / hh) * TMP;
485 }
486 }
487#endif
488
489 phifp = (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)Z * qfp;
490 phifp[0] += Ke;
491
492#ifndef NEW_METHOD1
493 {
494 Z.set_zero_same_struct();
495 residue_function.get_residue(
496 value, time + h, 0, esf,
497 b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref>::null, Z,
498 b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref>::null, 0, 0, 0);
499 tmp = (b2000::b2linalg::Matrix<T, b2000::b2linalg::Msym_compressed_col_st_constref>)Z
500 * eigenvector;
501 tmp -= Ke;
502 tmp *= 1 / h;
503 }
504#endif
505
506#ifdef NEW_METHOD1
507 {
508 tmp = (diff_coeff[0] / h) * Ke;
509 for (unsigned int i = 1; i < diff_coeff.size(); ++i) {
510 TMP.set_zero_same_struct();
511 residue_function.get_residue(
512 value, time + i * h, 0, esf, b2linalg::Vector<T, b2linalg::Vdense_ref>::null,
513 TMP, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, 0);
514 tmp += (diff_coeff[i] / h)
515 * (b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>)TMP
516 * eigenvector;
517 }
518 }
519#endif
520
521 phifp[1] += tmp;
522
523 phifp *= -1;
524 phifp[0] = inverse(K) * phifp[0];
525 phifp[1] = inverse(K) * phifp[1];
526
527 // computation of delta lambda with constraint := eigenvector * eigenvector - 1
528 double dc;
529 if (use_constraint_with_time) {
530 dc = -(c + 2 * time * (eigenvector * phifp[0]))
531 / (eigenvector_eigenvector + 2 * time * (eigenvector * phifp[1]));
532 } else {
533 dc = -(c + 2 * (eigenvector * phifp[0])) / (2 * (eigenvector * phifp[1]));
534 }
535
536 const double s = case_.get_double("NRS", 1);
537 qfp *= s;
538 phifp *= s;
539 time += s * dc;
540 qfp[0] += dc * qfp[1];
541 value += qfp[0];
542 phifp[0] += dc * phifp[1];
543 eigenvector += phifp[0];
544
545 // if (std::abs(dc) < eps2)
546 // break;
547 }
548}
549
554template <typename T, typename MATRIX_FORMAT>
555void b2000::solver::TaylorExpansionsBucklingSolver<T, MATRIX_FORMAT>::solve_on_case(Case& case_) {
556 logging::Logger& logger = logging::get_logger("solver.taylor_expansions_buckling");
557
558 if (model_ == nullptr) { Exception() << THROW; }
559
560 // Compute the static linear solution
561 logger << logging::info << "Start the Taylor expansions buckling solver for the case "
562 << case_.get_id() << "." << logging::LOGFLUSH;
563 model_->set_case(case_);
564
565 logger << logging::info << "Assemble the linear problem." << logging::LOGFLUSH;
566
567 SolutionTaylorExpansionsBuckling<T>& solution =
568 dynamic_cast<SolutionTaylorExpansionsBuckling<T>&>(model_->get_solution());
569
570 int nb_eigeinvalue;
571 double sigma;
572 char which[3];
573 solution.which_modes(nb_eigeinvalue, which, sigma);
574 const int path_approximation_order =
575 case_.get_int("APPROXIMATION_ORDER_PATH", case_.get_int("APPROXIMATION_ORDER", 2));
576 const int path_taylor_order =
577 case_.get_int("TAYLOR_ORDER_PATH", case_.get_int("TAYLOR_ORDER", 1));
578 const double path_h_scale = case_.get_double("H_SCALE_PATH", case_.get_double("H_SCALE", 1.0));
579 const int K_approximation_order =
580 case_.get_int("APPROXIMATION_ORDER_K", case_.get_int("APPROXIMATION_ORDER", 2));
581 const int K_taylor_order = case_.get_int("TAYLOR_ORDER_K", case_.get_int("TAYLOR_ORDER", 1));
582 const double K_h_scale = case_.get_double("H_SCALE_K", case_.get_double("H_SCALE", 1.0));
583 const double tol_eigenvalue_imag = case_.get_double("TOL_EIGENVALUE_IMAG", 1e-6);
584
585 DefaultStaticResidueFunction<T, MATRIX_FORMAT> residue_function;
586 residue_function.init(*model_, case_);
587
588 const size_t nb_dof = model_->get_domain().get_number_of_dof();
589 const size_t nb_rdof = residue_function.get_number_of_dof();
590 const size_t nb_lndof = residue_function.get_number_of_non_lagrange_dof();
591
592 b2linalg::Vector<T, b2linalg::Vdense> rdof(nb_rdof);
593 double time = 0;
594
595 logger << logging::info << "Compute the Taylor expansion of the path solution."
596 << logging::LOGFLUSH;
597
598 EquilibriumSolution es(true);
599
600 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible_ext> K(
601 nb_rdof, 1, model_->get_domain().get_dof_connectivity_type(), case_);
602 b2linalg::Matrix<T, b2linalg::Mrectangle> d_dof_and_time_d_path(nb_rdof + 1, path_taylor_order);
603 {
604 b2linalg::Vector<T, b2linalg::Vdense> q(nb_rdof);
605 b2linalg::Vector<T, b2linalg::Vdense> arc(nb_rdof + 1);
606 arc[nb_rdof] = 1.0;
607
608 residue_function.get_d_dof_d_path(
609 rdof, time, K, q, arc, d_dof_and_time_d_path, path_approximation_order, path_h_scale);
610 }
611
612 b2linalg::Matrix<T, b2linalg::Mcompressed_col> reduc_d_value_d_dof_trans;
613 double K_h_scale1;
614 {
615 b2linalg::Matrix<T, b2linalg::Mrectangle> path_te(nb_dof, d_dof_and_time_d_path.size2());
616 b2linalg::Vector<T, b2linalg::Vdense> reduc_r(nb_dof);
617 residue_function.mrbc_get_nonlinear_value(
618 rdof, time, false, b2linalg::Vector<T, b2linalg::Vdense_ref>::null,
619 reduc_d_value_d_dof_trans, reduc_r, 0);
620 path_te = transposed(reduc_d_value_d_dof_trans)
621 * d_dof_and_time_d_path(
622 b2linalg::Interval(0, nb_lndof),
623 b2linalg::Interval(0, d_dof_and_time_d_path.size2()));
624 b2linalg::Vector<T, b2linalg::Vdense> time_te(path_te.size2());
625 int s = 1;
626 for (size_t j = 0; j != path_te.size2(); ++j) {
627 s *= j + 1;
628 time_te[j] = d_dof_and_time_d_path(nb_rdof, j) / s;
629 path_te[j] *= 1.0 / s;
630 path_te[j] += time_te[j] * reduc_r;
631 }
632 solution.set_path(path_te, time_te);
633 K_h_scale1 = 1.0 / path_te[0].norm2();
634 }
635
636 logger << logging::info << "Compute the Taylor expansion of the stiffened matrix."
637 << logging::LOGFLUSH;
638
639 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> >
640 d2_residue_d_dof_d_path(K_taylor_order);
641 residue_function.get_d2_residue_d_dof_d_path(
642 rdof, time, K, d_dof_and_time_d_path, d2_residue_d_dof_d_path, K_approximation_order,
643 K_h_scale * K_h_scale1);
644
645 {
646 int s = 1;
647 for (size_t i = 0; i != d2_residue_d_dof_d_path.size(); ++i) {
648 s *= i + 1;
649 d2_residue_d_dof_d_path[i] *= -1.0 / s;
650 }
651 }
652
653 logger << logging::info << "Polynomial eigenvalue problem resolution" << logging::LOGFLUSH;
654
655 b2linalg::Vector<T, b2linalg::Vdense> eigvaluer(nb_eigeinvalue);
656 b2linalg::Vector<T, b2linalg::Vdense> eigvaluei(nb_eigeinvalue);
657 b2linalg::Matrix<T, b2linalg::Mrectangle> eigvector_red(nb_rdof, nb_eigeinvalue);
658
659 eigvector_red = eigenvector(
660 (const b2linalg::Matrix<T, b2linalg::Msym_compressed_col_update_inv>&)(K), 'P',
661 d2_residue_d_dof_d_path, 'N', eigvaluer, eigvaluei, which, sigma);
662
663 logger << logging::info << "Compute the buckling displacement and modes" << logging::LOGFLUSH;
664
665 const int nonlinear_refinement = case_.get_int("NONLINEAR_REFINEMENT", 0);
666
667 for (int j = 0; j != nb_eigeinvalue; ++j) {
668 if (std::abs(eigvaluei[j] / eigvaluer[j]) > tol_eigenvalue_imag) { continue; }
669 b2linalg::Vector<double> dr(d_dof_and_time_d_path.size1());
670 dr(b2linalg::Interval(0, nb_rdof)) = rdof;
671 dr[nb_rdof] = time;
672 {
673 b2linalg::Vector<double> p(d_dof_and_time_d_path.size2());
674 double s = 1;
675 for (size_t i = 0; i != p.size(); ++i) {
676 s *= eigvaluer[j] / (i + 1);
677 p[i] = s;
678 }
679 dr += d_dof_and_time_d_path * p;
680 }
681
682 if (j < nonlinear_refinement) {
683 logger << logging::info
684 << "Compute the nonlinear refinement on the buckling value and modes " << j
685 << "." << logging::LOGFLUSH;
686#ifdef NEW_METHOD2
687 solve_nonlinear_refinement_sha2(
688 case_, residue_function, dr(Interval(0, nb_rdof)), dr[nb_rdof], eigvector_red[j],
689 K);
690#else
691 solve_nonlinear_refinement(
692 case_, residue_function, dr(b2linalg::Interval(0, nb_rdof)), dr[nb_rdof],
693 eigvector_red[j], K);
694#endif
695 }
696
697 b2linalg::Vector<double> d(nb_dof);
698 residue_function.get_solution(dr(b2linalg::Interval(0, nb_rdof)), dr[nb_rdof], false, d);
699 b2linalg::Vector<double> ev(nb_dof);
700 ev = transposed(reduc_d_value_d_dof_trans) * eigvector_red[j];
701 ev *= 1.0 / ev.norm2();
702
703 // TODO: compute nbc, residue and gradient at dof
704 solution.set_solution(
705 j, eigvaluer[j], eigvaluei[j], d, dr[nb_rdof],
706 b2linalg::Vector<T, b2linalg::Vdense_constref>::null,
707 b2linalg::Vector<T, b2linalg::Vdense_constref>::null, ev);
708 }
709
710 RTable attributes;
711 solution.terminate_case(true, attributes);
712 case_.warn_on_non_used_key();
713 logger.LOG(logging::info, "End of Taylor expansions buckling solver");
714}
715
716#endif // B2LINEARIZE_PREBUCKLING_SOLVER_H_
THROW
#define THROW
Definition b2exception.H:198
b2solver.H
Interface to C++ representations of FE solvers.
b2000::Case::get_number_of_stage
virtual size_t get_number_of_stage() const =0
b2000::Dictionary::get_double
virtual double get_double(const std::string &key) const =0
b2000::Dictionary::get_int
virtual int get_int(const std::string &key) const =0
b2000::Model::get_case_list
virtual CaseList & get_case_list()=0
b2000::Solver::case_
Case * case_
This also.
Definition b2solver.H:93
b2000::Solver::model_
Model * model_
Definition b2solver.H:91
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829
b2000::get_forward_coefficient
void get_forward_coefficient(int l, int z, int p, std::vector< double > &a)
Definition b2numerical_differentiation.C:1982