Flow123d  master-f44eb46
assembly_elasticity.hh
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 assembly_elasticity.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_ELASTICITY_HH_
19 #define ASSEMBLY_ELASTICITY_HH_
20 
23 #include "mechanics/elasticity.hh"
24 #include "fem/fe_p.hh"
25 #include "fem/fe_values.hh"
27 #include "coupling/balance.hh"
29 
30 
31 /**
32  * Auxiliary container class for Finite element and related objects of given dimension.
33  */
34 template <unsigned int dim>
36 {
37 public:
38  typedef typename Elasticity::EqFields EqFields;
39  typedef typename Elasticity::EqData EqData;
40 
41  static constexpr const char * name() { return "StiffnessAssemblyElasticity"; }
42 
43  /// Constructor.
45  : AssemblyBase<dim>(1), eq_fields_(eq_fields), eq_data_(eq_data) {
53  }
54 
55  /// Destructor.
57 
58  /// Initialize auxiliary vectors and other data members
59  void initialize(ElementCacheMap *element_cache_map) {
60  //this->balance_ = eq_data_->balance_;
61  this->element_cache_map_ = element_cache_map;
62 
63  shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
64  shared_ptr<FE_P<dim-1>> fe_p_low = std::make_shared< FE_P<dim-1> >(1);
65  fe_ = std::make_shared<FESystem<dim>>(fe_p, FEVector, 3);
66  fe_low_ = std::make_shared<FESystem<dim-1>>(fe_p_low, FEVector, 3);
67  fe_values_.initialize(*this->quad_, *fe_,
71  fe_values_sub_.initialize(*this->quad_low_, *fe_low_,
73 
74  n_dofs_ = fe_->n_dofs();
75  n_dofs_sub_ = fe_low_->n_dofs();
77  dof_indices_.resize(n_dofs_);
78  side_dof_indices_.resize(2);
79  side_dof_indices_[0].resize(n_dofs_sub_); // index 0 = element with lower dimension,
80  side_dof_indices_[1].resize(n_dofs_); // index 1 = side of element with higher dimension
82  local_matrix_ngh_.resize(2);
83  for (uint m=0; m<2; ++m) {
84  local_matrix_ngh_[m].resize(2);
85  for (uint n=0; n<2; ++n)
86  local_matrix_ngh_[m][n].resize(n_dofs_*n_dofs_);
87  }
90  if (dim>1) vec_view_sub_ = &fe_values_sub_.vector_view(0);
91  }
92 
93 
94  /// Assemble integral over element
95  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
96  {
97  if (cell.dim() != dim) return;
98 
99  ElementAccessor<3> elm_acc = cell.elm();
100 
101  fe_values_.reinit(elm_acc);
103 
104  // assemble the local stiffness matrix
105  for (unsigned int i=0; i<n_dofs_; i++)
106  for (unsigned int j=0; j<n_dofs_; j++)
107  local_matrix_[i*n_dofs_+j] = 0;
108 
109  unsigned int k=0;
110  for (auto p : this->bulk_points(element_patch_idx) )
111  {
112  for (unsigned int i=0; i<n_dofs_; i++)
113  {
114  for (unsigned int j=0; j<n_dofs_; j++)
116  2*eq_fields_->lame_mu(p)*arma::dot(vec_view_->sym_grad(j,k), vec_view_->sym_grad(i,k))
118  )*fe_values_.JxW(k);
119  }
120  k++;
121  }
123  }
124 
125  /// Assembles boundary integral.
126  inline void boundary_side_integral(DHCellSide cell_side)
127  {
128  ASSERT_EQ(cell_side.dim(), dim).error("Dimension of element mismatch!");
129  if (!cell_side.cell().is_own()) return;
130 
131  Side side = cell_side.side();
132  const DHCellAccessor &dh_cell = cell_side.cell();
133  dh_cell.get_dof_indices(dof_indices_);
134  fe_values_side_.reinit(side);
135 
136  for (unsigned int i=0; i<n_dofs_; i++)
137  for (unsigned int j=0; j<n_dofs_; j++)
138  local_matrix_[i*n_dofs_+j] = 0;
139 
140  auto p_side = *( this->boundary_points(cell_side).begin() );
141  auto p_bdr = p_side.point_bdr( side.cond().element_accessor() );
142  unsigned int bc_type = eq_fields_->bc_type(p_bdr);
143  double side_measure = cell_side.measure();
144  if (bc_type == EqFields::bc_type_displacement)
145  {
146  unsigned int k=0;
147  for (auto p : this->boundary_points(cell_side) ) {
148  for (unsigned int i=0; i<n_dofs_; i++)
149  for (unsigned int j=0; j<n_dofs_; j++)
150  local_matrix_[i*n_dofs_+j] += (eq_fields_->dirichlet_penalty(p) / side_measure) *
151  arma::dot(vec_view_side_->value(i,k),vec_view_side_->value(j,k)) * fe_values_side_.JxW(k);
152  k++;
153  }
154  }
155  else if (bc_type == EqFields::bc_type_displacement_normal)
156  {
157  unsigned int k=0;
158  for (auto p : this->boundary_points(cell_side) ) {
159  for (unsigned int i=0; i<n_dofs_; i++)
160  for (unsigned int j=0; j<n_dofs_; j++)
161  local_matrix_[i*n_dofs_+j] += (eq_fields_->dirichlet_penalty(p) / side_measure) *
162  arma::dot(vec_view_side_->value(i,k), fe_values_side_.normal_vector(k)) *
164  k++;
165  }
166  }
167 
169  }
170 
171 
172  /// Assembles between elements of different dimensions.
173  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
174  if (dim == 1) return;
175  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
176 
177  cell_lower_dim.get_dof_indices(side_dof_indices_[0]);
178  ElementAccessor<3> cell_sub = cell_lower_dim.elm();
179  fe_values_sub_.reinit(cell_sub);
180 
181  DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element( neighb_side.element().idx() );
182  cell_higher_dim.get_dof_indices(side_dof_indices_[1]);
183  fe_values_side_.reinit(neighb_side.side());
184 
185  // Element id's for testing if they belong to local partition.
186  bool own_element_id[2];
187  own_element_id[0] = cell_lower_dim.is_own();
188  own_element_id[1] = cell_higher_dim.is_own();
189 
190  for (unsigned int n=0; n<2; ++n)
191  for (unsigned int i=0; i<n_dofs_; i++)
192  for (unsigned int m=0; m<2; ++m)
193  for (unsigned int j=0; j<n_dofs_; j++)
194  local_matrix_ngh_[n][m][i*(n_dofs_)+j] = 0;
195 
196  // set transmission conditions
197  unsigned int k=0;
198  for (auto p_high : this->coupling_points(neighb_side) )
199  {
200  auto p_low = p_high.lower_dim(cell_lower_dim);
202 
203  for (int n=0; n<2; n++)
204  {
205  if (!own_element_id[n]) continue;
206 
207  for (unsigned int i=0; i<n_dofs_ngh_[n]; i++)
208  {
209  arma::vec3 vi = (n==0) ? arma::zeros(3) : vec_view_side_->value(i,k);
210  arma::vec3 vf = (n==1) ? arma::zeros(3) : vec_view_sub_->value(i,k);
211  arma::mat33 gvft = (n==0) ? vec_view_sub_->grad(i,k) : arma::zeros(3,3);
212 
213  for (int m=0; m<2; m++)
214  for (unsigned int j=0; j<n_dofs_ngh_[m]; j++) {
215  arma::vec3 ui = (m==0) ? arma::zeros(3) : vec_view_side_->value(j,k);
216  arma::vec3 uf = (m==1) ? arma::zeros(3) : vec_view_sub_->value(j,k);
217  arma::mat33 guft = (m==0) ? mat_t(vec_view_sub_->grad(j,k),nv) : arma::zeros(3,3);
218  double divuft = (m==0) ? arma::trace(guft) : 0;
219 
220  local_matrix_ngh_[n][m][i*n_dofs_ngh_[m] + j] +=
221  eq_fields_->fracture_sigma(p_low)*(
222  arma::dot(vf-vi,
223  2/eq_fields_->cross_section(p_low)*(eq_fields_->lame_mu(p_low)*(uf-ui)+(eq_fields_->lame_mu(p_low)+eq_fields_->lame_lambda(p_low))*(arma::dot(uf-ui,nv)*nv))
224  + eq_fields_->lame_mu(p_low)*arma::trans(guft)*nv
225  + eq_fields_->lame_lambda(p_low)*divuft*nv
226  )
227  - arma::dot(gvft, eq_fields_->lame_mu(p_low)*arma::kron(nv,ui.t()) + eq_fields_->lame_lambda(p_low)*arma::dot(ui,nv)*arma::eye(3,3))
228  )*fe_values_sub_.JxW(k);
229  }
230 
231  }
232  }
233  k++;
234  }
235 
236  for (unsigned int n=0; n<2; ++n)
237  for (unsigned int m=0; m<2; ++m)
239  }
240 
241 
242 
243 private:
244  inline arma::mat33 mat_t(const arma::mat33 &m, const arma::vec3 &n)
245  {
246  arma::mat33 mt = m - m*arma::kron(n,n.t());
247  return mt;
248  }
249 
250 
251 
252  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
253  shared_ptr<FiniteElement<dim-1>> fe_low_; ///< Finite element for the solution of the advection-diffusion equation (dim-1).
254 
255  /// Data objects shared with Elasticity
258 
259  /// Sub field set contains fields used in calculation.
261 
262  unsigned int n_dofs_; ///< Number of dofs
263  unsigned int n_dofs_sub_; ///< Number of dofs (on lower dim element)
264  std::vector<unsigned int> n_dofs_ngh_; ///< Number of dofs on lower and higher dimension element (vector of 2 items)
265  FEValues<3> fe_values_; ///< FEValues of cell object (FESystem of P disc finite element type)
266  FEValues<3> fe_values_side_; ///< FEValues of side object
267  FEValues<3> fe_values_sub_; ///< FEValues of lower dimension cell object
268 
269  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
270  vector<vector<LongIdx> > side_dof_indices_; ///< 2 items vector of DOF indices in neighbour calculation.
271  vector<PetscScalar> local_matrix_; ///< Auxiliary vector for assemble methods
272  vector<vector<vector<PetscScalar>>> local_matrix_ngh_; ///< Auxiliary vectors for assemble ngh integral
273  const FEValuesViews::Vector<3> * vec_view_; ///< Vector view in cell integral calculation.
274  const FEValuesViews::Vector<3> * vec_view_side_; ///< Vector view in boundary / neighbour calculation.
275  const FEValuesViews::Vector<3> * vec_view_sub_; ///< Vector view of low dim element in neighbour calculation.
276 
277  template < template<IntDim...> class DimAssembly>
278  friend class GenericAssembly;
279 
280 };
281 
282 
283 template <unsigned int dim>
285 {
286 public:
287  typedef typename Elasticity::EqFields EqFields;
288  typedef typename Elasticity::EqData EqData;
289 
290  static constexpr const char * name() { return "RhsAssemblyElasticity"; }
291 
292  /// Constructor.
293  RhsAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
294  : AssemblyBase<dim>(1), eq_fields_(eq_fields), eq_data_(eq_data) {
297  this->used_fields_ += eq_fields_->load;
302  this->used_fields_ += eq_fields_->bc_type;
307  }
308 
309  /// Destructor.
311 
312  /// Initialize auxiliary vectors and other data members
313  void initialize(ElementCacheMap *element_cache_map) {
314  //this->balance_ = eq_data_->balance_;
315  this->element_cache_map_ = element_cache_map;
316 
317  shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
318  shared_ptr<FE_P<dim-1>> fe_p_low = std::make_shared< FE_P<dim-1> >(1);
319  fe_ = std::make_shared<FESystem<dim>>(fe_p, FEVector, 3);
320  fe_low_ = std::make_shared<FESystem<dim-1>>(fe_p_low, FEVector, 3);
321  fe_values_.initialize(*this->quad_, *fe_,
325  fe_values_side_.initialize(*this->quad_low_, *fe_,
327  fe_values_sub_.initialize(*this->quad_low_, *fe_low_,
329  n_dofs_ = fe_->n_dofs();
330  n_dofs_sub_ = fe_low_->n_dofs();
332  dof_indices_.resize(n_dofs_);
333  side_dof_indices_.resize(2);
334  side_dof_indices_[0].resize(n_dofs_sub_); // index 0 = element with lower dimension,
335  side_dof_indices_[1].resize(n_dofs_); // index 1 = side of element with higher dimension
336  local_rhs_.resize(n_dofs_);
337  local_rhs_ngh_.resize(2);
338  for (uint n=0; n<2; ++n) local_rhs_ngh_[n].resize(n_dofs_);
342  if (dim>1) vec_view_sub_ = &fe_values_sub_.vector_view(0);
343  }
344 
345 
346  /// Assemble integral over element
347  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
348  {
349  if (cell.dim() != dim) return;
350  if (!cell.is_own()) return;
351 
352  ElementAccessor<3> elm_acc = cell.elm();
353 
354  fe_values_.reinit(elm_acc);
356 
357  // assemble the local stiffness matrix
358  fill_n(&(local_rhs_[0]), n_dofs_, 0);
359  //local_source_balance_vector.assign(n_dofs_, 0);
360  //local_source_balance_rhs.assign(n_dofs_, 0);
361 
362  // compute sources
363  unsigned int k=0;
364  for (auto p : this->bulk_points(element_patch_idx) )
365  {
366  for (unsigned int i=0; i<n_dofs_; i++)
367  local_rhs_[i] += (
368  arma::dot(eq_fields_->load(p), vec_view_->value(i,k))
370  -arma::dot(eq_fields_->initial_stress(p), vec_view_->grad(i,k))
372  ++k;
373  }
375 
376 // for (unsigned int i=0; i<n_dofs_; i++)
377 // {
378 // for (unsigned int k=0; k<qsize_; k++) // point range
379 // local_source_balance_vector[i] -= 0;//sources_sigma[k]*fe_values_[vec_view_].value(i,k)*fe_values_.JxW(k);
380 //
381 // local_source_balance_rhs[i] += local_rhs_[i];
382 // }
383 // balance_->add_source_matrix_values(subst_idx, elm_acc.region().bulk_idx(), dof_indices_, local_source_balance_vector);
384 // balance_->add_source_vec_values(subst_idx, elm_acc.region().bulk_idx(), dof_indices_, local_source_balance_rhs);
385  }
386 
387  /// Assembles boundary integral.
388  inline void boundary_side_integral(DHCellSide cell_side)
389  {
390  ASSERT_EQ(cell_side.dim(), dim).error("Dimension of element mismatch!");
391  if (!cell_side.cell().is_own()) return;
392 
393  Side side = cell_side.side();
394  const DHCellAccessor &dh_cell = cell_side.cell();
395  dh_cell.get_dof_indices(dof_indices_);
397 
398  auto p_side = *( this->boundary_points(cell_side).begin() );
399  auto p_bdr = p_side.point_bdr( side.cond().element_accessor() );
400  unsigned int bc_type = eq_fields_->bc_type(p_bdr);
401 
402  fill_n(&(local_rhs_[0]), n_dofs_, 0);
403  // local_flux_balance_vector.assign(n_dofs_, 0);
404  // local_flux_balance_rhs = 0;
405 
406  unsigned int k = 0;
407 
408  // addtion from initial stress
409  for (auto p : this->boundary_points(cell_side) )
410  {
411  for (unsigned int i=0; i<n_dofs_; i++)
414  vec_view_bdr_->value(i,k)) *
416  ++k;
417  }
418 
419  k = 0;
420  if (bc_type == EqFields::bc_type_displacement)
421  {
422  double side_measure = cell_side.measure();
423  for (auto p : this->boundary_points(cell_side) )
424  {
425  auto p_bdr = p.point_bdr( side.cond().element_accessor() );
426  for (unsigned int i=0; i<n_dofs_; i++)
427  local_rhs_[i] += (eq_fields_->dirichlet_penalty(p) / side_measure) *
428  arma::dot(eq_fields_->bc_displacement(p_bdr), vec_view_bdr_->value(i,k)) *
430  ++k;
431  }
432  }
433  else if (bc_type == EqFields::bc_type_displacement_normal)
434  {
435  double side_measure = cell_side.measure();
436  for (auto p : this->boundary_points(cell_side) )
437  {
438  auto p_bdr = p.point_bdr( side.cond().element_accessor() );
439  for (unsigned int i=0; i<n_dofs_; i++)
440  local_rhs_[i] += (eq_fields_->dirichlet_penalty(p) / side_measure) *
444  ++k;
445  }
446  }
447  else if (bc_type == EqFields::bc_type_traction)
448  {
449  for (auto p : this->boundary_points(cell_side) )
450  {
451  auto p_bdr = p.point_bdr( side.cond().element_accessor() );
452  for (unsigned int i=0; i<n_dofs_; i++)
456  ++k;
457  }
458  }
459  else if (bc_type == EqFields::bc_type_stress)
460  {
461  for (auto p : this->boundary_points(cell_side) )
462  {
463  auto p_bdr = p.point_bdr( side.cond().element_accessor() );
464  for (unsigned int i=0; i<n_dofs_; i++)
465  // stress is multiplied by inward normal to obtain traction
470  ++k;
471  }
472  }
474 
475 
476 // balance_->add_flux_matrix_values(subst_idx, loc_b, side_dof_indices, local_flux_balance_vector);
477 // balance_->add_flux_vec_value(subst_idx, loc_b, local_flux_balance_rhs);
478  // ++loc_b;
479  }
480 
481 
482  /// Assembles between elements of different dimensions.
483  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
484  if (dim == 1) return;
485  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
486 
487  cell_lower_dim.get_dof_indices(side_dof_indices_[0]);
488  ElementAccessor<3> cell_sub = cell_lower_dim.elm();
489  fe_values_sub_.reinit(cell_sub);
490 
491  DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element( neighb_side.element().idx() );
492  cell_higher_dim.get_dof_indices(side_dof_indices_[1]);
493  fe_values_side_.reinit(neighb_side.side());
494 
495  // Element id's for testing if they belong to local partition.
496  bool own_element_id[2];
497  own_element_id[0] = cell_lower_dim.is_own();
498  own_element_id[1] = cell_higher_dim.is_own();
499 
500  for (unsigned int n=0; n<2; ++n)
501  for (unsigned int i=0; i<n_dofs_; i++)
502  local_rhs_ngh_[n][i] = 0;
503 
504  // set transmission conditions
505  unsigned int k=0;
506  for (auto p_high : this->coupling_points(neighb_side) )
507  {
508  auto p_low = p_high.lower_dim(cell_lower_dim);
510 
511  for (int n=0; n<2; n++)
512  {
513  if (!own_element_id[n]) continue;
514 
515  for (unsigned int i=0; i<n_dofs_ngh_[n]; i++)
516  {
517  arma::vec3 vi = (n==0) ? arma::zeros(3) : vec_view_side_->value(i,k);
518  arma::vec3 vf = (n==1) ? arma::zeros(3) : vec_view_sub_->value(i,k);
519 
520  local_rhs_ngh_[n][i] -= eq_fields_->fracture_sigma(p_low) * eq_fields_->cross_section(p_high) *
521  arma::dot(vf-vi, eq_fields_->potential_load(p_high) * nv) * fe_values_sub_.JxW(k);
522  }
523  }
524  ++k;
525  }
526 
527  for (unsigned int n=0; n<2; ++n)
529  }
530 
531 
532 
533 private:
534  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
535  shared_ptr<FiniteElement<dim-1>> fe_low_; ///< Finite element for the solution of the advection-diffusion equation (dim-1).
536 
537  /// Data objects shared with Elasticity
540 
541  /// Sub field set contains fields used in calculation.
543 
544  unsigned int n_dofs_; ///< Number of dofs
545  unsigned int n_dofs_sub_; ///< Number of dofs (on lower dim element)
546  std::vector<unsigned int> n_dofs_ngh_; ///< Number of dofs on lower and higher dimension element (vector of 2 items)
547  FEValues<3> fe_values_; ///< FEValues of cell object (FESystem of P disc finite element type)
548  FEValues<3> fe_values_bdr_side_; ///< FEValues of side (boundary integral) object
549  FEValues<3> fe_values_side_; ///< FEValues of side (neighbour integral) object
550  FEValues<3> fe_values_sub_; ///< FEValues of lower dimension cell object
551 
552  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
553  vector<vector<LongIdx> > side_dof_indices_; ///< 2 items vector of DOF indices in neighbour calculation.
554  vector<PetscScalar> local_rhs_; ///< Auxiliary vector for assemble methods
555  vector<vector<PetscScalar>> local_rhs_ngh_; ///< Auxiliary vectors for assemble ngh integral
556  const FEValuesViews::Vector<3> * vec_view_; ///< Vector view in cell integral calculation.
557  const FEValuesViews::Vector<3> * vec_view_bdr_; ///< Vector view in boundary calculation.
558  const FEValuesViews::Vector<3> * vec_view_side_; ///< Vector view in neighbour calculation.
559  const FEValuesViews::Vector<3> * vec_view_sub_; ///< Vector view of low dim element in neighbour calculation.
560 
561 
562  template < template<IntDim...> class DimAssembly>
563  friend class GenericAssembly;
564 
565 };
566 
567 template <unsigned int dim>
569 {
570 public:
571  typedef typename Elasticity::EqFields EqFields;
572  typedef typename Elasticity::EqData EqData;
573 
574  static constexpr const char * name() { return "OutpuFieldsAssemblyElasticity"; }
575 
576  /// Constructor.
578  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
581  this->used_fields_ += eq_fields_->lame_mu;
584  }
585 
586  /// Destructor.
588 
589  /// Initialize auxiliary vectors and other data members
590  void initialize(ElementCacheMap *element_cache_map) {
591  //this->balance_ = eq_data_->balance_;
592  this->element_cache_map_ = element_cache_map;
593 
594  shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
595  fe_ = std::make_shared<FESystem<dim>>(fe_p, FEVector, 3);
596  fv_.initialize(*this->quad_, *fe_,
598  fsv_.initialize(*this->quad_low_, *fe_,
600  n_dofs_ = fe_->n_dofs();
601  vec_view_ = &fv_.vector_view(0);
602  // if (dim>1) ??
604 
611  }
612 
613 
614  /// Assemble integral over element
615  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
616  {
617  if (cell.dim() != dim) return;
618  if (!cell.is_own()) return;
619  DHCellAccessor cell_tensor = cell.cell_with_other_dh(eq_data_->dh_tensor_.get());
620  DHCellAccessor cell_scalar = cell.cell_with_other_dh(eq_data_->dh_scalar_.get());
621 
622  auto elm = cell.elm();
623 
624  fv_.reinit(elm);
628 
629  auto p = *( this->bulk_points(element_patch_idx).begin() );
630 
632  double div = 0;
633  for (unsigned int i=0; i<n_dofs_; i++)
634  {
635  stress += (2*eq_fields_->lame_mu(p)*vec_view_->sym_grad(i,0) + eq_fields_->lame_lambda(p)*vec_view_->divergence(i,0)*arma::eye(3,3))
638  }
639 
640  arma::mat33 stress_dev = stress - arma::trace(stress)/3*arma::eye(3,3);
641  double von_mises_stress = sqrt(1.5*arma::dot(stress_dev, stress_dev));
642  double mean_stress = arma::trace(stress) / 3;
644 
645  for (unsigned int i=0; i<3; i++)
646  for (unsigned int j=0; j<3; j++)
647  output_stress_vec_.add( dof_indices_tensor_[i*3+j], stress(i,j) );
648  output_von_mises_stress_vec_.set( dof_indices_scalar_[0], von_mises_stress );
650 
652  }
653 
654 
655  /// Assembles between elements of different dimensions.
656  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
657  if (dim == 1) return;
658  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
659 
661  normal_stress_.zeros();
662 
663  DHCellAccessor cell_higher_dim = neighb_side.cell();
664  DHCellAccessor cell_tensor = cell_lower_dim.cell_with_other_dh(eq_data_->dh_tensor_.get());
665  DHCellAccessor cell_scalar = cell_lower_dim.cell_with_other_dh(eq_data_->dh_scalar_.get());
666  fsv_.reinit(neighb_side.side());
667 
668  dof_indices_ = cell_higher_dim.get_loc_dof_indices();
669  auto p_high = *( this->coupling_points(neighb_side).begin() );
670  auto p_low = p_high.lower_dim(cell_lower_dim);
671 
672  for (unsigned int i=0; i<n_dofs_; i++)
673  {
675  arma::mat33 grad = -arma::kron(vec_view_side_->value(i,0)*output_vec_.get(dof_indices_[i]), fsv_.normal_vector(0).t()) / eq_fields_->cross_section(p_low);
676  normal_stress_ += eq_fields_->lame_mu(p_low)*(grad+grad.t()) + eq_fields_->lame_lambda(p_low)*arma::trace(grad)*arma::eye(3,3);
677  }
678 
681  for (unsigned int i=0; i<3; i++)
682  for (unsigned int j=0; j<3; j++)
686  }
687 
688 
689 
690 private:
691  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
692 
693  /// Data objects shared with Elasticity
696 
697  /// Sub field set contains fields used in calculation.
699 
700  unsigned int n_dofs_; ///< Number of dofs
701  FEValues<3> fv_; ///< FEValues of cell object (FESystem of P disc finite element type)
702  FEValues<3> fsv_; ///< FEValues of side (neighbour integral) object
703 
704  LocDofVec dof_indices_; ///< Vector of local DOF indices of vector fields
705  LocDofVec dof_indices_scalar_; ///< Vector of local DOF indices of scalar fields
706  LocDofVec dof_indices_tensor_; ///< Vector of local DOF indices of tensor fields
707  const FEValuesViews::Vector<3> * vec_view_; ///< Vector view in cell integral calculation.
708  const FEValuesViews::Vector<3> * vec_view_side_; ///< Vector view in neighbour calculation.
709 
710  double normal_displacement_; ///< Holds constributions of normal displacement.
711  arma::mat33 normal_stress_; ///< Holds constributions of normal stress.
712 
713  /// Data vectors of output fields (FieldFE).
720 
721  template < template<IntDim...> class DimAssembly>
722  friend class GenericAssembly;
723 
724 };
725 
726 
727 
728 /**
729  * Container class for assembly of constraint matrix for contact condition.
730  */
731 template <unsigned int dim>
733 {
734 public:
735  typedef typename Elasticity::EqFields EqFields;
736  typedef typename Elasticity::EqData EqData;
737 
738  static constexpr const char * name() { return "ConstraintAssemblyElasticity"; }
739 
740  /// Constructor.
742  : AssemblyBase<dim>(1), eq_fields_(eq_fields), eq_data_(eq_data) {
746  }
747 
748  /// Destructor.
750 
751  /// Initialize auxiliary vectors and other data members
752  void initialize(ElementCacheMap *element_cache_map) {
753  this->element_cache_map_ = element_cache_map;
754 
755  shared_ptr<FE_P<dim>> fe_p = std::make_shared< FE_P<dim> >(1);
756  fe_ = std::make_shared<FESystem<dim>>(fe_p, FEVector, 3);
757  fe_values_side_.initialize(*this->quad_low_, *fe_,
759 
760  n_dofs_ = fe_->n_dofs();
761  dof_indices_.resize(n_dofs_);
762  local_matrix_.resize(n_dofs_*n_dofs_);
764  }
765 
766 
767  /// Assembles between elements of different dimensions.
768  inline void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side) {
769  if (dim == 1) return;
770  if (!cell_lower_dim.is_own()) return;
771 
772  ASSERT_EQ(cell_lower_dim.dim(), dim-1).error("Dimension of element mismatch!");
773 
774  DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element( neighb_side.element().idx() );
775  cell_higher_dim.get_dof_indices(dof_indices_);
776 
777  fe_values_side_.reinit(neighb_side.side());
778 
779  for (unsigned int i=0; i<n_dofs_; i++)
780  local_matrix_[i] = 0;
781 
782  // Assemble matrix and vector for contact conditions in the form B*x <= c,
783  // where B*x is the average jump of normal displacements and c is the average cross-section on element.
784  // Positive value means that the fracture closes.
785  unsigned int k=0;
786  double local_vector = 0;
787  for (auto p_high : this->coupling_points(neighb_side) )
788  {
789  auto p_low = p_high.lower_dim(cell_lower_dim);
791 
792  local_vector += (eq_fields_->cross_section(p_low) - eq_fields_->cross_section_min(p_low))*fe_values_side_.JxW(k) / cell_lower_dim.elm().measure() / cell_lower_dim.elm()->n_neighs_vb();
793 
794  for (unsigned int i=0; i<n_dofs_; i++)
795  {
796  local_matrix_[i] += eq_fields_->cross_section(p_high)*arma::dot(vec_view_side_->value(i,k), nv)*fe_values_side_.JxW(k) / cell_lower_dim.elm().measure();
797  }
798  k++;
799  }
800 
801  int arow[1] = { eq_data_->constraint_idx[cell_lower_dim.elm_idx()] };
802  MatSetValues(eq_data_->constraint_matrix, 1, arow, n_dofs_, dof_indices_.data(), &(local_matrix_[0]), ADD_VALUES);
803  VecSetValue(eq_data_->constraint_vec, arow[0], local_vector, ADD_VALUES);
804  }
805 
806 
807 
808 private:
809 
810 
811 
812  shared_ptr<FiniteElement<dim>> fe_; ///< Finite element for the solution of the advection-diffusion equation.
813 
814  /// Data objects shared with Elasticity
817 
818  /// Sub field set contains fields used in calculation.
820 
821  unsigned int n_dofs_; ///< Number of dofs
822  FEValues<3> fe_values_side_; ///< FEValues of side object
823 
824  vector<LongIdx> dof_indices_; ///< Vector of global DOF indices
825  vector<vector<LongIdx> > side_dof_indices_; ///< 2 items vector of DOF indices in neighbour calculation.
826  vector<PetscScalar> local_matrix_; ///< Auxiliary vector for assemble methods
827  const FEValuesViews::Vector<3> * vec_view_side_; ///< Vector view in boundary / neighbour calculation.
828 
829  template < template<IntDim...> class DimAssembly>
830  friend class GenericAssembly;
831 
832 };
833 
834 
835 #endif /* ASSEMBLY_ELASTICITY_HH_ */
836 
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
Quadrature * quad_
Quadrature used in assembling methods.
Quadrature * quad_low_
Quadrature used in assembling methods (dim-1).
Range< CouplingPoint > coupling_points(const DHCellSide &cell_side) const
Return CouplingPoint range of appropriate dimension.
int active_integrals_
Holds mask of active integrals.
Range< BoundaryPoint > boundary_points(const DHCellSide &cell_side) const
Return BoundaryPoint range of appropriate dimension.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
ElementAccessor< 3 > element_accessor()
FEValues< 3 > fe_values_side_
FEValues of side object.
EqFields * eq_fields_
Data objects shared with Elasticity.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
static constexpr const char * name()
ConstraintAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
unsigned int n_dofs_
Number of dofs.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in boundary / neighbour calculation.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
FieldSet used_fields_
Sub field set contains fields used in calculation.
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
bool is_own() const
Return true if accessor represents own element (false for ghost element)
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
DHCellAccessor cell_with_other_dh(const DOFHandlerMultiDim *dh) const
Create new accessor with same local idx and given DOF handler. Actual and given DOF handler must be c...
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
unsigned int dim() const
Return dimension of element appropriate to cell.
unsigned int elm_idx() const
Return serial idx to element of loc_ele_idx_.
Side accessor allows to iterate over sides of DOF handler cell.
Side side() const
Return Side of given cell and side_idx.
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
unsigned int dim() const
Return dimension of element appropriate to the side.
double measure() const
ElementAccessor< 3 > element() const
LinSys * ls
Linear algebraic system.
Definition: elasticity.hh:135
std::shared_ptr< DOFHandlerMultiDim > dh_scalar_
Definition: elasticity.hh:128
std::shared_ptr< DOFHandlerMultiDim > dh_tensor_
Definition: elasticity.hh:129
std::shared_ptr< DOFHandlerMultiDim > dh_
Objects for distribution of dofs.
Definition: elasticity.hh:127
std::map< LongIdx, LongIdx > constraint_idx
Definition: elasticity.hh:141
Field< 3, FieldValue< 3 >::Scalar > potential_load
Potential of an additional (external) load.
Definition: elasticity.hh:79
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::TensorFixed > > output_stress_ptr
Definition: elasticity.hh:101
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_von_mises_stress_ptr
Definition: elasticity.hh:102
Field< 3, FieldValue< 3 >::Scalar > cross_section_min
Definition: elasticity.hh:78
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_cross_section_ptr
Definition: elasticity.hh:104
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::VectorFixed > > output_field_ptr
Definition: elasticity.hh:100
BCField< 3, FieldValue< 3 >::Enum > bc_type
Definition: elasticity.hh:66
Field< 3, FieldValue< 3 >::Scalar > lame_lambda
Definition: elasticity.hh:95
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_div_ptr
Definition: elasticity.hh:105
BCField< 3, FieldValue< 3 >::VectorFixed > bc_displacement
Definition: elasticity.hh:67
Field< 3, FieldValue< 3 >::TensorFixed > initial_stress
Definition: elasticity.hh:74
BCField< 3, FieldValue< 3 >::VectorFixed > bc_traction
Definition: elasticity.hh:68
Field< 3, FieldValue< 3 >::Scalar > ref_potential_load
Potential of reference external load on boundary. TODO: Switch to BCField when possible.
Definition: elasticity.hh:80
Field< 3, FieldValue< 3 >::VectorFixed > load
Definition: elasticity.hh:70
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > output_mean_stress_ptr
Definition: elasticity.hh:103
Field< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures.
Definition: elasticity.hh:73
Field< 3, FieldValue< 3 >::Scalar > dirichlet_penalty
Definition: elasticity.hh:96
Field< 3, FieldValue< 3 >::Scalar > cross_section
Pointer to DarcyFlow field cross_section.
Definition: elasticity.hh:77
BCField< 3, FieldValue< 3 >::TensorFixed > bc_stress
Definition: elasticity.hh:69
Field< 3, FieldValue< 3 >::Scalar > lame_mu
Definition: elasticity.hh:94
double measure() const
Computes the measure of the element.
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
Directing class of FieldValueCache.
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
arma::mat::fixed< spacedim, spacedim > grad(unsigned int function_no, unsigned int point_no) const
Return gradient of vector-valued shape function.
arma::mat::fixed< spacedim, spacedim > sym_grad(unsigned int function_no, unsigned int point_no) const
Return symmetric gradient of vector-valued shape function.
arma::vec::fixed< spacedim > value(unsigned int function_no, unsigned int point_no) const
Return value of vector-valued shape function.
double divergence(unsigned int function_no, unsigned int point_no) const
Return divergence of vector-valued shape function.
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point.
Definition: fe_values.hh:223
const FEValuesViews::Vector< spacedim > & vector_view(unsigned int i) const
Accessor to vector values of multicomponent FE.
Definition: fe_values.hh:277
void initialize(Quadrature &_quadrature, FiniteElement< DIM > &_fe, UpdateFlags _flags)
Initialize structures and calculates cell-independent data.
Definition: fe_values.cc:137
void reinit(const ElementAccessor< spacedim > &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:537
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
Definition: fe_values.hh:257
Conforming Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:86
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Generic class of assemblation.
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
LocDofVec dof_indices_scalar_
Vector of local DOF indices of scalar fields.
VectorMPI output_vec_
Data vectors of output fields (FieldFE).
unsigned int n_dofs_
Number of dofs.
FieldSet used_fields_
Sub field set contains fields used in calculation.
FEValues< 3 > fv_
FEValues of cell object (FESystem of P disc finite element type)
static constexpr const char * name()
LocDofVec dof_indices_
Vector of local DOF indices of vector fields.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
LocDofVec dof_indices_tensor_
Vector of local DOF indices of tensor fields.
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in neighbour calculation.
EqFields * eq_fields_
Data objects shared with Elasticity.
arma::mat33 normal_stress_
Holds constributions of normal stress.
FEValues< 3 > fsv_
FEValues of side (neighbour integral) object.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
double normal_displacement_
Holds constributions of normal displacement.
OutpuFieldsAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
FEValues< 3 > fe_values_sub_
FEValues of lower dimension cell object.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
const FEValuesViews::Vector< 3 > * vec_view_bdr_
Vector view in boundary calculation.
vector< vector< PetscScalar > > local_rhs_ngh_
Auxiliary vectors for assemble ngh integral.
FEValues< 3 > fe_values_bdr_side_
FEValues of side (boundary integral) object.
const FEValuesViews::Vector< 3 > * vec_view_sub_
Vector view of low dim element in neighbour calculation.
EqFields * eq_fields_
Data objects shared with Elasticity.
Elasticity::EqFields EqFields
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
FEValues< 3 > fe_values_side_
FEValues of side (neighbour integral) object.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
std::vector< unsigned int > n_dofs_ngh_
Number of dofs on lower and higher dimension element (vector of 2 items)
Elasticity::EqData EqData
FEValues< 3 > fe_values_
FEValues of cell object (FESystem of P disc finite element type)
RhsAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
unsigned int n_dofs_sub_
Number of dofs (on lower dim element)
vector< PetscScalar > local_rhs_
Auxiliary vector for assemble methods.
unsigned int n_dofs_
Number of dofs.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
static constexpr const char * name()
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in neighbour calculation.
void boundary_side_integral(DHCellSide cell_side)
Assembles boundary integral.
FieldSet used_fields_
Sub field set contains fields used in calculation.
Boundary cond() const
FEValues< 3 > fe_values_side_
FEValues of side object.
void dimjoin_intergral(DHCellAccessor cell_lower_dim, DHCellSide neighb_side)
Assembles between elements of different dimensions.
arma::mat33 mat_t(const arma::mat33 &m, const arma::vec3 &n)
FieldSet used_fields_
Sub field set contains fields used in calculation.
const FEValuesViews::Vector< 3 > * vec_view_
Vector view in cell integral calculation.
shared_ptr< FiniteElement< dim > > fe_
Finite element for the solution of the advection-diffusion equation.
static constexpr const char * name()
const FEValuesViews::Vector< 3 > * vec_view_side_
Vector view in boundary / neighbour calculation.
shared_ptr< FiniteElement< dim-1 > > fe_low_
Finite element for the solution of the advection-diffusion equation (dim-1).
vector< PetscScalar > local_matrix_
Auxiliary vector for assemble methods.
vector< vector< LongIdx > > side_dof_indices_
2 items vector of DOF indices in neighbour calculation.
vector< LongIdx > dof_indices_
Vector of global DOF indices.
EqFields * eq_fields_
Data objects shared with Elasticity.
FEValues< 3 > fe_values_sub_
FEValues of lower dimension cell object.
const FEValuesViews::Vector< 3 > * vec_view_sub_
Vector view of low dim element in neighbour calculation.
unsigned int n_dofs_sub_
Number of dofs (on lower dim element)
unsigned int n_dofs_
Number of dofs.
FEValues< 3 > fe_values_
FEValues of cell object (FESystem of P disc finite element type)
vector< vector< vector< PetscScalar > > > local_matrix_ngh_
Auxiliary vectors for assemble ngh integral.
std::vector< unsigned int > n_dofs_ngh_
Number of dofs on lower and higher dimension element (vector of 2 items)
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
StiffnessAssemblyElasticity(EqFields *eq_fields, EqData *eq_data)
Constructor.
void boundary_side_integral(DHCellSide cell_side)
Assembles boundary integral.
double get(unsigned int pos) const
Return value on given position.
Definition: vector_mpi.hh:105
void set(unsigned int pos, double val)
Set value on given position.
Definition: vector_mpi.hh:111
void add(unsigned int pos, double val)
Add value to item on given position.
Definition: vector_mpi.hh:125
FEM for linear elasticity.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ FEVector
@ coupling
@ boundary
@ bulk
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Definitions of particular quadrature rules on simplices.
@ update_values
Shape function values.
Definition: update_flags.hh:87
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_side_JxW_values
Transformed quadrature weight for cell sides.
@ update_gradients
Shape function gradients.
Definition: update_flags.hh:95
@ update_quadrature_points
Transformed quadrature points.