5 #ifndef _hydrogels_three_fields_h_
6 #define _hydrogels_three_fields_h_
9 #include <deal2lkit/parsed_function.h>
30 template <
int dim,
int spacedim,
typename LAC>
44 template <
typename EnergyType,
typename Res
idualType>
46 FEValuesCache<dim,spacedim> &scratch,
47 std::vector<EnergyType> &energies,
48 std::vector<std::vector<ResidualType> > &local_residuals,
49 bool compute_only_system_terms)
const;
53 const std::vector<shared_ptr<LATrilinos::BlockMatrix> >,
75 mutable shared_ptr<TrilinosWrappers::PreconditionAMG> U_prec;
76 mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> p_prec;
77 mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> c_prec;
81 template <
int dim,
int spacedim,
typename LAC>
85 "FESystem[FE_Q(2)^d-FE_Q(1)-FE_Q(2)]",
89 template <
int dim,
int spacedim,
typename LAC>
94 couplings[0] =
"1,1,0;1,0,1;0,1,1";
95 couplings[1] =
"0,0,0;0,1,0;0,0,0";
98 template <
int dim,
int spacedim,
typename LAC>
99 template <
typename EnergyType,
typename Res
idualType>
103 FEValuesCache<dim,spacedim> &fe_cache,
104 std::vector<EnergyType> &energies,
105 std::vector<std::vector<ResidualType> > &,
106 bool compute_only_system_terms)
const
108 EnergyType alpha = this->alpha;
109 this->reinit(alpha, cell, fe_cache);
115 auto &Fs = fe_cache.get_deformation_gradients(
"solution",
"Fu", displacement, alpha);
117 auto &ps = fe_cache.get_values(
"solution",
"p", pressure, alpha);
119 auto &cs = fe_cache.get_values(
"solution",
"c", concentration, alpha);
121 const unsigned int n_q_points = ps.size();
123 auto &JxW = fe_cache.get_JxW_values();
125 for (
unsigned int q=0; q<n_q_points; ++q)
129 const EnergyType &c = cs[q];
130 const EnergyType &p = ps[q];
132 const EnergyType I = trace(C);
133 const EnergyType J = determinant(F);
136 EnergyType psi = ( 0.5*G*l0_3*(l02*I - dim)
138 + (l0_3*R*T/Omega)*((Omega*l03*c)*std::log((Omega*l03*c)/(1.+Omega*l03*c))
139 + chi*((Omega*l03*c)/(1.+Omega*l03*c)) )
141 - (mu0)*c - p*(J-l0_3-Omega*c)) ;
143 energies[0] += psi*JxW[q];
145 if (!compute_only_system_terms)
146 energies[1] += 0.5*(p*p)*JxW[q];
151 template <
int dim,
int spacedim,
typename LAC>
163 template <
int dim,
int spacedim,
typename LAC>
171 mu0 = R*T*(std::log((l03-1.)/l03) + l0_3 + chi*l0_6) + G*Omega/l0;
175 template <
int dim,
int spacedim,
typename LAC>
179 const std::vector<shared_ptr<LATrilinos::BlockMatrix> > matrices,
184 std::vector<std::vector<bool> > constant_modes;
201 U_prec->initialize (matrices[0]->block(0,0), Amg_data);
202 p_prec->initialize (matrices[1]->block(1,1));
203 c_prec->initialize (matrices[0]->block(2,2));
207 auto A = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,0) );
208 auto Bt = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,1) );
209 auto Z02 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,2) );
211 auto B = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,0) );
212 auto Z11 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,1) );
213 auto C = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,2) );
215 auto Z20 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,0) );
216 auto D = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,1) );
217 auto E = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,2) );
219 auto P0 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,0));
220 auto P1 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[1]->block(1,1));
221 auto P2 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,2));
236 const std::array<std::array<LinearOperator<LATrilinos::VectorType::BlockType>, 3 >, 3 > matrix_array = {{
243 system_op = block_operator<3, 3, LATrilinos::VectorType >(matrix_array);
245 const std::array<LinearOperator<LATrilinos::VectorType::BlockType>, 3 > diagonal_array = {{ P0_i, P1_i, P2_i }};
248 prec_op = block_diagonal_operator<3,LATrilinos::VectorType>(diagonal_array);
void declare_parameters(ParameterHandler &prm)
Definition: hydrogels_three_fields.h:152
void compute_system_operators(const DoFHandler< dim, spacedim > &, const std::vector< shared_ptr< LATrilinos::BlockMatrix > >, LinearOperator< LATrilinos::VectorType > &, LinearOperator< LATrilinos::VectorType > &) const
Definition: hydrogels_three_fields.h:178
Definition: lac_type.h:33
double aggregation_threshold
std::vector< std::vector< bool > > constant_modes
ActiveSelector::active_cell_iterator active_cell_iterator
This is the class that users should derive from.
Definition: pde_system_interface.h:46
unsigned int smoother_sweeps
void set_matrix_couplings(std::vector< std::string > &couplings) const
Definition: hydrogels_three_fields.h:92
void parse_parameters_call_back()
Definition: hydrogels_three_fields.h:164
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)
~HydroGelThreeFields()
Definition: hydrogels_three_fields.h:35
Definition: hydrogels_three_fields.h:31
void energies_and_residuals(const typename DoFHandler< dim, spacedim >::active_cell_iterator &cell, FEValuesCache< dim, spacedim > &scratch, std::vector< EnergyType > &energies, std::vector< std::vector< ResidualType > > &local_residuals, bool compute_only_system_terms) const
Definition: hydrogels_three_fields.h:102
HydroGelThreeFields()
Definition: hydrogels_three_fields.h:82
const FiniteElement< dim, spacedim > & get_fe() const
bool higher_order_elements