pi-DoMUS: Parallel Deal.II MUltiphysics Solver
compressible_neo_hookean.h
Go to the documentation of this file.
1 
5 #ifndef _compressible_neo_hookean_h_
6 #define _compressible_neo_hookean_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 CompressibleNeoHookeanInterface : public ConservativeInterface<dim,spacedim,dim, CompressibleNeoHookeanInterface<dim,spacedim> >
28 {
29  typedef FEValuesCache<dim,spacedim> Scratch;
30  typedef Assembly::CopyData::piDoMUSPreconditioner<dim,dim> CopyPreconditioner;
31  typedef Assembly::CopyData::piDoMUSSystem<dim,dim> CopySystem;
34 
35 public:
36 
37  /* specific and useful functions for this very problem */
39 
41 
44 
45 
46  /* these functions MUST have the follwowing names
47  * because they are called by the ConservativeInterface class
48  */
49 
50  template<typename Number>
52  Scratch &,
53  CopyPreconditioner &,
54  Number &energy) const;
55 
56  template<typename Number>
58  Scratch &,
59  CopySystem &,
60  Number &energy) const;
61 
65  const std::vector<shared_ptr<MAT> >,
67  LinearOperator<VEC> &) const;
68 
69 private:
70  double E;
71  double nu;
72  double mu;
73  double lambda;
74 
75 
76  mutable shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
77  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> Mp_preconditioner;
78  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> T_preconditioner;
79 
80 };
81 
82 template <int dim, int spacedim>
84  ConservativeInterface<dim,spacedim,dim,CompressibleNeoHookeanInterface<dim,spacedim> >("Compressible NeoHookean Interface",
85  "FESystem[FE_Q(1)^d]",
86  "u,u,u", "1", "1","1")
87 {};
88 
89 
90 template <int dim, int spacedim>
91 template<typename Number>
93  Scratch &fe_cache,
94  CopyPreconditioner &data,
95  Number &energy) const
96 {
97  Number alpha = this->alpha;
98  this->reinit(alpha, cell, fe_cache);
99 
100  const FEValuesExtractors::Vector displacement(0);
101  auto &us = fe_cache.get_values("solution", "u", displacement, alpha);
102 
103  const unsigned int n_q_points = us.size();
104 
105  auto &JxW = fe_cache.get_JxW_values();
106  energy = 0;
107  for (unsigned int q=0; q<n_q_points; ++q)
108  {
109  const Tensor <1, dim, Number> &u = us[q];
110 
111  energy += 0.5*(u*u)*JxW[q];
112  }
113 }
114 
115 template <int dim, int spacedim>
116 template<typename Number>
118  Scratch &fe_cache,
119  CopySystem &data,
120  Number &energy) const
121 {
122  Number alpha = this->alpha;
123  this->reinit(alpha, cell, fe_cache);
124 
125  const FEValuesExtractors::Vector displacement(0);
126  auto &us = fe_cache.get_values("solution", "u", displacement, alpha);
127  auto &us_dot = fe_cache.get_values("solution_dot", "u_dot", displacement, alpha);
128  auto &Fs = fe_cache.get_deformation_gradients("solution", "Fu", displacement, alpha);
129 
130  const unsigned int n_q_points = us.size();
131 
132  auto &JxW = fe_cache.get_JxW_values();
133  energy = 0;
134  for (unsigned int q=0; q<n_q_points; ++q)
135  {
136 
137  const Tensor <1, dim, Number> &u = us[q];
138  const Tensor <1, dim, Number> &u_dot = us_dot[q];
139  const Tensor <2, dim, Number> &F = Fs[q];
140  const Tensor<2, dim, Number> C = transpose(F)*F;
141 
142  Number Ic = trace(C);
143  Number J = determinant(F);
144  Number lnJ = std::log (J);
145 
146  Number psi = (mu/2.)*(Ic-dim) - mu*lnJ + (lambda/2.)*(lnJ)*(lnJ);
147 
148  energy += (psi)*JxW[q];
149 
150  }
151 
152 }
153 
154 
155 template <int dim, int spacedim>
157 {
159  this->add_parameter(prm, &E, "Young's modulus", "10.0", Patterns::Double(0.0));
160  this->add_parameter(prm, &nu, "Poisson's ratio", "0.3", Patterns::Double(0.0));
161 }
162 
163 template <int dim, int spacedim>
165 {
167  mu = E/(2.0*(1.+nu));
168  lambda = (E *nu)/((1.+nu)*(1.-2.*nu));
169 }
170 
171 template <int dim, int spacedim>
172 void
175  const TrilinosWrappers::BlockSparseMatrix &preconditioner_matrix,
176  const std::vector<shared_ptr<MAT> >,
177  LinearOperator<VEC> &system_op,
178  LinearOperator<VEC> &prec_op) const
179 {
180 
181  std::vector<std::vector<bool> > constant_modes;
182  FEValuesExtractors::Vector displacement(0);
183  DoFTools::extract_constant_modes (dh, dh.get_fe().component_mask(displacement),
184  constant_modes);
185 
186  Mp_preconditioner.reset (new TrilinosWrappers::PreconditionJacobi());
187  Amg_preconditioner.reset (new TrilinosWrappers::PreconditionAMG());
188 
190  Amg_data.constant_modes = constant_modes;
191  Amg_data.elliptic = true;
192  Amg_data.higher_order_elements = true;
193  Amg_data.smoother_sweeps = 2;
194  Amg_data.aggregation_threshold = 0.02;
195 
196 // Mp_preconditioner->initialize (preconditioner_matrix.block(0,0));
197  Amg_preconditioner->initialize (preconditioner_matrix.block(0,0),
198  Amg_data);
199 
200 
201  auto A = linear_operator< TrilinosWrappers::MPI::Vector >( matrix.block(0,0) );
202 
203  static ReductionControl solver_control_pre(5000, 1e-8);
204  static SolverCG<TrilinosWrappers::MPI::Vector> solver_CG(solver_control_pre);
205  auto A_inv = inverse_operator( A, solver_CG, *Amg_preconditioner);
206 
207  auto P00 = A_inv;
208 
209  // ASSEMBLE THE PROBLEM:
210  system_op = block_operator<1, 1, VEC >({{
211  {{ A }}
212  }
213  });
214 
215 
216  //const auto S = linear_operator<VEC>(matrix);
217 
218  prec_op = block_operator<1, 1, VEC >({{
219  {{ P00}} ,
220  }
221  });
222 }
223 
224 
226 
227 #endif
228 
Definition: conservative.h:20
std::vector< std::vector< bool > > constant_modes
void declare_parameters(ParameterHandler &prm)
Definition: compressible_neo_hookean.h:63
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: compressible_neo_hookean.h:14
CompressibleNeoHookeanInterface()
Definition: compressible_neo_hookean.h:54
void preconditioner_energy(const typename DoFHandler< dim, spacedim >::active_cell_iterator &, Scratch &, CopyPreconditioner &, Number &energy) const
Definition: compressible_neo_hookean.h:92
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)
void compute_system_operators(const DoFHandler< dim, spacedim > &, const std::vector< shared_ptr< LATrilinos::BlockMatrix > >, LinearOperator< LATrilinos::VectorType > &, LinearOperator< LATrilinos::VectorType > &) const
Definition: compressible_neo_hookean.h:130
BlockType & block(const unsigned int row, const unsigned int column)
~CompressibleNeoHookeanInterface()
Definition: compressible_neo_hookean.h:38
void parse_parameters_call_back()
Definition: compressible_neo_hookean.h:71
void extract_constant_modes(const DoFHandlerType &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
void system_energy(const typename DoFHandler< dim, spacedim >::active_cell_iterator &, Scratch &, CopySystem &, Number &energy) const
Definition: compressible_neo_hookean.h:117
const FiniteElement< dim, spacedim > & get_fe() const