b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2solver_utils.H
1//------------------------------------------------------------------------
2// b2solver_utils.H --
3//
4//
5// written by Mathias Doreille
6// Johannes Wendler <johannes.wendler@dlr.de>
7//
8// Copyright (c) 2008-2016,2017,2018 SMR Engineering & Development SA
9// 2502 Bienne, Switzerland
10//
11// Copyright (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
12// Linder Höhe, 51147 Köln
13//
14// All Rights Reserved. Proprietary source code. The contents of
15// this file may not be disclosed to third parties, copied or
16// duplicated in any form, in whole or in part, without the prior
17// written permission of SMR.
18//------------------------------------------------------------------------
19
20// #define CHECK_SECOND_VARIATION 1
21#define CHECK_USE_CENTRAL_DIFFERENCE 1
22// #define SYMMETRIC_AS_UNSYMMETRIC 1
23// #define SOLVER_UTILS_USE_RECURSIVE_TREE 1
24// #define LOG_CONGESTION 1
25// #define CHECK_MULTITHREADING 1
26
27#ifndef B2SOLVER_UTILS_H_
28#define B2SOLVER_UTILS_H_
29
30#include "b2ppconfig.h"
31#include "model/b2domain.H"
32#include "model/b2model.H"
33#include "model/b2solution.H"
34#include "model/b2solver.H"
35#include "utils/b2linear_algebra.H"
36#ifdef HAVE_LIB_TBB
37#include <tbb/global_control.h>
38#endif // HAVE_LIB_TBB
39#include <cmath>
40#include <ctime>
41#include <list>
42#include <string>
43
44#include "utils/b2logging.H"
45#include "utils/b2threading.H"
46
47namespace {
48
49#ifdef LOG_CONGESTION
50inline uint64_t get_time_stamp() {
51 struct timespec tp;
52 clock_gettime(CLOCK_REALTIME, &tp);
53 return tp.tv_sec * 1000000000ull + tp.tv_nsec;
54}
55
56std::list<std::string> locklog;
57#endif
58
59#ifdef CHECK_MULTITHREADING
60void foo(uint64_t& a) { a = uint64_t(pthread_self()); }
61
62class ApplyFoo {
63 uint64_t* const my_a;
64
65public:
66 void operator()(const tbb::blocked_range<size_t>& r) const {
67 uint64_t* a = my_a;
68 for (size_t i = r.begin(); i != r.end(); ++i) { foo(a[i]); }
69 }
70 ApplyFoo(uint64_t a[]) : my_a(a) {}
71};
72
73size_t check_parallel_for() {
74 std::vector<uint64_t> v(10000000, 0);
75
76 tbb::parallel_for(tbb::blocked_range<size_t>(0, v.size()), ApplyFoo(&v[0]));
77
78 std::set<uint64_t> ids;
79 for (auto i : v) { ids.insert(i); }
80 std::cerr << "solver check_parallel_for: approximate number of threads=" << ids.size()
81 << std::endl;
82 return ids.size();
83}
84#endif /* DEBUG_MULTITHREADING */
85
86} // namespace
87
88namespace b2000 {
89
90template <typename T>
91class DofFieldFlux : public Object {
92public:
93 virtual void add_flux(
94 b2linalg::Index& index, const b2linalg::Vector<T, b2linalg::Vdense>& v) = 0;
95};
96
97namespace solver {
98
99template <
100 typename T, typename MATRIX_FORMAT, typename Args, bool useRayleighDampening = false,
101 bool INERTIA = false>
102struct ElementsTaskLoop {
103 ElementsTaskLoop(Args& args_, b2Mutex& mutex_) : args(args_), mutex(mutex_) {}
104
105 Args& args;
106 b2Mutex& mutex;
107
108 void operator()(size_t begin, size_t end) const {
109 const bool compute_value_with_h = args.linear && !args.dof.is_null()
110 && !args.d_value_d_dof.is_null() && !args.value.is_null();
111 b2linalg::Index index;
112 b2linalg::Vector<T, b2linalg::Vdense> value;
113 b2linalg::Vector<T, b2linalg::Vdense> d_value_d_time;
114 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_value_d_dof(
115 args.d_value_d_dof_e_flags.size());
116 for (size_t i = begin; i != end; ++i) {
117 Element* element = args.list_elements[i];
118 GradientContainer* gc =
119 (args.gradient_container != 0
120 && args.gradient_container->set_current_element(*element)
121 ? args.gradient_container
122 : 0);
123 if (auto te = dynamic_cast<TypedElement<T>*>(element); te != nullptr) {
124 te->get_value(
125 args.model, args.linear, args.equilibrium_solution, args.time,
126 args.delta_time,
127 compute_value_with_h
128 ? b2linalg::Matrix<T, b2linalg::Mrectangle_constref>::null
129 : args.dof, // for rayleigh: b2linalg::Matrix<T,
130 // b2linalg::Mrectangle_constref>(args.dof[0]),
131 gc, args.solver_hints, index,
132 (args.value.is_null() ? b2linalg::Vector<T, b2linalg::Vdense>::null : value),
133 args.d_value_d_dof_e_flags, d_value_d_dof,
134 (args.d_value_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense>::null
135 : d_value_d_time));
136 }
137 if (index.empty()) { continue; }
138 if (compute_value_with_h) {
139 {
140 // TODO add the missing expression in b2linalg to avoid temporary
141 b2linalg::Vector<T> tmp = args.dof[0](index);
142 value += d_value_d_dof[0] * tmp;
143 }
144 if constexpr (!useRayleighDampening) {
145 for (size_t i = 1; i != args.dof.size2(); ++i) {
146 // TODO add the missing expression in b2linalg to avoid temporary
147 b2linalg::Vector<T> tmp = args.dof[i](index);
148 value += d_value_d_dof[i] * tmp;
149 }
150 }
151 }
152 if constexpr (useRayleighDampening) {
153 if (!value.empty()) {
154 // TODO add the missing expression in b2linalg to avoid temporary
155 b2linalg::Vector<T> tmp = args.dof[1](index);
156 value += args.d_fe_d_dofdot_elem[i] * tmp;
157 if (INERTIA) {
158 tmp = args.dof[2](index);
159 value += args.d_fe_d_dofdotdot_elem[i] * tmp;
160 }
161 }
162 }
163 if (args.logger) {
164 *args.logger << logging::data << "elementary index "
165 << logging::DBName("EDOF", element->get_object_name()) << index
166 << " of element id " << element->get_object_name()
167 << logging::LOGFLUSH;
168 for (size_t i = 0; i != args.dof.size2(); ++i) {
169 b2linalg::Vector<T> tmp = args.dof[i](index);
170 *(args.logger)
171 << logging::data << "elementary dof "
172 << logging::DBName("ELDOF", i, element->get_object_name()) << tmp
173 << " of element id " << element->get_object_name() << logging::LOGFLUSH;
174 }
175 if (!args.value.is_null()) {
176 *(args.logger)
177 << logging::data << "elementary value "
178 << logging::DBName("ELFV", element->get_object_name()) << value
179 << " of element id " << element->get_object_name() << logging::LOGFLUSH;
180 }
181 if (!args.d_value_d_dof.is_null()) {
182 for (size_t i = 0; i != d_value_d_dof.size(); ++i) {
183 b2linalg::Matrix<T> tmp_element(d_value_d_dof[i]);
184 *(args.logger) << logging::data << "elementary d_value_d_dof "
185 << logging::DBName("ELSV", i, element->get_object_name())
186 << tmp_element << " of element id "
187 << element->get_object_name() << logging::LOGFLUSH;
188 }
189 }
190 if (!args.d_value_d_time.is_null()) {
191 *(args.logger)
192 << logging::data << "elementary d_value_d_time "
193 << logging::DBName("ELTV", element->get_object_name()) << d_value_d_time
194 << " of element id " << element->get_object_name() << logging::LOGFLUSH;
195 }
196 }
197 if (!d_value_d_dof.empty()) {
198 if constexpr (useRayleighDampening) {
199 d_value_d_dof[0] += T(args.d_value_d_dof_coeff[1]) * args.d_fe_d_dofdot_elem[i];
200 if (INERTIA) {
201 d_value_d_dof[0] +=
202 T(args.d_value_d_dof_coeff[2]) * args.d_fe_d_dofdotdot_elem[i];
203 }
204 } else {
205 if (!args.d_value_d_dof_coeff.is_null() && args.d_value_d_dof_coeff[0] != 1.0) {
206 d_value_d_dof[0] *= args.d_value_d_dof_coeff[0];
207 }
208 for (size_t i = 1; i != d_value_d_dof.size(); ++i) {
209 d_value_d_dof[0] += T(args.d_value_d_dof_coeff[i]) * d_value_d_dof[i];
210 }
211 }
212 }
213 if (!args.d_value_d_time.is_null()) {
214 d_value_d_time += d_value_d_dof[0] * args.d_r_d_time(index);
215 }
216#ifdef LOG_CONGESTION
217 std::ostringstream os;
218 os << get_time_stamp() << " thread_id=" << pthread_self()
219 << " 2nd var=" << args.d_value_d_dof.is_null() << " acquiring" << std::endl;
220#endif
221 {
222 // #ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
223 b2Mutex::scoped_lock lock(mutex);
224 // #else
225 // b2Mutex::scoped_lock lock(args.dof_mutex);
226 // #endif
227
228#ifdef LOG_CONGESTION
229 os << get_time_stamp() << " thread_id=" << pthread_self() << " got it" << std::endl;
230#endif
231 if (!value.empty()) {
232 if (args.dof_field_flux) { args.dof_field_flux->add_flux(index, value); }
233 args.value(index) += value;
234 }
235 if (!args.d_value_d_time.is_null()) {
236 args.d_value_d_time(index) += d_value_d_time;
237 }
238 if (!args.d_value_d_dof.is_null()) {
239 args.d_value_d_dof += args.trans_d_r_d_dof
240 * StructuredMatrix(args.nb_dof, index, d_value_d_dof[0])
241 * transposed(args.trans_d_r_d_dof);
242 }
243#ifdef LOG_CONGESTION
244 os << get_time_stamp() << " thread_id=" << pthread_self() << " releasing"
245 << std::endl;
246 locklog.push_back(os.str());
247#endif
248 }
249 if constexpr (useRayleighDampening) {
250 if (gc) {
251 const double el_volume = gc->get_element_volume();
252 const GradientContainer::InternalElementPosition pos = {0, 0, 0, 0};
253 if (el_volume == 0) { gc->set_current_position(pos, el_volume); }
254 // this do not work and break the other sampling point fields
255 // else
256 // gc->set_current_volume(el_volume);
257 const double inv_el_volume = el_volume == 0 ? 0 : 1 / el_volume;
258
259 b2linalg::Vector<T> v = args.dof[1](index);
260 b2linalg::Vector<T> tmp;
261 tmp = args.d_fe_d_dofdot_elem[i] * v;
262 const double rate_of_dissipated_energie = v * tmp;
263 args.dissipated_energy[i] +=
264 0.5 * args.delta_time
265 * (rate_of_dissipated_energie + args.last_rate_dissipated_energy[i]);
266 args.last_rate_dissipated_energy[i] = rate_of_dissipated_energie;
267 if (INERTIA) { tmp = args.d_fe_d_dofdotdot_elem[i] * v; }
268 const double kinetic_energie = INERTIA ? 0.5 * v * tmp : 0;
269 if (args.energy) {
270 b2Mutex::scoped_lock lock(mutex);
271 args.energy[0] += kinetic_energie;
272 args.energy[1] += rate_of_dissipated_energie;
273 }
274
275 static const b2000::GradientContainer::FieldDescription descr_energy = {
276 "RAYLEIGH_ENERGY_DENSITY",
277 "Kinetic.RateDissipated.Dissipated",
278 "kinetic, rate of dissipated energy and dissipated energy, per unit "
279 "volume, in the undeformed configuration",
280 3,
281 3,
282 -1,
283 0,
284 false,
285 typeid(double)};
286 const double e[3] = {
287 kinetic_energie * inv_el_volume,
288 rate_of_dissipated_energie * inv_el_volume,
289 args.dissipated_energy[i] * inv_el_volume};
290 gc->set_field_value(descr_energy, e);
291 } else if (args.energy) {
292 b2linalg::Vector<T> v = args.dof[1](index);
293 b2linalg::Vector<T> tmp;
294 tmp = args.d_fe_d_dofdot_elem[i] * v;
295 const double rate_of_dissipated_energie = v * tmp;
296 args.dissipated_energy[i] +=
297 0.5 * args.delta_time
298 * (rate_of_dissipated_energie + args.last_rate_dissipated_energy[i]);
299 args.last_rate_dissipated_energy[i] = rate_of_dissipated_energie;
300 if (INERTIA) { tmp = args.d_fe_d_dofdotdot_elem[i] * v; }
301 const double kinetic_energie = INERTIA ? 0.5 * v * tmp : 0;
302 b2Mutex::scoped_lock lock(mutex);
303 args.energy[0] += kinetic_energie;
304 args.energy[1] += rate_of_dissipated_energie;
305 }
306 }
307 }
308 }
309};
310
311template <typename T, typename MATRIX_FORMAT>
312class SolverUtils {
313public:
314 void get_elements_values(
315 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
316 const double delta_time, const EquilibriumSolution equilibrium_solution,
317 const bool linear, const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof,
318 const b2linalg::Vector<T, b2linalg::Vdense_constref>& d_r_d_time,
319 const b2linalg::Vector<double, b2linalg::Vdense_constref>& d_value_d_dof_coeff,
320 Model& model, b2linalg::Vector<T, b2linalg::Vdense_ref> value,
321 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_value_d_dof,
322 b2linalg::Vector<T, b2linalg::Vdense_ref> d_value_d_time,
323 GradientContainer* gradient_container, SolverHints* solver_hints, logging::Logger& logger,
324 DofFieldFlux<T>* dof_field_flux = 0) {
325 Domain& domain = model.get_domain();
326 domain.set_dof(dof);
327 if (list_elements.empty()) { compute_tree(domain, trans_d_r_d_dof); }
328 std::vector<bool> d_value_d_dof_e_flags(
329 d_value_d_dof.is_null() && d_value_d_time.is_null()
330 ? 0
331 : (d_value_d_dof_coeff.is_null() ? 1 : d_value_d_dof_coeff.size()),
332 true);
333 Args args = {
334 dof,
335 time,
336 delta_time,
337 equilibrium_solution,
338 linear,
339 trans_d_r_d_dof,
340 d_r_d_time,
341 d_value_d_dof_coeff,
342 model,
343 domain.get_number_of_dof(),
344 value,
345 d_value_d_dof_e_flags,
346 d_value_d_dof,
347 d_value_d_time,
348 gradient_container,
349 solver_hints,
350 (logger.is_enabled_for(logging::data) ? &logger : 0),
351 dof_field_flux,
352 list_elements,
353 tree};
354
355#ifdef CHECK_MULTITHREADING
356 check_parallel_for();
357#endif
358
359#ifdef LOG_CONGESTION
360 locklog.clear();
361#endif
362 b2Mutex mutex;
363
364 b2RunPossiblyParallel<ElementsTaskLoop<T, MATRIX_FORMAT, Args, false> >(
365 args.tree[0].first, args.tree[0].second, args, mutex);
366
367#ifdef LOG_CONGESTION
368 for (auto s : locklog) { std::cout << s; }
369#endif
370 }
371
372private:
373 struct Args {
374 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof;
375 const double time;
376 const double delta_time;
377 const EquilibriumSolution equilibrium_solution;
378 const bool linear;
379 const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof;
380 const b2linalg::Vector<T, b2linalg::Vdense_constref>& d_r_d_time;
381 const b2linalg::Vector<double, b2linalg::Vdense_constref>& d_value_d_dof_coeff;
382 Model& model;
383 const size_t nb_dof;
384 b2linalg::Vector<T, b2linalg::Vdense_ref>& value;
385 std::vector<bool>& d_value_d_dof_e_flags;
386 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_value_d_dof;
387 b2linalg::Vector<T, b2linalg::Vdense_ref>& d_value_d_time;
388 GradientContainer* gradient_container;
389 SolverHints* solver_hints;
390 logging::Logger* logger;
391 DofFieldFlux<T>* dof_field_flux;
392 std::vector<Element*>& list_elements;
393 std::vector<std::pair<size_t, size_t> >& tree;
394#ifndef SOLVER_UTILS_USE_RECURSIVE_TREE
395 b2Mutex dof_mutex;
396#endif
397 };
398
399 void compute_tree(
400 Domain& domain, const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof) {
401 element_iterator = domain.get_element_iterator();
402 const size_t nb_element = domain.get_number_of_elements();
403 list_elements.reserve(nb_element);
404#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
405 size_t i_ptr = 0;
406 std::vector<std::vector<size_t> > dof(domain.get_number_of_dof());
407 std::vector<size_t> g_ptr(nb_element + 1);
408 g_ptr[0] = 0;
409 std::vector<size_t> g_ind;
410 g_ind.reserve(nb_element);
411 b2linalg::Index index;
412#endif
413 Domain::ElementIterator* i = element_iterator.get();
414 if (i == 0) { TypeError() << "Bad solver type." << THROW; }
415 for (Element* element = i->next(); element != 0; element = i->next()) {
416 list_elements.push_back(element);
417#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
418 element->get_dof_numbering(index);
419 std::set<size_t> s;
420 std::set<size_t> ind(index.begin(), index.end());
421 trans_d_r_d_dof.this_prod_index(index);
422 ind.insert(index.begin(), index.end());
423 for (std::set<size_t>::const_iterator j = ind.begin(); j != ind.end(); ++j) {
424 std::vector<size_t>& d = dof[*j];
425 s.insert(d.begin(), d.end());
426 d.push_back(i_ptr);
427 }
428 g_ind.insert(g_ind.end(), s.begin(), s.end());
429 g_ptr[++i_ptr] = g_ind.size();
430#endif
431 }
432#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
433 dof.clear();
434 std::vector<size_t> inv_perm(nb_element);
435 int nlevel =
436 int(std::log(tbb::global_control::active_value(
437 tbb::global_control::max_allowed_parallelism))
438 / std::log(2));
439 tree_nlevel = nlevel;
440 tree.resize(int(std::pow(2.0, nlevel + 1)));
441 mlevel_nested_dissection_with_separator(
442 nb_element, &g_ptr[0], &g_ind[0], nlevel, &inv_perm[0], &tree[0]);
443 std::vector<Element*> list_elements_tmp(nb_element);
444 for (size_t i = 0; i != nb_element; ++i) {
445 list_elements_tmp[inv_perm[i]] = list_elements[i];
446 }
447 list_elements.swap(list_elements_tmp);
448#else
449 tree.push_back(std::pair<size_t, size_t>(0, list_elements.size()));
450 tree_nlevel = 0;
451#endif
452 }
453
454 std::unique_ptr<Domain::ElementIterator> element_iterator;
455 std::vector<Element*> list_elements;
456 std::vector<std::pair<size_t, size_t> > tree;
457 size_t tree_nlevel;
458};
459
460template <typename T, typename MATRIX_FORMAT, bool INERTIA = true>
461class SolverUtilsOrderTwoRayleighDamping {
462public:
463 void init(
464 const double rayleigh_damping_alpha_, const double rayleigh_damping_beta_, Model& model,
465 logging::Logger& logger) {
466 rayleigh_damping_alpha = rayleigh_damping_alpha_;
467 rayleigh_damping_beta = rayleigh_damping_beta_;
468 last_rate_dissipated_energy.assign(last_rate_dissipated_energy.size(), 0);
469 dissipated_energy.assign(dissipated_energy.size(), 0);
470 global_last_rate_dissipated_energy = 0;
471 global_dissipated_energy = 0;
472 }
473
474 void scale_raylay_damping(double s) {
475 for (size_t i = 0; i != d_fe_d_dofdot_elem.size(); ++i) { d_fe_d_dofdot_elem[i] *= s; }
476 }
477
478 void get_elements_values(
479 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
480 const double delta_time, const EquilibriumSolution equilibrium_solution,
481 const bool linear, const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof,
482 const b2linalg::Vector<T, b2linalg::Vdense_constref> d_r_d_time,
483 const b2linalg::Vector<double, b2linalg::Vdense_constref>& d_value_d_dof_coeff,
484 Model& model, b2linalg::Vector<T, b2linalg::Vdense_ref> value,
485 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_value_d_dof,
486 b2linalg::Vector<T, b2linalg::Vdense_ref> d_value_d_time,
487 GradientContainer* gradient_container, SolverHints* solver_hints, logging::Logger& logger,
488 DofFieldFlux<T>* dof_field_flux = 0, double* energy = 0) {
489 Domain& domain = model.get_domain();
490
491 if (energy) { std::fill_n(energy, 3, 0.0); }
492
493 domain.set_dof(dof);
494 std::vector<bool> d_value_d_dof_e_flags(
495 d_value_d_dof.is_null() && d_value_d_time.is_null() ? 0 : 1, true);
496
497 if (list_elements.empty()) {
498 compute_tree(domain, trans_d_r_d_dof);
499
500 std::vector<bool> d_value_d_dof_flags(3);
501 d_value_d_dof_flags[0] = d_value_d_dof_flags[2] = true;
502 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_fe_d_dof_elem(3);
503 d_fe_d_dofdot_elem.resize(list_elements.size());
504 if (INERTIA) { d_fe_d_dofdotdot_elem.resize(list_elements.size()); }
505 b2linalg::Index index;
506 for (size_t i = 0; i != list_elements.size(); ++i) {
507 Element* e = list_elements[i];
508 if (auto te = dynamic_cast<TypedElement<T>*>(e); te != nullptr) {
509 te->get_value(
510 model,
511 false, // linear
512 false, // equilibrium_solution
513 1, // time
514 0, // delta_time
515 b2linalg::Matrix<T, b2linalg::Mrectangle_constref>::null,
516 0, // gradient_container
517 solver_hints, index, b2linalg::Vector<T, b2linalg::Vdense>::null,
518 d_value_d_dof_flags, d_fe_d_dof_elem,
519 b2linalg::Vector<T, b2linalg::Vdense>::null);
520 }
521 d_fe_d_dofdot_elem[i] = d_fe_d_dof_elem[0];
522 if (INERTIA) { d_fe_d_dofdotdot_elem[i] = d_fe_d_dof_elem[2]; }
523 d_fe_d_dofdot_elem[i] *= rayleigh_damping_beta;
524 d_fe_d_dofdot_elem[i] += rayleigh_damping_alpha * d_fe_d_dof_elem[2];
525 }
526 last_rate_dissipated_energy.assign(d_fe_d_dofdot_elem.size(), 0);
527 dissipated_energy.assign(d_fe_d_dofdot_elem.size(), 0);
528 }
529
530 Args args = {
531 dof,
532 time,
533 delta_time,
534 equilibrium_solution,
535 linear,
536 trans_d_r_d_dof,
537 d_r_d_time,
538 d_value_d_dof_coeff,
539 model,
540 domain.get_number_of_dof(),
541 value,
542 d_value_d_dof_e_flags,
543 d_value_d_dof,
544 d_value_d_time,
545 gradient_container,
546 solver_hints,
547 (logger.is_enabled_for(logging::data) ? &logger : 0),
548 dof_field_flux,
549 energy,
550 list_elements,
551 tree,
552 d_fe_d_dofdot_elem,
553 d_fe_d_dofdotdot_elem,
554 last_rate_dissipated_energy,
555 dissipated_energy,
556 global_last_rate_dissipated_energy,
557 global_dissipated_energy};
558
559 b2Mutex mutex;
560 b2runRecursively<ElementsTaskLoop<T, MATRIX_FORMAT, Args, true, INERTIA> >(
561 0, tree_nlevel, args, mutex);
562
563 if (energy) {
564 global_dissipated_energy +=
565 0.5 * delta_time * (global_last_rate_dissipated_energy + energy[1]);
566 global_last_rate_dissipated_energy = energy[1];
567 energy[2] = global_dissipated_energy;
568 }
569 }
570
571private:
572 struct Args {
573 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref> dof;
574 const double time;
575 const double delta_time;
576 const EquilibriumSolution equilibrium_solution;
577 const bool linear;
578 const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof;
579 const b2linalg::Vector<T, b2linalg::Vdense_constref>& d_r_d_time;
580 const b2linalg::Vector<double, b2linalg::Vdense_constref>& d_value_d_dof_coeff;
581 Model& model;
582 const size_t nb_dof;
583 b2linalg::Vector<T, b2linalg::Vdense_ref>& value;
584 const std::vector<bool>& d_value_d_dof_e_flags;
585 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_value_d_dof;
586 b2linalg::Vector<T, b2linalg::Vdense_ref>& d_value_d_time;
587 GradientContainer* gradient_container;
588 SolverHints* solver_hints;
589 logging::Logger* logger;
590 DofFieldFlux<T>* dof_field_flux;
591 double* energy;
592 std::vector<Element*>& list_elements;
593 std::vector<std::pair<size_t, size_t> >& tree;
594 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> >& d_fe_d_dofdot_elem;
595 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> >& d_fe_d_dofdotdot_elem;
596 std::vector<double>& last_rate_dissipated_energy;
597 std::vector<double>& dissipated_energy;
598 double& global_last_rate_dissipated_energy;
599 double& global_dissipated_energy;
600 };
601
602 void compute_tree(
603 Domain& domain, const b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_r_d_dof) {
604 element_iterator = domain.get_element_iterator();
605 const size_t nb_element = domain.get_number_of_elements();
606 list_elements.reserve(nb_element);
607#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
608 size_t i_ptr = 0;
609 std::vector<std::vector<size_t> > dof(domain.get_number_of_dof());
610 std::vector<size_t> g_ptr(nb_element + 1);
611 g_ptr[0] = 0;
612 std::vector<size_t> g_ind;
613 g_ind.reserve(nb_element);
614 b2linalg::Index index;
615#endif
616 Domain::ElementIterator* i = element_iterator.get();
617 if (i == 0) { TypeError() << "Bad solver type." << THROW; }
618 for (Element* element = i->next(); element != 0; element = i->next()) {
619 list_elements.push_back(element);
620#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
621 element->get_dof_numbering(index);
622 std::set<size_t> s;
623 std::set<size_t> ind(index.begin(), index.end());
624 trans_d_r_d_dof.this_prod_index(index);
625 ind.insert(index.begin(), index.end());
626 for (std::set<size_t>::const_iterator j = ind.begin(); j != ind.end(); ++j) {
627 std::vector<size_t>& d = dof[*j];
628 s.insert(d.begin(), d.end());
629 d.push_back(i_ptr);
630 }
631 g_ind.insert(g_ind.end(), s.begin(), s.end());
632 g_ptr[++i_ptr] = g_ind.size();
633#endif
634 }
635#ifdef SOLVER_UTILS_USE_RECURSIVE_TREE
636 dof.clear();
637 std::vector<size_t> inv_perm(nb_element);
638 int nlevel =
639 int(std::log(tbb::global_control::active_value(
640 tbb::global_control::max_allowed_parallelism))
641 / std::log(2));
642 tree_nlevel = nlevel;
643 tree.resize(int(std::pow(2.0, nlevel + 1)));
644 mlevel_nested_dissection_with_separator(
645 nb_element, &g_ptr[0], &g_ind[0], nlevel, &inv_perm[0], &tree[0]);
646 std::vector<Element*> list_elements_tmp(nb_element);
647 for (size_t i = 0; i != nb_element; ++i) {
648 list_elements_tmp[inv_perm[i]] = list_elements[i];
649 }
650 list_elements.swap(list_elements_tmp);
651#else
652 tree.push_back(std::pair<size_t, size_t>(0, list_elements.size()));
653 tree_nlevel = 0;
654#endif
655 }
656
657 std::unique_ptr<Domain::ElementIterator> element_iterator;
658 std::vector<Element*> list_elements;
659 std::vector<std::pair<size_t, size_t> > tree;
660 size_t tree_nlevel;
661
662 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_fe_d_dofdot_elem;
663 std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_fe_d_dofdotdot_elem;
664 std::vector<double> last_rate_dissipated_energy;
665 std::vector<double> dissipated_energy;
666 double global_last_rate_dissipated_energy;
667 double global_dissipated_energy;
668 double rayleigh_damping_alpha;
669 double rayleigh_damping_beta;
670};
671
672} // namespace solver
673
674} // namespace b2000
675#endif /*B2SOLVER_UTILS_H_*/
THROW
#define THROW
Definition b2exception.H:198
b2solver.H
Interface to C++ representations of FE solvers.
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::TypeError
GenericException< TypeError_name > TypeError
Definition b2exception.H:325