49 std::string equation_name = std::string(name_) +
"_FE";
51 std::string(equation_name),
52 "FEM for linear elasticity.")
56 "Settings for computing balance.")
58 "Parameters of output stream.")
60 "Linear solver for elasticity.")
63 .make_field_descriptor_type(equation_name)),
65 "Input fields of the equation.")
67 EqFields().output_fields.make_output_type(equation_name,
""),
68 IT::Default(
"{ \"fields\": [ \"displacement\" ] }"),
69 "Setting of the field output.")
75 Input::register_class< Elasticity, Mesh &, const Input::Record>(std::string(name_) +
"_FE") +
80 double lame_mu(
double young,
double poisson)
82 return young*0.5/(poisson+1.);
88 return young*poisson/((poisson+1.)*(1.-2.*poisson));
93 inline double operator() (
double young,
double poisson) {
94 return young * 0.5 / (poisson+1.);
100 inline double operator() (
double young,
double poisson) {
101 return young * poisson / ((poisson+1.)*(1.-2.*poisson));
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.")
124 "Prescribed traction.")
126 "Prescribed stress tensor.")
136 "Type of boundary condition.")
138 .input_default(
"\"traction\"")
139 .input_selection( get_bc_type_selection() )
142 *
this+=bc_displacement
143 .name(
"bc_displacement")
144 .description(
"Prescribed displacement on boundary.")
146 .input_default(
"0.0")
151 .description(
"Prescribed traction on boundary.")
153 .input_default(
"0.0")
158 .description(
"Prescribed stress on boundary.")
160 .input_default(
"0.0")
165 .description(
"Prescribed bulk load.")
166 .units(
UnitSI().N().m(-3) )
167 .input_default(
"0.0")
171 .name(
"young_modulus")
172 .description(
"Young's modulus.")
174 .input_default(
"0.0")
175 .flags_add(in_main_matrix & in_rhs);
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);
184 *
this+=fracture_sigma
185 .name(
"fracture_sigma")
187 "Coefficient of transfer of forces through fractures.")
189 .input_default(
"1.0")
190 .flags_add(in_main_matrix & in_rhs);
192 *
this+=initial_stress
193 .name(
"initial_stress")
194 .description(
"Initial stress tensor.")
196 .input_default(
"0.0")
199 *
this += region_id.name(
"region_id")
203 *
this += subdomain.name(
"subdomain")
208 .name(
"cross_section")
209 .units(
UnitSI().m(3).md() )
210 .flags(input_copy & in_time_term & in_main_matrix & in_rhs);
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");
218 *
this+=potential_load
219 .name(
"potential_load")
221 .flags(input_copy & in_rhs);
224 .name(
"displacement")
225 .description(
"Displacement vector field output.")
227 .flags(equation_result);
229 *
this += output_stress
231 .description(
"Stress tensor output.")
233 .flags(equation_result);
235 *
this += output_von_mises_stress
236 .name(
"von_mises_stress")
237 .description(
"von Mises stress output.")
239 .flags(equation_result);
241 *
this += output_mean_stress
243 .description(
"mean stress output.")
245 .flags(equation_result);
247 *
this += output_cross_section
248 .name(
"cross_section_updated")
249 .description(
"Cross-section after deformation - output.")
251 .flags(equation_result);
253 *
this += output_divergence
254 .name(
"displacement_divergence")
255 .description(
"Displacement divergence output.")
256 .units(
UnitSI().dimensionless() )
257 .flags(equation_result);
259 *
this +=
lame_mu.name(
"lame_mu")
260 .description(
"Field lame_mu.")
261 .input_default(
"0.0")
265 .description(
"Field lame_lambda.")
266 .input_default(
"0.0")
269 *
this += dirichlet_penalty.name(
"dirichlet_penalty")
270 .description(
"Field dirichlet_penalty.")
271 .input_default(
"0.0")
275 output_fields += *
this;
277 this->add_coords_field();
278 this->set_default_fieldset();
284 ASSERT_EQ(fe_order, 1)(fe_order).error(
"Unsupported polynomial order for finite elements in Elasticity");
288 std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>(mesh, fe);
289 dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh);
291 dh_->distribute_dofs(ds);
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);
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);
310 stiffness_assembly_(nullptr),
311 rhs_assembly_(nullptr),
312 constraint_assembly_(nullptr),
313 output_fields_assembly_(nullptr)
319 eq_data_ = std::make_shared<EqData>();
331 ASSERT( time_from_rec.
is_default() ).error(
"Duplicate key 'time', time in elasticity is already initialized from parent class!");
344 DebugOut().fmt(
"Mechanics: solution size {}\n",
eq_data_->dh_->n_global_dofs());
367 eq_fields_->output_field_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(
eq_data_->dh_);
371 eq_fields_->output_stress_ptr = create_field_fe<3, FieldValue<3>::TensorFixed>(
eq_data_->dh_tensor_);
375 eq_fields_->output_von_mises_stress_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_scalar_);
379 eq_fields_->output_mean_stress_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_scalar_);
383 eq_fields_->output_cross_section_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_scalar_);
387 eq_fields_->output_div_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_scalar_);
408 std::string petsc_default_opts;
409 petsc_default_opts =
"-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg";
415 #ifndef FLOW123D_HAVE_PERMON
416 ASSERT(
false).error(
"Flow123d was not built with PERMON library, therefore contact conditions are unsupported.");
421 unsigned int n_own_constraints = 0;
422 for (
auto cell :
eq_data_->dh_->own_range())
423 if (cell.elm()->n_neighs_vb() > 0)
425 unsigned int n_constraints = 0;
427 if (elm->n_neighs_vb() > 0)
428 eq_data_->constraint_idx[elm.idx()] = n_constraints++;
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);
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();
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();
506 std::stringstream ss;
520 MatSetOption(*
eq_data_->ls->get_matrix(), MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
527 MessageOut().fmt(
"[mech solver] lin. it: {}, reason: {}, residual: {}\n",
578 if (
eq_data_->ls->get_matrix() == NULL
581 DebugOut() <<
"Mechanics: Assembling matrix.\n";
592 DebugOut() <<
"Mechanics: Assembling right hand side.\n";
601 MessageOut().fmt(
"[mech solver] lin. it: {}, reason: {}, residual: {}\n",
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);
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
static const Input::Type::Record & get_input_type()
Main balance input record type.
void create_dh(Mesh *mesh, unsigned int fe_order)
Create DOF handler objects.
static const Input::Type::Selection & get_bc_type_selection()
std::shared_ptr< EqData > eq_data_
Data for model parameters.
static const int registrar
Registrar of class to factory.
void update_solution() override
Computes the solution in one time instant.
GenericAssembly< RhsAssemblyElasticity > * rhs_assembly_
void update_output_fields()
void assemble_constraint_matrix()
GenericAssembly< OutpuFieldsAssemblyElasticity > * output_fields_assembly_
static constexpr const char * name_
GenericAssembly< StiffnessAssemblyElasticity > * stiffness_assembly_
general assembly objects, hold assembly objects of appropriate dimension
void solve_linear_system()
Solve without updating time step and without output.
std::shared_ptr< OutputTime > output_stream_
void initialize() override
Input::Record input_rec
Record with input specification.
void calculate_cumulative_balance()
void output_data()
Postprocesses the solution and writes to output file.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
bool has_contact_
Indicator of contact conditions on fractures.
GenericAssembly< ConstraintAssemblyElasticity > * constraint_assembly_
Elasticity(Mesh &init_mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
Constructor.
void zero_time_step() override
Initialize solution in the zero time.
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
void next_time()
Pass to next time and update equation data.
std::shared_ptr< Balance > balance() const
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
std::shared_ptr< FieldSet > eq_fieldset_
void init_user_fields(Input::Array user_fields, FieldSet &output_fields)
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
static bool print_message_table(ostream &stream, std::string equation_name)
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
static auto region_id(Mesh &mesh) -> IndexField
static auto subdomain(Mesh &mesh) -> IndexField
static const Input::Type::Record & get_input_type()
void set_solution(Vec sol_vec)
virtual void set_from_input(const Input::Record in_rec)
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
unsigned int max_edge_sides(unsigned int dim) const
static const Input::Type::Record & get_input_type()
The specification of 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.
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.
static UnitSI & dimensionless()
Returns dimensionless unit.
double lame_mu(double young, double poisson)
double lame_lambda(double young, double poisson)
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)
Class FEValues calculates finite element data on the actual cells such as shape function values,...
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.
#define MessageOut()
Macro defining 'message' record of log.
#define DebugOut()
Macro defining 'debug' record of log.
Definitions of particular quadrature rules on simplices.
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.