Reduced Lagrange Multipliers
 
Loading...
Searching...
No Matches
coupled_elasticity.h
Go to the documentation of this file.
1/* ---------------------------------------------------------------------
2 */
3
4#ifndef dealii_distributed_lagrange_multiplier_coupled_elasticity_h
5#define dealii_distributed_lagrange_multiplier_coupled_elasticity_h
6
7#ifdef ENABLE_COUPLED_PROBLEMS
8
9# include <deal.II/base/function.h>
10# include <deal.II/base/parsed_convergence_table.h>
11# include <deal.II/base/quadrature_lib.h>
12# include <deal.II/base/timer.h>
13
14# include <deal.II/lac/block_linear_operator.h>
15# include <deal.II/lac/generic_linear_algebra.h>
16# include <deal.II/lac/linear_operator.h>
17# include <deal.II/lac/linear_operator_tools.h>
18# define FORCE_USE_OF_TRILINOS
19namespace LA
20{
21# if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
22 !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
23 using namespace dealii::LinearAlgebraPETSc;
24# define USE_PETSC_LA
25# elif defined(DEAL_II_WITH_TRILINOS)
26 using namespace dealii::LinearAlgebraTrilinos;
27# else
28# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
29# endif
30} // namespace LA
31# include <deal.II/base/conditional_ostream.h>
32# include <deal.II/base/index_set.h>
33# include <deal.II/base/parameter_acceptor.h>
34# include <deal.II/base/parsed_function.h>
35# include <deal.II/base/utilities.h>
36# include <deal.II/base/work_stream.h>
37
38# include <deal.II/distributed/grid_refinement.h>
39# include <deal.II/distributed/solution_transfer.h>
40# include <deal.II/distributed/tria.h>
41
42# include <deal.II/dofs/dof_handler.h>
43# include <deal.II/dofs/dof_renumbering.h>
44# include <deal.II/dofs/dof_tools.h>
45
46# include <deal.II/fe/fe_nothing.h>
47# include <deal.II/fe/fe_q.h>
48# include <deal.II/fe/fe_system.h>
49# include <deal.II/fe/fe_values.h>
50# include <deal.II/fe/mapping_fe_field.h>
51# include <deal.II/fe/mapping_q.h>
52
53# include <deal.II/grid/grid_generator.h>
54# include <deal.II/grid/grid_in.h>
55# include <deal.II/grid/grid_out.h>
56# include <deal.II/grid/grid_refinement.h>
57# include <deal.II/grid/manifold_lib.h>
58# include <deal.II/grid/tria.h>
59# include <deal.II/grid/tria_accessor.h>
60# include <deal.II/grid/tria_iterator.h>
61
62# include <deal.II/lac/affine_constraints.h>
63# include <deal.II/lac/dynamic_sparsity_pattern.h>
64# include <deal.II/lac/full_matrix.h>
65# include <deal.II/lac/petsc_precondition.h>
66# include <deal.II/lac/petsc_solver.h>
67# include <deal.II/lac/petsc_sparse_matrix.h>
68# include <deal.II/lac/petsc_vector.h>
69# include <deal.II/lac/solver_cg.h>
70# include <deal.II/lac/solver_gmres.h>
71# include <deal.II/lac/solver_minres.h>
72# include <deal.II/lac/sparsity_tools.h>
73# include <deal.II/lac/vector.h>
74
75# include <deal.II/meshworker/dof_info.h>
76# include <deal.II/meshworker/integration_info.h>
77# include <deal.II/meshworker/loop.h>
78# include <deal.II/meshworker/scratch_data.h>
79# include <deal.II/meshworker/simple.h>
80
81# include <deal.II/numerics/data_out.h>
82# include <deal.II/numerics/error_estimator.h>
83# include <deal.II/numerics/vector_tools.h>
84
85# include <deal.II/opencascade/manifold_lib.h>
86# include <deal.II/opencascade/utilities.h>
87
88# include "inclusions.h"
89
90
91# ifdef DEAL_II_WITH_OPENCASCADE
92# include <TopoDS.hxx>
93# endif
94# include <deal.II/base/hdf5.h>
95
96# include <cmath>
97# include <fstream>
98# include <iomanip>
99# include <iostream>
100# include <memory>
101
102
103
104template <int dim, int spacedim = dim>
105class CoupledElasticityProblemParameters : public ParameterAcceptor
106{
107public:
108 CoupledElasticityProblemParameters();
109
110 std::string output_directory = ".";
111 std::string output_name = "solution";
112 unsigned int fe_degree = 1;
113 unsigned int initial_refinement = 5;
114 std::list<types::boundary_id> dirichlet_ids{0};
115 std::list<types::boundary_id> neumann_ids{};
116 std::set<types::boundary_id> normal_flux_ids{};
117 std::string domain_type = "";
118 std::string name_of_grid = "hyper_cube";
119 std::string arguments_for_grid = "-1: 1: false";
120 std::string refinement_strategy = "fixed_fraction";
121 double coarsening_fraction = 0.0;
122 double refinement_fraction = 0.3;
123 unsigned int n_refinement_cycles = 1;
124 unsigned int max_cells = 20000;
125
126 double Lame_mu = 1;
127 double Lame_lambda = 1;
128
129 mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>> rhs;
130 mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
131 exact_solution;
132 mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>> bc;
133 mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
134 Neumann_bc;
135
136 std::string weight_expression = "1.";
137
138 mutable ParameterAcceptorProxy<ReductionControl> inner_control;
139 mutable ParameterAcceptorProxy<ReductionControl> outer_control;
140
141 bool output_results = false;
142
143 mutable ParsedConvergenceTable convergence_table;
144
145 // Time dependency.
146 double initial_time = 0.0;
147 double final_time = 0.0;
148 double dt = 5e-3;
149};
150
151
152template <int dim, int spacedim = dim>
153class CoupledElasticityProblem : public Subscriptor
154{
155public:
156 CoupledElasticityProblem(
157 const CoupledElasticityProblemParameters<dim, spacedim> &par);
158
159 void
160 run();
161 void
162 run_timestep0();
163 void
164 run_timestep();
165 void
166 compute_coupling_pressure();
167
168 void
169 update_inclusions_data(std::vector<double> new_data);
170 void
171 update_inclusions_data(std::vector<double> new_data,
172 std::vector<int> cells_per_vessel);
173
174 std::vector<std::vector<double>>
175 split_pressure_over_inclusions(std::vector<int>, Vector<double>) const;
176
177 unsigned int
178 n_vessels() const
179 {
180 return inclusions.get_n_vessels();
181 };
182
183 TrilinosWrappers::MPI::Vector coupling_pressure;
184 TrilinosWrappers::MPI::Vector
185 // Vector<double>
186 coupling_pressure_at_inclusions;
187
188private:
189 void
190 make_grid();
191 void
192 setup_fe();
193 void
194 setup_dofs();
195 void
196 assemble_elasticity_system();
197 void
198 assemble_coupling();
199 void
200 reassemble_coupling_rhs();
201
202 void
203 check_boundary_ids();
204
208 IndexSet
209 assemble_coupling_sparsity(DynamicSparsityPattern &dsp);
210
211 void
212 solve();
213
214 void
215 refine_and_transfer();
216
217 void
218 refine_and_transfer_around_inclusions();
219
220 void
221 execute_actual_refine_and_transfer();
222
223 std::string
224 output_solution() const;
225
226 void
227 output_results() const;
228
229 void
230 print_parameters() const;
231
232 void
233 compute_internal_and_boundary_stress(bool openfilefirsttime) const;
234
235 void
236 compute_face_stress(bool /* openfilefirsttime */){};
237
238 // TrilinosWrappers::MPI::Vector
239 // output_pressure(bool openfilefirsttime) /*const*/;
240 void
241 output_coupling_pressure(bool) const;
242
243 // std::vector<types::global_dof_index>
244 // get_inclusions_of_vessel() const;
245
246 const CoupledElasticityProblemParameters<dim, spacedim> &par;
247 MPI_Comm mpi_communicator;
248 ConditionalOStream pcout;
249 mutable TimerOutput computing_timer;
250 parallel::distributed::Triangulation<spacedim> tria;
251 std::unique_ptr<FiniteElement<spacedim>> fe;
252 Inclusions<spacedim> inclusions;
253 std::unique_ptr<Quadrature<spacedim>> quadrature;
254 std::unique_ptr<Quadrature<spacedim - 1>> face_quadrature_formula;
255 DoFHandler<spacedim> dh;
256 std::vector<IndexSet> owned_dofs;
257 std::vector<IndexSet> relevant_dofs;
258
259 AffineConstraints<double> constraints;
260 AffineConstraints<double> inclusion_constraints;
261 AffineConstraints<double> mean_value_constraints;
262
263 LA::MPI::SparseMatrix stiffness_matrix;
264 LA::MPI::SparseMatrix coupling_matrix;
265 LA::MPI::SparseMatrix inclusion_matrix;
266 LA::MPI::BlockVector solution;
267 LA::MPI::BlockVector locally_relevant_solution;
268 LA::MPI::BlockVector system_rhs;
269 std::vector<std::vector<BoundingBox<spacedim>>> global_bounding_boxes;
270 unsigned int cycle = 0;
271
272 FEValuesExtractors::Vector displacement;
273
274 // Postprocessing values
275 std::map<types::boundary_id, Tensor<1, spacedim>> forces;
276 std::map<types::boundary_id, Tensor<1, spacedim>> average_displacements;
277 std::map<types::boundary_id, Tensor<1, spacedim>> average_normals;
278 std::map<types::boundary_id, double> areas;
279 // std::vector<BaseClass::BlockType> pressure_records;
280
281 // Time dependency.
282 double current_time = 0.0;
283
284 // mutable std::unique_ptr<HDF5::File> pressure_file;
285 // std::ofstream pressure_file;
286 // std::ofstream forces_file;
287};
288
289#endif
290#endif
affine_constraints.h
block_linear_operator.h
conditional_ostream.h
grid_refinement.h
solution_transfer.h
dof_handler.h
dof_renumbering.h
dof_tools.h
dynamic_sparsity_pattern.h
Subscriptor
EnableObserverPointer Subscriptor
error_estimator.h
fe_values.h
fe_nothing.h
fe_q.h
fe_system.h
full_matrix.h
generic_linear_algebra.h
grid_refinement.h
manifold_lib.h
grid_generator.h
grid_in.h
grid_out.h
index_set.h
integration_info.h
linear_operator.h
linear_operator_tools.h
mapping_fe_field.h
mapping_q.h
dof_info.h
e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
WorkStream::run
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
data_out.h
manifold_lib.h
parameter_acceptor.h
parsed_convergence_table.h
parsed_function.h
petsc_precondition.h
petsc_solver.h
petsc_sparse_matrix.h
petsc_vector.h
quadrature_lib.h
scratch_data.h
solver_cg.h
solver_gmres.h
solver_minres.h
sparsity_tools.h
tria_accessor.h
tria_iterator.h
vector_tools.h
work_stream.h