Flow123d  release_3.0.0-1106-ga3b2e4c
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 
22 #include "io/output_time.hh"
24 #include "fem/mapping_p1.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 "coupling/balance.hh"
32 #include "mesh/neighbours.h"
33 
34 #include "fields/multi_field.hh"
35 #include "fields/generic_field.hh"
36 #include "input/factory.hh"
37 
38 
39 
40 
41 using namespace Input::Type;
42 
43 
45  std::string equation_name = std::string(Elasticity::EqData::name()) + "_FE";
46  return IT::Record(
47  std::string(equation_name),
48  "FEM for linear elasticity.")
50  "Time governor setting for the secondary equation.")
51  .declare_key("balance", Balance::get_input_type(), Default("{}"),
52  "Settings for computing balance.")
54  "Parameters of output stream.")
56  "Linear solver for elasticity.")
57  .declare_key("input_fields", Array(
59  .make_field_descriptor_type(equation_name)),
61  "Input fields of the equation.")
62  .declare_key("output",
63  EqData().output_fields.make_output_type(equation_name, ""),
64  IT::Default("{ \"fields\": [ "+Elasticity::EqData::default_output_field()+" ] }"),
65  "Setting of the field output.")
66  .close();
67 }
68 
69 const int Elasticity::registrar =
70  Input::register_class< Elasticity, Mesh &, const Input::Record>(std::string(Elasticity::EqData::name()) + "_FE") +
72 
73 
74 
75 namespace Mechanics {
76 
77 FEObjects::FEObjects(Mesh *mesh_, unsigned int fe_order)
78 {
79  unsigned int q_order;
80 
81  switch (fe_order)
82  {
83  case 1:
84  q_order = 2;
85  fe0_ = new FESystem<0>(std::make_shared<FE_P<0> >(1), FEVector, 3);
86  fe1_ = new FESystem<1>(std::make_shared<FE_P<1> >(1), FEVector, 3);
87  fe2_ = new FESystem<2>(std::make_shared<FE_P<2> >(1), FEVector, 3);
88  fe3_ = new FESystem<3>(std::make_shared<FE_P<3> >(1), FEVector, 3);
89  break;
90 
91  default:
92  q_order=0;
93  xprintf(PrgErr, "Unsupported polynomial order %d for finite elements in Elasticity", fe_order);
94  break;
95  }
96 
97  q0_ = new QGauss<0>(q_order);
98  q1_ = new QGauss<1>(q_order);
99  q2_ = new QGauss<2>(q_order);
100  q3_ = new QGauss<3>(q_order);
101 
102  map1_ = new MappingP1<1,3>;
103  map2_ = new MappingP1<2,3>;
104  map3_ = new MappingP1<3,3>;
105 
106  ds_ = std::make_shared<EqualOrderDiscreteSpace>(mesh_, fe0_, fe1_, fe2_, fe3_);
107  dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
108 
109  dh_->distribute_dofs(ds_);
110 }
111 
112 
114 {
115  delete fe1_;
116  delete fe2_;
117  delete fe3_;
118  delete q0_;
119  delete q1_;
120  delete q2_;
121  delete q3_;
122  delete map1_;
123  delete map2_;
124  delete map3_;
125 }
126 
127 template<> FiniteElement<0> *FEObjects::fe<0>() { return 0; }
128 template<> FiniteElement<1> *FEObjects::fe<1>() { return fe1_; }
129 template<> FiniteElement<2> *FEObjects::fe<2>() { return fe2_; }
130 template<> FiniteElement<3> *FEObjects::fe<3>() { return fe3_; }
131 
132 template<> Quadrature<0> *FEObjects::q<0>() { return q0_; }
133 template<> Quadrature<1> *FEObjects::q<1>() { return q1_; }
134 template<> Quadrature<2> *FEObjects::q<2>() { return q2_; }
135 template<> Quadrature<3> *FEObjects::q<3>() { return q3_; }
136 
137 template<> MappingP1<1,3> *FEObjects::mapping<1>() { return map1_; }
138 template<> MappingP1<2,3> *FEObjects::mapping<2>() { return map2_; }
139 template<> MappingP1<3,3> *FEObjects::mapping<3>() { return map3_; }
140 
141 std::shared_ptr<DOFHandlerMultiDim> FEObjects::dh() { return dh_; }
142 
143 
144 
145 } // namespace Mechanics
146 
147 
148 
149 
150 
151 
152 
154  return Selection("Elasticity_BC_Type", "Types of boundary conditions for heat transfer model.")
155  .add_value(bc_type_displacement, "displacement",
156  "Prescribed displacement.")
157  .add_value(bc_type_traction, "traction",
158  "Prescribed traction.")
159  .close();
160 }
161 
162 
164 {
165  *this+=bc_type
166  .name("bc_type")
167  .description(
168  "Type of boundary condition.")
169  .units( UnitSI::dimensionless() )
170  .input_default("\"traction\"")
171  .input_selection( get_bc_type_selection() )
173 
174  *this+=bc_displacement
175  .name("bc_displacement")
176  .description("Prescribed displacement.")
177  .units( UnitSI().m() )
178  .input_default("0.0")
179  .flags_add(in_rhs);
180 
181  *this+=bc_traction
182  .name("bc_traction")
183  .description("Prescribed traction.")
184  .units( UnitSI().Pa() )
185  .input_default("0.0")
186  .flags_add(in_rhs);
187 
188  *this+=load
189  .name("load")
190  .description("Prescribed load.")
191  .units( UnitSI().N().m(-3) )
192  .input_default("0.0")
193  .flags_add(in_rhs);
194 
195  *this+=young_modulus
196  .name("young_modulus")
197  .description("Young modulus.")
198  .units( UnitSI().Pa() )
199  .input_default("0.0")
200  .flags_add(in_main_matrix);
201 
202  *this+=poisson_ratio
203  .name("poisson_ratio")
204  .description("Poisson ratio.")
205  .units( UnitSI().dimensionless() )
206  .input_default("0.0")
207  .flags_add(in_main_matrix);
208 
209  *this+=fracture_sigma
210  .name("fracture_sigma")
211  .description(
212  "Coefficient of diffusive transfer through fractures (for each substance).")
213  .units( UnitSI::dimensionless() )
214  .input_default("1.0")
215  .flags_add(FieldFlag::in_main_matrix);
216 
217  *this += region_id.name("region_id")
218  .units( UnitSI::dimensionless())
220 
221  *this += subdomain.name("subdomain")
222  .units( UnitSI::dimensionless() )
224 
225  *this+=cross_section
226  .name("cross_section")
227  .units( UnitSI().m(3).md() )
228  .flags(input_copy & in_time_term & in_main_matrix);
229 
230  *this+=output_field
231  .name("displacement")
232  .units( UnitSI().m() )
233  .flags(equation_result);
234 
235 
236  // add all input fields to the output list
237  output_fields += *this;
238 
239 }
240 
241 Elasticity::Elasticity(Mesh & init_mesh, const Input::Record in_rec)
242  : EquationBase(init_mesh, in_rec),
243  input_rec(in_rec),
244  allocation_done(false)
245 {
246  // Can not use name() + "constructor" here, since START_TIMER only accepts const char *
247  // due to constexpr optimization.
249 
250  this->eq_data_ = &data_;
251  time_ = new TimeGovernor(in_rec.val<Input::Record>("time"));
252 
253 
254  // Set up physical parameters.
255  data_.set_mesh(init_mesh);
258 
259  // create finite element structures and distribute DOFs
260  feo = new Mechanics::FEObjects(mesh_, 1);
261  DebugOut().fmt("Mechanics: solution size {}\n", feo->dh()->n_global_dofs());
262 
263 }
264 
265 
267 {
268  output_stream_ = OutputTime::create_output_stream("mechanics", input_rec.val<Input::Record>("output_stream"), time().get_unit_string());
269 
270  data_.set_components({"displacement"});
271  data_.set_input_list( input_rec.val<Input::Array>("input_fields"), time() );
272 
273 // balance_ = std::make_shared<Balance>("mechanics", mesh_);
274 // balance_->init_from_input(input_rec.val<Input::Record>("balance"), *time_);
275  // initialization of balance object
276 // data_.balance_idx_ = {
277 // balance_->add_quantity("force_x"),
278 // balance_->add_quantity("force_y"),
279 // balance_->add_quantity("force_z")
280 // };
281 // balance_->units(UnitSI().kg().m().s(-2));
282 
283  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
284  std::shared_ptr<FieldFE<3, FieldValue<3>::VectorFixed> > output_field_ptr(new FieldFE<3, FieldValue<3>::VectorFixed>);
285  output_vec = output_field_ptr->set_fe_data(feo->dh());
286  data_.output_field.set_field(mesh_->region_db().get_region_set("ALL"), output_field_ptr, 0.);
288 
289  // set time marks for writing the output
291 
292  // equation default PETSc solver options
293  std::string petsc_default_opts;
294  petsc_default_opts = "-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg";
295 
296  // allocate matrix and vector structures
297  ls = new LinSys_PETSC(feo->dh()->distr().get(), petsc_default_opts);
300 
301  // initialization of balance object
302 // balance_->allocate(feo->dh()->distr()->lsize(),
303 // max(feo->fe<1>()->n_dofs(), max(feo->fe<2>()->n_dofs(), feo->fe<3>()->n_dofs())));
304 
305 }
306 
307 
309 {
310 // delete time_;
311 
312  MatDestroy(&stiffness_matrix);
313  VecDestroy(&rhs);
314  delete ls;
315  delete feo;
316 
317 }
318 
319 
320 
322 {
323  VecScatter output_scatter;
324  VecScatterCreateToZero(ls->get_solution(), &output_scatter, PETSC_NULL);
325  // gather solution to output_vec
326  VecScatterBegin(output_scatter, ls->get_solution(), output_vec.petsc_vec(), INSERT_VALUES, SCATTER_FORWARD);
327  VecScatterEnd(output_scatter, ls->get_solution(), output_vec.petsc_vec(), INSERT_VALUES, SCATTER_FORWARD);
328  VecScatterDestroy(&(output_scatter));
329 }
330 
331 
332 
333 
335 {
339  std::stringstream ss; // print warning message with table of uninitialized fields
340  if ( FieldCommon::print_message_table(ss, "mechanics") ) {
341  WarningOut() << ss.str();
342  }
343 
344  ( (LinSys_PETSC *)ls )->set_initial_guess_nonzero();
345 
346  // check first time assembly - needs preallocation
348 
349 
350  // after preallocation we assemble the matrices and vectors required for balance of forces
351 // for (auto subst_idx : data_.balance_idx_)
352 // balance_->calculate_instant(subst_idx, ls->get_solution());
353 
354 // update_solution();
356  MatSetOption(*ls->get_matrix(), MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
357  ls->mat_zero_entries();
359  set_sources();
361  ls->finish_assembly();
362  ls->apply_constrains(1.0);
363  ls->solve();
364  output_data();
365 }
366 
367 
368 
370 {
371  // preallocate system matrix
372  ls->start_allocation();
373  stiffness_matrix = NULL;
374  rhs = NULL;
375 
377  set_sources();
379 
380  allocation_done = true;
381 }
382 
383 
384 
385 
387 {
388  START_TIMER("DG-ONE STEP");
389 
390  time_->next_time();
391  time_->view("MECH");
392 
393  START_TIMER("data reinit");
395  END_TIMER("data reinit");
396 
397  // assemble stiffness matrix
398  if (stiffness_matrix == NULL
400  || flux_changed)
401  {
403  ls->mat_zero_entries();
405  ls->finish_assembly();
406 
407  if (stiffness_matrix == NULL)
408  MatConvert(*( ls->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix);
409  else
410  MatCopy(*( ls->get_matrix() ), stiffness_matrix, DIFFERENT_NONZERO_PATTERN);
411  }
412 
413  // assemble right hand side (due to sources and boundary conditions)
414  if (rhs == NULL
415  || data_.subset(FieldFlag::in_rhs).changed()
416  || flux_changed)
417  {
419  ls->rhs_zero_entries();
420  set_sources();
422  ls->finish_assembly();
423  ls->apply_constrains(1.0);
424 
425  if (rhs == nullptr) VecDuplicate(*( ls->get_rhs() ), &rhs);
426  VecCopy(*( ls->get_rhs() ), rhs);
427  }
428 
429  flux_changed = false;
430 
431 
432  START_TIMER("solve");
433  ls->solve();
434  END_TIMER("solve");
435 
437 
438  output_data();
439 
440  END_TIMER("DG-ONE STEP");
441 }
442 
443 
444 
445 
447 {
448  START_TIMER("MECH-OUTPUT");
449 
450  // gather the solution from all processors
451  data_.output_fields.set_time( this->time().step(), LimitSide::left);
452  //if (data_.output_fields.is_field_output_time(data_.output_field, this->time().step()) )
454  data_.output_fields.output(this->time().step());
455  output_stream_->write_time_frame();
456 
457 // START_TIMER("MECH-balance");
458 // balance_->calculate_instant(subst_idx, ls->get_solution());
459 // balance_->output();
460 // END_TIMER("MECH-balance");
461 
462  END_TIMER("MECH-OUTPUT");
463 }
464 
465 
466 
468 {
469 // if (balance_->cumulative())
470 // {
471 // balance_->calculate_cumulative(subst_idx, ls->get_solution());
472 // }
473 }
474 
475 
476 
477 
478 
479 
480 
481 
483 {
484  START_TIMER("assemble_stiffness");
485  START_TIMER("assemble_volume_integrals");
486  assemble_volume_integrals<1>();
487  assemble_volume_integrals<2>();
488  assemble_volume_integrals<3>();
489  END_TIMER("assemble_volume_integrals");
490 
491  START_TIMER("assemble_fluxes_boundary");
492  assemble_fluxes_boundary<1>();
493  assemble_fluxes_boundary<2>();
494  assemble_fluxes_boundary<3>();
495  END_TIMER("assemble_fluxes_boundary");
496 
497  START_TIMER("assemble_fluxes_elem_side");
498  assemble_fluxes_element_side<1>();
499  assemble_fluxes_element_side<2>();
500  assemble_fluxes_element_side<3>();
501  END_TIMER("assemble_fluxes_elem_side");
502  END_TIMER("assemble_stiffness");
503 }
504 
505 
506 double lame_mu(double young, double poisson)
507 {
508  return young*0.5/(poisson+1.);
509 }
510 
511 
512 double lame_lambda(double young, double poisson)
513 {
514  return young*poisson/((poisson+1.)*(1.-2.*poisson));
515 }
516 
517 
518 template<unsigned int dim>
520 {
521  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
523  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
524  vector<int> dof_indices(ndofs);
525  vector<arma::vec3> velocity(qsize);
526  vector<double> young(qsize), poisson(qsize), csection(qsize);
527  PetscScalar local_matrix[ndofs*ndofs];
528  auto vec = fe_values.vector_view(0);
529 
530  // assemble integral over elements
531  for (auto cell : feo->dh()->own_range())
532  {
533  if (cell.dim() != dim) continue;
534  ElementAccessor<3> elm_acc = cell.elm();
535 
536  fe_values.reinit(elm_acc);
537  cell.get_dof_indices(dof_indices);
538 
539  data_.cross_section.value_list(fe_values.point_list(), elm_acc, csection);
540  data_.young_modulus.value_list(fe_values.point_list(), elm_acc, young);
541  data_.poisson_ratio.value_list(fe_values.point_list(), elm_acc, poisson);
542 
543  // assemble the local stiffness matrix
544  for (unsigned int i=0; i<ndofs; i++)
545  for (unsigned int j=0; j<ndofs; j++)
546  local_matrix[i*ndofs+j] = 0;
547 
548  for (unsigned int k=0; k<qsize; k++)
549  {
550  double mu = lame_mu(young[k], poisson[k]);
551  double lambda = lame_lambda(young[k], poisson[k]);
552  for (unsigned int i=0; i<ndofs; i++)
553  {
554  for (unsigned int j=0; j<ndofs; j++)
555  local_matrix[i*ndofs+j] += csection[k]*(
556  2*mu*arma::dot(vec.sym_grad(j,k), vec.sym_grad(i,k))
557  + lambda*vec.divergence(j,k)*vec.divergence(i,k)
558  )*fe_values.JxW(k);
559  }
560  }
561  ls->mat_set_values(ndofs, dof_indices.data(), ndofs, dof_indices.data(), local_matrix);
562  }
563 }
564 
565 
566 
568 {
569  START_TIMER("assemble_sources");
570 // balance_->start_source_assembly(subst_idx);
571  set_sources<1>();
572  set_sources<2>();
573  set_sources<3>();
574 // balance_->finish_source_assembly(subst_idx);
575  END_TIMER("assemble_sources");
576 }
577 
578 
579 template<unsigned int dim>
581 {
582  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
584  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
585  vector<arma::vec3> load(qsize);
586  vector<double> csection(qsize);
587  vector<int> dof_indices(ndofs);
588  PetscScalar local_rhs[ndofs];
589  vector<PetscScalar> local_source_balance_vector(ndofs), local_source_balance_rhs(ndofs);
590  auto vec = fe_values.vector_view(0);
591 
592  // assemble integral over elements
593  for (auto cell : feo->dh()->own_range())
594  {
595  if (cell.dim() != dim) continue;
596  ElementAccessor<3> elm_acc = cell.elm();
597 
598  fe_values.reinit(elm_acc);
599  cell.get_dof_indices(dof_indices);
600 
601  // assemble the local stiffness matrix
602  fill_n(local_rhs, ndofs, 0);
603  local_source_balance_vector.assign(ndofs, 0);
604  local_source_balance_rhs.assign(ndofs, 0);
605 
606  data_.cross_section.value_list(fe_values.point_list(), elm_acc, csection);
607  data_.load.value_list(fe_values.point_list(), elm_acc, load);
608 
609  // compute sources
610  for (unsigned int k=0; k<qsize; k++)
611  {
612  for (unsigned int i=0; i<ndofs; i++)
613  local_rhs[i] += arma::dot(load[k], vec.value(i,k))*csection[k]*fe_values.JxW(k);
614  }
615  ls->rhs_set_values(ndofs, &(dof_indices[0]), local_rhs);
616 
617 // for (unsigned int i=0; i<ndofs; i++)
618 // {
619 // for (unsigned int k=0; k<qsize; k++)
620 // local_source_balance_vector[i] -= 0;//sources_sigma[k]*fe_values[vec].value(i,k)*fe_values.JxW(k);
621 //
622 // local_source_balance_rhs[i] += local_rhs[i];
623 // }
624 // balance_->add_source_matrix_values(subst_idx, elm_acc.region().bulk_idx(), dof_indices, local_source_balance_vector);
625 // balance_->add_source_vec_values(subst_idx, elm_acc.region().bulk_idx(), dof_indices, local_source_balance_rhs);
626  }
627 }
628 
629 
630 
631 
632 
633 
634 
635 template<unsigned int dim>
637 {
638  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
640  const unsigned int ndofs = feo->fe<dim>()->n_dofs();
641  vector<int> side_dof_indices(ndofs);
642  PetscScalar local_matrix[ndofs*ndofs];
643 
644  // assemble boundary integral
645  for (unsigned int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
646  {
647  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
648  if (edg->n_sides > 1) continue;
649  // check spatial dimension
650  if (edg->side(0)->dim() != dim-1) continue;
651  // skip edges lying not on the boundary
652  if (edg->side(0)->cond() == NULL) continue;
653 
654  SideIter side = edg->side(0);
655  ElementAccessor<3> cell = side->element();
656  feo->dh()->cell_accessor_from_element(cell.idx()).get_dof_indices(side_dof_indices);
657  fe_values_side.reinit(cell, side->side_idx());
658 // unsigned int bc_type = data_.bc_type.value(side->centre(), side->cond()->element_accessor());
659 
660 // for (unsigned int i=0; i<ndofs; i++)
661 // for (unsigned int j=0; j<ndofs; j++)
662 // local_matrix[i*ndofs+j] = 0;
663 //
664 // ls->mat_set_values(ndofs, side_dof_indices.data(), ndofs, side_dof_indices.data(), local_matrix);
665  }
666 }
667 
668 
670 {
671  arma::mat33 mt = m - m*arma::kron(n,n.t());
672  return mt;
673 }
674 
675 
676 template<unsigned int dim>
678 {
679  if (dim == 1) return;
680  FEValues<dim-1,3> fe_values_sub(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe<dim-1>(),
682  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
684 
685  vector<FEValuesSpaceBase<3>*> fv_sb(2);
686  const unsigned int ndofs_side = feo->fe<dim>()->n_dofs(); // number of local dofs
687  const unsigned int ndofs_sub = feo->fe<dim-1>()->n_dofs();
688  const unsigned int qsize = feo->q<dim-1>()->size(); // number of quadrature points
689  vector<vector<int> > side_dof_indices(2);
690  vector<unsigned int> n_dofs = { ndofs_sub, ndofs_side };
691  vector<double> frac_sigma(qsize);
692  vector<double> csection_lower(qsize), csection_higher(qsize), young(qsize), poisson(qsize), alpha(qsize);
693  PetscScalar local_matrix[2][2][(ndofs_side)*(ndofs_side)];
694  PetscScalar local_rhs[2][ndofs_side];
695  auto vec_side = fe_values_side.vector_view(0);
696  auto vec_sub = fe_values_sub.vector_view(0);
697 
698  // index 0 = element with lower dimension,
699  // index 1 = side of element with higher dimension
700  side_dof_indices[0].resize(ndofs_sub);
701  side_dof_indices[1].resize(ndofs_side);
702  fv_sb[0] = &fe_values_sub;
703  fv_sb[1] = &fe_values_side;
704 
705  // assemble integral over sides
706  for (unsigned int inb=0; inb<feo->dh()->n_loc_nb(); inb++)
707  {
708  Neighbour *nb = &mesh_->vb_neighbours_[feo->dh()->nb_index(inb)];
709  // skip neighbours of different dimension
710  if (nb->element()->dim() != dim-1) continue;
711 
712  ElementAccessor<3> cell_sub = nb->element();
713  feo->dh()->cell_accessor_from_element(cell_sub.idx()).get_dof_indices(side_dof_indices[0]);
714  fe_values_sub.reinit(cell_sub);
715 
716  ElementAccessor<3> cell = nb->side()->element();
717  feo->dh()->cell_accessor_from_element(cell.idx()).get_dof_indices(side_dof_indices[1]);
718  fe_values_side.reinit(cell, nb->side()->side_idx());
719 
720  // Element id's for testing if they belong to local partition.
721  bool own_element_id[2];
722  own_element_id[0] = feo->dh()->cell_accessor_from_element(cell_sub.idx()).is_own();
723  own_element_id[1] = feo->dh()->cell_accessor_from_element(cell.idx()).is_own();
724 
725  data_.cross_section.value_list(fe_values_sub.point_list(), cell_sub, csection_lower);
726  data_.cross_section.value_list(fe_values_sub.point_list(), cell, csection_higher);
727  data_.fracture_sigma.value_list(fe_values_sub.point_list(), cell_sub, frac_sigma);
728  data_.young_modulus.value_list(fe_values_sub.point_list(), cell_sub, young);
729  data_.poisson_ratio.value_list(fe_values_sub.point_list(), cell_sub, poisson);
730 
731  for (unsigned int n=0; n<2; ++n)
732  {
733  for (unsigned int i=0; i<ndofs_side; i++)
734  {
735  for (unsigned int m=0; m<2; ++m)
736  {
737  for (unsigned int j=0; j<ndofs_side; j++)
738  local_matrix[n][m][i*(ndofs_side)+j] = 0;
739  }
740  local_rhs[n][i] = 0;
741  }
742  }
743 
744  // set transmission conditions
745  for (unsigned int k=0; k<qsize; k++)
746  {
747  arma::vec3 nv = fe_values_side.normal_vector(k);
748  double mu = lame_mu(young[k], poisson[k]);
749  double lambda = lame_lambda(young[k], poisson[k]);
750 
751  for (int n=0; n<2; n++)
752  {
753  if (!own_element_id[n]) continue;
754 
755  for (unsigned int i=0; i<n_dofs[n]; i++)
756  {
757  arma::vec3 vi = (n==0)?arma::zeros(3):vec_side.value(i,k);
758  arma::vec3 vf = (n==1)?arma::zeros(3):vec_sub.value(i,k);
759  arma::mat33 gvft = (n==0)?vec_sub.grad(i,k):arma::zeros(3,3);
760 
761  for (int m=0; m<2; m++)
762  for (unsigned int j=0; j<n_dofs[m]; j++) {
763  arma::vec3 ui = (m==0)?arma::zeros(3):vec_side.value(j,k);
764  arma::vec3 uf = (m==1)?arma::zeros(3):vec_sub.value(j,k);
765  arma::mat33 guit = (m==1)?mat_t(vec_side.grad(j,k),nv):arma::zeros(3,3);
766  double divuit = (m==1)?arma::trace(guit):0;
767 
768  local_matrix[n][m][i*n_dofs[m] + j] +=
769  (
770  frac_sigma[k]*arma::dot(vf-vi,
771  2/csection_lower[k]*(mu*(uf-ui)+(mu+lambda)*(arma::dot(uf-ui,nv)*nv))
772  + mu*arma::trans(guit)*nv
773  + lambda*divuit*nv
774  )
775  - arma::dot(gvft, mu*arma::kron(nv,ui.t()) + lambda*arma::dot(ui,nv)*arma::eye(3,3))
776  )*fv_sb[0]->JxW(k);
777  }
778  }
779  }
780  }
781 
782  for (unsigned int n=0; n<2; ++n)
783  {
784  for (unsigned int m=0; m<2; ++m)
785  ls->mat_set_values(n_dofs[n], side_dof_indices[n].data(), n_dofs[m], side_dof_indices[m].data(), local_matrix[n][m]);
786 
787  ls->rhs_set_values(n_dofs[n], side_dof_indices[n].data(), local_rhs[n]);
788  }
789  }
790 }
791 
792 
793 
794 
795 
796 
797 
799 {
800  START_TIMER("assemble_bc");
801 // balance_->start_flux_assembly(subst_idx);
802  set_boundary_conditions<1>();
803  set_boundary_conditions<2>();
804  set_boundary_conditions<3>();
805 // balance_->finish_flux_assembly(subst_idx);
806  END_TIMER("assemble_bc");
807 }
808 
809 
810 
811 template<unsigned int dim>
813 {
814  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
816  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
817  vector<int> side_dof_indices(ndofs);
818  unsigned int loc_b=0;
819  double local_rhs[ndofs];
820  vector<PetscScalar> local_flux_balance_vector(ndofs);
821  PetscScalar local_flux_balance_rhs;
822  vector<arma::vec3> bc_values(qsize), bc_traction(qsize);
823  vector<double> csection(qsize);
824  auto vec = fe_values_side.vector_view(0);
825 
826  for (auto cell : feo->dh()->own_range())
827  {
828  ElementAccessor<3> elm = cell.elm();
829  if (elm->boundary_idx_ == nullptr) continue;
830 
831  for (unsigned int si=0; si<elm->n_sides(); ++si)
832  {
833  const Edge *edg = elm.side(si)->edge();
834  if (edg->n_sides > 1) continue;
835  // skip edges lying not on the boundary
836  if (edg->side(0)->cond() == NULL) continue;
837 
838  if (edg->side(0)->dim() != dim-1)
839  {
840  if (edg->side(0)->cond() != nullptr) ++loc_b;
841  continue;
842  }
843 
844  SideIter side = edg->side(0);
845  ElementAccessor<3> cell = side->element();
846  ElementAccessor<3> bc_cell = side->cond()->element_accessor();
847 
848  unsigned int bc_type = data_.bc_type.value(side->centre(), bc_cell);
849 
850  fe_values_side.reinit(cell, side->side_idx());
851 
852  data_.cross_section.value_list(fe_values_side.point_list(), cell, csection);
853  // The b.c. data are fetched for all possible b.c. types since we allow
854  // different bc_type for each substance.
855  data_.bc_displacement.value_list(fe_values_side.point_list(), bc_cell, bc_values);
856  data_.bc_traction.value_list(fe_values_side.point_list(), bc_cell, bc_traction);
857 
858  feo->dh()->cell_accessor_from_element(cell.idx()).get_dof_indices(side_dof_indices);
859 
860  fill_n(local_rhs, ndofs, 0);
861  local_flux_balance_vector.assign(ndofs, 0);
862  local_flux_balance_rhs = 0;
863 
864  if (bc_type == EqData::bc_type_displacement)
865  {
866  // We solve a local problem to determine which local dof
867  // to assign the boundary value to.
868  // TODO: Should be possibly optimized by inverting the matrix only once within reference element.
869  arma::mat d_mat(ndofs,ndofs);
870  arma::vec d_vec(ndofs);
871  for (unsigned int i=0; i<ndofs; i++)
872  {
873  d_vec(i) = 0;
874  for (unsigned int k=0; k<qsize; k++)
875  d_vec(i) += arma::dot(vec.value(i,k),bc_values[k])*fe_values_side.JxW(k);
876  for (unsigned int j=i; j<ndofs; j++)
877  {
878  d_mat(i,j) = 0;
879  for (unsigned int k=0; k<qsize; k++)
880  d_mat(i,j) += arma::dot(vec.value(i,k),vec.value(j,k))*fe_values_side.JxW(k);
881  d_mat(j,i) = d_mat(i,j);
882  }
883  }
884  arma::vec d_val = pinv(d_mat)*d_vec;
885 
886  for (unsigned int i=0; i<ndofs; i++)
887  if (norm(d_mat.row(i)) > 0)
888  ls->add_constraint(side_dof_indices[i], d_val(i));
889  }
890  else if (bc_type == EqData::bc_type_traction)
891  {
892  for (unsigned int k=0; k<qsize; k++)
893  {
894  for (unsigned int i=0; i<ndofs; i++)
895  local_rhs[i] += csection[k]*arma::dot(vec.value(i,k),bc_traction[k])*fe_values_side.JxW(k);
896  }
897  }
898  ls->rhs_set_values(ndofs, &(side_dof_indices[0]), local_rhs);
899 
900 
901 
902 // balance_->add_flux_matrix_values(subst_idx, loc_b, side_dof_indices, local_flux_balance_vector);
903 // balance_->add_flux_vec_value(subst_idx, loc_b, local_flux_balance_rhs);
904  ++loc_b;
905  }
906  }
907 }
908 
909 
910 
911 
TimeGovernor & time()
Definition: equation.hh:148
BCField< 3, FieldValue< 3 >::Enum > bc_type
Definition: elasticity.hh:124
void output_type(OutputTime::DiscreteSpace rt)
Definition: field_set.hh:211
const Edge * edge() const
Definition: side_impl.hh:66
BCField< 3, FieldValue< 3 >::VectorFixed > bc_displacement
Definition: elasticity.hh:125
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell...
Shape function values.
Definition: update_flags.hh:87
Field< 3, FieldValue< 3 >::Scalar > young_modulus
Definition: elasticity.hh:128
FieldSet * eq_data_
Definition: equation.hh:232
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:48
static auto subdomain(Mesh &mesh) -> IndexField
static constexpr const char * name()
Definition: elasticity.hh:115
Transformed quadrature weight for cell sides.
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:329
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
unsigned int * boundary_idx_
Definition: elements.h:83
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:598
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Definition: field_flag.hh:49
void reinit(ElementAccessor< 3 > &cell, unsigned int sid)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:610
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
void preallocate()
Definition: elasticity.cc:369
Field< 3, FieldValue< 3 >::Scalar > region_id
Definition: elasticity.hh:134
unsigned int side_idx() const
Definition: side_impl.hh:82
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
void set_from_input(const Input::Record in_rec) override
virtual void start_add_assembly()
Definition: linsys.hh:341
void output(TimeStep step)
Boundary * cond() const
Definition: side_impl.hh:71
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
VectorMPI output_vec
Vector of solution data.
Definition: elasticity.hh:299
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
void next_time()
Proceed to the next time according to current estimated time step.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: elasticity.cc:44
arma::vec3 centre() const
Centre of side.
Definition: sides.cc:121
FiniteElement< dim > * fe()
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
void set_boundary_conditions()
Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions.
Definition: elasticity.cc:798
static std::shared_ptr< OutputTime > create_output_stream(const std::string &equation_name, const Input::Record &in_rec, std::string unit_str)
This method delete all object instances of class OutputTime stored in output_streams vector...
Definition: output_time.cc:185
LinSys * ls
Linear algebra system for the transport equation.
Definition: elasticity.hh:290
Definition: mesh.h:76
Fields computed from the mesh data.
virtual void start_allocation()
Definition: linsys.hh:333
Class FEValues calculates finite element data on the actual cells such as shape function values...
virtual void finish_assembly()=0
Field< 3, FieldValue< 3 >::Scalar > cross_section
Pointer to DarcyFlow field cross_section.
Definition: elasticity.hh:133
Class FESystem for compound finite elements.
int n_sides
Definition: edges.h:36
~Elasticity()
Destructor.
Definition: elasticity.cc:308
Definition: edges.h:26
SideIter side(const unsigned int loc_index)
Definition: accessors.hh:137
double lame_mu(double young, double poisson)
Definition: elasticity.cc:506
const RegionDB & region_db() const
Definition: mesh.h:143
std::shared_ptr< DOFHandlerMultiDim > dh()
const TimeStep & step(int index=-1) const
void add_constraint(int row, double value)
Definition: linsys.hh:494
Elasticity(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
Definition: elasticity.cc:241
void update_solution() override
Computes the solution in one time instant.
Definition: elasticity.cc:386
arma::mat33 mat_t(const arma::mat33 &m, const arma::vec3 &n)
Definition: elasticity.cc:669
Basic time management functionality for unsteady (and steady) solvers (class Equation).
ElementAccessor< 3 > element() const
Definition: side_impl.hh:53
virtual void value_list(const std::vector< Point > &point_list, const ElementAccessor< spacedim > &elm, std::vector< typename Value::return_type > &value_list) const
Definition: field.hh:403
bool flux_changed
Indicator of change in advection vector field.
Definition: elasticity.hh:319
Input::Record input_rec
Record with input specification.
Definition: elasticity.hh:304
Class for declaration of inputs sequences.
Definition: type_base.hh:346
virtual void apply_constrains(double scalar)=0
Quadrature< dim > * q()
Vec rhs
Vector of right hand side.
Definition: elasticity.hh:284
void view(const char *name="") const
Symmetric Gauss-Legendre quadrature formulae on simplices.
FEM for linear elasticity.
bool allocation_done
Indicates whether matrices have been preallocated.
Definition: elasticity.hh:316
unsigned int dim() const
Definition: elements.h:124
ElementAccessor< 3 > element()
Definition: neighbours.h:163
Transformed quadrature points.
unsigned int dim() const
Definition: side_impl.hh:38
BCField< 3, FieldValue< 3 >::VectorFixed > bc_traction
Definition: elasticity.hh:126
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:101
const Vec & get_solution()
Definition: linsys.hh:282
double JxW(const unsigned int point_no) override
Return the product of Jacobian determinant and the quadrature weight at given quadrature point...
Definition: fe_values.hh:336
Definitions of basic Lagrangean finite elements with polynomial shape functions.
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
SideIter side()
Definition: neighbours.h:147
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
const Ret val(const string &key) const
unsigned int n_sides() const
Definition: elements.h:135
static auto region_id(Mesh &mesh) -> IndexField
#define xprintf(...)
Definition: system.hh:92
#define START_TIMER(tag)
Starts a timer with specified tag.
std::shared_ptr< DOFHandlerMultiDim > dh()
Definition: elasticity.cc:141
void set_sources()
Assembles the right hand side due to volume sources.
Definition: elasticity.cc:567
Mesh * mesh_
Definition: equation.hh:223
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.
virtual Value::return_type const & value(const Point &p, const ElementAccessor< spacedim > &elm) const
Definition: field.hh:389
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:501
void zero_time_step() override
Initialize solution in the zero time.
Definition: elasticity.cc:334
Vec & petsc_vec()
Getter for PETSC vector of output data (e.g. can be used by scatters).
Definition: vector_mpi.hh:120
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
Shape function gradients.
Definition: update_flags.hh:95
void output_data()
Postprocesses the solution and writes to output file.
Definition: elasticity.cc:446
Mat stiffness_matrix
The stiffness matrix.
Definition: elasticity.hh:287
const FEValuesViews::Vector< dim, spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
Definition: fe_values.hh:388
EqData & data()
Definition: elasticity.hh:187
virtual const Vec * get_rhs()
Definition: linsys.hh:203
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no) override
Returns the normal vector to a side at given quadrature point.
Definition: fe_values.hh:368
Normal vectors.
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
Field< 3, FieldValue< 3 >::VectorFixed > load
Definition: elasticity.hh:127
void set_solution(Vec sol_vec)
Definition: linsys.hh:290
Field< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures.
Definition: elasticity.hh:130
void assemble_fluxes_boundary()
Assembles the fluxes on the boundary.
Definition: elasticity.cc:636
void mark_input_times(const TimeGovernor &tg)
Definition: field_set.hh:218
vector< Neighbour > vb_neighbours_
Definition: mesh.h:273
void set_input_list(Input::Array input_list, const TimeGovernor &tg)
Definition: field_set.hh:190
Field< 3, FieldValue< 3 >::Scalar > poisson_ratio
Definition: elasticity.hh:129
void initialize(std::shared_ptr< OutputTime > stream, Mesh *mesh, Input::Record in_rec, const TimeGovernor &tg)
const Selection & close() const
Close the Selection, no more values can be added.
void set_components(const std::vector< string > &names)
Definition: field_set.hh:177
void assemble_fluxes_element_side()
Assembles the fluxes between elements of different dimensions.
Definition: elasticity.cc:677
EquationOutput output_fields
Definition: elasticity.hh:139
Definition: system.hh:64
static const int registrar
Registrar of class to factory.
Definition: elasticity.hh:197
void set_field(const RegionSet &domain, FieldBasePtr field, double time=0.0)
Definition: field.impl.hh:198
bool set_time(const TimeStep &time, LimitSide limit_side)
Definition: field_set.cc:157
std::vector< Edge > edges
Vector of MH edges, this should not be part of the geometrical mesh.
Definition: mesh.h:252
void initialize() override
Definition: elasticity.cc:266
EqData data_
Field data for model parameters.
Definition: elasticity.hh:265
double lame_lambda(double young, double poisson)
Definition: elasticity.cc:512
static const Input::Type::Selection & get_bc_type_selection()
Definition: elasticity.cc:153
static const Input::Type::Record & get_input_type()
Definitions of particular quadrature rules on simplices.
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:246
virtual SolveInfo solve()=0
#define END_TIMER(tag)
Ends a timer with specified tag.
static string default_output_field()
Definition: elasticity.hh:117
vector< arma::vec::fixed< spacedim > > & point_list()
Return coordinates of all quadrature points in the actual cell system.
Definition: fe_values.hh:357
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_mesh(const Mesh &mesh)
Definition: field_set.hh:183
void assemble_stiffness_matrix()
Assembles the stiffness matrix.
Definition: elasticity.cc:482
Record type proxy class.
Definition: type_record.hh:182
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
Class for representation SI units of Fields.
Definition: unit_si.hh:40
void assemble_volume_integrals()
Assembles the volume integrals into the stiffness matrix.
Definition: elasticity.cc:519
virtual const Mat * get_matrix()
Definition: linsys.hh:187
void calculate_cumulative_balance()
Definition: elasticity.cc:467
MappingP1< dim, 3 > * mapping()
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
#define DebugOut()
Macro defining &#39;debug&#39; record of log.
Definition: logger.hh:252
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:96
unsigned int idx() const
Return local idx of element in boundary / bulk part of element vector.
Definition: accessors.hh:111
Field< 3, FieldValue< 3 >::VectorFixed > output_field
Definition: elasticity.hh:137
SideIter side(const unsigned int i) const
Definition: edges.h:31
void output_vector_gather()
Definition: elasticity.cc:321
bool changed() const
Definition: field_set.cc:165
Template for classes storing finite set of named values.
Field< 3, FieldValue< 3 >::Scalar > subdomain
Definition: elasticity.hh:135
Definitions of Raviart-Thomas finite elements.
ElementAccessor< 3 > element_accessor()
Definition: boundaries.cc:48
std::shared_ptr< OutputTime > output_stream_
Definition: elasticity.hh:301
FEObjects(Mesh *mesh_, unsigned int fe_order)
TimeGovernor * time_
Definition: equation.hh:224
Definition: field.hh:56
Mechanics::FEObjects * feo
Finite element objects.
Definition: elasticity.hh:274
void reinit(ElementAccessor< 3 > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:539
Transformed quadrature weights.
Calculates finite element data on a side.
Definition: fe_values.hh:582