b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2dynamic_residue_function.H
1//------------------------------------------------------------------------
2// b2dynamic_rsidue_function.H --
3//
4// Collection of dynamic residue functions.
5//
6// written by Mathias Doreille
7// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
8// Harald Klimach <harald.klimach@dlr.de>
9//
10// (c) 2004-2012,2016 SMR Engineering & Development SA
11// 2502 Bienne, Switzerland
12//
13// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
14// Linder Höhe, 51147 Köln
15//
16// All Rights Reserved. Proprietary source code. The contents of
17// this file may not be disclosed to third parties, copied or
18// duplicated in any form, in whole or in part, without the prior
19// written permission of SMR.
20//------------------------------------------------------------------------
21
22#ifndef B2DYNAMIC_RESIDUE_FUNCTION_H_
23#define B2DYNAMIC_RESIDUE_FUNCTION_H_
24
25#include "b2ppconfig.h"
26#include "model/b2boundary_condition.H"
27#include "model/b2case.H"
28#include "model/b2domain.H"
29#include "model/b2model.H"
30#include "model/b2solution.H"
31#include "model/b2solver.H"
32#include "solvers/b2solver_utils.H"
33#include "solvers/b2static_nonlinear_utile.H"
34
35namespace b2000::solver {
36
37template <typename T, typename MATRIX_FORMAT>
38class DynamicResidueFunction;
39
40template <typename T, typename MATRIX_FORMAT>
41class DynamicResidueFunction : public StaticResidueFunctionForTerminationTest<T> {
42public:
43 DynamicResidueFunction()
44 : logger(logging::get_logger("solver.static_nonlinear.ResidueFunction")) {}
45
46 virtual void init(Model& model, const AttributeList& attribute) = 0;
47
48 virtual void set_case(Case* case_) = 0;
49
50 virtual int get_order() const = 0;
51
52 virtual b2linalg::SparseMatrixConnectivityType get_dof_connectivity_type() const = 0;
53
54 virtual void save_solution(
55 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof, const double time, const double dtime,
56 const RTable& attributes, bool end_of_stage = false,
57 bool unconverged_subcycle = false) = 0;
58
59 virtual void get_residue(
60 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
61 const double delta_time, const EquilibriumSolution equilibrium_solution,
62 b2linalg::Vector<T, b2linalg::Vdense>& residue,
63 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
64 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
65 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
66 CorrectionTerminationTest<T>* termination_test = 0) = 0;
67
68 virtual const b2linalg::Index& get_dof_field() = 0;
69
70 typedef ObjectTypeIncomplete<DynamicResidueFunction> type_t;
71 static type_t type;
72
73protected:
74 logging::Logger& logger;
75};
76
77template <typename T, typename MATRIX_FORMAT>
78typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t
79 DynamicResidueFunction<T, MATRIX_FORMAT>::type(
80 "DynamicResidueFunction", type_name<T, MATRIX_FORMAT>(),
81 StringList("DYNAMIC_RESIDUE_FUNCTION"), b2000::solver::module);
82
83template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
84class DynamicOrderNResidueFunction : public DynamicResidueFunction<T, MATRIX_FORMAT> {
85public:
86 typedef TypedNaturalBoundaryCondition<T, typename MATRIX_FORMAT::sparse_inversible> nbc_t;
87 typedef TypedEssentialBoundaryCondition<T> ebc_t;
88 typedef TypedModelReductionBoundaryCondition<T> mrbc_t;
89
90 DynamicOrderNResidueFunction()
91 : logger(logging::get_logger("solver.dynamic_nonlinear.residue")),
92 number_of_dof(0),
93 number_of_dof_red(0),
94 model(nullptr),
95 domain(nullptr),
96 mrbc(0),
97 ebc(0),
98 nbc(0),
99 solution(0),
100 constraint_has_penalty(false),
101 penalty_factor(0),
102 lagrange_mult_scale(1),
103 rcfo_only_spc(false) {}
104
105 void set_case(Case* case_) {}
106
107 void init(Model& model_, const AttributeList& attribute) {
108 model = &model_;
109 domain = &model->get_domain();
110 ebc = &dynamic_cast<ebc_t&>(model->get_essential_boundary_condition(ebc_t::type.get_subtype(
111 "TypedEssentialBoundaryConditionListComponent" + type_name<T>())));
112 mrbc = &dynamic_cast<mrbc_t&>(
113 model->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
114 "TypedModelReductionBoundaryConditionListComponent" + type_name<T>())));
115 nbc = &dynamic_cast<nbc_t&>(model->get_natural_boundary_condition(nbc_t::type.get_subtype(
116 "TypedNaturalBoundaryConditionListComponent"
117 + type_name<T, typename MATRIX_FORMAT::sparse_inversible>())));
118 solution = &dynamic_cast<SolutionDynamicNonlinear<T>&>(model->get_solution());
119 number_of_dof = domain->get_number_of_dof();
120 nbc_sol.resize(number_of_dof);
121 fi.resize(number_of_dof);
122
123 constraint_has_penalty =
124 attribute.get_string("CONSTRAINT_METHOD", "REDUCTION") == "PENALTY";
125
126 constraint_has_penalty |=
127 attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "PENALTY";
128
129 number_of_dof_red = mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1)
130 + (constraint_has_penalty ? 0 : ebc->get_size(false));
131
132 if (constraint_has_penalty) {
133 penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1e10);
134 } else {
135 if (attribute.get_bool("REMOVE_DEPENDENT_CONSTRAINTS", false)) {
136 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
137 ebc->get_nonlinear_value(
138 fi, 0, false, b2linalg::Vector<T, b2linalg::Vdense_ref>::null,
139 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
140 trans_d_l_d_dof.get_dep_columns(ebc_to_remove, 3e-13, 1.5);
141 number_of_dof_red -= ebc_to_remove.size();
142 if (!ebc_to_remove.empty()) {
143 DynamicResidueFunction<T, MATRIX_FORMAT>::logger
144 << logging::info << "Remove " << ebc_to_remove.size()
145 << " dependent constraints." << logging::LOGFLUSH;
146 }
147 }
148 penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1);
149 if (attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "LAGRANGE") {
150 penalty_factor = 0;
151 }
152 lagrange_mult_scale = attribute.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
153 }
154
155 rcfo_only_spc = attribute.get_bool("RCFO_ONLY_SPC", false);
156 if (rcfo_only_spc
157 && attribute.get_string("CONSTRAINT_METHOD", "REDUCTION") == "REDUCTION") {
158 Exception() << "It is not possible to compute the reaction "
159 "forces with rcfo_only_spc=yes and "
160 "constraint_method=reduction (which is the default). "
161 "Instead, set constraint_method either to lagrange, "
162 "augmented_lagrange, or penalty."
163 << THROW;
164 }
165 }
166
167 void mrbc_get_nonlinear_value(
168 const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, double time,
169 EquilibriumSolution equilibrium_solution, b2linalg::Vector<T, b2linalg::Vdense_ref> value,
170 b2linalg::Matrix<T, b2linalg::Mcompressed_col>& d_value_d_dof_trans,
171 b2linalg::Vector<T, b2linalg::Vdense_ref> d_value_d_time, SolverHints* solver_hints) {
172 mrbc->get_nonlinear_value(
173 dof, time, equilibrium_solution, value, d_value_d_dof_trans,
174 b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(d_value_d_time), solver_hints);
175 }
176
177 size_t get_number_of_dof() {
178 return mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1)
179 + (constraint_has_penalty ? 0 : ebc->get_size(false) - ebc_to_remove.size());
180 }
181
182 size_t get_number_of_non_lagrange_dof() {
183 return mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
184 }
185
186 b2linalg::SparseMatrixConnectivityType get_dof_connectivity_type() const {
187 return domain->get_dof_connectivity_type();
188 }
189
190 int get_order() const { return domain->get_value_info().size() - 1; }
191
192 void save_solution(
193 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
194 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
195 const size_t mrbc_size =
196 mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
197 b2linalg::Interval disp_range(0, mrbc_size);
198 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(number_of_dof, dof_red.size2());
199 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
200 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(number_of_dof, dof_red.size2() - 1);
201 mrbc->get_nonlinear_value(
202 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
203 for (size_t j = 1; j < dof.size2(); ++j) {
204 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
205 dof[j] += d_r_d_time[j - 1];
206 }
207 if (rcfo_only_spc) {
208 // recompute fi using only the spc contribution
209
210 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
211 b2linalg::Vector<T> l(ebc->get_size(false));
212 ebc->get_nonlinear_value(
213 dof[0], time, false, l, trans_d_l_d_dof,
214 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
215 if (!ebc_to_remove.empty()) {
216 l.remove(ebc_to_remove);
217 trans_d_l_d_dof.remove_col(ebc_to_remove);
218 }
219
220 if (constraint_has_penalty) {
221 l *= -penalty_factor;
222 } else {
223 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
224 if (MATRIX_FORMAT::symmetric && penalty_factor != 0) {
225 // We apply penalty method with a small factor to obtain a definite positive
226 // matrix.
227 l *= -penalty_factor;
228 l += -lagrange_mult_scale * dof_red[0](lagrange_mult_range);
229 } else {
230 l = -lagrange_mult_scale * dof_red[0](lagrange_mult_range);
231 }
232 }
233 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
234 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
235 }
236 fi = trans_d_l_d_dof * l;
237 }
238 if (unconverged_subcycle) {
239 solution->set_solution(dof, time, nbc_sol, fi, unconverged_subcycle);
240 return;
241 }
242 solution->new_cycle(end_of_stage);
243 solution->set_solution(dof, time, nbc_sol, fi);
244
245 typename SolutionDynamicNonlinear<T>::gradient_container gc =
246 solution->get_gradient_container();
247 domain->set_dof(dof);
248 EquilibriumSolution es(1, 0);
249 solver_utile.get_elements_values(
250 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
251 b2linalg::Vector<double, b2linalg::Vdense>::null, *model,
252 b2linalg::Vector<T, b2linalg::Vdense>::null,
253 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
254 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, logger, 0);
255 std::string domain_state_id = solution->get_domain_state_id();
256 if (!domain_state_id.empty()) { domain->save_state(domain_state_id); }
257 solution->terminate_cycle(time, dtime, attributes);
258 if (end_of_stage) { solution->terminate_stage(true); }
259 }
260
261 void get_residue(
262 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
263 const double delta_time, const EquilibriumSolution equilibrium_solution,
264 b2linalg::Vector<T, b2linalg::Vdense>& residue,
265 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
266 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
267 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
268 CorrectionTerminationTest<T>* termination_test = 0) {
269 b2linalg::Interval disp_range(
270 0, mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
271 b2linalg::Interval lagrange_mult_range(
272 mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
273 number_of_dof_red);
274 if (!residue.is_null()) { residue.set_zero(); }
275 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
276
277 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
278 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(number_of_dof, dof_red.size2() - 1);
279 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(number_of_dof, dof_red.size2());
280 mrbc->get_nonlinear_value(
281 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
282 d_r_d_time, solver_hints);
283 for (size_t j = 1; j < dof.size2(); ++j) {
284 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
285 dof[j] += d_r_d_time[j - 1];
286 }
287 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
288 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(number_of_dof);
289 nbc_sol.set_zero();
290 nbc->get_nonlinear_value(
291 dof[0], time, equilibrium_solution,
292 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
293 : b2linalg::Vector<T, b2linalg::Vdense_ref>(nbc_sol),
294 d_residue_d_dofc.is_null()
295 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
296 : d_fe_d_dof,
297 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
298 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
299 solver_hints);
300 if (termination_test && !residue.is_null()) { termination_test->add_flux(nbc_sol); }
301 fi = -nbc_sol;
302 d_fe_d_time = -d_fe_d_time;
303 if (d_fe_d_dof.size1() != 0) {
304 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
305 }
306
307 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
308 b2linalg::Vector<T> l;
309
310 if (constraint_has_penalty) {
311 if (!residue.is_null()) { l.resize(ebc->get_size(false)); }
312 ebc->get_nonlinear_value(
313 dof[0], time, equilibrium_solution,
314 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
315 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
316 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
317 } else {
318 if (ebc_to_remove.empty()) {
319 ebc->get_nonlinear_value(
320 dof[0], time, equilibrium_solution,
321 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
322 : residue(lagrange_mult_range)),
323 trans_d_l_d_dof,
324 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
325 : d_residue_d_time(lagrange_mult_range)),
326 solver_hints);
327 } else {
328 b2linalg::Vector<T> l(residue.is_null() ? 0 : ebc->get_size(false));
329 ebc->get_nonlinear_value(
330 dof[0], time, equilibrium_solution,
331 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
332 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
333 trans_d_l_d_dof,
334 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
335 : d_residue_d_time(lagrange_mult_range)),
336 solver_hints);
337 if (!residue.is_null()) {
338 l.remove(ebc_to_remove);
339 residue(lagrange_mult_range) = l;
340 }
341 trans_d_l_d_dof.remove_col(ebc_to_remove);
342 }
343 if (!d_residue_d_time.is_null()) {
344 d_residue_d_time(lagrange_mult_range) +=
345 transposed(trans_d_l_d_dof) * d_r_d_time[0];
346 }
347 }
348
349 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
350 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
351 }
352
353 solver_utile.get_elements_values(
354 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
355 d_residue_d_dof_coeff, *model, fi, d_residue_d_dofc,
356 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
357 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
358 0, solver_hints, logger, termination_test);
359
360 if (!residue.is_null()) {
361 b2linalg::Vector<T> tmp = fi;
362 if (constraint_has_penalty) {
363 tmp += (trans_d_l_d_dof * l) * penalty_factor;
364 } else {
365 if (MATRIX_FORMAT::symmetric && penalty_factor != 0) {
366 // We apply penalty method with a small factor to obtain a definite positive
367 // matrix.
368 b2linalg::Vector<T> tmp1;
369 tmp1 = residue(lagrange_mult_range) * penalty_factor;
370 tmp1 += dof_red[0](lagrange_mult_range) * lagrange_mult_scale;
371 tmp += trans_d_l_d_dof * tmp1;
372 } else {
373 tmp +=
374 lagrange_mult_scale * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
375 }
376 residue(lagrange_mult_range) *= lagrange_mult_scale;
377 }
378 residue(disp_range) = trans_d_r_d_dof * tmp;
379 }
380 if (!d_residue_d_dofc.is_null()) {
381 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
382 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
383 if (penalty_factor != 0 && (constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
384 // We apply penalty method but with a small factor to
385 // obtain a definite positive matrix when Lagrange multipliers method is used.
386 d_residue_d_dofc(disp_range, disp_range) +=
387 (trans_LR * transposed(trans_LR)) * penalty_factor;
388 }
389 if (!constraint_has_penalty) {
390 trans_LR *= lagrange_mult_scale;
391 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
392 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
393 if (!MATRIX_FORMAT::symmetric) {
394 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
395 }
396 }
397 }
398 if (!d_residue_d_time.is_null()) {
399 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
400 }
401 }
402
403 const b2linalg::Index& get_dof_field() {
404 if (dof_field.empty()
405 && mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1) != 0) {
406 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
407 mrbc->get_nonlinear_value(
408 b2linalg::Vector<double, b2linalg::Vdense_constref>::null, 0, false,
409 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, trans_d_r_d_dof,
410 b2linalg::Matrix<T, b2linalg::Mrectangle_ref>::null, 0);
411 dof_field.resize(trans_d_r_d_dof.size1());
412 const b2linalg::Index& df = model->get_domain().get_dof_field().second;
413 size_t* colind;
414 size_t* rowind;
415 T* value;
416 const size_t s = trans_d_r_d_dof.size2();
417 trans_d_r_d_dof.get_values(colind, rowind, value);
418 const size_t* rowind_ptr = rowind;
419 for (size_t j = 0; j != s; ++j) {
420 const size_t* rowind_end = rowind + colind[j + 1];
421 for (; rowind_ptr != rowind_end; ++rowind_ptr) {
422 dof_field[*rowind_ptr] = std::max(dof_field[*rowind_ptr], df[j]);
423 }
424 }
425 }
426 return dof_field;
427 }
428
429 typedef ObjectTypeComplete<
430 DynamicOrderNResidueFunction, typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
431 type_t;
432 static type_t type;
433
434protected:
435 logging::Logger& logger;
436 size_t number_of_dof;
437 size_t number_of_dof_red;
438 Model* model;
439 Domain* domain;
440 mrbc_t* mrbc;
441 ebc_t* ebc;
442 nbc_t* nbc;
443 SolutionDynamicNonlinear<T>* solution;
444 b2linalg::Vector<T, b2linalg::Vdense> fi;
445 b2linalg::Vector<T, b2linalg::Vdense> nbc_sol;
446 b2linalg::Index ebc_to_remove;
447 SolverUtils<T, MATRIX_FORMAT> solver_utile;
448 bool constraint_has_penalty;
449 double penalty_factor;
450 double lagrange_mult_scale;
451 bool rcfo_only_spc;
452 b2linalg::Index dof_field;
453};
454
455template <typename T, typename MATRIX_FORMAT>
456typename DynamicOrderNResidueFunction<T, MATRIX_FORMAT>::type_t
457 DynamicOrderNResidueFunction<T, MATRIX_FORMAT>::type(
458 "DynamicOrderNResidueFunction", type_name<T, MATRIX_FORMAT>(),
459 StringList(
460 "ORDER_ONE_DYNAMIC_RESIDUE_FUNCTION", "ORDER_TWO_DYNAMIC_RESIDUE_FUNCTION",
461 "ORDER_N_DYNAMIC_RESIDUE_FUNCTION"),
462 b2000::solver::module, &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
463
464template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
465class DynamicOrderTwoRayleighDampingResidueFunction
466 : public DynamicOrderNResidueFunction<T, MATRIX_FORMAT> {
467public:
468 typedef DynamicOrderNResidueFunction<T, MATRIX_FORMAT> parent;
469
470 void init(Model& model, const AttributeList& attribute) {
471 parent::init(model, attribute);
472 if (parent::get_order() != 2) {
473 Exception() << "Rayleigh damping can only be used for models "
474 "of order 2"
475 << THROW;
476 }
477 solver_utile.init(
478 rayleigh_damping_alpha, rayleigh_damping_beta, *parent::model, parent::logger);
479 }
480
481 void set_case(Case* case_) {
482 rayleigh_damping_alpha = case_->get_double("RAYLEIGH_ALPHA");
483 rayleigh_damping_beta = case_->get_double("RAYLEIGH_BETA");
484 }
485
486 void save_solution(
487 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
488 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
489 const size_t mrbc_size =
490 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
491 b2linalg::Interval disp_range(0, mrbc_size);
492 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
493 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
494 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
495 parent::number_of_dof, dof_red.size2() - 1);
496 parent::mrbc->get_nonlinear_value(
497 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
498 for (size_t j = 1; j < dof.size2(); ++j) {
499 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
500 dof[j] += d_r_d_time[j - 1];
501 }
502 if (parent::rcfo_only_spc) {
503 // recompute fi using only the spc contribution
504
505 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
506 b2linalg::Vector<T> l(parent::ebc->get_size(false));
507 parent::ebc->get_nonlinear_value(
508 dof[0], time, false, l, trans_d_l_d_dof,
509 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
510 if (!parent::ebc_to_remove.empty()) {
511 l.remove(parent::ebc_to_remove);
512 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
513 }
514
515 if (parent::constraint_has_penalty) {
516 l *= -parent::penalty_factor;
517 } else {
518 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
519 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
520 // We apply penalty method with a small factor to obtain a definite positive
521 // matrix.
522 l *= -parent::penalty_factor;
523 l += -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
524 } else {
525 l = -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
526 }
527 }
528 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
529 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
530 }
531 parent::fi = trans_d_l_d_dof * l;
532 }
533 if (unconverged_subcycle) {
534 parent::solution->set_solution(
535 dof, time, parent::nbc_sol, parent::fi, unconverged_subcycle);
536 return;
537 }
538 parent::solution->new_cycle(end_of_stage);
539 parent::solution->set_solution(dof, time, parent::nbc_sol, parent::fi);
540
541 typename SolutionDynamicNonlinear<T>::gradient_container gc =
542 parent::solution->get_gradient_container();
543 parent::domain->set_dof(dof);
544 EquilibriumSolution es(1, 0);
545 double energy[3];
546 solver_utile.get_elements_values(
547 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
548 b2linalg::Vector<double, b2linalg::Vdense>::null, *parent::model,
549 b2linalg::Vector<T, b2linalg::Vdense>::null,
550 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
551 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, parent::logger, 0, energy);
552 parent::solution->set_value("RAYLEIGH_KINEMATIC_ENERGY", energy[0]);
553 parent::solution->set_value("RAYLEIGH_RATE_DAMPING_DISSIPATED_ENERGY", energy[1]);
554 parent::solution->set_value("RAYLEIGH_DAMPING_DISSIPATED_ENERGY", energy[2]);
555 std::string domain_state_id = parent::solution->get_domain_state_id();
556 if (!domain_state_id.empty()) { parent::domain->save_state(domain_state_id); }
557 parent::solution->terminate_cycle(time, dtime, attributes);
558 if (end_of_stage) { parent::solution->terminate_stage(true); }
559 }
560
561 void get_residue(
562 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
563 const double delta_time, const EquilibriumSolution equilibrium_solution,
564 b2linalg::Vector<T, b2linalg::Vdense>& residue,
565 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
566 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
567 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
568 CorrectionTerminationTest<T>* termination_test = 0) {
569 // logging::get_logger("solver.static_nonlinear.residue_function")
570 b2linalg::Interval disp_range(
571 0, parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
572 b2linalg::Interval lagrange_mult_range(
573 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
574 parent::number_of_dof_red);
575
576 if (!residue.is_null()) { residue.set_zero(); }
577 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
578
579 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
580 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
581 parent::domain->get_number_of_dof(), dof_red.size2() - 1);
582 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
583 parent::mrbc->get_nonlinear_value(
584 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
585 d_r_d_time, solver_hints);
586 for (size_t j = 1; j < dof.size2(); ++j) {
587 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
588 dof[j] += d_r_d_time[j - 1];
589 }
590 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
591 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(parent::domain->get_number_of_dof());
592 parent::nbc_sol.set_zero();
593 parent::nbc->get_nonlinear_value(
594 dof[0], time, equilibrium_solution,
595 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
596 : b2linalg::Vector<T, b2linalg::Vdense_ref>(parent::nbc_sol),
597 d_residue_d_dofc.is_null()
598 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
599 : d_fe_d_dof,
600 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
601 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
602 solver_hints);
603 if (termination_test && !residue.is_null()) { termination_test->add_flux(parent::nbc_sol); }
604 parent::fi = -parent::nbc_sol;
605 d_fe_d_time = -d_fe_d_time;
606 if (d_fe_d_dof.size1() != 0) {
607 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
608 }
609
610 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
611 b2linalg::Vector<T> l;
612
613 if (parent::constraint_has_penalty) {
614 if (!residue.is_null()) { l.resize(parent::ebc->get_size(false)); }
615 parent::ebc->get_nonlinear_value(
616 dof[0], time, equilibrium_solution,
617 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
618 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
619 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
620 } else {
621 if (parent::ebc_to_remove.empty()) {
622 parent::ebc->get_nonlinear_value(
623 dof[0], time, equilibrium_solution,
624 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
625 : residue(lagrange_mult_range),
626 trans_d_l_d_dof,
627 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
628 : d_residue_d_time(lagrange_mult_range),
629 solver_hints);
630 } else {
631 b2linalg::Vector<T> l(residue.is_null() ? 0 : parent::ebc->get_size(false));
632 parent::ebc->get_nonlinear_value(
633 dof[0], time, equilibrium_solution,
634 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
635 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
636 trans_d_l_d_dof,
637 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
638 : d_residue_d_time(lagrange_mult_range)),
639 solver_hints);
640 if (!residue.is_null()) {
641 l.remove(parent::ebc_to_remove);
642 residue(lagrange_mult_range) = l;
643 }
644 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
645 }
646 if (!d_residue_d_time.is_null()) {
647 d_residue_d_time(lagrange_mult_range) +=
648 transposed(trans_d_l_d_dof) * d_r_d_time[0];
649 }
650 }
651
652 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
653 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
654 }
655
656 solver_utile.get_elements_values(
657 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
658 d_residue_d_dof_coeff, *parent::model, parent::fi, d_residue_d_dofc,
659 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
660 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
661 0, solver_hints, parent::logger, termination_test);
662
663 if (!residue.is_null()) {
664 b2linalg::Vector<T> tmp = parent::fi;
665 if (parent::constraint_has_penalty) {
666 tmp += (trans_d_l_d_dof * l) * parent::penalty_factor;
667 } else {
668 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
669 // We apply penalty method with a small factor to obtain a definite positive
670 // matrix.
671 b2linalg::Vector<T> tmp1;
672 tmp1 = residue(lagrange_mult_range) * parent::penalty_factor;
673 tmp1 += dof_red[0](lagrange_mult_range) * parent::lagrange_mult_scale;
674 tmp += trans_d_l_d_dof * tmp1;
675 } else {
676 tmp += parent::lagrange_mult_scale
677 * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
678 }
679 residue(lagrange_mult_range) *= parent::lagrange_mult_scale;
680 }
681 residue(disp_range) = trans_d_r_d_dof * tmp;
682 }
683 if (!d_residue_d_dofc.is_null()) {
684 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
685 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
686 if (parent::penalty_factor != 0
687 && (parent::constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
688 // We apply penalty method but with a small factor to
689 // obtain a definite positive matrix when Lagrange multipliers method is used.
690 d_residue_d_dofc(disp_range, disp_range) +=
691 (trans_LR * transposed(trans_LR)) * parent::penalty_factor;
692 }
693 if (!parent::constraint_has_penalty) {
694 trans_LR *= parent::lagrange_mult_scale;
695 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
696 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
697 if (!MATRIX_FORMAT::symmetric) {
698 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
699 }
700 }
701 }
702 if (!d_residue_d_time.is_null()) {
703 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
704 }
705 }
706
707 typedef ObjectTypeComplete<
708 DynamicOrderTwoRayleighDampingResidueFunction,
709 typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
710 type_t;
711 static type_t type;
712
713private:
714 double rayleigh_damping_alpha;
715 double rayleigh_damping_beta;
716 SolverUtilsOrderTwoRayleighDamping<T, MATRIX_FORMAT> solver_utile;
717};
718
719template <typename T, typename MATRIX_FORMAT>
720typename DynamicOrderTwoRayleighDampingResidueFunction<T, MATRIX_FORMAT>::type_t
721 DynamicOrderTwoRayleighDampingResidueFunction<T, MATRIX_FORMAT>::type(
722 "DynamicOrderTwoRayleighDampingResidueFunction", type_name<T, MATRIX_FORMAT>(),
723 StringList("RAYLEIGH_DAMPING_DYNAMIC_RESIDUE_FUNCTION"), b2000::solver::module,
724 &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
725
726template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
727class QuasistaticViscousDampingResidueFunction
728 : public DynamicOrderNResidueFunction<T, MATRIX_FORMAT> {
729public:
730 typedef DynamicOrderNResidueFunction<T, MATRIX_FORMAT> parent;
731
732 void init(Model& model, const AttributeList& attribute) {
733 parent::init(model, attribute);
734 if (parent::get_order() != 2) {
735 Exception() << "Quasi-static viscous damping can only be "
736 "used for models of order 2"
737 << THROW;
738 }
739 solver_utile.init(
740 rayleigh_damping_alpha, rayleigh_damping_beta, *parent::model, parent::logger);
741 }
742
743 void set_case(Case* case_) {
744 const char* alpha_s = "RAYLEIGH_ALPHA";
745 const char* beta_s = "RAYLEIGH_BETA";
746 if (case_->has_key(alpha_s) && case_->has_key(beta_s)) {
747 rayleigh_damping_alpha = case_->get_double("RAYLEIGH_ALPHA");
748 rayleigh_damping_beta = case_->get_double("RAYLEIGH_BETA");
749 dissipated_energy_fraction = 0;
750 } else {
751 dissipated_energy_fraction = case_->get_double("DISSIPATED_ENERGY_FRACTION", 1.e-4);
752 rayleigh_damping_alpha = 1e-10; // we set a very small dissipation for the first step
753 rayleigh_damping_beta = 0;
754 }
755 }
756
757 int get_order() const { return 1; }
758
759 void save_solution(
760 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
761 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
762 const size_t mrbc_size =
763 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
764 b2linalg::Interval disp_range(0, mrbc_size);
765 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
766 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
767 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
768 parent::number_of_dof, dof_red.size2() - 1);
769 parent::mrbc->get_nonlinear_value(
770 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
771 for (size_t j = 1; j < dof.size2(); ++j) {
772 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
773 dof[j] += d_r_d_time[j - 1];
774 }
775 if (parent::rcfo_only_spc) {
776 // recompute fi using only the spc contribution
777
778 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
779 b2linalg::Vector<T> l(parent::ebc->get_size(false));
780 parent::ebc->get_nonlinear_value(
781 dof[0], time, false, l, trans_d_l_d_dof,
782 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
783 if (!parent::ebc_to_remove.empty()) {
784 l.remove(parent::ebc_to_remove);
785 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
786 }
787
788 if (parent::constraint_has_penalty) {
789 l *= -parent::penalty_factor;
790 } else {
791 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
792 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
793 // We apply penalty method with a small factor to obtain a definite positive
794 // matrix.
795 l *= -parent::penalty_factor;
796 l += -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
797 } else {
798 l = -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
799 }
800 }
801 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
802 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
803 }
804 parent::fi = trans_d_l_d_dof * l;
805 }
806 if (unconverged_subcycle) {
807 parent::solution->set_solution(
808 dof, time, parent::nbc_sol, parent::fi, unconverged_subcycle);
809 return;
810 }
811 parent::solution->new_cycle(end_of_stage);
812 parent::solution->set_solution(dof, time, parent::nbc_sol, parent::fi);
813
814 typename SolutionDynamicNonlinear<T>::gradient_container gc =
815 parent::solution->get_gradient_container();
816 parent::domain->set_dof(dof);
817 EquilibriumSolution es(1, 0);
818 double energy[3];
819 solver_utile.get_elements_values(
820 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
821 b2linalg::Vector<double, b2linalg::Vdense>::null, *parent::model,
822 b2linalg::Vector<T, b2linalg::Vdense>::null,
823 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
824 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, parent::logger, 0, energy);
825 parent::solution->set_value("RAYLEIGH_KINEMATIC_ENERGY", energy[0]);
826 parent::solution->set_value("RAYLEIGH_RATE_DAMPING_DISSIPATED_ENERGY", energy[1]);
827 parent::solution->set_value("RAYLEIGH_DAMPING_DISSIPATED_ENERGY", energy[2]);
828 if (dissipated_energy_fraction != 0) {
829 double res = 0;
830 for (size_t i = 0; i != dof.size1(); ++i) {
831 // res += std::max(norm(dof(i, 0) * parent::nbc_sol[i]), norm(dof(i, 0) *
832 // parent::fi[i]));
833 res += dof(i, 0) * (parent::nbc_sol[i] + parent::fi[i]);
834 }
835 double sfactor = dissipated_energy_fraction * res / (energy[1] * dtime);
836 parent::logger << logging::debug << "scale factor for the stabilised method "
837 << sfactor * 1e-10 << logging::LOGFLUSH;
838 solver_utile.scale_raylay_damping(sfactor);
839 dissipated_energy_fraction = 0;
840 }
841 std::string domain_state_id = parent::solution->get_domain_state_id();
842 if (!domain_state_id.empty()) { parent::domain->save_state(domain_state_id); }
843 parent::solution->terminate_cycle(time, dtime, attributes);
844 if (end_of_stage) { parent::solution->terminate_stage(true); }
845 }
846
847 void get_residue(
848 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
849 const double delta_time, const EquilibriumSolution equilibrium_solution,
850 b2linalg::Vector<T, b2linalg::Vdense>& residue,
851 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
852 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
853 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
854 CorrectionTerminationTest<T>* termination_test = 0) {
855 // logging::get_logger("solver.static_nonlinear.residue_function")
856 b2linalg::Interval disp_range(
857 0, parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
858 b2linalg::Interval lagrange_mult_range(
859 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
860 parent::number_of_dof_red);
861
862 if (!residue.is_null()) { residue.set_zero(); }
863 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
864
865 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
866 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
867 parent::domain->get_number_of_dof(), dof_red.size2() - 1);
868 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
869 parent::mrbc->get_nonlinear_value(
870 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
871 d_r_d_time, solver_hints);
872 for (size_t j = 1; j < dof.size2(); ++j) {
873 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
874 dof[j] += d_r_d_time[j - 1];
875 }
876 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
877 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(parent::domain->get_number_of_dof());
878 parent::nbc_sol.set_zero();
879 parent::nbc->get_nonlinear_value(
880 dof[0], time, equilibrium_solution,
881 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
882 : b2linalg::Vector<T, b2linalg::Vdense_ref>(parent::nbc_sol),
883 d_residue_d_dofc.is_null()
884 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
885 : d_fe_d_dof,
886 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
887 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
888 solver_hints);
889 if (termination_test && !residue.is_null()) { termination_test->add_flux(parent::nbc_sol); }
890 parent::fi = -parent::nbc_sol;
891 d_fe_d_time = -d_fe_d_time;
892 if (d_fe_d_dof.size1() != 0) {
893 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
894 }
895
896 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
897 b2linalg::Vector<T> l;
898
899 if (parent::constraint_has_penalty) {
900 if (!residue.is_null()) { l.resize(parent::ebc->get_size(false)); }
901 parent::ebc->get_nonlinear_value(
902 dof[0], time, equilibrium_solution,
903 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
904 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
905 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
906 } else {
907 if (parent::ebc_to_remove.empty()) {
908 parent::ebc->get_nonlinear_value(
909 dof[0], time, equilibrium_solution,
910 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
911 : residue(lagrange_mult_range),
912 trans_d_l_d_dof,
913 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
914 : d_residue_d_time(lagrange_mult_range),
915 solver_hints);
916 } else {
917 b2linalg::Vector<T> l(residue.is_null() ? 0 : parent::ebc->get_size(false));
918 parent::ebc->get_nonlinear_value(
919 dof[0], time, equilibrium_solution,
920 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
921 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
922 trans_d_l_d_dof,
923 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
924 : d_residue_d_time(lagrange_mult_range)),
925 solver_hints);
926 if (!residue.is_null()) {
927 l.remove(parent::ebc_to_remove);
928 residue(lagrange_mult_range) = l;
929 }
930 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
931 }
932 if (!d_residue_d_time.is_null()) {
933 d_residue_d_time(lagrange_mult_range) +=
934 transposed(trans_d_l_d_dof) * d_r_d_time[0];
935 }
936 }
937
938 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
939 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
940 }
941
942 solver_utile.get_elements_values(
943 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
944 d_residue_d_dof_coeff, *parent::model, parent::fi, d_residue_d_dofc,
945 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
946 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
947 0, solver_hints, parent::logger, termination_test);
948
949 if (!residue.is_null()) {
950 b2linalg::Vector<T> tmp = parent::fi;
951 if (parent::constraint_has_penalty) {
952 tmp += (trans_d_l_d_dof * l) * parent::penalty_factor;
953 } else {
954 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
955 // We apply penalty method with a small factor to obtain a definite positive
956 // matrix.
957 b2linalg::Vector<T> tmp1;
958 tmp1 = residue(lagrange_mult_range) * parent::penalty_factor;
959 tmp1 += dof_red[0](lagrange_mult_range) * parent::lagrange_mult_scale;
960 tmp += trans_d_l_d_dof * tmp1;
961 } else {
962 tmp += parent::lagrange_mult_scale
963 * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
964 }
965 residue(lagrange_mult_range) *= parent::lagrange_mult_scale;
966 }
967 residue(disp_range) = trans_d_r_d_dof * tmp;
968 }
969 if (!d_residue_d_dofc.is_null()) {
970 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
971 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
972 if (parent::penalty_factor != 0
973 && (parent::constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
974 // We apply penalty method but with a small factor to
975 // obtain a definite positive matrix when Lagrange multipliers method is used.
976 d_residue_d_dofc(disp_range, disp_range) +=
977 (trans_LR * transposed(trans_LR)) * parent::penalty_factor;
978 }
979 if (!parent::constraint_has_penalty) {
980 trans_LR *= parent::lagrange_mult_scale;
981 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
982 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
983 if (!MATRIX_FORMAT::symmetric) {
984 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
985 }
986 }
987 }
988 if (!d_residue_d_time.is_null()) {
989 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
990 }
991 }
992
993 typedef ObjectTypeComplete<
994 QuasistaticViscousDampingResidueFunction,
995 typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
996 type_t;
997 static type_t type;
998
999private:
1000 double rayleigh_damping_alpha;
1001 double rayleigh_damping_beta;
1002 double dissipated_energy_fraction;
1003 SolverUtilsOrderTwoRayleighDamping<T, MATRIX_FORMAT, false> solver_utile;
1004};
1005
1006template <typename T, typename MATRIX_FORMAT>
1007typename QuasistaticViscousDampingResidueFunction<T, MATRIX_FORMAT>::type_t
1008 QuasistaticViscousDampingResidueFunction<T, MATRIX_FORMAT>::type(
1009 "QuasistaticViscousDampingResidueFunction", type_name<T, MATRIX_FORMAT>(),
1010 StringList(
1011 "QUASISTATIC_VISCOUS_DAMPING_DYNAMIC_RESIDUE_FUNCTION",
1012 "ARTIFICIAL_DAMPING_DYNAMIC_RESIDUE_FUNCTION"),
1013 b2000::solver::module, &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
1014
1015} // namespace b2000::solver
1016
1017#endif // B2DYNAMIC_RESIDUE_FUNCTION_H_
b2boundary_condition.H
THROW
#define THROW
Definition b2exception.H:198
b2solver.H
Interface to C++ representations of FE solvers.
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829