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