pi-DoMUS: Parallel Deal.II MUltiphysics Solver
neo_hookean_two_fields.h
Go to the documentation of this file.
1 
5 #ifndef _neo_hookean_two_fields_interface_h_
6 #define _neo_hookean_two_fields_interface_h_
7 
9 #include <deal2lkit/parsed_function.h>
10 
11 
12 #include <deal.II/fe/fe_values.h>
17 
21 
22 #include <deal.II/lac/solver_cg.h>
24 
25 
26 template <int dim, int spacedim>
27 class NeoHookeanTwoFieldsInterface : public ConservativeInterface<dim,spacedim,dim+1, NeoHookeanTwoFieldsInterface<dim,spacedim> >
28 {
29 public:
30 
31  typedef FEValuesCache<dim,spacedim> Scratch;
32  typedef Assembly::CopyData::piDoMUSPreconditioner<dim,dim> CopyPreconditioner;
33  typedef Assembly::CopyData::piDoMUSSystem<dim,dim> CopySystem;
36 
37  /* specific and useful functions for this very problem */
39 
40  virtual void declare_parameters (ParameterHandler &prm);
41  virtual void parse_parameters_call_back ();
42 
43  /* these functions MUST have the follwowing names
44  * because they are called by the ConservativeInterface class
45  */
46  template<typename Number>
48  Scratch &,
49  CopyPreconditioner &,
50  Number &energy) const;
51 
52  template<typename Number>
54  Scratch &,
55  CopySystem &,
56  Number &energy) const;
57 
61  const std::vector<shared_ptr<MAT> >,
63  BlockLinearOperator<VEC> &) const;
64 
65 private:
66  double E;
67  double nu;
68  double mu;
69  double lambda;
70 
71  mutable shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
72  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> Mp_preconditioner;
73  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> T_preconditioner;
74 
75 };
76 
77 template <int dim, int spacedim>
79  ConservativeInterface<dim,spacedim,dim+1,NeoHookeanTwoFieldsInterface<dim,spacedim> >("NeoHookean Interface",
80  "FESystem[FE_Q(2)^d-FE_DGP(1)]",
81  "u,u,u,p", "1,1; 1,0", "1,0; 0,1",
82  "1,0")
83 {};
84 
85 
86 
87 template <int dim, int spacedim>
88 template<typename Number>
90  Scratch &fe_cache,
91  CopyPreconditioner &data,
92  Number &energy) const
93 {
94  Number alpha = this->alpha;
95  this->reinit(alpha, cell, fe_cache);
96 
97  auto &JxW = fe_cache.get_JxW_values();
98 
99  const FEValuesExtractors::Vector displacement(0);
100  const FEValuesExtractors::Scalar pressure(dim);
101  auto &ps = fe_cache.get_values("solution","p", pressure, alpha);
102  auto &grad_us = fe_cache.get_gradients("solution","grad_u",displacement, alpha);
103 
104  const unsigned int n_q_points = ps.size();
105 
106  energy = 0;
107  for (unsigned int q=0; q<n_q_points; ++q)
108  {
109  const Number &p = ps[q];
110  const Tensor <2, dim, Number> &grad_u = grad_us[q];
111 
112  energy += (scalar_product(grad_u,grad_u) +
113  0.5*p*p)*JxW[q];
114  }
115 }
116 
117 template <int dim, int spacedim>
118 template<typename Number>
120  Scratch &fe_cache,
121  CopySystem &data,
122  Number &energy) const
123 {
124  Number alpha = this->alpha;
125  this->reinit(alpha, cell, fe_cache);
126 
127  const FEValuesExtractors::Vector displacement(0);
128  auto &us = fe_cache.get_values("solution", "u", displacement, alpha);
129  auto &us_dot = fe_cache.get_values("solution_dot", "u_dot", displacement, alpha);
130  auto &Fs = fe_cache.get_deformation_gradients("solution", "Fu", displacement, alpha);
131 
132  const FEValuesExtractors::Scalar pressure(dim);
133  auto &ps = fe_cache.get_values("solution","p", pressure, alpha);
134 
135  auto &JxW = fe_cache.get_JxW_values();
136 
137  const unsigned int n_q_points = ps.size();
138 
139  energy = 0;
140  for (unsigned int q=0; q<n_q_points; ++q)
141  {
143 
144  const Tensor <1, dim, Number> &u = us[q];
145  const Tensor <1, dim, Number> &u_dot = us_dot[q];
146  const Number &p = ps[q];
147  const Tensor <2, dim, Number> &F = Fs[q];
148  const Tensor<2, dim, Number> C = transpose(F)*F;
149 
150  Number Ic = trace(C);
151  Number J = determinant(F);
152 
153  Number psi = (mu/2.)*(Ic-dim) +p*(J-1.);
154  energy += (u*u_dot + psi)*JxW[q];
155  }
156 }
157 
158 
159 template <int dim, int spacedim>
161 {
163  this->add_parameter(prm, &E, "Young's modulus", "10.0", Patterns::Double(0.0));
164  this->add_parameter(prm, &nu, "Poisson's ratio", "0.3", Patterns::Double(0.0));
165 }
166 
167 template <int dim, int spacedim>
169 {
171  mu = E/(2.0*(1.+nu));
172  lambda = (E *nu)/((1.+nu)*(1.-2.*nu));
173 }
174 
175 template <int dim, int spacedim>
176 void
179  const TrilinosWrappers::BlockSparseMatrix &preconditioner_matrix,
180  const std::vector<shared_ptr<MAT> >,
181  BlockLinearOperator<VEC> &system_op,
182  BlockLinearOperator<VEC> &prec_op) const
183 {
184 
185  std::vector<std::vector<bool> > constant_modes;
186  FEValuesExtractors::Vector displacement(0);
187  DoFTools::extract_constant_modes (dh, dh.get_fe().component_mask(displacement),
188  constant_modes);
189 
190  Mp_preconditioner.reset (new TrilinosWrappers::PreconditionJacobi());
191  Amg_preconditioner.reset (new TrilinosWrappers::PreconditionAMG());
192 
194  Amg_data.constant_modes = constant_modes;
195  Amg_data.elliptic = true;
196  Amg_data.higher_order_elements = true;
197  Amg_data.smoother_sweeps = 2;
198  Amg_data.aggregation_threshold = 0.02;
199 
200  Mp_preconditioner->initialize (preconditioner_matrix.block(1,1));
201  Amg_preconditioner->initialize (preconditioner_matrix.block(0,0),
202  Amg_data);
203 
204 
205  // SYSTEM MATRIX:
206  auto A = linear_operator< TrilinosWrappers::MPI::Vector >( matrix.block(0,0) );
207  auto Bt = linear_operator< TrilinosWrappers::MPI::Vector >( matrix.block(0,1) );
208  // auto B = transpose_operator(Bt);
209  auto B = linear_operator< TrilinosWrappers::MPI::Vector >( matrix.block(1,0) );
210  auto ZeroP = 0*linear_operator< TrilinosWrappers::MPI::Vector >( matrix.block(1,1) );
211 
212  auto Mp = linear_operator< TrilinosWrappers::MPI::Vector >( preconditioner_matrix.block(1,1) );
213 
214  static ReductionControl solver_control_pre(5000, 1e-8);
215  static SolverCG<TrilinosWrappers::MPI::Vector> solver_CG(solver_control_pre);
216  auto A_inv = inverse_operator( A, solver_CG, *Amg_preconditioner);
217  auto Schur_inv = inverse_operator( Mp, solver_CG, *Mp_preconditioner);
218 
219  auto P00 = A_inv;
220  auto P01 = null_operator(Bt);
221  auto P10 = Schur_inv * B * A_inv;
222  auto P11 = -1 * Schur_inv;
223 
224  // ASSEMBLE THE PROBLEM:
225  system_op = block_operator<2, 2, VEC >({{
226  {{ A, Bt }} ,
227  {{ B, ZeroP }}
228  }
229  });
230 
231 
232  //const auto S = linear_operator<VEC>(matrix);
233 
234  prec_op = block_operator<2, 2, VEC >({{
235  {{ P00, P01 }} ,
236  {{ P10, P11 }}
237  }
238  });
239 }
240 
241 
243 
244 #endif
245 
Definition: conservative.h:20
std::vector< std::vector< bool > > constant_modes
ActiveSelector::active_cell_iterator active_cell_iterator
TrilinosWrappers::MPI::BlockVector VEC
Definition: neo_hookean_two_fields.h:34
virtual void parse_parameters_call_back()
Definition: neo_hookean_two_fields.h:168
NeoHookeanTwoFieldsInterface()
Definition: neo_hookean_two_fields.h:62
virtual void declare_parameters(ParameterHandler &prm)
Definition: neo_hookean_two_fields.h:116
void preconditioner_energy(const typename DoFHandler< dim, spacedim >::active_cell_iterator &, Scratch &, CopyPreconditioner &, Number &energy) const
Definition: neo_hookean_two_fields.h:89
LinearOperator< Range, Domain > null_operator(const LinearOperator< Range, Domain > &op)
void compute_system_operators(const DoFHandler< dim, spacedim > &, const std::vector< shared_ptr< LATrilinos::BlockMatrix > >, LinearOperator< LATrilinos::VectorType > &, LinearOperator< LATrilinos::VectorType > &) const
Definition: neo_hookean_two_fields.h:131
Assembly::CopyData::piDoMUSSystem< dim, dim > CopySystem
Definition: neo_hookean_two_fields.h:33
LinearOperator< typename Solver::vector_type, typename Solver::vector_type > inverse_operator(const LinearOperator< typename Solver::vector_type, typename Solver::vector_type > &op, Solver &solver, const Preconditioner &preconditioner)
FEValuesCache< dim, spacedim > Scratch
Definition: neo_hookean_two_fields.h:31
TrilinosWrappers::BlockSparseMatrix MAT
Definition: neo_hookean_two_fields.h:35
Definition: neo_hookean_two_fields.h:27
BlockType & block(const unsigned int row, const unsigned int column)
Assembly::CopyData::piDoMUSPreconditioner< dim, dim > CopyPreconditioner
Definition: neo_hookean_two_fields.h:32
void system_energy(const typename DoFHandler< dim, spacedim >::active_cell_iterator &, Scratch &, CopySystem &, Number &energy) const
Definition: neo_hookean_two_fields.h:119
void extract_constant_modes(const DoFHandlerType &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
const FiniteElement< dim, spacedim > & get_fe() const