b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2solution_database.H
1//------------------------------------------------------------------------
2// b2solution_database.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) 2004-2012,2016,2017 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 B2SOLUTION_DATABASE_H_
22#define B2SOLUTION_DATABASE_H_
23
24#include <complex>
25#include <limits>
26#include <map>
27#include <string>
28#include <type_traits>
29
30#include "b2case_database.H"
31#include "b2domain_database.H"
32#include "b2ppconfig.h"
33#include "model/b2boundary_condition.H"
34#include "model/b2solution.H"
35#include "utils/b2constants.H"
36#include "utils/b2database.H"
37#include "utils/b2object.H"
38#include "utils/b2threading.H"
39
40#ifdef HAVE_LIB_TBB
41#include "tbb/enumerable_thread_specific.h"
42#endif
43
44namespace b2000::b2dbv3 {
45
46struct sol_type_properties {
47 bool is_complex = false;
48 bool is_csda = false;
49 sol_type_properties(bool is_complex, bool is_csda) : is_complex(is_complex), is_csda(is_csda) {}
50};
51
52std::string sol_type_name(std::string basename, sol_type_properties);
53
54class GradientContainerOffload {
55public:
56 virtual void offload_data() = 0;
57};
58
59class GradientContainer;
60
61class GradientSolution : virtual public b2000::Solution {
62public:
63 void set_model(b2000::Model& model_) override {
64 model = &model_;
65 database = dynamic_cast<DataBaseFieldSet*>(&model_);
66 if (database == nullptr) { TypeError() << THROW; }
67 domain = dynamic_cast<Domain*>(&model_.get_domain());
68 if (domain == nullptr) { TypeError() << THROW; }
69 list_branch.clear();
70 for (size_t i = 0; i != domain->list_branch.size(); ++i) {
71 list_branch.push_back(Branch(this));
72 }
73 }
74
75 void set_case(b2000::Case& case__) override {
76 case_ = &dynamic_cast<b2000::b2dbv3::Case&>(case__);
77 need_refinement = false;
78 // is the "container" name that holds all SP fields: The name is arbitrary,
79 // and SP gradient fields are saved under their respective names.
80 db_grad_name = "GRADIENTS";
81 if (db_grad_name.case_undef()) { db_grad_name.case_(case__.get_id()); }
82 for (size_t i = 0; i != list_branch.size(); ++i) {
83 list_branch[i].field_list_i.clear();
84 list_branch[i].field_list_d.clear();
85 list_branch[i].field_list_cs.clear();
86 list_branch[i].field_list_c.clear();
87 }
88 if (case__.has_key("COMPUTE_FIELDS")) {
89 std::string s = case__.get_string("COMPUTE_FIELDS");
90 for (size_t i = 0; i < s.size();) {
91 const size_t j = s.find('.', i);
92 list_fields_to_compute.insert(s.substr(i, j - i));
93 if (j == std::string::npos) { break; }
94 i = j + 1;
95 }
96 list_fields_to_compute_neg = false;
97 }
98 }
99
100 void set_subcase_id(int subcase_id_) override {
101 if (subcase_id == subcase_id_) { return; }
102 terminate_gradient();
103 subcase_id = subcase_id_;
104 }
105
106 void set_system_null_space(
107 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& nks, bool normalize = true,
108 const std::string suffix = "NS") override {
109 SetName sol = case_->get_string("DOF_SOL", std::string("DISP")) + "_" + suffix;
110 if (sol.case_undef()) { sol.case_(case_->get_id()); }
111 database->add_field_entry(
112 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
113 database->add_field_entry(
114 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
115 b2linalg::Vector<double> tmp;
116 for (size_t j = 0; j != nks.size2(); ++j) {
117 sol.subcase(j + 1);
118 tmp = nks[j];
119 if (normalize) { tmp.normalize(); }
120 domain->save_global_dof(sol, b2linalg::Vector<double, b2linalg::Vdense_constref>(tmp));
121 }
122 database->sync();
123 }
124
125 void set_system_null_space(
126 const b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>& nks,
127 bool normalize = true, const std::string suffix = "NS") override {
128 SetName sol = case_->get_string("DOF_SOL", std::string("DISP")) + "_" + suffix;
129 if (sol.case_undef()) { sol.case_(case_->get_id()); }
130 database->add_field_entry(
131 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
132 database->add_field_entry(
133 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
134 b2linalg::Vector<b2000::csda<double>> tmp;
135 for (size_t j = 0; j != nks.size2(); ++j) {
136 sol.subcase(j + 1);
137 tmp = nks[j];
138 if (normalize) { tmp.normalize(); }
139 domain->save_global_dof(
140 sol, b2linalg::Vector<b2000::csda<double>, b2linalg::Vdense_constref>(tmp));
141 }
142 database->sync();
143 }
144
145 void set_system_null_space(
146 const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& nks,
147 bool normalize = true, const std::string suffix = "NS") override {
148 SetName sol = case_->get_string("DOF_SOL", std::string("DISP")) + "_" + suffix;
149 if (sol.case_undef()) { sol.case_(case_->get_id()); }
150 database->add_field_entry(
151 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
152 database->add_field_entry(
153 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
154 b2linalg::Vector<std::complex<double>> tmp;
155 for (size_t j = 0; j != nks.size2(); ++j) {
156 sol.subcase(j + 1);
157 tmp = nks[j];
158 if (normalize) { tmp.normalize(); }
159 domain->save_global_dof(
160 sol, b2linalg::Vector<std::complex<double>, b2linalg::Vdense_constref>(tmp));
161 }
162 database->sync();
163 }
164
165 bool is_need_refinement() const override { return need_refinement; }
166
167 void set_need_refinement() override { need_refinement = true; }
168
169 b2000::Solution::gradient_container get_gradient_container() override;
170
171 void add_sub_solution(Solution& s) { sub_solution.push_back(&s); }
172
173 void remove_sub_solution(Solution& s) {
174 sub_solution_t::iterator i = std::find(sub_solution.begin(), sub_solution.end(), &s);
175 if (i == sub_solution.end()) { Exception() << THROW; }
176 sub_solution.erase(i);
177 }
178
179 void set_value(const std::string& key, bool value) override { key_bool_value[key] = value; }
180 void set_value(const std::string& key, int value) override { key_int_value[key] = value; }
181 void set_value(const std::string& key, double value) override { key_double_value[key] = value; }
182 void set_value(const std::string& key, std::string& value) override {
183 key_string_value[key] = value;
184 }
185
186private:
187 using key_bool_value_t = std::map<std::string, bool>;
188 key_bool_value_t key_bool_value;
189 using key_int_value_t = std::map<std::string, int>;
190 key_int_value_t key_int_value;
191 using key_double_value_t = std::map<std::string, double>;
192 key_double_value_t key_double_value;
193 using key_string_value_t = std::map<std::string, std::string>;
194 key_string_value_t key_string_value;
195
196protected:
197 void store_value_in_rtable_and_clear(RTable& rt) {
198 for (key_bool_value_t::iterator i = key_bool_value.begin(); i != key_bool_value.end();
199 ++i) {
200 rt.set(i->first, i->second);
201 }
202 key_bool_value.clear();
203 for (key_int_value_t::iterator i = key_int_value.begin(); i != key_int_value.end(); ++i) {
204 rt.set(i->first, i->second);
205 }
206 key_int_value.clear();
207 for (key_double_value_t::iterator i = key_double_value.begin(); i != key_double_value.end();
208 ++i) {
209 rt.set(i->first, i->second);
210 }
211 key_double_value.clear();
212 for (key_string_value_t::iterator i = key_string_value.begin(); i != key_string_value.end();
213 ++i) {
214 rt.set(i->first, i->second);
215 }
216 key_string_value.clear();
217 }
218
219 DataBaseFieldSet* database{};
220 Domain* domain{};
221 b2000::Model* model{};
222 b2000::b2dbv3::Case* case_{};
223 int subcase_id{};
224 bool need_refinement{};
225
226 using sub_solution_t = std::vector<b2000::Solution*>;
227 sub_solution_t sub_solution{};
228
229 void set_cycle(int cycle_, bool end_of_stage_) {
230 cycle = cycle_;
231 end_of_stage = end_of_stage_;
232 db_grad_name.cycle(cycle);
233 }
234
235 void terminate_gradient() {
236 if (have_gradient) {
237 for (size_t i = 0; i != list_branch.size(); ++i) {
238 SetName grad_id = db_grad_name;
239 grad_id.branch(domain->list_branch[i].id);
240 if (subcase_id != 0) { grad_id.subcase(subcase_id); }
241 list_branch[i].save_fields(grad_id, *database);
242 }
243 }
244 }
245
246 // Collects all SP field names of SP fields created in any branch and adds
247 // the names to the list_sampling_fields set.
248 void add_gradient_solution_case_information(RTable& rtable) {
249 std::set<std::string> list_sampling_fields;
250 for (size_t i = 0; i != list_branch.size(); ++i) {
251 for (Branch::field_list_i_t::const_iterator j = list_branch[i].field_list_i.begin();
252 j != list_branch[i].field_list_i.end(); ++j) {
253 list_sampling_fields.insert(j->first.name);
254 }
255 for (Branch::field_list_d_t::const_iterator j = list_branch[i].field_list_d.begin();
256 j != list_branch[i].field_list_d.end(); ++j) {
257 list_sampling_fields.insert(j->first.name);
258 }
259 for (Branch::field_list_cs_t::const_iterator j = list_branch[i].field_list_cs.begin();
260 j != list_branch[i].field_list_cs.end(); ++j) {
261 list_sampling_fields.insert(j->first.name);
262 }
263 for (Branch::field_list_c_t::const_iterator j = list_branch[i].field_list_c.begin();
264 j != list_branch[i].field_list_c.end(); ++j) {
265 list_sampling_fields.insert(j->first.name);
266 }
267 }
268 if (!list_sampling_fields.empty()) {
269 std::string list_sampling_fields_string;
270 for (std::set<std::string>::const_iterator i = list_sampling_fields.begin();;) {
271 list_sampling_fields_string += *i;
272 if (++i == list_sampling_fields.end()) { break; }
273 list_sampling_fields_string += ',';
274 }
275 rtable.set("SP_SOL", list_sampling_fields_string);
276 }
277 }
278
279private:
280 bool list_fields_to_compute_neg{true};
281 std::set<std::string> list_fields_to_compute;
282
283 friend class GradientContainer;
284 friend class Branch;
285
286 SetName db_grad_name;
287 int cycle{-1};
288 bool end_of_stage{};
289 bool have_gradient{};
290 std::set<GradientContainerOffload*> gradient_container_set;
291
292 b2Mutex mutex;
293
294 struct Branch {
295 Branch(GradientSolution* parent_)
296 : parent(parent_), old_ip_list_size(0), old_elem_ip_list_size(0) {}
297
298 GradientSolution* parent;
299
300 struct CompInternalElementPosition
301 : public std::function<bool(
302 b2000::GradientContainer::InternalElementPosition,
303 b2000::GradientContainer::InternalElementPosition)> {
304 constexpr bool operator()(
305 const b2000::GradientContainer::InternalElementPosition& a,
306 const b2000::GradientContainer::InternalElementPosition& b) const {
307 constexpr double eps = 1e-7;
308 if (a.layer != b.layer) { return a.layer < b.layer; }
309 if (std::abs(a.t - b.t) > eps) { return a.t < b.t; }
310 if (std::abs(a.s - b.s) > eps) { return a.s < b.s; }
311 if (std::abs(a.r - b.r) > eps) { return a.r < b.r; }
312 return false;
313 // if (a.t != b.t)
314 // return a.t < b.t;
315 // if (a.s != b.s)
316 // return a.s < b.s;
317 // return a.r < b.r;
318 }
319 };
320
321 using ip_list_t = std::map<
322 b2000::GradientContainer::InternalElementPosition, mcInt32,
323 CompInternalElementPosition>;
324 ip_list_t ip_list;
325 size_t old_ip_list_size;
326 std::string old_ip_list_name;
327
328 struct ElemIPptr {
329 mcInt32 branch_elem_id;
330 ip_list_t::iterator ip;
331 bool operator<(const ElemIPptr& b) const {
332 if (branch_elem_id != b.branch_elem_id) {
333 return branch_elem_id < b.branch_elem_id;
334 }
335 return ip == b.ip ? false : CompInternalElementPosition()(ip->first, b.ip->first);
336 }
337 };
338
339 using elem_ip_list_t = std::map<ElemIPptr, mcInt32>;
340 elem_ip_list_t elem_ip_list;
341 size_t old_elem_ip_list_size;
342 std::string old_elem_ip_list_name;
343
344 struct ElemIP {
345 mcInt32 branch_elem_id;
346 mcInt32 ip;
347 ElemIP(const ElemIPptr& ei) : branch_elem_id(ei.branch_elem_id), ip(ei.ip->second) {}
348 };
349
350 SetName old_coor_ip_name;
351 SetName old_mbase_ip_name;
352
353 struct CompFieldDescription : public std::function<bool(
354 b2000::GradientContainer::FieldDescription,
355 b2000::GradientContainer::FieldDescription)> {
356 bool operator()(
357 const b2000::GradientContainer::FieldDescription& a,
358 const b2000::GradientContainer::FieldDescription& b) const {
359 return a.name < b.name;
360 }
361 };
362
363 template <typename T>
364 struct FieldValueReduction {
365 FieldValueReduction()
366 : max(0),
367 min(0),
368 max_set(false),
369 min_set(false),
370 max_pos(),
371 min_pos(),
372 integration(0) {}
373
374 void reset() {
375 max = 0;
376 min = 0;
377 max_set = false;
378 min_set = false;
379 integration = 0;
380 }
381
382 void add(const T v, const elem_ip_list_t::iterator pos, const double area) {
383 if (!max_set || v > max) {
384 max = v;
385 max_pos = pos;
386 max_set = true;
387 }
388 if (!min_set || v < min) {
389 min = v;
390 min_pos = pos;
391 min_set = true;
392 }
393 integration += area * v;
394 }
395
396 T max;
397 T min;
398 bool max_set;
399 bool min_set;
400 elem_ip_list_t::iterator max_pos;
401 elem_ip_list_t::iterator min_pos;
402 T integration;
403 };
404
405 template <typename T>
406 struct Field {
407 std::vector<elem_ip_list_t::iterator> elem_ip;
408 std::vector<T> values;
409 std::vector<FieldValueReduction<T>> vr_value;
410 };
411 using field_list_i_t = std::map<
412 b2000::GradientContainer::FieldDescription, Field<int>, CompFieldDescription>;
413 field_list_i_t field_list_i;
414 field_list_i_t::iterator field_list_i_cur;
415
416 using field_list_d_t = std::map<
417 b2000::GradientContainer::FieldDescription, Field<double>, CompFieldDescription>;
418 field_list_d_t field_list_d;
419 field_list_d_t::iterator field_list_d_cur;
420
421 using field_list_cs_t = std::map<
422 b2000::GradientContainer::FieldDescription, Field<b2000::csda<double>>,
423 CompFieldDescription>;
424 field_list_cs_t field_list_cs;
425 field_list_cs_t::iterator field_list_cs_cur;
426
427 using field_list_c_t = std::map<
428 b2000::GradientContainer::FieldDescription, Field<std::complex<double>>,
429 CompFieldDescription>;
430 field_list_c_t field_list_c;
431 field_list_c_t::iterator field_list_c_cur;
432
433 // for old gradient api (Fortran element)
434 b2linalg::Matrix<mcOff, b2linalg::Mrectangle> strs;
435
436 void save_fields(const SetName& grad_id, DataBaseFieldSet& database) {
437 if (strs.size2() != 0) {
438 database.set(grad_id, strs);
439 strs.resize(0, 0);
440 RTable desc;
441 desc.set("TYPE", "B2000-GRAD");
442 database.set_desc(grad_id, desc);
443 }
444
445 // Collect any remaining data from any gradient
446 // container object.
447 for (std::set<GradientContainerOffload*>::iterator i =
448 parent->gradient_container_set.begin();
449 i != parent->gradient_container_set.end(); ++i) {
450 (*i)->offload_data();
451 }
452
453 bool new_ip = false;
454 size_t ip_id = 0;
455 for (Branch::ip_list_t::iterator i = ip_list.begin(); i != ip_list.end();) {
456 if (i->second != 0) {
457 ip_list.erase(i++);
458 new_ip = true;
459 } else {
460 i->second = ++ip_id;
461 ++i;
462 }
463 }
464 if (ip_list.empty()) { return; }
465 if (new_ip || ip_list.size() != old_ip_list_size) {
466 std::vector<b2000::GradientContainer::InternalElementPosition> ip_list_v;
467 ip_list_v.reserve(ip_list.size());
468 for (Branch::ip_list_t::iterator i = ip_list.begin(); i != ip_list.end(); ++i) {
469 ip_list_v.push_back(i->first);
470 ++ip_list_v.back().layer;
471 }
472 SetName ip_name(grad_id);
473 ip_name.gname("IP");
474 ip_name.subcase(0);
475 std::vector<DataBase::ArrayTableCol> col_att;
476 col_att.push_back(DataBase::ArrayTableCol(
477 "rst", 3, ip_list_v.size(), &ip_list_v.front().r,
478 sizeof(b2000::GradientContainer::InternalElementPosition)));
479 col_att.push_back(DataBase::ArrayTableCol(
480 "layer", 1, ip_list_v.size(), &ip_list_v.front().layer,
481 sizeof(b2000::GradientContainer::InternalElementPosition)));
482 database.set(ip_name, col_att);
483 new_ip = true;
484 old_ip_list_name = ip_name;
485 old_ip_list_size = ip_list.size();
486 }
487
488 bool new_elem_ip = false;
489 size_t elem_ip_id = 0;
490 for (Branch::elem_ip_list_t::iterator i = elem_ip_list.begin();
491 i != elem_ip_list.end();) {
492 if (i->second != 0) {
493 elem_ip_list.erase(i++);
494 new_elem_ip = true;
495 } else {
496 i->second = ++elem_ip_id;
497 ++i;
498 }
499 }
500 if (new_ip || new_elem_ip || elem_ip_list.size() != old_elem_ip_list_size) {
501 std::vector<ElemIP> elem_ip_list_v;
502 elem_ip_list_v.reserve(elem_ip_list.size());
503 for (Branch::elem_ip_list_t::iterator i = elem_ip_list.begin();
504 i != elem_ip_list.end(); ++i) {
505 elem_ip_list_v.push_back(i->first);
506 }
507 SetName elem_ip_name(grad_id);
508 elem_ip_name.gname("ELEM_IP");
509 elem_ip_name.subcase(0);
510
511 std::vector<DataBase::ArrayTableCol> col_att;
512 col_att.push_back(DataBase::ArrayTableCol(
513 "elem", 1, elem_ip_list_v.size(), &elem_ip_list_v.front().branch_elem_id,
514 sizeof(ElemIP)));
515 col_att.push_back(DataBase::ArrayTableCol(
516 "ip", 1, elem_ip_list_v.size(), &elem_ip_list_v.front().ip, sizeof(ElemIP)));
517 database.set(elem_ip_name, col_att);
518 {
519 RTable desc;
520 std::ostringstream o;
521 SetName etab_name("ETAB");
522 etab_name.branch(grad_id.get_branch());
523 o << 1 << "->" << etab_name << "," << 2 << "->" << old_ip_list_name;
524 desc.set("JOIN", o.str());
525 database.set_desc(elem_ip_name, desc);
526 }
527 new_elem_ip = true;
528 old_elem_ip_list_name = elem_ip_name;
529 old_elem_ip_list_size = elem_ip_list.size();
530 }
531
532 save_fields_typed<int>(grad_id, database, new_elem_ip, field_list_i);
533 save_fields_typed<double>(grad_id, database, new_elem_ip, field_list_d);
534 save_fields_typed<b2000::csda<double>>(grad_id, database, new_elem_ip, field_list_cs);
535 save_fields_typed<std::complex<double>>(grad_id, database, new_elem_ip, field_list_c);
536 }
537
538 template <typename T>
539 void save_fields_typed(
540 const SetName& grad_id, DataBaseFieldSet& database, const bool new_elem_ip,
541 std::map<b2000::GradientContainer::FieldDescription, Field<T>, CompFieldDescription>&
542 field_list) {
543 // Collect any remaining data from any gradient
544 // container object.
545 for (std::set<GradientContainerOffload*>::iterator i =
546 parent->gradient_container_set.begin();
547 i != parent->gradient_container_set.end(); ++i) {
548 (*i)->offload_data();
549 }
550
551 for (typename std::map<
552 b2000::GradientContainer::FieldDescription, Field<T>,
553 CompFieldDescription>::iterator f = field_list.begin();
554 f != field_list.end(); ++f) {
555 if (f->second.elem_ip.empty()) { continue; }
556
557 SetName field_name(grad_id);
558 field_name.gname(f->first.name);
559
560 // remove the COOR field if the elem_ip is not new
561 if (f->first.name == "COOR_IP") {
562 if (!new_elem_ip && field_name.get_case() == old_coor_ip_name.get_case()
563 && field_name.get_subcase() == old_coor_ip_name.get_subcase()) {
564 f->second.elem_ip.clear();
565 f->second.values.clear();
566 continue;
567 } else {
568 old_coor_ip_name = field_name;
569 }
570 }
571 if (f->first.name == "MBASE_IP") {
572 if (!new_elem_ip && field_name.get_case() == old_mbase_ip_name.get_case()
573 && field_name.get_subcase() == old_mbase_ip_name.get_subcase()) {
574 f->second.elem_ip.clear();
575 f->second.values.clear();
576 continue;
577 } else {
578 old_mbase_ip_name = field_name;
579 }
580 }
581
582 database.add_field_entry(
583 f->first.description, field_name.get_gname(), "FIELD-SAMPLING", "BRANCH",
584 false, false, false);
585 std::vector<mcInt32> elem_ip_tmp;
586 std::vector<DataBase::ArrayTableCol> col_att;
587 bool field_pos_sorted = true;
588 if (f->second.elem_ip.size() > 1) {
589 const std::vector<elem_ip_list_t::iterator>::const_iterator i_end =
590 f->second.elem_ip.end() - 1;
591 for (std::vector<elem_ip_list_t::iterator>::const_iterator i =
592 f->second.elem_ip.begin();
593 i != i_end;) {
594 std::vector<elem_ip_list_t::iterator>::const_iterator ii = i++;
595 if (!((*i)->second > (*ii)->second)) {
596 field_pos_sorted = false;
597 break;
598 }
599 }
600 }
601 if (f->second.elem_ip.size() == old_elem_ip_list_size) {
602 if (!field_pos_sorted) {
603 const int field_nb_comp = f->first.nb_comp;
604 std::vector<std::pair<mcInt32, size_t>> elem_ip_tmp1(
605 f->second.elem_ip.size());
606 for (size_t i = 0; i != elem_ip_tmp1.size(); ++i) {
607 elem_ip_tmp1[i].first = f->second.elem_ip[i]->second;
608 elem_ip_tmp1[i].second = i * field_nb_comp;
609 }
610 std::sort(elem_ip_tmp1.begin(), elem_ip_tmp1.end(), CompareFirstOfPair());
611 std::vector<T> values_tmp;
612 values_tmp.reserve(elem_ip_tmp1.size());
613 const typename std::vector<T>::const_iterator fi = f->second.values.begin();
614 for (size_t i = 0; i != elem_ip_tmp1.size(); ++i) {
615 values_tmp.insert(
616 values_tmp.end(), fi + elem_ip_tmp1[i].second,
617 fi + elem_ip_tmp1[i].second + field_nb_comp);
618 }
619 f->second.values.swap(values_tmp);
620 }
621 } else {
622 elem_ip_tmp.resize(f->second.elem_ip.size());
623 if (field_pos_sorted) {
624 for (size_t i = 0; i != elem_ip_tmp.size(); ++i) {
625 elem_ip_tmp[i] = f->second.elem_ip[i]->second;
626 }
627 } else {
628 const int field_nb_comp = f->first.nb_comp;
629 std::vector<std::pair<mcInt32, size_t>> elem_ip_tmp1(
630 f->second.elem_ip.size());
631 for (size_t i = 0; i != elem_ip_tmp1.size(); ++i) {
632 elem_ip_tmp1[i].first = f->second.elem_ip[i]->second;
633 elem_ip_tmp1[i].second = i * field_nb_comp;
634 }
635 std::sort(elem_ip_tmp1.begin(), elem_ip_tmp1.end(), CompareFirstOfPair());
636 std::vector<T> values_tmp;
637 values_tmp.reserve(elem_ip_tmp1.size());
638 const typename std::vector<T>::const_iterator fi = f->second.values.begin();
639 for (size_t i = 0; i != elem_ip_tmp1.size(); ++i) {
640 values_tmp.insert(
641 values_tmp.end(), fi + elem_ip_tmp1[i].second,
642 fi + elem_ip_tmp1[i].second + field_nb_comp);
643 elem_ip_tmp[i] = elem_ip_tmp1[i].first;
644 }
645 f->second.values.swap(values_tmp);
646 }
647 col_att.push_back(DataBase::ArrayTableCol(
648 "pos", 1, elem_ip_tmp.size(), &elem_ip_tmp.front(), sizeof(mcInt32)));
649 }
650 f->second.elem_ip.clear();
651 col_att.push_back(DataBase::ArrayTableCol(
652 "value", f->first.nb_comp, f->second.values.size() / f->first.nb_comp,
653 &f->second.values.front(), sizeof(T) * f->first.nb_comp));
654 database.set(field_name, col_att);
655 f->second.values.clear();
656 {
657 RTable desc;
658 std::ostringstream o;
659 o << col_att.size() - 1 << "->" << old_elem_ip_list_name;
660 desc.set("JOIN", o.str());
661 desc.set("DOMAIN", "FIELD-SAMPLING");
662 desc.set("COMP_NAMES", f->first.components_name);
663 desc.set("DESCR", f->first.description);
664 desc.set("DOMAIN_NDIM", f->first.domain_ndim);
665 desc.set("NB_COMP", f->first.nb_comp);
666 if (f->first.tensor_order) {
667 desc.set("TENSOR_ORDER", f->first.tensor_order);
668 desc.set("TENSOR_SIZE", f->first.tensor_size);
669 desc.set("TENSOR_SYM", (f->first.tensor_symmetric ? 1 : 0));
670 }
671
672 // storage of the reduced value
673 if (f->first.tensor_order < 0) {
674 std::vector<T> max;
675 std::vector<T> min;
676 std::vector<mcInt32> max_elem_ip;
677 std::vector<mcInt32> min_elem_ip;
678 std::vector<T> integral;
679 assert((int)f->second.vr_value.size() == f->first.nb_comp);
680 for (int i = 0; i != f->first.nb_comp; ++i) {
681 assert(f->second.vr_value[i].max_pos != elem_ip_list.end());
682 assert(f->second.vr_value[i].min_pos != elem_ip_list.end());
683 max.push_back(f->second.vr_value[i].max);
684 min.push_back(f->second.vr_value[i].min);
685 max_elem_ip.push_back(f->second.vr_value[i].max_pos->second);
686 min_elem_ip.push_back(f->second.vr_value[i].min_pos->second);
687 integral.push_back(f->second.vr_value[i].integration);
688 f->second.vr_value[i].reset();
689 }
690 desc.set("MAX", max);
691 desc.set("MAX_ELEM_IP", max_elem_ip);
692 desc.set("MIN", min);
693 desc.set("MIN_ELEM_IP", min_elem_ip);
694 desc.set("INTEGRAL", integral);
695 }
696 database.set_desc(field_name, desc);
697 }
698 }
699 }
700 };
701
702 std::vector<Branch> list_branch;
703};
704
705class GradientContainer : public b2000::GradientContainer, public GradientContainerOffload {
706public:
707 GradientContainer(GradientSolution& solution_) : solution(solution_) {}
708
709 ~GradientContainer() override {
710 b2Mutex::scoped_lock lock(solution.mutex);
711
712 offload_data();
713 std::set<GradientContainerOffload*>::iterator i =
714 solution.gradient_container_set.find(this);
715 if (i != solution.gradient_container_set.end()) {
716 solution.gradient_container_set.erase(i);
717 }
718 }
719
720 bool set_current_element(const Element& e) override {
721 Context& c = get_context();
722 const std::pair<ssize_t, size_t> new_eid =
723 solution.domain->get_internal_branch_and_element_id(e);
724 if (new_eid != c.eid) {
725 b2Mutex::scoped_lock lock(solution.mutex);
726
727 offload_data_for_context(c);
728 c.eid = new_eid;
729 c.branch = &solution.list_branch[c.eid.first];
730 c.element_volume = 0;
731 }
732 assert(c.branch != nullptr);
733 return true;
734 }
735
736 void set_current_position(const InternalElementPosition& pos, const double volume) override {
737 Context& c = get_context();
738 assert(c.branch != nullptr);
739
740 c.point = c.point_data_map.insert(
741 c.point, PointDataMap::value_type(PointKey(c.eid, pos), PointData()));
742 (*c.point).second.volume = volume;
743 c.element_volume += volume;
744 }
745
746 void set_current_volume(const double volume) override {
747 Context& c = get_context();
748 assert(c.point != c.point_data_map.end());
749 (*c.point).second.volume = volume;
750 c.element_volume += volume;
751 }
752
753 double get_element_volume() override { return get_context().element_volume; }
754
755 void set_field_value(const FieldDescription& descr, const int* value) override {
756 Context& c = get_context();
757 assert(c.point != c.point_data_map.end());
758 FieldDescriptionSet::const_iterator i = c.descr_set.insert(descr).first;
759 PointData& p = (*c.point).second;
760 p.idata.descr.push_back(i);
761 p.idata.data.insert(p.idata.data.end(), value, value + descr.nb_comp);
762 }
763
764 void set_field_value(const FieldDescription& descr, const double* value) override {
765 Context& c = get_context();
766 assert(c.point != c.point_data_map.end());
767 FieldDescriptionSet::const_iterator i = c.descr_set.insert(descr).first;
768 PointData& p = (*c.point).second;
769 p.ddata.descr.push_back(i);
770 p.ddata.data.insert(p.ddata.data.end(), value, value + descr.nb_comp);
771 }
772
773 void set_field_value(const FieldDescription& descr, const b2000::csda<double>* value) override {
774 Context& c = get_context();
775 assert(c.point != c.point_data_map.end());
776 FieldDescriptionSet::const_iterator i = c.descr_set.insert(descr).first;
777 PointData& p = (*c.point).second;
778 p.csdata.descr.push_back(i);
779 p.csdata.data.insert(p.csdata.data.end(), value, value + descr.nb_comp);
780 }
781
782 void set_field_value(
783 const FieldDescription& descr, const std::complex<double>* value) override {
784 Context& c = get_context();
785 assert(c.point != c.point_data_map.end());
786 FieldDescriptionSet::const_iterator i = c.descr_set.insert(descr).first;
787 PointData& p = (*c.point).second;
788 p.cdata.descr.push_back(i);
789 p.cdata.data.insert(p.cdata.data.end(), value, value + descr.nb_comp);
790 }
791
792 void set_failure_index(const double failure_index) override {
793 Context& c = get_context();
794 assert(c.point != c.point_data_map.end());
795 (*c.point).second.failure_index = std::max((*c.point).second.failure_index, failure_index);
796 }
797
798 void set_current_discrete_position(const DiscreteInternalElementPosition& pos) override {
799 UnimplementedError() << THROW;
800 }
801
802 void set_discrete_value(const FieldDescription& descr, const int* value) override {
803 UnimplementedError() << THROW;
804 }
805
806 void set_discrete_value(const FieldDescription& descr, const double* value) override {
807 UnimplementedError() << THROW;
808 }
809
810 void set_discrete_value(
811 const FieldDescription& descr, const b2000::csda<double>* value) override {
812 UnimplementedError() << THROW;
813 }
814
815 void set_discrete_value(
816 const FieldDescription& descr, const std::complex<double>* value) override {
817 UnimplementedError() << THROW;
818 }
819
820 bool compute_field_value(const std::string& name) const override {
821 std::set<std::string>::const_iterator i = solution.list_fields_to_compute.find(name);
822 if (i == solution.list_fields_to_compute.end()) {
823 return solution.list_fields_to_compute_neg;
824 }
825 return true;
826 }
827
828private:
829 GradientSolution& solution;
830
831 // Identifies an integration point in a domain.
832 struct PointKey {
833 std::pair<ssize_t, size_t> eid;
834 InternalElementPosition ip;
835
836 PointKey(const std::pair<ssize_t, size_t> eid_, InternalElementPosition ip_)
837 : eid(eid_), ip(ip_) {}
838
839 bool operator<(const PointKey& o) const {
840 static constexpr double eps = 1e-7;
841 if (eid != o.eid) { return eid < o.eid; }
842 if (ip.layer != o.ip.layer) { return ip.layer < o.ip.layer; }
843 if (std::abs(ip.t - o.ip.t) > eps) { return ip.t < o.ip.t; }
844 if (std::abs(ip.s - o.ip.s) > eps) { return ip.s < o.ip.s; }
845 if (std::abs(ip.r - o.ip.r) > eps) { return ip.r < o.ip.r; }
846 return false;
847 }
848 };
849
850 struct CompFieldDescription : public std::function<bool(FieldDescription, FieldDescription)> {
851 bool operator()(const FieldDescription& a, const FieldDescription& b) const {
852 return a.name < b.name;
853 }
854 };
855
856 using FieldDescriptionSet = std::set<FieldDescription, CompFieldDescription>;
857
858 template <typename T>
859 struct Data {
860 std::vector<FieldDescriptionSet::const_iterator> descr;
861 std::vector<T> data;
862 };
863 using IData = Data<int>;
864 using DData = Data<double>;
865 using CSData = Data<b2000::csda<double>>;
866 using CData = Data<std::complex<double>>;
867
868 // Contains all data that is defined at an integration point.
869 struct PointData {
870 double volume;
871 double failure_index;
872 IData idata;
873 DData ddata;
874 CSData csdata;
875 CData cdata;
876
877 PointData() : volume(0), failure_index(0) {}
878 };
879 using PointDataMap = std::map<PointKey, PointData>;
880
881 struct Context {
882 GradientSolution::Branch* branch;
883 std::pair<ssize_t, size_t> eid;
884 FieldDescriptionSet descr_set;
885 PointDataMap point_data_map;
886 PointDataMap::iterator point;
887 double element_volume;
888
889 Context() : branch(nullptr), eid(-1, 0), point(point_data_map.begin()), element_volume(0) {}
890
891 void reset() {
892 branch = nullptr;
893 eid = std::pair<ssize_t, size_t>(-1, 0);
894 point_data_map.clear();
895 point = point_data_map.begin();
896 }
897 };
898
899#ifdef HAVE_LIB_TBB
900 tbb::enumerable_thread_specific<Context> contexts;
901#else
902 Context the_context;
903#endif
904
905 Context& get_context() {
906#ifdef HAVE_LIB_TBB
907 return contexts.local();
908#else
909 return the_context;
910#endif
911 }
912
913 void offload_data() override {
914#ifdef HAVE_LIB_TBB
915 for (tbb::enumerable_thread_specific<Context>::iterator i = contexts.begin();
916 i != contexts.end(); ++i) {
917 Context& c = *i;
918 offload_data_for_context(c);
919 }
920#else
921 offload_data_for_context(the_context);
922#endif
923 }
924
925 void offload_idata_for_context_and_point(
926 Context& c, const PointData& point_data,
927 GradientSolution::Branch::elem_ip_list_t::iterator& current_elem_ip) {
928 size_t o = 0;
929 for (size_t j = 0; j != point_data.idata.descr.size(); ++j) {
930 const FieldDescription& descr = *(point_data.idata.descr[j]);
931 const int* value = &point_data.idata.data[o];
932 o += descr.nb_comp;
933
934 GradientSolution::Branch::field_list_i_t::iterator k =
935 c.branch->field_list_i.find(descr);
936 if (k == c.branch->field_list_i.end()) {
937 k = c.branch->field_list_i
938 .insert(GradientSolution::Branch::field_list_i_t::value_type(
939 descr, GradientSolution::Branch::Field<int>()))
940 .first;
941 }
942 GradientSolution::Branch::Field<int>& f = (*k).second;
943 f.elem_ip.push_back(current_elem_ip);
944 f.values.insert(f.values.end(), value, value + descr.nb_comp);
945 }
946 }
947
948 void offload_ddata_for_context_and_point(
949 Context& c, const PointData& point_data,
950 GradientSolution::Branch::elem_ip_list_t::iterator& current_elem_ip) {
951 size_t o = 0;
952 for (size_t j = 0; j != point_data.ddata.descr.size(); ++j) {
953 const FieldDescription& descr = *(point_data.ddata.descr[j]);
954 const double* value = &point_data.ddata.data[o];
955 o += descr.nb_comp;
956
957 GradientSolution::Branch::field_list_d_t::iterator k =
958 c.branch->field_list_d.find(descr);
959 if (k == c.branch->field_list_d.end()) {
960 k = c.branch->field_list_d
961 .insert(GradientSolution::Branch::field_list_d_t::value_type(
962 descr, GradientSolution::Branch::Field<double>()))
963 .first;
964 }
965 GradientSolution::Branch::Field<double>& f = (*k).second;
966 f.elem_ip.push_back(current_elem_ip);
967 f.values.insert(f.values.end(), value, value + descr.nb_comp);
968 if (descr.tensor_order < 0) {
969 // we compute reduced value for non tensorial field
970 if (f.vr_value.empty()) { f.vr_value.resize(descr.nb_comp); }
971 for (int l = 0; l != descr.nb_comp; ++l) {
972 f.vr_value[l].add(value[l], current_elem_ip, point_data.volume);
973 }
974 }
975 }
976 }
977
978 void offload_csdata_for_context_and_point(
979 Context& c, const PointData& point_data,
980 GradientSolution::Branch::elem_ip_list_t::iterator& current_elem_ip) {
981 size_t o = 0;
982 for (size_t j = 0; j != point_data.csdata.descr.size(); ++j) {
983 const FieldDescription& descr = *(point_data.csdata.descr[j]);
984 const b2000::csda<double>* value = &point_data.csdata.data[o];
985 o += descr.nb_comp;
986
987 GradientSolution::Branch::field_list_cs_t::iterator k =
988 c.branch->field_list_cs.find(descr);
989 if (k == c.branch->field_list_cs.end()) {
990 k = c.branch->field_list_cs
991 .insert(GradientSolution::Branch::field_list_cs_t::value_type(
992 descr, GradientSolution::Branch::Field<b2000::csda<double>>()))
993 .first;
994 }
995 GradientSolution::Branch::Field<b2000::csda<double>>& f = (*k).second;
996 f.elem_ip.push_back(current_elem_ip);
997 f.values.insert(f.values.end(), value, value + descr.nb_comp);
998 if (descr.tensor_order < 0) {
999 // we compute reduced value for non tensorial field
1000 if (f.vr_value.empty()) { f.vr_value.resize(descr.nb_comp); }
1001 for (int l = 0; l != descr.nb_comp; ++l) {
1002 f.vr_value[l].add(value[l], current_elem_ip, point_data.volume);
1003 }
1004 }
1005 }
1006 }
1007
1008 void offload_cdata_for_context_and_point(
1009 Context& c, const PointData& point_data,
1010 GradientSolution::Branch::elem_ip_list_t::iterator& current_elem_ip) {
1011 size_t o = 0;
1012 for (size_t j = 0; j != point_data.cdata.descr.size(); ++j) {
1013 const FieldDescription& descr = *(point_data.cdata.descr[j]);
1014 const std::complex<double>* value = &point_data.cdata.data[o];
1015 o += descr.nb_comp;
1016
1017 GradientSolution::Branch::field_list_c_t::iterator k =
1018 c.branch->field_list_c.find(descr);
1019 if (k == c.branch->field_list_c.end()) {
1020 k = c.branch->field_list_c
1021 .insert(GradientSolution::Branch::field_list_c_t::value_type(
1022 descr, GradientSolution::Branch::Field<std::complex<double>>()))
1023 .first;
1024 }
1025 GradientSolution::Branch::Field<std::complex<double>>& f = (*k).second;
1026 f.elem_ip.push_back(current_elem_ip);
1027 f.values.insert(f.values.end(), value, value + descr.nb_comp);
1028 }
1029 }
1030
1031 void offload_data_for_context(Context& c) {
1032 if (c.point_data_map.empty()) { return; }
1033 assert(c.branch != nullptr);
1034
1035 // Get integration point with the maximum failure index.
1036 PointDataMap::const_iterator imax = c.point_data_map.begin();
1037 for (PointDataMap::const_iterator i = c.point_data_map.begin(); i != c.point_data_map.end();
1038 ++i) {
1039 if ((*i).second.failure_index > (*imax).second.failure_index) { imax = i; }
1040 }
1041
1042 GradientSolution::Branch::ip_list_t::iterator current_ip = c.branch->ip_list.end();
1043 GradientSolution::Branch::elem_ip_list_t::iterator current_elem_ip =
1044 c.branch->elem_ip_list.end();
1045
1046 // Loop over all integration points and add field data.
1047 for (PointDataMap::const_iterator i = c.point_data_map.begin(); i != c.point_data_map.end();
1048 ++i) {
1049 const PointKey& key = (*i).first;
1050 const PointData& data = (*i).second;
1051
1052 // Filter by failure index.
1053 if ((*imax).second.failure_index > data.failure_index) { continue; }
1054
1055 // Add integration point.
1056 {
1057 current_ip = c.branch->ip_list.insert(
1058 current_ip, GradientSolution::Branch::ip_list_t::value_type(key.ip, 0));
1059 current_ip->second = 0;
1060 const GradientSolution::Branch::ElemIPptr elem_ip_id = {
1061 mcInt32(c.eid.second) + 1, current_ip};
1062 current_elem_ip = c.branch->elem_ip_list.insert(
1063 current_elem_ip,
1064 GradientSolution::Branch::elem_ip_list_t::value_type(elem_ip_id, 0));
1065 current_elem_ip->second = 0;
1066 }
1067
1068 offload_idata_for_context_and_point(c, data, current_elem_ip);
1069 offload_ddata_for_context_and_point(c, data, current_elem_ip);
1070 offload_csdata_for_context_and_point(c, data, current_elem_ip);
1071 offload_cdata_for_context_and_point(c, data, current_elem_ip);
1072 }
1073
1074 // Reset.
1075 c.reset();
1076 }
1077};
1078
1079inline b2000::Solution::gradient_container GradientSolution::get_gradient_container() {
1080 int gi = case_->get_int("GRADIENTS", 0);
1081 have_gradient = false;
1082 if (cycle == -1) {
1083 if (gi == 0) { return b2000::Solution::gradient_container(nullptr); }
1084 } else if (!((gi > 0 && (cycle % gi == 0 || end_of_stage)) || (gi == -1 && end_of_stage))) {
1085 return b2000::Solution::gradient_container(nullptr);
1086 }
1087 have_gradient = true;
1088
1089 b2dbv3::GradientContainer* gc = new b2dbv3::GradientContainer(*this);
1090 gradient_container_set.insert(gc);
1091 return b2000::Solution::gradient_container(gc);
1092}
1093
1094template <typename T>
1095class SolutionStaticLinear : public b2000::b2dbv3::GradientSolution,
1096 public b2000::SolutionStaticLinear<T> {
1097public:
1098 bool compute_dof() const override { return true; }
1099
1100 void set_subcase_id(int subcase_id_) override {
1101 b2000::b2dbv3::GradientSolution::set_subcase_id(subcase_id_);
1102 }
1103
1104 void set_dof(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1105 SetName sol = case_->get_string("DOF_SOL", std::string("DISP"));
1106 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1107 if (subcase_id != 0) { sol.subcase(subcase_id); }
1108 database->add_field_entry(
1109 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
1110 database->add_field_entry(
1111 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
1112 domain->save_global_dof(sol, dof);
1113 }
1114
1115 void get_dof(b2linalg::Vector<T, b2linalg::Vdense_ref> dof) override {
1116 SetName sol = case_->get_string("DOF_SOL", std::string("DISP"));
1117 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1118 if (subcase_id != 0) { sol.subcase(subcase_id); }
1119 domain->load_global_dof(sol, dof);
1120 }
1121
1122 bool compute_natural_boundary_condition() const override { return true; }
1123
1124 void set_natural_boundary_condition(
1125 const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1126 SetName sol = case_->get_string("NBC_SOL", "FORC");
1127 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1128 if (subcase_id != 0) { sol.subcase(subcase_id); }
1129 database->add_field_entry(
1130 "NBC Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
1131 database->add_field_entry(
1132 "NBC Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1133 domain->save_global_dof(sol, dof);
1134 }
1135
1136 bool compute_constraint_value() const override { return case_->get_bool("COMPUTE_RCFO", true); }
1137
1138 void set_constraint_value(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1139 SetName sol = case_->get_string("RESIDUUM_SOL", "RCFO");
1140 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1141 if (subcase_id != 0) { sol.subcase(subcase_id); }
1142 database->add_field_entry(
1143 "Residuum Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
1144 database->add_field_entry(
1145 "Residuum Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1146 domain->save_global_dof(sol, dof);
1147 }
1148
1149 void terminate_case(bool success, const RTable& attributes) override {
1150 terminate_gradient();
1151 RTable& rtable = case_->get_solution_rtable();
1152 rtable.update(attributes);
1153 store_value_in_rtable_and_clear(rtable);
1154 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1155 rtable.set(
1156 "DOF_SOL", SetName(case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
1157 rtable.set(
1158 "NBC_SOL", SetName(case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
1159 if (compute_constraint_value()) {
1160 rtable.set(
1161 "RESIDUUM_SOL",
1162 SetName(case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
1163 }
1164 add_gradient_solution_case_information(rtable);
1165 SetName name("SOLUTION");
1166 name.case_(case_->get_id());
1167 database->set(name, rtable);
1168 database->sync();
1169 }
1170
1171 using type_t = ObjectTypeComplete<SolutionStaticLinear, Solution::type_t>;
1172 static type_t type;
1173};
1174
1175template <typename T>
1176typename SolutionStaticLinear<T>::type_t SolutionStaticLinear<T>::type(
1177 "SolutionStaticLinear", type_name<T>(),
1178 StringList("STATIC_LINEAR_SOLUTION", "LINEAR", "LINEAR_P_REFINEMENT"), module,
1179 &b2000::Solution::type);
1180
1181template <typename T>
1182class SolutionFreeVibration : public b2000::b2dbv3::GradientSolution,
1183 public b2000::SolutionFreeVibration<T> {
1184public:
1185 void set_case(b2000::Case& case__) override {
1186 GradientSolution::set_case(case__);
1187 nmodes = 0;
1188 }
1189
1190 void which_modes(
1191 int& number_of_mode, char which[3], T& sigma_min, T& sigma_max) const override {
1192 sigma_min = case_->get_csda_double("SHIFT", 0.0);
1193 sigma_min *= (M_PIl * 2);
1194 sigma_min *= sigma_min;
1195 sigma_max = sigma_min;
1196 number_of_mode = case_->get_int("NMODES", 1);
1197 std::string whichs = case_->get_string("WHICH_EIGENVALUES", "SM");
1198 strcpy(which, whichs.c_str());
1199 }
1200
1201 void set_mode(
1202 int mode, T eigenvalue,
1203 const b2linalg::Vector<T, b2linalg::Vdense_constref>& eigenvector) override {
1204 database->add_field_entry("Eigenmodes", "MODE", "NODE", "BRANCH", false, false, false);
1205 database->add_field_entry("Eigenmodes", "MODE_E", "CELL", "BRANCH", false, false, false);
1206 SetName name("MODE");
1207 if (name.case_undef()) { name.case_(case_->get_id()); }
1208 nmodes = ++mode;
1209 name.subcase(mode);
1210 RTable rtable;
1211 rtable.set("MODE", mode);
1212 rtable.set("EIGVAL", eigenvalue);
1213 bool pos_eigenvalue = true;
1214 if constexpr (std::is_same<T, double>::value) { pos_eigenvalue = eigenvalue >= 0; }
1215 if (pos_eigenvalue) {
1216 T w = sqrt(eigenvalue);
1217 rtable.set("OMEGA", w);
1218 T freq = w / static_cast<T>(2 * M_PIl);
1219 rtable.set("FREQ", freq);
1220 }
1221 rtable.set("TYPE", "FREQ");
1222 if (case_->get_bool("M_ORTHONORMAL", true) == true) {
1223 // If "M_ORTHONORMAL" is absent, true is supposed.
1224 // Eigenvectors M-orthogonalized: MODK=omega**2, MODM=1
1225 rtable.set("MODK", eigenvalue);
1226 rtable.set("MODM", 1.0);
1227 domain->save_global_dof(name, eigenvector, rtable);
1228 } else {
1229 // Eigenvectors normalized to 1: scaled by 1/inv_norm_inf
1230 T inv_norm_inf = 1.0 / eigenvector.norm_inf();
1231 rtable.set("MODK", inv_norm_inf * inv_norm_inf * eigenvalue);
1232 rtable.set("MODM", inv_norm_inf * inv_norm_inf);
1233 b2linalg::Vector<T, b2linalg::Vdense> tmp;
1234 tmp = inv_norm_inf * eigenvector;
1235 domain->save_global_dof(
1236 name, (const b2linalg::Vector<T, b2linalg::Vdense_constref>)(tmp), rtable);
1237 }
1238 }
1239
1240 void terminate_case(bool success, const RTable& attributes) override {
1241 terminate_gradient();
1242 RTable& rtable = case_->get_solution_rtable();
1243 store_value_in_rtable_and_clear(rtable);
1244 rtable.update(attributes);
1245 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1246 rtable.set("NMODES", nmodes);
1247 rtable.set("DOF_SOL", "MODE");
1248 add_gradient_solution_case_information(rtable);
1249 SetName name("SOLUTION");
1250 name.case_(case_->get_id());
1251 database->set(name, rtable);
1252 database->sync();
1253 }
1254
1255 using type_t = ObjectTypeComplete<SolutionFreeVibration, Solution::type_t>;
1256 static type_t type;
1257
1258private:
1259 int nmodes;
1260};
1261
1262template <typename T>
1263typename SolutionFreeVibration<T>::type_t SolutionFreeVibration<T>::type(
1264 "SolutionFreeVibration", type_name<T>(),
1265 StringList(
1266 "FREE_VIBRATION_SOLUTION", "VIBRATION", "FREE_VIBRATION",
1267 "FREQUENCY_DEPENDENT_FREE_VIBRATION"),
1268 module, &b2000::Solution::type);
1269
1270template <typename T>
1271class SolutionLinearisedPrebuckling : public b2000::b2dbv3::GradientSolution,
1272 public b2000::SolutionLinearisedPrebuckling<T> {
1273public:
1274 void set_case(b2000::Case& case__) override {
1275 GradientSolution::set_case(case__);
1276 nmodes = 0;
1277 }
1278
1279 void set_subcase_id(int subcase_id_) override {
1280 b2000::b2dbv3::GradientSolution::set_subcase_id(subcase_id_);
1281 }
1282
1283 bool compute_dof() const override { return true; }
1284
1285 void set_dof(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1286 SetName sol = case_->get_string("DOF_SOL", std::string("DISP"));
1287 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1288 database->add_field_entry(
1289 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
1290 database->add_field_entry(
1291 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
1292 domain->save_global_dof(sol, dof);
1293 }
1294
1295 void get_dof(b2linalg::Vector<T, b2linalg::Vdense_ref> dof) override {
1296 SetName sol = case_->get_string("DOF_SOL", std::string("DISP"));
1297 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1298 domain->load_global_dof(sol, dof);
1299 }
1300
1301 bool compute_natural_boundary_condition() const override { return true; }
1302
1303 void set_natural_boundary_condition(
1304 const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1305 SetName sol = case_->get_string("NBC_SOL", "FORC");
1306 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1307 database->add_field_entry(
1308 "NBC Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
1309 database->add_field_entry(
1310 "NBC Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1311 domain->save_global_dof(sol, dof);
1312 }
1313
1314 bool compute_constraint_value() const override { return case_->get_bool("COMPUTE_RCFO", true); }
1315
1316 void set_constraint_value(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
1317 SetName sol = case_->get_string("RESIDUUM_SOL", "REFO");
1318 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1319 database->add_field_entry(
1320 "Residuum Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
1321 database->add_field_entry(
1322 "Residuum Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1323 domain->save_global_dof(sol, dof);
1324 }
1325 void which_modes(int& number_of_mode, char which[3], T& sigma) const override {
1326 sigma = case_->get_csda_double("SHIFT", 0.0);
1327 number_of_mode = case_->get_int("NMODES", 1);
1328 std::string whichs = case_->get_string("WHICH_EIGENVALUES", "SM");
1329 strcpy(which, whichs.c_str());
1330 }
1331
1332 void set_mode(
1333 int mode, T eigenvalue,
1334 const b2linalg::Vector<T, b2linalg::Vdense_constref>& eigenvector) override {
1335 database->add_field_entry("Buckling modes", "BMODE", "NODE", "BRANCH", false, false, false);
1336 database->add_field_entry(
1337 "Buckling modes", "BMODE_E", "CELL", "BRANCH", false, false, false);
1338 SetName name("BMODE");
1339 if (name.case_undef()) { name.case_(case_->get_id()); }
1340 RTable rtable;
1341 rtable.set("EIGVAL", eigenvalue);
1342 rtable.set("TYPE", "Crit load factor");
1343 {
1344 b2linalg::Vector<T, b2linalg::Vdense> tmp;
1345 tmp = (static_cast<T>(1.0) / eigenvector.norm_inf()) * eigenvector;
1346 domain->save_global_dof(
1347 name.subcase(mode + 1),
1348 (const b2linalg::Vector<T, b2linalg::Vdense_constref>)(tmp), rtable);
1349 }
1350 nmodes = mode + 1;
1351 }
1352
1353 void terminate_case(bool success, const RTable& attributes) override {
1354 terminate_gradient();
1355 RTable& rtable = case_->get_solution_rtable();
1356 store_value_in_rtable_and_clear(rtable);
1357 rtable.update(attributes);
1358 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1359 rtable.set("NMODES", nmodes);
1360 rtable.set(
1361 "DOF_SOL", SetName(case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
1362 rtable.set(
1363 "NBC_SOL", SetName(case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
1364 if (compute_constraint_value()) {
1365 rtable.set(
1366 "RESIDUUM_SOL",
1367 SetName(case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
1368 }
1369 add_gradient_solution_case_information(rtable);
1370 SetName name("SOLUTION");
1371 name.case_(case_->get_id());
1372 database->set(name, rtable);
1373 database->sync();
1374 }
1375
1376 using type_t = ObjectTypeComplete<SolutionLinearisedPrebuckling, Solution::type_t>;
1377 static type_t type;
1378
1379private:
1380 int nmodes;
1381};
1382template <typename T>
1383typename SolutionLinearisedPrebuckling<T>::type_t SolutionLinearisedPrebuckling<T>::type(
1384 "SolutionLinearisedPrebuckling", type_name<T>(),
1385 StringList("LINEARISED_PREBUCKLING_SOLUTION", "L_BIFURCATION", "LINEARISED_PREBUCKLING"),
1386 module, &b2000::Solution::type);
1387
1388template <typename T>
1389class SolutionIterative : public GradientSolution, virtual public b2000::SolutionIterative {
1390public:
1391 SolutionIterative()
1392 : cycle(0),
1393 subcycle(0),
1394 end_of_stage(false),
1395 sync_db_interval(-1),
1396 current_time(0),
1397 time_at_start_stage(0) {}
1398
1399 void set_case(b2000::Case& case__) override {
1400 GradientSolution::set_case(case__);
1401 // std::cerr << "SolutionIterative set_case for case__=" << case__.get_id()
1402 // << " case_=" << case_->get_id()
1403 // << " stage_id=" << case__.get_stage_id()
1404 // << " get_stage_number=" << case__.get_stage_number()
1405 // << " cycle=" << cycle
1406 // << " current_time=" << current_time
1407 // << std::endl;
1408 if (case_ != &case__) {
1409 // Is this a new stage???
1410 cycle = case__.get_int("STEP_INIT_NUM", 0);
1411 time_at_start_stage = 0;
1412 current_time = 0;
1413 // std::cerr << " case_ != &case__, cycle=" << cycle << std::endl;
1414 } else {
1415 // std::cerr << " case_ == &case__ cycle=" << cycle
1416 // << " time_at_start_stage=" << time_at_start_stage
1417 // << " current_time=" << current_time
1418 // << " STEP_INIT_NUM=" << case__.get_int("STEP_INIT_NUM", 0)
1419 // << std::endl;
1420 int c = case__.get_int("STEP_INIT_NUM", 0);
1421 if (c > 0) {
1422 logging::Logger& logger = logging::get_logger("solver");
1423 logger << logging::info << "Restart at step=" << c << logging::LOGFLUSH;
1424 cycle = c;
1425 }
1426 }
1427 sync_db_interval = case__.get_int("SYNC_DB_INTERVAL", -1);
1428 }
1429
1430 void new_cycle(bool end_of_stage_ = false) override {
1431 rtable.clear();
1432 ++cycle;
1433 subcycle = 0;
1434 end_of_stage = end_of_stage_;
1435 GradientSolution::set_cycle(cycle, end_of_stage_);
1436 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1437 s != GradientSolution::sub_solution.end(); ++s) {
1438 dynamic_cast<b2000::SolutionIterative&>(**s).new_cycle(end_of_stage_);
1439 }
1440 }
1441
1442 void set_system_null_space(
1443 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& nks, bool normalize = true,
1444 const std::string suffix = "NS") override {
1445 SetName sol = case_->get_string("DOF_SOL", std::string("DISP")) + "_" + suffix;
1446 if (sol.case_undef()) { sol.case_(case_->get_id()); }
1447 sol.cycle(cycle + 1);
1448 database->add_field_entry(
1449 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
1450 database->add_field_entry(
1451 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
1452 b2linalg::Vector<T> tmp;
1453 for (size_t j = 0; j != nks.size2(); ++j) {
1454 sol.subcase(j + 1);
1455 tmp = nks[j];
1456 if (normalize) { tmp.normalize(); }
1457 domain->save_global_dof(sol, b2linalg::Vector<T, b2linalg::Vdense_constref>(tmp));
1458 }
1459 }
1460
1461 void terminate_cycle(double time, double dtime, const RTable& attributes) override {
1462 current_time = time;
1463 terminate_gradient();
1464 store_value_in_rtable_and_clear(rtable);
1465 rtable.set("STAGE_ID", case_->get_stage_id());
1466 rtable.set("STAGE_NUMBER", case_->get_stage_number() + 1);
1467 rtable.set("TIME", time + time_at_start_stage);
1468 rtable.set("STAGE_TIME", time);
1469 rtable.set("DELTA_TIME", dtime);
1470 rtable.update(attributes);
1471 SetName name("SOLUTION-STEP");
1472 name.cycle(cycle);
1473 name.case_(case_->get_id());
1474 database->set(name, rtable);
1475 if (sync_db_interval == -1) {
1476 database->sync_delayed();
1477 } else if (cycle % std::max(1, sync_db_interval) == 0) {
1478 database->sync();
1479 }
1480 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1481 s != GradientSolution::sub_solution.end(); ++s) {
1482 dynamic_cast<b2000::SolutionIterative&>(**s).terminate_cycle(time, dtime, attributes);
1483 }
1484 }
1485
1486 void terminate_stage(bool success = true) override {
1487 RTable& rtable = case_->get_solution_rtable();
1488 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1489 rtable.set("LAST_STEP", cycle);
1490 add_gradient_solution_case_information(rtable);
1491 SetName name("SOLUTION-STAGE");
1492 name.case_(case_->get_id());
1493 name.subcase(case_->get_stage_number() + 1);
1494 database->set(name, rtable);
1495 database->sync();
1496 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1497 s != GradientSolution::sub_solution.end(); ++s) {
1498 dynamic_cast<b2000::SolutionIterative&>(**s).terminate_stage(success);
1499 }
1500 time_at_start_stage = current_time;
1501 }
1502
1503 size_t get_cycle_id() override { return cycle; }
1504
1505 size_t get_subcycle_id() override { return subcycle; }
1506
1507 std::string get_domain_state_id() override {
1508 return SetName("STATE.GLOB").cycle(cycle).case_(case_->get_id());
1509 }
1510
1511 using nbc_t = TypedNaturalBoundaryCondition<T, b2linalg::Msym_compressed_col_update_inv>;
1512
1513protected:
1514 RTable& get_cycle_rtable() { return rtable; }
1515 int cycle;
1516 int subcycle;
1517 bool end_of_stage;
1518 int sync_db_interval;
1519 double current_time;
1520 double time_at_start_stage;
1521 RTable rtable;
1522};
1523
1524template <typename T>
1525class SolutionStaticNonlinear : virtual public SolutionIterative<T>,
1526 public b2000::SolutionStaticNonlinear<T> {
1527public:
1528 void terminate_case(bool success, const RTable& attributes) override {
1529 RTable rtable;
1530 rtable.update(attributes);
1531 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1532 rtable.set("NSTEPS", this->cycle);
1533 rtable.set("NSTAGES", this->case_->get_stage_number() + 1);
1534 rtable.set(
1535 "DOF_SOL",
1536 SetName(this->case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
1537 rtable.set(
1538 "NBC_SOL",
1539 SetName(this->case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
1540 rtable.set(
1541 "RESIDUUM_SOL",
1542 SetName(this->case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
1543 this->add_gradient_solution_case_information(rtable);
1544 SetName name("SOLUTION");
1545 name.case_(this->case_->get_id());
1546 this->database->set(name, rtable);
1547 this->database->sync();
1548 this->cycle = 0;
1549 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1550 s != GradientSolution::sub_solution.end(); ++s) {
1551 dynamic_cast<b2000::SolutionIterative&>(**s).terminate_case(success, attributes);
1552 }
1553 }
1554
1555 void set_dof(
1556 double time, const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
1557 const b2linalg::Vector<T, b2linalg::Vdense_constref>& nbc,
1558 const b2linalg::Vector<T, b2linalg::Vdense_constref>& residue,
1559 bool unconverged_subcycle = false) {
1560 RTable rtable;
1561 rtable.set("STAGE_ID", this->case_->get_stage_id());
1562 rtable.set("STAGE", this->case_->get_stage_number());
1563 rtable.set("LAMBDA", time + this->time_at_start_stage);
1564 rtable.set("STAGE_LAMBDA", time);
1565 if (!dof.is_null()) {
1566 SetName name = this->case_->get_string("DOF_SOL", "DISP");
1567 if (unconverged_subcycle) {
1568 name.cycle(this->cycle + 1);
1569 name.subcycle(this->subcycle + 1);
1570 } else {
1571 name.cycle(this->cycle);
1572 }
1573 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1574 this->database->add_field_entry(
1575 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
1576 this->database->add_field_entry(
1577 "DOF Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
1578 this->domain->save_global_dof(name, dof, rtable);
1579 }
1580 if (!nbc.is_null()) {
1581 SetName name = this->case_->get_string("NBC_SOL", "FORC");
1582 if (unconverged_subcycle) {
1583 name.cycle(this->cycle + 1);
1584 name.subcycle(this->subcycle + 1);
1585 } else {
1586 name.cycle(this->cycle);
1587 }
1588 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1589 this->database->add_field_entry(
1590 "NBC Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
1591 this->database->add_field_entry(
1592 "NBC Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1593 this->domain->save_global_dof(name, nbc, rtable);
1594 }
1595 if (!residue.is_null()) {
1596 SetName name = this->case_->get_string("RESIDUUM_SOL", "REFO");
1597 if (unconverged_subcycle) {
1598 name.cycle(this->cycle + 1);
1599 name.subcycle(this->subcycle + 1);
1600 } else {
1601 name.cycle(this->cycle);
1602 }
1603 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1604 this->database->add_field_entry(
1605 "Residuum Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
1606 this->database->add_field_entry(
1607 "Residuum Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true,
1608 true);
1609 this->domain->save_global_dof(name, residue, rtable);
1610 }
1611 if (this->case_->get_int("SAVE_NBC_COMP", 0) == 1) {
1612 SetName name;
1613 name.cycle(this->cycle);
1614 name.case_(this->case_->get_id());
1615 dynamic_cast<typename SolutionIterative<T>::nbc_t&>(
1616 this->model->get_natural_boundary_condition(
1617 SolutionIterative<T>::nbc_t::type.get_subtype(
1618 "TypedNaturalBoundaryConditionListComponent"
1619 + type_name<T, b2linalg::Msym_compressed_col_update_inv>())))
1620 .save_nonlinear_value(name, dof, time);
1621 }
1622 if (unconverged_subcycle) {
1623 ++this->subcycle;
1624 this->database->sync();
1625 }
1626 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1627 s != GradientSolution::sub_solution.end(); ++s) {
1628 dynamic_cast<b2000::SolutionStaticNonlinear<T>&>(**s).set_dof(
1629 time, dof, nbc, residue, unconverged_subcycle);
1630 }
1631 }
1632
1633 void get_dof(double& time, b2linalg::Vector<T, b2linalg::Vdense_ref> dof) {
1634 SetName name = this->case_->get_string("DOF_SOL", "DISP");
1635 name.cycle(this->cycle);
1636 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1637 RTable desc;
1638 this->domain->load_global_dof(name, dof, &desc);
1639 time = desc.get_double("LAMBDA");
1640 }
1641
1642 void set_equilibrium(
1643 double time, const b2linalg::Vector<T, b2linalg::Vdense_constref>& residue) {
1644 RTable rtable;
1645 rtable.set("STAGE_ID", this->case_->get_stage_id());
1646 rtable.set("STAGE", this->case_->get_stage_number());
1647 rtable.set("LAMBDA", time + this->time_at_start_stage);
1648 rtable.set("STAGE_LAMBDA", time);
1649 SetName name = "NR_EQUILIBRIUM";
1650 name.cycle(this->cycle + 1);
1651 name.subcycle(this->subcycle + 1);
1652 name.case_(this->case_->get_id());
1653 this->database->add_field_entry(
1654 "Residue", name.get_gname(), "NODE", "BRANCH", false, true, true);
1655 this->database->add_field_entry(
1656 "Residue", name.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1657 this->domain->save_global_dof(name, residue, rtable);
1658 for (GradientSolution::sub_solution_t::iterator s = GradientSolution::sub_solution.begin();
1659 s != GradientSolution::sub_solution.end(); ++s) {
1660 dynamic_cast<b2000::SolutionStaticNonlinear<T>&>(**s).set_equilibrium(time, residue);
1661 }
1662 }
1663
1664 using type_t = ObjectTypeComplete<SolutionStaticNonlinear, Solution::type_t>;
1665 static type_t type;
1666};
1667
1668template <typename T>
1669typename SolutionStaticNonlinear<T>::type_t SolutionStaticNonlinear<T>::type(
1670 "SolutionStaticNonlinear", type_name<T>(),
1671 StringList("STATIC_NONLINEAR_SOLUTION", "NONLINEAR", "COLLAPSE"), module,
1672 &b2000::Solution::type);
1673
1674template <typename T>
1675class SolutionDynamicNonlinear : virtual public SolutionIterative<T>,
1676 public b2000::SolutionDynamicNonlinear<T> {
1677public:
1678 SolutionDynamicNonlinear() {
1679 last_rate_bc_work = 0;
1680 last_bc_work = 0;
1681 }
1682
1683 void terminate_case(bool success, const RTable& attributes) override {
1684 RTable& rtable = this->case_->get_solution_rtable();
1685 rtable.update(attributes);
1686 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1687 rtable.set("NSTEPS", this->cycle);
1688 rtable.set("NSTAGES", this->case_->get_stage_number() + 1);
1689 rtable.set(
1690 "DOF_SOL",
1691 SetName(this->case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
1692 rtable.set(
1693 "DOF_DOT_SOL",
1694 SetName(this->case_->get_string("DOF_SOL", std::string("DISP")) + "D").get_gname());
1695 rtable.set(
1696 "DOF_DOTDOT_SOL",
1697 SetName(this->case_->get_string("DOF_SOL", std::string("DISP")) + "DD").get_gname());
1698 rtable.set(
1699 "NBC_SOL",
1700 SetName(this->case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
1701 rtable.set(
1702 "RESIDUUM_SOL",
1703 SetName(this->case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
1704 this->add_gradient_solution_case_information(rtable);
1705 SetName name("SOLUTION");
1706 name.case_(this->case_->get_id());
1707 this->database->set(name, rtable);
1708 this->database->sync();
1709 this->cycle = 0;
1710 }
1711
1712 void set_solution(
1713 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
1714 const b2linalg::Vector<T, b2linalg::Vdense_constref>& nbc,
1715 const b2linalg::Vector<T, b2linalg::Vdense_constref>& residue,
1716 bool unconverged_subcycle = false) override {
1717 RTable rtable;
1718 rtable.set("STAGE_ID", this->case_->get_stage_id());
1719 rtable.set("STAGE", this->case_->get_stage_number());
1720 rtable.set("TIME", time + this->time_at_start_stage);
1721 rtable.set("STAGE_TIME", time);
1722 if (!dof.is_null()) {
1723 std::string dname = "";
1724 for (size_t i = 0; i != dof.size2(); ++i) {
1725 SetName name = this->case_->get_string("DOF_SOL", "DISP") + dname;
1726 if (unconverged_subcycle) {
1727 name.cycle(this->cycle + 1);
1728 name.subcycle(this->subcycle + 1);
1729 } else {
1730 name.cycle(this->cycle);
1731 }
1732 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1733 this->database->add_field_entry(
1734 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
1735 this->database->add_field_entry(
1736 "DOF Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, false,
1737 false);
1738 this->domain->save_global_dof(name, dof[i], rtable);
1739 dname += "D";
1740 }
1741 }
1742 if (!nbc.is_null()) {
1743 SetName name = this->case_->get_string("NBC_SOL", "FORC");
1744 if (unconverged_subcycle) {
1745 name.cycle(this->cycle + 1);
1746 name.subcycle(this->subcycle + 1);
1747 } else {
1748 name.cycle(this->cycle);
1749 }
1750 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1751 this->database->add_field_entry(
1752 "NBC Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
1753 this->database->add_field_entry(
1754 "NBC Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
1755 this->domain->save_global_dof(name, nbc, rtable);
1756 }
1757 if (!residue.is_null()) {
1758 SetName name = this->case_->get_string("RESIDUUM_SOL", "REFO");
1759 if (unconverged_subcycle) {
1760 name.cycle(this->cycle + 1);
1761 name.subcycle(this->subcycle + 1);
1762 } else {
1763 name.cycle(this->cycle);
1764 }
1765 if (name.case_undef()) { name.case_(this->case_->get_id()); }
1766 this->database->add_field_entry(
1767 "Residuum Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
1768 this->database->add_field_entry(
1769 "Residuum Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true,
1770 true);
1771 this->domain->save_global_dof(name, residue, rtable);
1772 }
1773 if (this->case_->get_int("SAVE_NBC_COMP", 0) == 1) {
1774 SetName name;
1775 if (unconverged_subcycle) {
1776 name.cycle(this->cycle + 1);
1777 name.subcycle(this->subcycle + 1);
1778 } else {
1779 name.cycle(this->cycle);
1780 }
1781 name.case_(this->case_->get_id());
1782 dynamic_cast<typename SolutionIterative<T>::nbc_t&>(
1783 this->model->get_natural_boundary_condition(
1784 SolutionIterative<T>::nbc_t::type.get_subtype(
1785 "TypedNaturalBoundaryConditionListComponent"
1786 + type_name<T, b2linalg::Msym_compressed_col_update_inv>())))
1787 .save_nonlinear_value(name, dof[0], time);
1788 }
1789 if (unconverged_subcycle) {
1790 ++this->subcycle;
1791 this->database->sync();
1792 } else if (!dof.is_null()) {
1793 // This branch is only tested for T = double
1794 T tmp = 0;
1795 if (!nbc.is_null()) { tmp += nbc * dof[1]; }
1796 if (!residue.is_null()) { tmp += residue * dof[1]; }
1797 last_bc_work +=
1798 (last_rate_bc_work + tmp) * static_cast<T>((time - this->current_time) / 2);
1799 last_rate_bc_work = tmp;
1800 RTable& rtable = this->get_cycle_rtable();
1801 rtable.set("RATE_BC_WORK", last_rate_bc_work);
1802 rtable.set("BC_WORK", last_bc_work);
1803 }
1804 }
1805
1806 using type_t = ObjectTypeComplete<SolutionDynamicNonlinear, Solution::type_t>;
1807 static type_t type;
1808
1809private:
1810 T last_rate_bc_work;
1811 T last_bc_work;
1812};
1813
1814template <typename T>
1815typename SolutionDynamicNonlinear<T>::type_t SolutionDynamicNonlinear<T>::type(
1816 "SolutionDynamicNonlinear", type_name<T>(),
1817 StringList(
1818 "DYNAMIC_NONLINEAR_SOLUTION", "DYNAMIC_NONLINEAR", "TRANSIENT_NONLINEAR",
1819 "DYNAMIC_LINEAR"),
1820 module, &b2000::Solution::type);
1821
1822class SolutionModelReduction
1823 : public b2000::b2dbv3::GradientSolution,
1824 public b2000::SolutionModelReduction<double, b2linalg::Mpacked> /*,
1825 public b2000::SolutionModelReduction<double, b2linalg::Mrectangle>,
1826 public b2000::SolutionModelReduction<std::complex<double>, b2linalg::Mpacked>,
1827 public b2000::SolutionModelReduction<std::complex<double>, b2linalg::Mrectangle> */
1828{
1829public:
1830 void get_interface_dof(b2linalg::Index& interface_dof) override {
1831 SetName name;
1832 name.case_(case_->get_string("SET_INTERFACE_ID", "interface"));
1833
1834 name.gname("FACESET");
1835 {
1836 std::set<size_t> interface_dof_set;
1837 for (b2dbv3::Domain::list_branch_t::const_iterator i = domain->list_branch.begin();
1838 i != domain->list_branch.end(); ++i) {
1839 name.branch(i->id);
1840 if (database->has_key(name)) {
1841 b2linalg::Matrix<int, b2linalg::Mrectangle> faceset(2, 0);
1842 database->get(name, faceset);
1843 b2linalg::Vector<double, b2linalg::Vdense> internal_coor(2);
1844 for (size_t j = 0; j != faceset.size2(); ++j) {
1845 b2linalg::Index dof_numbering;
1846 b2linalg::Matrix<double, b2linalg::Mrectangle> d_value_d_dof;
1847 if (auto tve = dynamic_cast<TypedElement<double>*>(
1848 i->list_all_elements[faceset(0, j) - 1]);
1849 tve != nullptr) {
1850 tve->face_field_value(
1851 faceset(1, j), "?", internal_coor,
1852 b2linalg::Matrix<double, b2linalg::Mrectangle_constref>::null, 0,
1853 b2linalg::Vector<double, b2linalg::Vdense>::null,
1854 b2linalg::Matrix<double, b2linalg::Mrectangle>::null,
1855 dof_numbering, d_value_d_dof, b2linalg::Index::null);
1856 }
1857 interface_dof_set.insert(dof_numbering.begin(), dof_numbering.end());
1858 }
1859 }
1860 }
1861
1862 if (!interface_dof_set.empty()) {
1863 interface_dof.assign(interface_dof_set.begin(), interface_dof_set.end());
1864 return;
1865 }
1866 }
1867
1868 name.gname("NODESET");
1869 domain->load_node_set(name, nodes_set);
1870 if (!nodes_set.empty()) {
1871 size_t s = 0;
1872 for (size_t n = 0; n != nodes_set.size(); ++n) {
1873 s += nodes_set[n]->get_number_of_dof();
1874 }
1875 interface_dof.assign(s, 0);
1876 size_t* i = &interface_dof[0];
1877 for (size_t n = 0; n != nodes_set.size(); ++n) {
1878 i = nodes_set[n]->get_global_dof_numbering(i);
1879 }
1880 return;
1881 }
1882
1883 DBError() << "The database don't have NODESET of FACESET interface set" << THROW;
1884 }
1885
1886 void set_reduction(
1887 const size_t nb_interface_dof,
1888 const b2linalg::Matrix<double, b2linalg::Mrectangle>& basis,
1889 std::vector<b2linalg::Matrix<double, b2linalg::Mpacked>>& reduced_matrix) override {
1890 database->add_field_entry(
1891 "ReducedBasis", "REDUCED_BASIS", "NODE", "BRANCH", false, false, false);
1892 database->add_field_entry(
1893 "ReducedBasis", "REDUCED_BASIS_E", "CELL", "BRANCH", false, false, false);
1894 SetName name("REDUCED_BASIS");
1895 name.case_(case_->get_id());
1896 {
1897 RTable desc;
1898 size_t j = 0;
1899 if (!nodes_set.empty()) {
1900 RTable desc_i;
1901 for (size_t n = 0; n != nodes_set.size(); ++n) {
1902 std::pair<size_t, size_t> ibni =
1903 domain->get_internal_branch_and_node_id(*nodes_set[n]);
1904 desc_i.set("NI_BRANCH", int(ibni.first) + 1);
1905 desc_i.set("NI_NODE", int(ibni.second) + 1);
1906 const int nd = nodes_set[n]->get_number_of_dof();
1907 for (int i = 0; i != nd; ++i, ++j) {
1908 desc_i.set("NI_DOF", i + 1);
1909 name.cycle(j + 1);
1910 domain->save_global_dof(name, basis[j], desc_i);
1911 }
1912 }
1913 }
1914 for (; j != basis.size2(); ++j) {
1915 name.cycle(j + 1);
1916 domain->save_global_dof(name, basis[j], desc);
1917 }
1918 }
1919 name.gname("REDUCED_SYSTEM");
1920 for (size_t i = 0; i != reduced_matrix.size(); ++i) {
1921 if (reduced_matrix[i].size1() != 0) {
1922 name.cycle(i);
1923 database->set(name, reduced_matrix[i]);
1924 }
1925 }
1926 current_case_nb_interface_dof = nb_interface_dof;
1927 current_case_nb_internal_dof = basis.size2() - nb_interface_dof;
1928 }
1929
1930 void terminate_case(bool success, const RTable& attributes) override {
1931 terminate_gradient();
1932 RTable& rtable = case_->get_solution_rtable();
1933 store_value_in_rtable_and_clear(rtable);
1934 rtable.update(attributes);
1935 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1936 rtable.set("NB_INTERFACE_DOF", int(current_case_nb_interface_dof));
1937 rtable.set("NB_INTERNAL_DOF", int(current_case_nb_internal_dof));
1938 SetName name("SOLUTION");
1939 name.case_(case_->get_id());
1940 database->set(name, rtable);
1941 database->sync();
1942 }
1943
1944 using type_t = ObjectTypeComplete<SolutionModelReduction, Solution::type_t>;
1945 static type_t type;
1946
1947private:
1948 std::vector<Node*> nodes_set;
1949 size_t current_case_nb_interface_dof;
1950 size_t current_case_nb_internal_dof;
1951};
1952
1953class SolutionBeamCrossSection : public b2000::SolutionBeamCrossSection,
1954 public b2000::b2dbv3::GradientSolution {
1955public:
1956 void terminate_case(bool success, const RTable& attributes) override {
1957 terminate_gradient();
1958 RTable& rtable = case_->get_solution_rtable();
1959 store_value_in_rtable_and_clear(rtable);
1960 rtable.update(attributes);
1961 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
1962 SetName name("SOLUTION");
1963 name.case_(case_->get_id());
1964 database->set(name, rtable);
1965 database->sync();
1966 }
1967
1968 void save_field(
1969 SetName& name, const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& m) {
1970 database->add_field_entry(
1971 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
1972
1973 RTable desc;
1974 desc.set("SYSTEM", "BRANCH");
1975 desc.set("TYPE", "NODE");
1976 size_t pos = 0;
1977
1978 for (size_t branch = 0; branch != domain->list_branch.size(); ++branch) {
1979 name.branch(domain->list_branch[branch].id);
1980 if (domain->list_branch[branch].max_number_of_dof_by_node == 0) { continue; }
1981 database->set(
1982 name,
1983 m(b2linalg::Interval(pos, pos + domain->list_branch[branch].list_nodes.size())));
1984 pos += domain->list_branch[branch].list_nodes.size();
1985 database->set_desc(name, desc);
1986 }
1987 }
1988
1989 void set_solution(
1990 const double neutral_point[2], const double inertia_mass, const double inertia_point[2],
1991 const double inertia_moment[3],
1992 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
1993 const b2linalg::Matrix<double> strain[7], const b2linalg::Matrix<double> stress[7],
1994 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& constitutive_matrix,
1995 const b2linalg::Vector<double, b2linalg::Vdense_constref>& constant_load) override {
1996 SetName name("BCS_CONSTITUTIVE_MATRIX");
1997 if (name.case_undef()) { name.case_(case_->get_id()); }
1998 database->set(name, constitutive_matrix);
1999 {
2000 RTable rt;
2001 std::vector<double> tmp(neutral_point, neutral_point + 2);
2002 rt.set("NEUTRAL_AXIS", tmp);
2003 rt.set("MASS", inertia_mass);
2004 tmp.assign(inertia_point, inertia_point + 2);
2005 rt.set("MASS_CENTER", tmp);
2006 tmp.assign(inertia_moment, inertia_moment + 3);
2007 rt.set("MASS_INERTIA_MOMENTS", tmp);
2008 {
2009 std::vector<double> tmp(&constant_load[0], &constant_load[6]);
2010 rt.set("CONSTANT_FORCE_AND_MOMENT", tmp);
2011 }
2012 database->set_desc(name, rt);
2013 }
2014 for (size_t i = 0; i != 7; ++i) {
2015 name.gname(case_->get_string("DOF_SOL", "DISP"));
2016 name.subcase(i + 1);
2017 database->add_field_entry(
2018 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
2019 domain->save_global_dof(name, dof[i]);
2020
2021 name.gname("STRAIN");
2022 save_field(name, strain[i]);
2023 name.gname("STRESS");
2024 save_field(name, stress[i]);
2025 }
2026 }
2027
2028 using type_t = ObjectTypeComplete<SolutionBeamCrossSection, Solution::type_t>;
2029 static type_t type;
2030};
2031
2032class SolutionTaylorExpansionsBuckling : public b2000::b2dbv3::GradientSolution,
2033 public b2000::SolutionTaylorExpansionsBuckling<double> {
2034public:
2035 SolutionTaylorExpansionsBuckling() : nbmodes(0) {}
2036
2037 void which_modes(int& number_of_mode, char which[3], double& sigma) const override {
2038 sigma = case_->get_double("SHIFT", 0);
2039 number_of_mode = case_->get_int("NMODES", 1);
2040 std::string whichs = case_->get_string("WHICH_EIGENVALUES", "RA");
2041 strcpy(which, whichs.c_str());
2042 }
2043
2044 void terminate_case(bool success, const RTable& attributes) override {
2045 RTable& rtable = case_->get_solution_rtable();
2046 store_value_in_rtable_and_clear(rtable);
2047 rtable.update(attributes);
2048 rtable.set("NMODE", nbmodes);
2049 rtable.set(
2050 "DOF_SOL", SetName(case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
2051 rtable.set(
2052 "NBC_SOL", SetName(case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
2053 rtable.set(
2054 "RESIDUUM_SOL",
2055 SetName(case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
2056 add_gradient_solution_case_information(rtable);
2057 SetName name("SOLUTION");
2058 name.case_(case_->get_id());
2059 database->set(name, rtable);
2060 database->sync();
2061 nbmodes = 0;
2062 }
2063
2064 void set_path(
2065 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& path,
2066 const b2linalg::Vector<double, b2linalg::Vdense_constref>& time) override {
2067 RTable rtable;
2068 SetName name = case_->get_string("PDOF_SOL", "PDISP");
2069 if (name.case_undef()) { name.case_(case_->get_id()); }
2070 database->add_field_entry(
2071 "DOF Path Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
2072 database->add_field_entry(
2073 "DOF Path Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
2074 for (size_t j = 0; j != path.size2(); ++j) {
2075 rtable.set("LAMBDA", time[j]);
2076 name.subcase(j + 1);
2077 domain->save_global_dof(name, path[j], rtable);
2078 }
2079 database->sync();
2080 }
2081
2082 void set_solution(
2083 const int mode, const double cc, const double cci,
2084 const b2linalg::Vector<double, b2linalg::Vdense_constref>& dof, const double time,
2085 const b2linalg::Vector<double, b2linalg::Vdense_constref>& nbc,
2086 const b2linalg::Vector<double, b2linalg::Vdense_constref>& residue,
2087 const b2linalg::Vector<double, b2linalg::Vdense_constref>& eigenvector) override {
2088 ++nbmodes;
2089 RTable rtable;
2090 rtable.set("SOLUTION_PATH_COORDINATE", cc);
2091 rtable.set("SOLUTION_PATH_COORDINATE_i", cci);
2092 rtable.set("LAMBDA", time);
2093 if (!dof.is_null()) {
2094 SetName name = case_->get_string("DOF_SOL", "DISP");
2095 name.subcase(mode + 1);
2096 if (name.case_undef()) { name.case_(case_->get_id()); }
2097 database->add_field_entry(
2098 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
2099 database->add_field_entry(
2100 "DOF Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
2101 domain->save_global_dof(name, dof, rtable);
2102 }
2103 if (!nbc.is_null()) {
2104 SetName name = case_->get_string("NBC_SOL", "FORC");
2105 name.subcase(mode + 1);
2106 if (name.case_undef()) { name.case_(case_->get_id()); }
2107 database->add_field_entry(
2108 "NBC Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
2109 database->add_field_entry(
2110 "NBC Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
2111 domain->save_global_dof(name, nbc, rtable);
2112 }
2113 if (!residue.is_null()) {
2114 SetName name = case_->get_string("RESIDUUM_SOL", "REFO");
2115 name.subcase(mode + 1);
2116 if (name.case_undef()) { name.case_(case_->get_id()); }
2117 database->add_field_entry(
2118 "Residuum Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
2119 database->add_field_entry(
2120 "Residuum Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true,
2121 true);
2122 domain->save_global_dof(name, residue, rtable);
2123 }
2124 {
2125 SetName name = case_->get_string("BMODE", "BMODE");
2126 name.subcase(mode + 1);
2127 if (name.case_undef()) { name.case_(case_->get_id()); }
2128 database->add_field_entry(
2129 "Buckling Mode", name.get_gname(), "NODE", "BRANCH", false, false, false);
2130 database->add_field_entry(
2131 "Buckling Mode", name.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
2132 domain->save_global_dof(name, eigenvector, rtable);
2133 }
2134 }
2135
2136 using type_t = ObjectTypeComplete<SolutionTaylorExpansionsBuckling, Solution::type_t>;
2137 static type_t type;
2138
2139private:
2140 int nbmodes;
2141};
2142
2143// Monolithic
2144template <typename T>
2145class SolutionMonolithic : virtual public SolutionIterative<T>,
2146 public b2000::SolutionMonolithic<T> {
2147public:
2148 // bool compute_dof() const override { return true; }
2149
2150 // void set_subcase_id(int subcase_id_) override {
2151 // b2000::b2dbv3::GradientSolution::set_subcase_id(subcase_id_);
2152 // }
2153
2154 void set_dof(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
2155 SetName sol = this->case_->get_string("DOF_SOL", std::string("DISP"));
2156 if (sol.case_undef()) { sol.case_(this->case_->get_id()); }
2157 if (this->subcase_id != 0) { sol.subcase(this->subcase_id); }
2158 this->database->add_field_entry(
2159 "DOF Solution", sol.get_gname(), "NODE", "BRANCH", false, false, false);
2160 this->database->add_field_entry(
2161 "DOF Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, false, false);
2162 this->domain->save_global_dof(sol, dof);
2163 }
2164
2165 // void get_dof(b2linalg::Vector<T, b2linalg::Vdense_ref> dof) override {
2166 // SetName sol = case_->get_string("DOF_SOL", std::string("DISP"));
2167 // if (sol.case_undef()) { sol.case_(case_->get_id()); }
2168 // if (subcase_id != 0) { sol.subcase(subcase_id); }
2169 // domain->load_global_dof(sol, dof);
2170 // }
2171
2172 // bool compute_natural_boundary_condition() const override { return true; }
2173
2174 void set_natural_boundary_condition(
2175 const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override {
2176 SetName sol = this->case_->get_string("NBC_SOL", "FORC");
2177 if (sol.case_undef()) { sol.case_(this->case_->get_id()); }
2178 if (this->subcase_id != 0) { sol.subcase(this->subcase_id); }
2179 this->database->add_field_entry(
2180 "NBC Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
2181 this->database->add_field_entry(
2182 "NBC Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
2183 this->domain->save_global_dof(sol, dof);
2184 }
2185
2186 bool compute_constraint_value() const override {
2187 return (this->case_->get_bool("COMPUTE_RCFO", false));
2188 }
2189
2190 // // TODO change dof input name to F_reac or similar.
2191 // void set_constraint_value(const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof) override
2192 // {
2193 // SetName sol = case_->get_string("RESIDUUM_SOL", "RCFO");
2194 // if (sol.case_undef()) { sol.case_(case_->get_id()); }
2195 // if (subcase_id != 0) { sol.subcase(subcase_id); }
2196 // database->add_field_entry(
2197 // "Residuum Solution", sol.get_gname(), "NODE", "BRANCH", false, true, true);
2198 // database->add_field_entry(
2199 // "Residuum Solution", sol.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
2200 // domain->save_global_dof(sol, dof);
2201 // }
2202
2203 // void terminate_case(bool success, const RTable& attributes) override {
2204 // this->terminate_gradient();
2205 // RTable& rtable = this->case_->get_solution_rtable();
2206 // rtable.update(attributes);
2207 // this->store_value_in_rtable_and_clear(rtable);
2208 // rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
2209 // rtable.set(
2210 // "DOF_SOL",
2211 // SetName(this->case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
2212 // rtable.set(
2213 // "NBC_SOL",
2214 // SetName(this->case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
2215 // if (compute_constraint_value()) {
2216 // rtable.set(
2217 // "RESIDUUM_SOL",
2218 // SetName(this->case_->get_string("RESIDUUM_SOL", std::string("RCFO")))
2219 // .get_gname());
2220 // }
2221 // this->add_gradient_solution_case_information(rtable);
2222 // SetName name("SOLUTION");
2223 // name.case_(this->case_->get_id());
2224 // this->database->set(name, rtable);
2225 // this->database->sync();
2226 // }
2227
2228 void terminate_case(bool success, const RTable& attributes) override {
2229 RTable& rtable = this->case_->get_solution_rtable();
2230 rtable.update(attributes);
2231 rtable.set("TERMINATION", success ? "NORMAL" : "FAILURE");
2232 rtable.set("NSTEPS", this->cycle);
2233 rtable.set("NSTAGES", this->case_->get_stage_number() + 1);
2234 rtable.set(
2235 "DOF_SOL",
2236 SetName(this->case_->get_string("DOF_SOL", std::string("DISP"))).get_gname());
2237 rtable.set(
2238 "DOF_DOT_SOL",
2239 SetName(this->case_->get_string("DOF_SOL", std::string("DISP")) + "D").get_gname());
2240 rtable.set(
2241 "DOF_DOTDOT_SOL",
2242 SetName(this->case_->get_string("DOF_SOL", std::string("DISP")) + "DD").get_gname());
2243 rtable.set(
2244 "NBC_SOL",
2245 SetName(this->case_->get_string("NBC_SOL", std::string("FORC"))).get_gname());
2246 rtable.set(
2247 "RESIDUUM_SOL",
2248 SetName(this->case_->get_string("RESIDUUM_SOL", std::string("RCFO"))).get_gname());
2249 this->add_gradient_solution_case_information(rtable);
2250 SetName name("SOLUTION");
2251 name.case_(this->case_->get_id());
2252 this->database->set(name, rtable);
2253 this->database->sync();
2254 this->cycle = 0;
2255 }
2256
2257 void set_solution(
2258 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
2259 const b2linalg::Vector<T, b2linalg::Vdense_constref>& nbc,
2260 const b2linalg::Vector<T, b2linalg::Vdense_constref>& residue,
2261 bool unconverged_subcycle = false) override {
2262 RTable rtable;
2263 rtable.set("STAGE_ID", this->case_->get_stage_id());
2264 rtable.set("STAGE", this->case_->get_stage_number());
2265 rtable.set("TIME", time + this->time_at_start_stage);
2266 rtable.set("STAGE_TIME", time);
2267 if (!dof.is_null()) {
2268 std::string dname = "";
2269 for (size_t i = 0; i != dof.size2(); ++i) {
2270 SetName name = this->case_->get_string("DOF_SOL", "DISP") + dname;
2271 if (unconverged_subcycle) {
2272 name.cycle(this->cycle + 1);
2273 name.subcycle(this->subcycle + 1);
2274 } else {
2275 name.cycle(this->cycle);
2276 }
2277 if (name.case_undef()) { name.case_(this->case_->get_id()); }
2278 this->database->add_field_entry(
2279 "DOF Solution", name.get_gname(), "NODE", "BRANCH", false, false, false);
2280 this->database->add_field_entry(
2281 "DOF Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, false,
2282 false);
2283 this->domain->save_global_dof(name, dof[i], rtable);
2284 dname += "D";
2285 }
2286 }
2287 if (!nbc.is_null()) {
2288 SetName name = this->case_->get_string("NBC_SOL", "FORC");
2289 if (unconverged_subcycle) {
2290 name.cycle(this->cycle + 1);
2291 name.subcycle(this->subcycle + 1);
2292 } else {
2293 name.cycle(this->cycle);
2294 }
2295 if (name.case_undef()) { name.case_(this->case_->get_id()); }
2296 this->database->add_field_entry(
2297 "NBC Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
2298 this->database->add_field_entry(
2299 "NBC Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true, true);
2300 this->domain->save_global_dof(name, nbc, rtable);
2301 }
2302 if (!residue.is_null()) {
2303 SetName name = this->case_->get_string("RESIDUUM_SOL", "REFO");
2304 if (unconverged_subcycle) {
2305 name.cycle(this->cycle + 1);
2306 name.subcycle(this->subcycle + 1);
2307 } else {
2308 name.cycle(this->cycle);
2309 }
2310 if (name.case_undef()) { name.case_(this->case_->get_id()); }
2311 this->database->add_field_entry(
2312 "Residuum Solution", name.get_gname(), "NODE", "BRANCH", false, true, true);
2313 this->database->add_field_entry(
2314 "Residuum Solution", name.get_gname() + "_E", "CELL", "BRANCH", false, true,
2315 true);
2316 this->domain->save_global_dof(name, residue, rtable);
2317 }
2318 if (this->case_->get_int("SAVE_NBC_COMP", 0) == 1) {
2319 SetName name;
2320 if (unconverged_subcycle) {
2321 name.cycle(this->cycle + 1);
2322 name.subcycle(this->subcycle + 1);
2323 } else {
2324 name.cycle(this->cycle);
2325 }
2326 name.case_(this->case_->get_id());
2327 dynamic_cast<typename SolutionIterative<T>::nbc_t&>(
2328 this->model->get_natural_boundary_condition(
2329 SolutionIterative<T>::nbc_t::type.get_subtype(
2330 "TypedNaturalBoundaryConditionListComponent"
2331 + type_name<T, b2linalg::Msym_compressed_col_update_inv>())))
2332 .save_nonlinear_value(name, dof[0], time);
2333 }
2334 if (unconverged_subcycle) {
2335 ++this->subcycle;
2336 this->database->sync();
2337 } else if (!dof.is_null()) {
2338 // This branch is only tested for T = double
2339 T tmp = 0;
2340 if (!nbc.is_null()) { tmp += nbc * dof[1]; }
2341 if (!residue.is_null()) { tmp += residue * dof[1]; }
2342 last_bc_work +=
2343 (last_rate_bc_work + tmp) * static_cast<T>((time - this->current_time) / 2);
2344 last_rate_bc_work = tmp;
2345 RTable& rtable = this->get_cycle_rtable();
2346 rtable.set("RATE_BC_WORK", last_rate_bc_work);
2347 rtable.set("BC_WORK", last_bc_work);
2348 }
2349 }
2350
2351 using type_t = ObjectTypeComplete<SolutionMonolithic, Solution::type_t>;
2352 static type_t type;
2353
2354private:
2355 T last_rate_bc_work{};
2356 T last_bc_work{};
2357};
2358
2359template <typename T>
2360typename SolutionMonolithic<T>::type_t SolutionMonolithic<T>::type(
2361 "SolutionMonolithic", type_name<T>(),
2362 StringList(
2363 "MLINEAR", "MSTATIC_LINEAR_SOLUTION", "MDYNAMIC_LINEAR", "MTRANSIENT_LINEAR",
2364 "MNONLINEAR", "MSTATIC_NONLINEAR_SOLUTION", "MDYNAMIC_NONLINEAR",
2365 "MTRANSIENT_NONLINEAR"),
2366 module, &b2000::Solution::type);
2367
2368} // namespace b2000::b2dbv3
2369
2370#endif
b2boundary_condition.H
b2000::details::sqrt
csda< T > sqrt(const csda< T > &a)
Square root of a csda number.
Definition b2csda.H:396
b2000::details::operator<
bool operator<(const csda< T1 > &a, const csda< T2 > &b)
Comparison of two csda numbers is performed on the real part only.
Definition b2csda.H:262
b2domain_database.H
THROW
#define THROW
Definition b2exception.H:198
b2000::Case
Definition b2case.H:110
b2000::Case::get_id
virtual std::string get_id() const =0
b2000::Dictionary::get_string
virtual std::string get_string(const std::string &key) const =0
b2000::Dictionary::get_int
virtual int get_int(const std::string &key) const =0
b2000::Dictionary::has_key
virtual bool has_key(const std::string &key) const =0
b2000::GradientContainer
Definition b2solution.H:54
b2000::Model
Definition b2model.H:69
b2000::Model::get_domain
virtual Domain & get_domain()=0
b2000::Solution
Definition b2solution.H:227
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829
b2000::TypeError
GenericException< TypeError_name > TypeError
Definition b2exception.H:325
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314
b2000::DBError
GenericException< DBError_name > DBError
Definition b2exception.H:340