Flow123d  master-f44eb46
elasticity.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file transport_dg.cc
15  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
16  * @author Jan Stebel
17  */
18 
19 #include "system/sys_profiler.hh"
20 #include "mechanics/elasticity.hh"
22 
23 #include "io/output_time.hh"
25 #include "fem/fe_values.hh"
26 #include "fem/fe_p.hh"
27 #include "fem/fe_rt.hh"
28 #include "fem/fe_system.hh"
29 #include "fields/field_fe.hh"
30 #include "la/linsys_PETSC.hh"
31 #include "la/linsys_PERMON.hh"
32 #include "coupling/balance.hh"
33 #include "mesh/neighbours.h"
35 
36 #include "fields/multi_field.hh"
37 #include "fields/generic_field.hh"
38 #include "fields/field_model.hh"
39 #include "input/factory.hh"
40 
41 
42 
43 
44 using namespace Input::Type;
45 
46 
47 
49  std::string equation_name = std::string(name_) + "_FE";
50  return IT::Record(
51  std::string(equation_name),
52  "FEM for linear elasticity.")
55  .declare_key("balance", Balance::get_input_type(), Default("{}"),
56  "Settings for computing balance.")
58  "Parameters of output stream.")
60  "Linear solver for elasticity.")
61  .declare_key("input_fields", Array(
63  .make_field_descriptor_type(equation_name)),
65  "Input fields of the equation.")
66  .declare_key("output",
67  EqFields().output_fields.make_output_type(equation_name, ""),
68  IT::Default("{ \"fields\": [ \"displacement\" ] }"),
69  "Setting of the field output.")
70  .declare_key("contact", Bool(), IT::Default("false"), "Indicates the use of contact conditions on fractures.")
71  .close();
72 }
73 
74 const int Elasticity::registrar =
75  Input::register_class< Elasticity, Mesh &, const Input::Record>(std::string(name_) + "_FE") +
77 
78 
79 
80 double lame_mu(double young, double poisson)
81 {
82  return young*0.5/(poisson+1.);
83 }
84 
85 
86 double lame_lambda(double young, double poisson)
87 {
88  return young*poisson/((poisson+1.)*(1.-2.*poisson));
89 }
90 
91 // Functor computing lame_mu
92 struct fn_lame_mu {
93  inline double operator() (double young, double poisson) {
94  return young * 0.5 / (poisson+1.);
95  }
96 };
97 
98 // Functor computing lame_lambda
100  inline double operator() (double young, double poisson) {
101  return young * poisson / ((poisson+1.)*(1.-2.*poisson));
102  }
103 };
104 
105 // Functor computing base of dirichlet_penalty (without dividing by side meassure)
107  inline double operator() (double lame_mu, double lame_lambda) {
108  return 1e3 * (2 * lame_mu + lame_lambda);
109  }
110 };
111 
112 
113 
114 
115 
116 
118  return Selection("Elasticity_BC_Type", "Types of boundary conditions for mechanics.")
119  .add_value(bc_type_displacement, "displacement",
120  "Prescribed displacement.")
121  .add_value(bc_type_displacement_normal, "displacement_n",
122  "Prescribed displacement in the normal direction to the boundary.")
123  .add_value(bc_type_traction, "traction",
124  "Prescribed traction.")
125  .add_value(bc_type_stress, "stress",
126  "Prescribed stress tensor.")
127  .close();
128 }
129 
130 
132 {
133  *this+=bc_type
134  .name("bc_type")
135  .description(
136  "Type of boundary condition.")
137  .units( UnitSI::dimensionless() )
138  .input_default("\"traction\"")
139  .input_selection( get_bc_type_selection() )
141 
142  *this+=bc_displacement
143  .name("bc_displacement")
144  .description("Prescribed displacement on boundary.")
145  .units( UnitSI().m() )
146  .input_default("0.0")
147  .flags_add(in_rhs);
148 
149  *this+=bc_traction
150  .name("bc_traction")
151  .description("Prescribed traction on boundary.")
152  .units( UnitSI().Pa() )
153  .input_default("0.0")
154  .flags_add(in_rhs);
155 
156  *this+=bc_stress
157  .name("bc_stress")
158  .description("Prescribed stress on boundary.")
159  .units( UnitSI().Pa() )
160  .input_default("0.0")
161  .flags_add(in_rhs);
162 
163  *this+=load
164  .name("load")
165  .description("Prescribed bulk load.")
166  .units( UnitSI().N().m(-3) )
167  .input_default("0.0")
168  .flags_add(in_rhs);
169 
170  *this+=young_modulus
171  .name("young_modulus")
172  .description("Young's modulus.")
173  .units( UnitSI().Pa() )
174  .input_default("0.0")
175  .flags_add(in_main_matrix & in_rhs);
176 
177  *this+=poisson_ratio
178  .name("poisson_ratio")
179  .description("Poisson's ratio.")
180  .units( UnitSI().dimensionless() )
181  .input_default("0.0")
182  .flags_add(in_main_matrix & in_rhs);
183 
184  *this+=fracture_sigma
185  .name("fracture_sigma")
186  .description(
187  "Coefficient of transfer of forces through fractures.")
188  .units( UnitSI::dimensionless() )
189  .input_default("1.0")
190  .flags_add(in_main_matrix & in_rhs);
191 
192  *this+=initial_stress
193  .name("initial_stress")
194  .description("Initial stress tensor.")
195  .units( UnitSI().Pa() )
196  .input_default("0.0")
197  .flags_add(in_rhs);
198 
199  *this += region_id.name("region_id")
200  .units( UnitSI::dimensionless())
202 
203  *this += subdomain.name("subdomain")
204  .units( UnitSI::dimensionless() )
206 
207  *this+=cross_section
208  .name("cross_section")
209  .units( UnitSI().m(3).md() )
210  .flags(input_copy & in_time_term & in_main_matrix & in_rhs);
211 
212  *this+=cross_section_min
213  .name("cross_section_min")
214  .description("Minimal cross-section of fractures.")
215  .units( UnitSI().m(3).md() )
216  .input_default("0.0");
217 
218  *this+=potential_load
219  .name("potential_load")
220  .units( UnitSI().m() )
221  .flags(input_copy & in_rhs);
222 
223  *this+=output_field
224  .name("displacement")
225  .description("Displacement vector field output.")
226  .units( UnitSI().m() )
227  .flags(equation_result);
228 
229  *this += output_stress
230  .name("stress")
231  .description("Stress tensor output.")
232  .units( UnitSI().Pa() )
233  .flags(equation_result);
234 
235  *this += output_von_mises_stress
236  .name("von_mises_stress")
237  .description("von Mises stress output.")
238  .units( UnitSI().Pa() )
239  .flags(equation_result);
240 
241  *this += output_mean_stress
242  .name("mean_stress")
243  .description("mean stress output.")
244  .units( UnitSI().Pa() )
245  .flags(equation_result);
246 
247  *this += output_cross_section
248  .name("cross_section_updated")
249  .description("Cross-section after deformation - output.")
250  .units( UnitSI().m() )
251  .flags(equation_result);
252 
253  *this += output_divergence
254  .name("displacement_divergence")
255  .description("Displacement divergence output.")
256  .units( UnitSI().dimensionless() )
257  .flags(equation_result);
258 
259  *this += lame_mu.name("lame_mu")
260  .description("Field lame_mu.")
261  .input_default("0.0")
262  .units( UnitSI().Pa() );
263 
264  *this += lame_lambda.name("lame_lambda")
265  .description("Field lame_lambda.")
266  .input_default("0.0")
267  .units( UnitSI().Pa() );
268 
269  *this += dirichlet_penalty.name("dirichlet_penalty")
270  .description("Field dirichlet_penalty.")
271  .input_default("0.0")
272  .units( UnitSI().Pa() );
273 
274  // add all input fields to the output list
275  output_fields += *this;
276 
277  this->add_coords_field();
278  this->set_default_fieldset();
279 
280 }
281 
282 void Elasticity::EqData::create_dh(Mesh * mesh, unsigned int fe_order)
283 {
284  ASSERT_EQ(fe_order, 1)(fe_order).error("Unsupported polynomial order for finite elements in Elasticity");
285  MixedPtr<FE_P> fe_p(1);
287 
288  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>(mesh, fe);
289  dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh);
290 
291  dh_->distribute_dofs(ds);
292 
293 
294  MixedPtr<FE_P_disc> fe_p_disc(0);
295  dh_scalar_ = make_shared<DOFHandlerMultiDim>(*mesh);
296  std::shared_ptr<DiscreteSpace> ds_scalar = std::make_shared<EqualOrderDiscreteSpace>( mesh, fe_p_disc);
297  dh_scalar_->distribute_dofs(ds_scalar);
298 
299 
301  dh_tensor_ = make_shared<DOFHandlerMultiDim>(*mesh);
302  std::shared_ptr<DiscreteSpace> dst = std::make_shared<EqualOrderDiscreteSpace>( mesh, fe_t);
303  dh_tensor_->distribute_dofs(dst);
304 }
305 
306 
307 Elasticity::Elasticity(Mesh & init_mesh, const Input::Record in_rec, TimeGovernor *tm)
308  : EquationBase(init_mesh, in_rec),
309  input_rec(in_rec),
310  stiffness_assembly_(nullptr),
311  rhs_assembly_(nullptr),
312  constraint_assembly_(nullptr),
313  output_fields_assembly_(nullptr)
314 {
315  // Can not use name() + "constructor" here, since START_TIMER only accepts const char *
316  // due to constexpr optimization.
318 
319  eq_data_ = std::make_shared<EqData>();
320  eq_fields_ = std::make_shared<EqFields>();
321  this->eq_fieldset_ = eq_fields_;
322 
323  auto time_rec = in_rec.val<Input::Record>("time");
324  if (tm == nullptr)
325  {
326  time_ = new TimeGovernor(time_rec);
327  }
328  else
329  {
330  TimeGovernor time_from_rec(time_rec);
331  ASSERT( time_from_rec.is_default() ).error("Duplicate key 'time', time in elasticity is already initialized from parent class!");
332  time_ = tm;
333  }
334 
335 
336  // Set up physical parameters.
337  eq_fields_->set_mesh(init_mesh);
340  eq_data_->balance_ = this->balance();
341 
342  // create finite element structures and distribute DOFs
343  eq_data_->create_dh(mesh_, 1);
344  DebugOut().fmt("Mechanics: solution size {}\n", eq_data_->dh_->n_global_dofs());
345 
346 }
347 
348 
350 {
351  output_stream_ = OutputTime::create_output_stream("mechanics", input_rec.val<Input::Record>("output_stream"), time().get_unit_conversion());
352 
353  eq_fields_->set_components({"displacement"});
354  eq_fields_->set_input_list( input_rec.val<Input::Array>("input_fields"), time() );
355 
356 // balance_ = std::make_shared<Balance>("mechanics", mesh_);
357 // balance_->init_from_input(input_rec.val<Input::Record>("balance"), *time_);
358  // initialization of balance object
359 // eq_data_->balance_idx_ = {
360 // balance_->add_quantity("force_x"),
361 // balance_->add_quantity("force_y"),
362 // balance_->add_quantity("force_z")
363 // };
364 // balance_->units(UnitSI().kg().m().s(-2));
365 
366  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
367  eq_fields_->output_field_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(eq_data_->dh_);
368  eq_fields_->output_field.set(eq_fields_->output_field_ptr, 0.);
369 
370  // setup output stress
371  eq_fields_->output_stress_ptr = create_field_fe<3, FieldValue<3>::TensorFixed>(eq_data_->dh_tensor_);
372  eq_fields_->output_stress.set(eq_fields_->output_stress_ptr, 0.);
373 
374  // setup output von Mises stress
375  eq_fields_->output_von_mises_stress_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_scalar_);
376  eq_fields_->output_von_mises_stress.set(eq_fields_->output_von_mises_stress_ptr, 0.);
377 
378  // setup output mean stress
379  eq_fields_->output_mean_stress_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_scalar_);
380  eq_fields_->output_mean_stress.set(eq_fields_->output_mean_stress_ptr, 0.);
381 
382  // setup output cross-section
383  eq_fields_->output_cross_section_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_scalar_);
384  eq_fields_->output_cross_section.set(eq_fields_->output_cross_section_ptr, 0.);
385 
386  // setup output divergence
387  eq_fields_->output_div_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_scalar_);
388  eq_fields_->output_divergence.set(eq_fields_->output_div_ptr, 0.);
389 
390  // read optional user fields
391  Input::Array user_fields_arr;
392  if (input_rec.opt_val("user_fields", user_fields_arr)) {
393  this->init_user_fields(user_fields_arr, eq_fields_->output_fields);
394  }
395 
396  eq_fields_->output_fields.set_mesh(*mesh_);
397  eq_fields_->output_field.output_type(OutputTime::CORNER_DATA);
398 
399  // set time marks for writing the output
400  eq_fields_->output_fields.initialize(output_stream_, mesh_, input_rec.val<Input::Record>("output"), this->time());
401 
402  // set instances of FieldModel
403  eq_fields_->lame_mu.set(Model<3, FieldValue<3>::Scalar>::create(fn_lame_mu(), eq_fields_->young_modulus, eq_fields_->poisson_ratio), 0.0);
404  eq_fields_->lame_lambda.set(Model<3, FieldValue<3>::Scalar>::create(fn_lame_lambda(), eq_fields_->young_modulus, eq_fields_->poisson_ratio), 0.0);
405  eq_fields_->dirichlet_penalty.set(Model<3, FieldValue<3>::Scalar>::create(fn_dirichlet_penalty(), eq_fields_->lame_mu, eq_fields_->lame_lambda), 0.0);
406 
407  // equation default PETSc solver options
408  std::string petsc_default_opts;
409  petsc_default_opts = "-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg";
410 
411  // allocate matrix and vector structures
412  LinSys *ls;
413  has_contact_ = input_rec.val<bool>("contact");
414  if (has_contact_) {
415 #ifndef FLOW123D_HAVE_PERMON
416  ASSERT(false).error("Flow123d was not built with PERMON library, therefore contact conditions are unsupported.");
417 #endif //FLOW123D_HAVE_PERMON
418  ls = new LinSys_PERMON(eq_data_->dh_->distr().get(), petsc_default_opts);
419 
420  // allocate constraint matrix and vector
421  unsigned int n_own_constraints = 0; // count locally owned cells with neighbours
422  for (auto cell : eq_data_->dh_->own_range())
423  if (cell.elm()->n_neighs_vb() > 0)
424  n_own_constraints++;
425  unsigned int n_constraints = 0; // count all cells with neighbours
426  for (auto elm : mesh_->elements_range())
427  if (elm->n_neighs_vb() > 0)
428  eq_data_->constraint_idx[elm.idx()] = n_constraints++;
429  unsigned int nnz = eq_data_->dh_->ds()->fe()[1_d]->n_dofs()*mesh_->max_edge_sides(1) +
430  eq_data_->dh_->ds()->fe()[2_d]->n_dofs()*mesh_->max_edge_sides(2) +
431  eq_data_->dh_->ds()->fe()[3_d]->n_dofs()*mesh_->max_edge_sides(3);
432  MatCreateAIJ(PETSC_COMM_WORLD, n_own_constraints, eq_data_->dh_->lsize(), PETSC_DECIDE, PETSC_DECIDE, nnz, 0, nnz, 0, &eq_data_->constraint_matrix);
433  VecCreateMPI(PETSC_COMM_WORLD, n_own_constraints, PETSC_DECIDE, &eq_data_->constraint_vec);
434  ((LinSys_PERMON*)ls)->set_inequality(eq_data_->constraint_matrix,eq_data_->constraint_vec);
435 
437  } else {
438  ls = new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
439  ((LinSys_PETSC*)ls)->set_initial_guess_nonzero();
440  }
441  ls->set_from_input( input_rec.val<Input::Record>("solver") );
442  ls->set_solution(eq_fields_->output_field_ptr->vec().petsc_vec());
443  eq_data_->ls = ls;
444 
448 
449  // initialization of balance object
450 // balance_->allocate(eq_data_->dh_->distr()->lsize(),
451 // max(feo->fe<1>()->n_dofs(), max(feo->fe<2>()->n_dofs(), feo->fe<3>()->n_dofs())));
452 
453 }
454 
455 
457 {
458 // delete time_;
459 
460  if (stiffness_assembly_!=nullptr) delete stiffness_assembly_;
461  if (rhs_assembly_!=nullptr) delete rhs_assembly_;
462  if (constraint_assembly_ != nullptr) delete constraint_assembly_;
464 
465  eq_data_.reset();
466  eq_fields_.reset();
467 }
468 
469 
470 
472 {
473  eq_fields_->set_time(time_->step(), LimitSide::right);
474 
475  // update ghost values of solution vector and prepare dependent fields
476  eq_fields_->output_field_ptr->vec().local_to_ghost_begin();
477  eq_fields_->output_stress_ptr->vec().zero_entries();
478  eq_fields_->output_cross_section_ptr->vec().zero_entries();
479  eq_fields_->output_div_ptr->vec().zero_entries();
480  eq_fields_->output_field_ptr->vec().local_to_ghost_end();
481 
482  // compute new output fields depending on solution (stress, divergence etc.)
484 
485  // update ghost values of computed fields
486  eq_fields_->output_stress_ptr->vec().local_to_ghost_begin();
487  eq_fields_->output_von_mises_stress_ptr->vec().local_to_ghost_begin();
488  eq_fields_->output_mean_stress_ptr->vec().local_to_ghost_begin();
489  eq_fields_->output_cross_section_ptr->vec().local_to_ghost_begin();
490  eq_fields_->output_div_ptr->vec().local_to_ghost_begin();
491  eq_fields_->output_stress_ptr->vec().local_to_ghost_end();
492  eq_fields_->output_von_mises_stress_ptr->vec().local_to_ghost_end();
493  eq_fields_->output_mean_stress_ptr->vec().local_to_ghost_end();
494  eq_fields_->output_cross_section_ptr->vec().local_to_ghost_end();
495  eq_fields_->output_div_ptr->vec().local_to_ghost_end();
496 }
497 
498 
499 
500 
502 {
504  eq_fields_->mark_input_times( *time_ );
505  eq_fields_->set_time(time_->step(), LimitSide::right);
506  std::stringstream ss; // print warning message with table of uninitialized fields
507  if ( FieldCommon::print_message_table(ss, "mechanics") ) {
508  WarningOut() << ss.str();
509  }
510 
511  preallocate();
512 
513 
514  // after preallocation we assemble the matrices and vectors required for balance of forces
515 // for (auto subst_idx : eq_data_->balance_idx_)
516 // balance_->calculate_instant(subst_idx, eq_data_->ls->get_solution());
517 
518 // update_solution();
519  eq_data_->ls->start_add_assembly();
520  MatSetOption(*eq_data_->ls->get_matrix(), MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
521  eq_data_->ls->mat_zero_entries();
522  eq_data_->ls->rhs_zero_entries();
525  eq_data_->ls->finish_assembly();
526  LinSys::SolveInfo si = eq_data_->ls->solve();
527  MessageOut().fmt("[mech solver] lin. it: {}, reason: {}, residual: {}\n",
528  si.n_iterations, si.converged_reason, eq_data_->ls->compute_residual());
529  output_data();
530 }
531 
532 
533 
535 {
536  // preallocate system matrix
537  eq_data_->ls->start_allocation();
540 
541  if (has_contact_)
543 }
544 
545 
546 
547 
549 {
550  START_TIMER("DG-ONE STEP");
551 
552  next_time();
554 
556 
557  output_data();
558 
559  END_TIMER("DG-ONE STEP");
560 }
561 
563 {
564  time_->next_time();
565  time_->view("MECH");
566 
567 }
568 
569 
570 
572 {
573  START_TIMER("data reinit");
574  eq_fields_->set_time(time_->step(), LimitSide::right);
575  END_TIMER("data reinit");
576 
577  // assemble stiffness matrix
578  if (eq_data_->ls->get_matrix() == NULL
579  || eq_fields_->subset(FieldFlag::in_main_matrix).changed())
580  {
581  DebugOut() << "Mechanics: Assembling matrix.\n";
582  eq_data_->ls->start_add_assembly();
583  eq_data_->ls->mat_zero_entries();
585  eq_data_->ls->finish_assembly();
586  }
587 
588  // assemble right hand side (due to sources and boundary conditions)
589  if (eq_data_->ls->get_rhs() == NULL
590  || eq_fields_->subset(FieldFlag::in_rhs).changed())
591  {
592  DebugOut() << "Mechanics: Assembling right hand side.\n";
593  eq_data_->ls->start_add_assembly();
594  eq_data_->ls->rhs_zero_entries();
596  eq_data_->ls->finish_assembly();
597  }
598 
599  START_TIMER("solve");
600  LinSys::SolveInfo si = eq_data_->ls->solve();
601  MessageOut().fmt("[mech solver] lin. it: {}, reason: {}, residual: {}\n",
602  si.n_iterations, si.converged_reason, eq_data_->ls->compute_residual());
603  END_TIMER("solve");
604 }
605 
606 
607 
608 
609 
610 
612 {
613  START_TIMER("MECH-OUTPUT");
614 
615  // gather the solution from all processors
616  eq_fields_->output_fields.set_time( this->time().step(), LimitSide::left);
617  //if (eq_fields_->output_fields.is_field_output_time(eq_fields_->output_field, this->time().step()) )
619  eq_fields_->output_fields.output(this->time().step());
620 
621 // START_TIMER("MECH-balance");
622 // balance_->calculate_instant(subst_idx, eq_data_->ls->get_solution());
623 // balance_->output();
624 // END_TIMER("MECH-balance");
625 
626  END_TIMER("MECH-OUTPUT");
627 }
628 
629 
630 
631 
633 {
634 // if (balance_->cumulative())
635 // {
636 // balance_->calculate_cumulative(subst_idx, eq_data_->ls->get_solution());
637 // }
638 }
639 
640 
642 {
643  MatZeroEntries(eq_data_->constraint_matrix);
644  VecZeroEntries(eq_data_->constraint_vec);
646  MatAssemblyBegin(eq_data_->constraint_matrix, MAT_FINAL_ASSEMBLY);
647  MatAssemblyEnd(eq_data_->constraint_matrix, MAT_FINAL_ASSEMBLY);
648  VecAssemblyBegin(eq_data_->constraint_vec);
649  VecAssemblyEnd(eq_data_->constraint_vec);
650 }
651 
652 
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:53
void create_dh(Mesh *mesh, unsigned int fe_order)
Create DOF handler objects.
Definition: elasticity.cc:282
static const Input::Type::Selection & get_bc_type_selection()
Definition: elasticity.cc:117
std::shared_ptr< EqData > eq_data_
Data for model parameters.
Definition: elasticity.hh:235
static const int registrar
Registrar of class to factory.
Definition: elasticity.hh:220
void update_solution() override
Computes the solution in one time instant.
Definition: elasticity.cc:548
GenericAssembly< RhsAssemblyElasticity > * rhs_assembly_
Definition: elasticity.hh:262
void update_output_fields()
Definition: elasticity.cc:471
void assemble_constraint_matrix()
Definition: elasticity.cc:641
GenericAssembly< OutpuFieldsAssemblyElasticity > * output_fields_assembly_
Definition: elasticity.hh:264
static constexpr const char * name_
Definition: elasticity.hh:257
GenericAssembly< StiffnessAssemblyElasticity > * stiffness_assembly_
general assembly objects, hold assembly objects of appropriate dimension
Definition: elasticity.hh:261
void solve_linear_system()
Solve without updating time step and without output.
Definition: elasticity.cc:571
std::shared_ptr< OutputTime > output_stream_
Definition: elasticity.hh:248
void initialize() override
Definition: elasticity.cc:349
Input::Record input_rec
Record with input specification.
Definition: elasticity.hh:251
void calculate_cumulative_balance()
Definition: elasticity.cc:632
void output_data()
Postprocesses the solution and writes to output file.
Definition: elasticity.cc:611
~Elasticity()
Destructor.
Definition: elasticity.cc:456
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: elasticity.cc:48
bool has_contact_
Indicator of contact conditions on fractures.
Definition: elasticity.hh:238
GenericAssembly< ConstraintAssemblyElasticity > * constraint_assembly_
Definition: elasticity.hh:263
Elasticity(Mesh &init_mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
Constructor.
Definition: elasticity.cc:307
void zero_time_step() override
Initialize solution in the zero time.
Definition: elasticity.cc:501
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
Definition: elasticity.hh:232
void preallocate()
Definition: elasticity.cc:534
void next_time()
Pass to next time and update equation data.
Definition: elasticity.cc:562
std::shared_ptr< Balance > balance() const
Definition: equation.hh:189
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:39
std::shared_ptr< FieldSet > eq_fieldset_
Definition: equation.hh:249
void init_user_fields(Input::Array user_fields, FieldSet &output_fields)
Definition: equation.cc:84
TimeGovernor * time_
Definition: equation.hh:241
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
Definition: equation.cc:46
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:95
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Definition: field_flag.hh:49
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
static auto region_id(Mesh &mesh) -> IndexField
static auto subdomain(Mesh &mesh) -> IndexField
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
bool opt_val(const string &key, Ret &value) const
const Ret val(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_solution(Vec sol_vec)
Definition: linsys.hh:290
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:641
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
unsigned int max_edge_sides(unsigned int dim) const
Definition: mesh.h:145
Definition: mesh.h:362
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
static std::shared_ptr< OutputTime > create_output_stream(const std::string &equation_name, const Input::Record &in_rec, const std::shared_ptr< TimeUnitConversion > &time_unit_conv)
This method delete all object instances of class OutputTime stored in output_streams vector.
Definition: output_time.cc:187
Basic time management functionality for unsteady (and steady) solvers (class Equation).
void view(const char *name="") const
const TimeStep & step(int index=-1) const
void next_time()
Proceed to the next time according to current estimated time step.
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
double lame_mu(double young, double poisson)
Definition: elasticity.cc:80
double lame_lambda(double young, double poisson)
Definition: elasticity.cc:86
FEM for linear elasticity.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas finite elements.
Class FESystem for compound finite elements.
MixedPtr< FESystem > mixed_fe_system(MixedPtr< FiniteElement > fe, Args &&... args)
Definition: fe_system.hh:176
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ FETensor
@ FEVector
PERMON QP solvers and FETI.
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
Definitions of particular quadrature rules on simplices.
int converged_reason
Definition: linsys.hh:108
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.