Reduced Lagrange Multipliers
 
Loading...
Searching...
No Matches
coupled_elasticity.cc
Go to the documentation of this file.
1/* ---------------------------------------------------------------------
2 *
3 * Copyright (C) 2000 - 2020 by the deal.II authors
4 *
5 * This file is part of the deal.II library.
6 *
7 * The deal.II library is free software; you can use it, redistribute
8 * it, and/or modify it under the terms of the GNU Lesser General
9 * Public License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 * The full text of the license can be found in the file LICENSE.md at
12 * the top level directory of deal.II.
13 *
14 * ---------------------------------------------------------------------
15 *
16 * Author: Wolfgang Bangerth, University of Heidelberg, 2000
17 * Modified by: Luca Heltai, 2020
18 */
19
20
21#include "coupled_elasticity.h"
22
23#include <deal.II/grid/grid_tools.h>
24
25#ifdef FALSE
26
27template <int dim, int spacedim>
28CoupledElasticityProblemParameters<dim, spacedim>::
29 CoupledElasticityProblemParameters()
30 : ParameterAcceptor("/Immersed Problem/")
31 , rhs("/Immersed Problem/Right hand side", spacedim)
32 , exact_solution("/Immersed Problem/Exact solution", spacedim)
33 , bc("/Immersed Problem/Dirichlet boundary conditions", spacedim)
34 , Neumann_bc("/Immersed Problem/Neumann boundary conditions", spacedim)
35 , inner_control("/Immersed Problem/Solver/Inner control")
36 , outer_control("/Immersed Problem/Solver/Outer control")
37 , convergence_table(std::vector<std::string>(spacedim, "u"))
38{
39 add_parameter("FE degree", fe_degree, "", this->prm, Patterns::Integer(1));
40 add_parameter("Output directory", output_directory);
41 add_parameter("Output name", output_name);
42 add_parameter("Output results", output_results);
43 add_parameter("Initial refinement", initial_refinement);
44 add_parameter("Dirichlet boundary ids", dirichlet_ids);
45 add_parameter("Neumann boundary ids", neumann_ids);
46 add_parameter("Normal flux boundary ids", normal_flux_ids);
47 enter_subsection("Grid generation");
48 {
49 add_parameter("Domain type", domain_type);
50 add_parameter("Grid generator", name_of_grid);
51 add_parameter("Grid generator arguments", arguments_for_grid);
52 }
53 leave_subsection();
54 enter_subsection("Refinement and remeshing");
55 {
56 add_parameter("Strategy",
57 refinement_strategy,
58 "",
59 this->prm,
60 Patterns::Selection("fixed_fraction|fixed_number|global"));
61 add_parameter("Coarsening fraction", coarsening_fraction);
62 add_parameter("Refinement fraction", refinement_fraction);
63 add_parameter("Maximum number of cells", max_cells);
64 add_parameter("Number of refinement cycles", n_refinement_cycles);
65 }
66 leave_subsection();
67 enter_subsection("Physical constants");
68 {
69 add_parameter("Lame mu", Lame_mu);
70 add_parameter("Lame lambda", Lame_lambda);
71 }
72 leave_subsection();
73 enter_subsection("Exact solution");
74 {
75 add_parameter("Weight expression", weight_expression);
76 }
77 leave_subsection();
78 enter_subsection("Time dependency");
79 {
80 add_parameter("Initial time", initial_time);
81 add_parameter("Final time", final_time);
82 add_parameter("Time step", dt);
83 }
84 leave_subsection();
85
86 this->prm.enter_subsection("Error");
87 convergence_table.add_parameters(this->prm);
88 this->prm.leave_subsection();
89}
90
91
92template <int dim, int spacedim>
93CoupledElasticityProblem<dim, spacedim>::CoupledElasticityProblem(
94 const CoupledElasticityProblemParameters<dim, spacedim> &par)
95 : par(par)
96 , mpi_communicator(MPI_COMM_WORLD)
97 , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
98 , computing_timer(mpi_communicator,
99 pcout,
100 TimerOutput::summary,
101 TimerOutput::wall_times)
102 , tria(mpi_communicator,
103 typename Triangulation<spacedim>::MeshSmoothing(
104 Triangulation<spacedim>::smoothing_on_refinement |
105 Triangulation<spacedim>::smoothing_on_coarsening))
106 , inclusions(spacedim)
107 , dh(tria)
108 , displacement(0)
109{}
110
111
112template <int dim, int spacedim>
113void
114read_grid_and_cad_files(const std::string &grid_file_name,
115 const std::string &ids_and_cad_file_names,
116 Triangulation<dim, spacedim> &tria)
117{
118 GridIn<dim, spacedim> grid_in;
119 grid_in.attach_triangulation(tria);
120 grid_in.read(grid_file_name);
121# ifdef DEAL_II_WITH_OPENCASCADE
122 using map_type = std::map<types::manifold_id, std::string>;
123 using Converter = Patterns::Tools::Convert<map_type>;
124 for (const auto &pair : Converter::to_value(ids_and_cad_file_names))
125 {
126 const auto &manifold_id = pair.first;
127 const auto &cad_file_name = pair.second;
128 const auto extension = boost::algorithm::to_lower_copy(
129 cad_file_name.substr(cad_file_name.find_last_of('.') + 1));
130 TopoDS_Shape shape;
131 if (extension == "iges" || extension == "igs")
132 shape = OpenCASCADE::read_IGES(cad_file_name);
133 else if (extension == "step" || extension == "stp")
134 shape = OpenCASCADE::read_STEP(cad_file_name);
135 else
136 AssertThrow(false,
137 ExcNotImplemented("We found an extension that we "
138 "do not recognize as a CAD file "
139 "extension. Bailing out."));
140 const auto n_elements = OpenCASCADE::count_elements(shape);
141 if ((std::get<0>(n_elements) == 0))
142 tria.set_manifold(
143 manifold_id,
144 OpenCASCADE::ArclengthProjectionLineManifold<dim, spacedim>(shape));
145 else if (spacedim == 3)
146 {
147 const auto t = reinterpret_cast<Triangulation<dim, 3> *>(&tria);
148 t->set_manifold(manifold_id,
149 OpenCASCADE::NormalToMeshProjectionManifold<dim, 3>(
150 shape));
151 }
152 else
153 tria.set_manifold(manifold_id,
154 OpenCASCADE::NURBSPatchManifold<dim, spacedim>(
155 TopoDS::Face(shape)));
156 }
157# else
158 (void)ids_and_cad_file_names;
159 AssertThrow(false, ExcNotImplemented("Generation of the grid failed."));
160# endif
161}
162
163
164
165template <int dim, int spacedim>
166void
167CoupledElasticityProblem<dim, spacedim>::make_grid()
168{
169 if (par.domain_type == "generate")
170 {
171 try
172 {
173 GridGenerator::generate_from_name_and_arguments(
174 tria, par.name_of_grid, par.arguments_for_grid);
175 }
176 catch (...)
177 {
178 pcout << "Generating from name and argument failed." << std::endl
179 << "Trying to read from file name." << std::endl;
180 read_grid_and_cad_files(par.name_of_grid,
181 par.arguments_for_grid,
182 tria);
183 }
184 }
185 else if (par.domain_type == "cylinder")
186 {
187 Assert(spacedim == 2, ExcInternalError());
188 GridGenerator::hyper_ball(tria, Point<spacedim>(), 1.);
189 std::cout << " ATTENTION: GRID: cirle of radius 1." << std::endl;
190 }
191 else if (par.domain_type == "cheese")
192 {
193 Assert(spacedim == 2, ExcInternalError());
194 GridGenerator::cheese(tria, std::vector<unsigned int>(2, 2));
195 }
196 else if (par.domain_type == "file")
197 {
198 GridIn<spacedim> gi;
199 gi.attach_triangulation(tria);
200 // std::ifstream infile(par.name_of_grid);
201 const std::string infile(par.name_of_grid);
202 Assert(!infile.empty(), ExcIO());
203 try
204 {
205 gi.read_msh(infile);
206 // gi.read_vtk(infile);
207 }
208 catch (...)
209 {
210 Assert(false, ExcInternalError());
211 }
212 }
213
214 tria.refine_global(par.initial_refinement);
215}
216
217
218
219template <int dim, int spacedim>
220void
221CoupledElasticityProblem<dim, spacedim>::setup_fe()
222{
223 TimerOutput::Scope t(computing_timer, "Initial setup");
224 fe = std::make_unique<FESystem<spacedim>>(FE_Q<spacedim>(par.fe_degree),
225 spacedim);
226 quadrature = std::make_unique<QGauss<spacedim>>(par.fe_degree + 1);
227 face_quadrature_formula =
228 std::make_unique<QGauss<spacedim - 1>>(par.fe_degree + 1);
229}
230
231
232template <int dim, int spacedim>
233void
234CoupledElasticityProblem<dim, spacedim>::setup_dofs()
235{
236 TimerOutput::Scope t(computing_timer, "Setup dofs");
237 dh.distribute_dofs(*fe);
238
239 owned_dofs.resize(2);
240 owned_dofs[0] = dh.locally_owned_dofs();
241 relevant_dofs.resize(2);
242 DoFTools::extract_locally_relevant_dofs(dh, relevant_dofs[0]);
243
244 FEFaceValues<spacedim> fe_face_values(*fe,
245 *face_quadrature_formula,
246 update_values | update_JxW_values |
247 update_quadrature_points |
248 update_normal_vectors);
249
250 {
251 constraints.reinit(relevant_dofs[0]);
252 DoFTools::make_hanging_node_constraints(dh, constraints);
253 for (const auto id : par.dirichlet_ids)
254 {
255 VectorTools::interpolate_boundary_values(dh, id, par.bc, constraints);
256 }
257 std::map<types::boundary_id, const Function<spacedim, double> *>
258 function_map;
259 for (const auto id : par.normal_flux_ids)
260 {
261 function_map.insert(
262 std::pair<types::boundary_id, const Function<spacedim, double> *>(
263 id, &par.Neumann_bc));
264 }
265 VectorTools::compute_nonzero_normal_flux_constraints(
266 dh, 0, par.normal_flux_ids, function_map, constraints);
267 constraints.close();
268
269 /*{
270 mean_value_constraints.clear();
271 mean_value_constraints.reinit(relevant_dofs[0]);
272
273 for (const auto id : par.normal_flux_ids)
274 {
275 const std::set<types::boundary_id > &boundary_ids={id};
276 const ComponentMask &component_mask=ComponentMask();
277 const IndexSet boundary_dofs = DoFTools::extract_boundary_dofs(dh,
278 component_mask, boundary_ids);
279
280 const types::global_dof_index first_boundary_dof =
281 boundary_dofs.nth_index_in_set(0);
282
283 mean_value_constraints.add_line(first_boundary_dof);
284 for (types::global_dof_index i : boundary_dofs)
285 if (i != first_boundary_dof)
286 mean_value_constraints.add_entry(first_boundary_dof, i, -1);
287 }
288 mean_value_constraints.close();
289
290 constraints.merge(mean_value_constraints);
291 }*/
292 }
293 {
294 DynamicSparsityPattern dsp(relevant_dofs[0]);
295 DoFTools::make_sparsity_pattern(dh, dsp, constraints, false);
296 SparsityTools::distribute_sparsity_pattern(dsp,
297 owned_dofs[0],
298 mpi_communicator,
299 relevant_dofs[0]);
300 stiffness_matrix.clear();
301 stiffness_matrix.reinit(owned_dofs[0],
302 owned_dofs[0],
303 dsp,
304 mpi_communicator);
305 }
306 inclusion_constraints.close();
307
308 if (inclusions.n_dofs() > 0)
309 {
310 auto inclusions_set =
311 Utilities::MPI::create_evenly_distributed_partitioning(
312 mpi_communicator, inclusions.n_inclusions());
313
314 owned_dofs[1] = inclusions_set.tensor_product(
315 complete_index_set(inclusions.n_dofs_per_inclusion()));
316
317 DynamicSparsityPattern dsp(dh.n_dofs(),
318 inclusions.n_dofs(),
319 relevant_dofs[0]);
320
321 relevant_dofs[1] = assemble_coupling_sparsity(dsp);
322 relevant_dofs[1].add_indices(owned_dofs[1]);
323 SparsityTools::distribute_sparsity_pattern(dsp,
324 owned_dofs[0],
325 mpi_communicator,
326 relevant_dofs[0]);
327 coupling_matrix.clear();
328 coupling_matrix.reinit(owned_dofs[0],
329 owned_dofs[1],
330 dsp,
331 mpi_communicator);
332
333 DynamicSparsityPattern idsp(relevant_dofs[1]);
334 for (const auto i : relevant_dofs[1])
335 idsp.add(i, i);
336
337 SparsityTools::distribute_sparsity_pattern(idsp,
338 owned_dofs[1],
339 mpi_communicator,
340 relevant_dofs[1]);
341 inclusion_matrix.clear();
342 inclusion_matrix.reinit(owned_dofs[1],
343 owned_dofs[1],
344 idsp,
345 mpi_communicator);
346 }
347
348 locally_relevant_solution.reinit(owned_dofs, relevant_dofs, mpi_communicator);
349 system_rhs.reinit(owned_dofs, mpi_communicator);
350 solution.reinit(owned_dofs, mpi_communicator);
351
352 if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
353 {
354 pcout << " Number of degrees of freedom: " << owned_dofs[0].size()
355 << " + " << owned_dofs[1].size()
356 << " (locally owned: " << owned_dofs[0].n_elements() << " + "
357 << owned_dofs[1].n_elements() << ")" << std::endl;
358 }
359}
360
361template <int dim, int spacedim>
362void
363CoupledElasticityProblem<dim, spacedim>::assemble_elasticity_system()
364{
365 stiffness_matrix = 0;
366 coupling_matrix = 0;
367 system_rhs = 0;
368 TimerOutput::Scope t(computing_timer, "Assemble Stiffness and rhs");
369 FEValues<spacedim> fe_values(*fe,
370 *quadrature,
371 update_values | update_gradients |
372 update_quadrature_points | update_JxW_values);
373 FEFaceValues<spacedim> fe_face_values(*fe,
374 *face_quadrature_formula,
375 update_values | update_JxW_values |
376 update_quadrature_points |
377 update_normal_vectors);
378
379 const unsigned int dofs_per_cell = fe->n_dofs_per_cell();
380 const unsigned int n_q_points = quadrature->size();
381 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
382 Vector<double> cell_rhs(dofs_per_cell);
383 std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(spacedim));
384 std::vector<Tensor<2, spacedim>> grad_phi_u(dofs_per_cell);
385 std::vector<double> div_phi_u(dofs_per_cell);
386 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
387
388 for (const auto &cell : dh.active_cell_iterators())
389 if (cell->is_locally_owned())
390 {
391 cell_matrix = 0;
392 cell_rhs = 0;
393 fe_values.reinit(cell);
394 par.rhs.vector_value_list(fe_values.get_quadrature_points(),
395 rhs_values);
396 for (unsigned int q = 0; q < n_q_points; ++q)
397 {
398 for (unsigned int k = 0; k < dofs_per_cell; ++k)
399 {
400 grad_phi_u[k] =
401 fe_values[displacement].symmetric_gradient(k, q);
402 div_phi_u[k] = fe_values[displacement].divergence(k, q);
403 }
404 for (unsigned int i = 0; i < dofs_per_cell; ++i)
405 {
406 for (unsigned int j = 0; j < dofs_per_cell; ++j)
407 {
408 cell_matrix(i, j) +=
409 (2 * par.Lame_mu *
410 scalar_product(grad_phi_u[i], grad_phi_u[j]) +
411 par.Lame_lambda * div_phi_u[i] * div_phi_u[j]) *
412 fe_values.JxW(q);
413 }
414 const auto comp_i = fe->system_to_component_index(i).first;
415 cell_rhs(i) += fe_values.shape_value(i, q) *
416 rhs_values[q][comp_i] * fe_values.JxW(q);
417 }
418 }
419
420
421 // Neumann boundary conditions
422 // for (const auto &f : cell->face_iterators()) ////
423 for (unsigned int f = 0; f < GeometryInfo<spacedim>::faces_per_cell;
424 ++f)
425 if (cell->face(f)->at_boundary())
426 {
427 // auto it = par.neumann_ids.find(cell->face(f)->boundary_id());
428 // if (it != par.neumann_ids.end())
429 if (std::find(par.neumann_ids.begin(),
430 par.neumann_ids.end(),
431 cell->face(f)->boundary_id()) !=
432 par.neumann_ids.end())
433 {
434 fe_face_values.reinit(cell, f);
435 for (unsigned int q = 0;
436 q < fe_face_values.n_quadrature_points;
437 ++q)
438 {
439 double neumann_value = 0;
440 for (int d = 0; d < spacedim; ++d)
441 neumann_value +=
442 par.Neumann_bc.value(
443 fe_face_values.quadrature_point(q), d) *
444 fe_face_values.normal_vector(q)[d];
445 neumann_value /= spacedim;
446 for (unsigned int i = 0; i < dofs_per_cell; ++i)
447 {
448 cell_rhs(i) += -neumann_value *
449 fe_face_values.shape_value(i, q) *
450 fe_face_values.JxW(q);
451 }
452 }
453 }
454 }
455 cell->get_dof_indices(local_dof_indices);
456 constraints.distribute_local_to_global(cell_matrix,
457 cell_rhs,
458 local_dof_indices,
459 stiffness_matrix,
460 system_rhs.block(0));
461 }
462 stiffness_matrix.compress(VectorOperation::add);
463 system_rhs.compress(VectorOperation::add);
464}
465
466
467
468template <int dim, int spacedim>
469IndexSet
470CoupledElasticityProblem<dim, spacedim>::assemble_coupling_sparsity(
471 DynamicSparsityPattern &dsp)
472{
473 TimerOutput::Scope t(computing_timer,
474 "Setup dofs: Assemble Coupling sparsity");
475
476 IndexSet relevant(inclusions.n_dofs());
477
478 std::vector<types::global_dof_index> dof_indices(fe->n_dofs_per_cell());
479 std::vector<types::global_dof_index> inclusion_dof_indices;
480
481 auto particle = inclusions.inclusions_as_particles.begin();
482 while (particle != inclusions.inclusions_as_particles.end())
483 {
484 const auto &cell = particle->get_surrounding_cell();
485 const auto dh_cell =
486 typename DoFHandler<spacedim>::cell_iterator(*cell, &dh);
487 dh_cell->get_dof_indices(dof_indices);
488 const auto pic =
489 inclusions.inclusions_as_particles.particles_in_cell(cell);
490 Assert(pic.begin() == particle, ExcInternalError());
491 std::set<types::global_dof_index> inclusion_dof_indices_set;
492 for (const auto &p : pic)
493 {
494 const auto ids = inclusions.get_dof_indices(p.get_id());
495 inclusion_dof_indices_set.insert(ids.begin(), ids.end());
496 }
497 inclusion_dof_indices.resize(0);
498 inclusion_dof_indices.insert(inclusion_dof_indices.begin(),
499 inclusion_dof_indices_set.begin(),
500 inclusion_dof_indices_set.end());
501 constraints.add_entries_local_to_global(dof_indices,
502 inclusion_constraints,
503 inclusion_dof_indices,
504 dsp);
505 relevant.add_indices(inclusion_dof_indices.begin(),
506 inclusion_dof_indices.end());
507 particle = pic.end();
508 }
509 return relevant;
510}
511
512
513
514template <int dim, int spacedim>
515void
516CoupledElasticityProblem<dim, spacedim>::assemble_coupling()
517{
518 TimerOutput::Scope t(computing_timer, "Assemble Coupling matrix");
519
520 // system_rhs.block(1) = 0;
521
522 // const FEValuesExtractors::Scalar scalar(0);
523 std::vector<types::global_dof_index> fe_dof_indices(fe->n_dofs_per_cell());
524 std::vector<types::global_dof_index> inclusion_dof_indices(
525 inclusions.n_dofs_per_inclusion());
526
527 FullMatrix<double> local_coupling_matrix(fe->n_dofs_per_cell(),
528 inclusions.n_dofs_per_inclusion());
529
530 FullMatrix<double> local_inclusion_matrix(inclusions.n_dofs_per_inclusion(),
531 inclusions.n_dofs_per_inclusion());
532
533 Vector<double> local_rhs(inclusions.n_dofs_per_inclusion());
534
535 auto particle = inclusions.inclusions_as_particles.begin();
536 while (particle != inclusions.inclusions_as_particles.end())
537 {
538 const auto &cell = particle->get_surrounding_cell();
539 const auto dh_cell =
540 typename DoFHandler<spacedim>::cell_iterator(*cell, &dh);
541 dh_cell->get_dof_indices(fe_dof_indices);
542 const auto pic =
543 inclusions.inclusions_as_particles.particles_in_cell(cell);
544 Assert(pic.begin() == particle, ExcInternalError());
545
546 auto p = pic.begin();
547 auto next_p = pic.begin();
548 while (p != pic.end())
549 {
550 const auto inclusion_id = inclusions.get_inclusion_id(p->get_id());
551 inclusion_dof_indices = inclusions.get_dof_indices(p->get_id());
552 local_coupling_matrix = 0;
553 local_inclusion_matrix = 0;
554 local_rhs = 0;
555
556 // Extract all points that refer to the same inclusion
557 std::vector<Point<spacedim>> ref_q_points;
558 for (; next_p != pic.end() &&
559 inclusions.get_inclusion_id(next_p->get_id()) == inclusion_id;
560 ++next_p)
561 ref_q_points.push_back(next_p->get_reference_location());
562 FEValues<spacedim, spacedim> fev(*fe,
563 ref_q_points,
564 update_values | update_gradients);
565 fev.reinit(dh_cell);
566 // double temp = 0;
567 for (unsigned int q = 0; q < ref_q_points.size(); ++q)
568 {
569 const auto id = p->get_id();
570 const auto &inclusion_fe_values = inclusions.get_fe_values(id);
571 const auto &real_q = p->get_location();
572 const auto ds =
573 inclusions.get_JxW(id); // /inclusions.get_radius(inclusion_id);
574
575 // Coupling and inclusions matrix
576 for (unsigned int j = 0; j < inclusions.n_dofs_per_inclusion();
577 ++j)
578 {
579 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
580 {
581 const auto comp_i =
582 fe->system_to_component_index(i).first;
583 if (comp_i == inclusions.get_component(j))
584 {
585 local_coupling_matrix(i, j) +=
586 (fev.shape_value(i, q)) * inclusion_fe_values[j] *
587 ds;
588 }
589 }
590 if (inclusions.data_file != "")
591 {
592 if (inclusions.inclusions_data[inclusion_id].size() > j)
593 {
594 auto temp =
595 inclusion_fe_values[j] * inclusion_fe_values[j] *
596 inclusions.get_rotated_inclusion_data(
597 inclusion_id)[j] /
598 // inclusions.inclusions_data[inclusion_id][j] / //
599 // data is always prescribed in relative coordinates
600 inclusions.get_radius(inclusion_id) * ds;
601 if (par.initial_time != par.final_time)
602 temp *= inclusions.inclusions_rhs.value(
603 real_q, inclusions.get_component(j));
604 local_rhs(j) += temp;
605 }
606 }
607 else
608 {
609 local_rhs(j) += inclusion_fe_values[j] *
610 inclusions.inclusions_rhs.value(
611 real_q, inclusions.get_component(j)) /
612 inclusions.get_radius(inclusion_id) * ds;
613 }
614 local_inclusion_matrix(j, j) +=
615 (inclusion_fe_values[j] * inclusion_fe_values[j] * ds);
616 }
617 ++p;
618 }
619 // I expect p and next_p to be the same now.
620 Assert(p == next_p, ExcInternalError());
621 // Add local matrices to global ones
622 constraints.distribute_local_to_global(local_coupling_matrix,
623 fe_dof_indices,
624 inclusion_constraints,
625 inclusion_dof_indices,
626 coupling_matrix);
627 inclusion_constraints.distribute_local_to_global(
628 local_rhs, inclusion_dof_indices, system_rhs.block(1));
629
630 // inclusion_constraints.distribute_local_to_global(
631 // local_inclusion_matrix, inclusion_dof_indices, inclusion_matrix);
632 }
633 particle = pic.end();
634 }
635 coupling_matrix.compress(VectorOperation::add);
636 inclusion_matrix.compress(VectorOperation::add);
637 system_rhs.compress(VectorOperation::add);
638}
639
640
641template <int dim, int spacedim>
642void
643CoupledElasticityProblem<dim, spacedim>::reassemble_coupling_rhs()
644{
645 TimerOutput::Scope t(computing_timer, "updating coupling rhs");
646
647 system_rhs.block(1) = 0;
648
649 std::vector<types::global_dof_index> fe_dof_indices(fe->n_dofs_per_cell());
650 std::vector<types::global_dof_index> inclusion_dof_indices(
651 inclusions.n_dofs_per_inclusion());
652
653 Vector<double> local_rhs(inclusions.n_dofs_per_inclusion());
654
655 auto particle = inclusions.inclusions_as_particles.begin();
656 while (particle != inclusions.inclusions_as_particles.end())
657 {
658 const auto &cell = particle->get_surrounding_cell();
659 const auto dh_cell =
660 typename DoFHandler<spacedim>::cell_iterator(*cell, &dh);
661 dh_cell->get_dof_indices(fe_dof_indices);
662 const auto pic =
663 inclusions.inclusions_as_particles.particles_in_cell(cell);
664 Assert(pic.begin() == particle, ExcInternalError());
665
666 auto p = pic.begin();
667 auto next_p = pic.begin();
668 while (p != pic.end())
669 {
670 const auto inclusion_id = inclusions.get_inclusion_id(p->get_id());
671 inclusion_dof_indices = inclusions.get_dof_indices(p->get_id());
672 local_rhs = 0;
673
674 // Extract all points that refer to the same inclusion
675 std::vector<Point<spacedim>> ref_q_points;
676 for (; next_p != pic.end() &&
677 inclusions.get_inclusion_id(next_p->get_id()) == inclusion_id;
678 ++next_p)
679 ref_q_points.push_back(next_p->get_reference_location());
680 FEValues<spacedim, spacedim> fev(*fe,
681 ref_q_points,
682 update_values | update_gradients);
683 fev.reinit(dh_cell);
684 // double temp = 0;
685 for (unsigned int q = 0; q < ref_q_points.size(); ++q)
686 {
687 const auto id = p->get_id();
688 const auto &inclusion_fe_values = inclusions.get_fe_values(id);
689 const auto &real_q = p->get_location();
690 const auto ds =
691 inclusions.get_JxW(id); // /inclusions.get_radius(inclusion_id);
692
693 // Coupling and inclusions matrix
694 for (unsigned int j = 0; j < inclusions.n_dofs_per_inclusion();
695 ++j)
696 {
697 if (inclusions.data_file != "")
698 {
699 if (inclusions.inclusions_data[inclusion_id].size() > j)
700 {
701 auto temp =
702 inclusion_fe_values[j] * inclusion_fe_values[j] *
703 inclusions.get_rotated_inclusion_data(
704 inclusion_id)[j] /
705 // inclusions.inclusions_data[inclusion_id][j] / //
706 // data is always prescribed in relative coordinates
707 inclusions.get_radius(inclusion_id) * ds;
708 if (par.initial_time != par.final_time)
709 temp *= inclusions.inclusions_rhs.value(
710 real_q, inclusions.get_component(j));
711 local_rhs(j) += temp;
712 }
713 }
714 else
715 {
716 local_rhs(j) += inclusion_fe_values[j] *
717 inclusions.inclusions_rhs.value(
718 real_q, inclusions.get_component(j)) /
719 inclusions.get_radius(inclusion_id) * ds;
720 }
721 }
722 ++p;
723 }
724 // I expect p and next_p to be the same now.
725 Assert(p == next_p, ExcInternalError());
726 // Add local matrices to global ones
727 inclusion_constraints.distribute_local_to_global(
728 local_rhs, inclusion_dof_indices, system_rhs.block(1));
729 }
730 particle = pic.end();
731 }
732 system_rhs.compress(VectorOperation::add);
733}
734
735
736
737template <int dim, int spacedim>
738void
739CoupledElasticityProblem<dim, spacedim>::solve()
740{
741 TimerOutput::Scope t(computing_timer, "Solve");
742 LA::MPI::PreconditionAMG prec_A;
743 {
744 // LA::MPI::PreconditionAMG::AdditionalData data;
745 TrilinosWrappers::PreconditionAMG::AdditionalData data;
746# ifdef USE_PETSC_LA
747 data.symmetric_operator = true;
748# endif
749 // informo il precondizionatore dei modi costanti del problema elastico
750 std::vector<std::vector<bool>> constant_modes;
751 const FEValuesExtractors::Vector displacement_components(0); // gia in .h
752 DoFTools::extract_constant_modes(
753 dh, fe->component_mask(displacement_components), constant_modes);
754 data.constant_modes = constant_modes;
755
756 prec_A.initialize(stiffness_matrix, data);
757 }
758
759 const auto A = linear_operator<LA::MPI::Vector>(stiffness_matrix);
760 auto invA = A;
761
762 const auto amgA = linear_operator(A, prec_A);
763
764 // for small radius you might need
765 // SolverFGMRES<LA::MPI::Vector>
766 SolverCG<LA::MPI::Vector> cg_stiffness(par.inner_control);
767 invA = inverse_operator(A, cg_stiffness, amgA);
768
769 // Some aliases
770 auto &u = solution.block(0);
771 auto &lambda = solution.block(1);
772
773 const auto &f = system_rhs.block(0);
774 auto &g = system_rhs.block(1);
775 // MPI_Barrier(mpi_communicator);
776 // g.print(std::cout);
777 // MPI_Barrier(mpi_communicator);
778 if (inclusions.n_dofs() == 0)
779 {
780 u = invA * f;
781 }
782 else
783 {
784 // std::cout << "matrix B" << std::endl;
785 // coupling_matrix.print(std::cout);
786 // std::cout << std::endl << " end matrix B" << std::endl;
787 const auto Bt = linear_operator<LA::MPI::Vector>(coupling_matrix);
788 const auto B = transpose_operator(Bt);
789 const auto M = linear_operator<LA::MPI::Vector>(inclusion_matrix);
790
791 // auto interp_g = g;
792 // interp_g = 0.1;
793 // g = C * interp_g;
794
795 // Schur complement
796 const auto S = B * invA * Bt;
797
798 // Schur complement preconditioner
799 // VERSION 1
800 // auto invS = S;
801 // SolverFGMRES<LA::MPI::Vector> cg_schur(par.outer_control);
802 SolverMinRes<LA::MPI::Vector> cg_schur(par.outer_control);
803 // invS = inverse_operator(S, cg_schur);
804 // VERSION2
805 auto invS = S;
806 auto S_inv_prec = B * invA * Bt + M;
807 // SolverCG<Vector<double>> cg_schur(par.outer_control);
808 // PrimitiveVectorMemory<Vector<double>> mem;
809 // SolverGMRES<Vector<double>> solver_gmres(
810 // par.outer_control, mem,
811 // SolverGMRES<Vector<double>>::AdditionalData(20));
812 invS = inverse_operator(S, cg_schur, S_inv_prec);
813
814 pcout << " f norm: " << f.l2_norm() << ", g norm: " << g.l2_norm()
815 << std::endl;
816 // pcout << " g: ";
817 // g.print(std::cout);
818
819 // Compute Lambda first
820 lambda = invS * (B * invA * f - g);
821 pcout << " Solved for lambda in " << par.outer_control.last_step()
822 << " iterations." << std::endl;
823
824 // Then compute u
825 u = invA * (f - Bt * lambda);
826 pcout << " u norm: " << u.l2_norm()
827 << ", lambda norm: " << lambda.l2_norm() << std::endl;
828 // std::cout << " lambda: ";
829 // lambda.print(std::cout);
830 }
831
832 pcout << " Solved for u in " << par.inner_control.last_step()
833 << " iterations." << std::endl;
834 constraints.distribute(u);
835 inclusion_constraints.distribute(lambda);
836 locally_relevant_solution = solution;
837}
838
839template <int dim, int spacedim>
840void
841CoupledElasticityProblem<dim, spacedim>::refine_and_transfer_around_inclusions()
842{
843 TimerOutput::Scope t(computing_timer, "Refine");
844 Vector<float> error_per_cell(tria.n_active_cells());
845 KellyErrorEstimator<spacedim>::estimate(dh,
846 QGauss<spacedim - 1>(par.fe_degree +
847 1),
848 {},
849 locally_relevant_solution.block(0),
850 error_per_cell);
851
852 const int material_id = 1;
853
854 auto particle = inclusions.inclusions_as_particles.begin();
855 while (particle != inclusions.inclusions_as_particles.end())
856 {
857 const auto &cell = particle->get_surrounding_cell();
858 const auto pic =
859 inclusions.inclusions_as_particles.particles_in_cell(cell);
860 Assert(pic.begin() == particle, ExcInternalError());
861
862 cell->set_refine_flag();
863 cell->set_material_id(material_id);
864
865 // for (auto f : cell->face_indices())
866 // for (unsigned int face_no = 0; face_no
867 // <GeometryInfo<spacedim>::faces_per_cell; ++face_no)
868 for (auto vertex : cell->vertex_indices())
869 {
870 for (const auto &neigh_i :
871 GridTools::find_cells_adjacent_to_vertex(dh, vertex))
872 // for (auto neigh_i = cell->neighbor(face_no)->begin_face();
873 // neigh_i < cell->neighbor(face_no)->end_face(); ++neigh_i)
874 {
875 if (!neigh_i->refine_flag_set())
876 {
877 neigh_i->set_refine_flag();
878 neigh_i->set_material_id(material_id);
879 for (auto vertey : neigh_i->vertex_indices())
880 for (const auto &neigh_j :
881 GridTools::find_cells_adjacent_to_vertex(dh, vertey))
882 // for (auto neigh_j =
883 // cell->neighbor(neigh_i)->begin_face(); neigh_j <
884 // cell->neighbor(neigh_i)->end_face(); ++neigh_j)
885 {
886 neigh_j->set_refine_flag();
887 neigh_j->set_material_id(material_id);
888 }
889 }
890 }
891 }
892
893 particle = pic.end();
894 }
895 execute_actual_refine_and_transfer();
896
897 for (unsigned int ref_cycle = 0; ref_cycle < par.n_refinement_cycles - 1;
898 ++ref_cycle)
899 {
900 for (const auto &cell : tria.active_cell_iterators())
901 {
902 if (cell->material_id() == material_id)
903 cell->set_refine_flag();
904 }
905 execute_actual_refine_and_transfer();
906 }
907}
908
909template <int dim, int spacedim>
910void
911CoupledElasticityProblem<dim, spacedim>::refine_and_transfer()
912{
913 TimerOutput::Scope t(computing_timer, "Refine");
914 Vector<float> error_per_cell(tria.n_active_cells());
915 KellyErrorEstimator<spacedim>::estimate(dh,
916 QGauss<spacedim - 1>(par.fe_degree +
917 1),
918 {},
919 locally_relevant_solution.block(0),
920 error_per_cell);
921 if (par.refinement_strategy == "fixed_fraction")
922 {
923 parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
924 tria, error_per_cell, par.refinement_fraction, par.coarsening_fraction);
925 }
926 else if (par.refinement_strategy == "fixed_number")
927 {
928 parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
929 tria,
930 error_per_cell,
931 par.refinement_fraction,
932 par.coarsening_fraction,
933 par.max_cells);
934 }
935 else if (par.refinement_strategy == "global")
936 for (const auto &cell : tria.active_cell_iterators())
937 cell->set_refine_flag();
938
939 execute_actual_refine_and_transfer();
940}
941template <int dim, int spacedim>
942void
943CoupledElasticityProblem<dim, spacedim>::execute_actual_refine_and_transfer()
944{
945 parallel::distributed::SolutionTransfer<spacedim, LA::MPI::Vector> transfer(
946 dh);
947 tria.prepare_coarsening_and_refinement();
948 inclusions.inclusions_as_particles.prepare_for_coarsening_and_refinement();
949 transfer.prepare_for_coarsening_and_refinement(
950 locally_relevant_solution.block(0));
951 tria.execute_coarsening_and_refinement();
952 inclusions.inclusions_as_particles.unpack_after_coarsening_and_refinement();
953 setup_dofs();
954 transfer.interpolate(solution.block(0));
955 constraints.distribute(solution.block(0));
956 locally_relevant_solution.block(0) = solution.block(0);
957}
958
959
960
961template <int dim, int spacedim>
962std::string
963CoupledElasticityProblem<dim, spacedim>::output_solution() const
964{
965 std::vector<std::string> solution_names(spacedim, "displacement");
966 std::vector<std::string> exact_solution_names(spacedim, "exact_displacement");
967
968
969 auto exact_vec(solution.block(0));
970 // VectorTools::interpolate(dh, par.bc, exact_vec);
971 VectorTools::interpolate(dh, par.exact_solution, exact_vec);
972 auto exact_vec_locally_relevant(locally_relevant_solution.block(0));
973 exact_vec_locally_relevant = exact_vec;
974
975 std::vector<DataComponentInterpretation::DataComponentInterpretation>
976 data_component_interpretation(
977 spacedim, DataComponentInterpretation::component_is_part_of_vector);
978
979 DataOut<spacedim> data_out;
980 data_out.attach_dof_handler(dh);
981 data_out.add_data_vector(locally_relevant_solution.block(0),
982 solution_names,
983 DataOut<spacedim>::type_dof_data,
984 data_component_interpretation);
985
986 data_out.add_data_vector(exact_vec_locally_relevant,
987 exact_solution_names,
988 DataOut<spacedim>::type_dof_data,
989 data_component_interpretation);
990
991 Vector<float> subdomain(tria.n_active_cells());
992 for (unsigned int i = 0; i < subdomain.size(); ++i)
993 subdomain(i) = tria.locally_owned_subdomain();
994 data_out.add_data_vector(subdomain, "subdomain");
995 data_out.build_patches();
996 const std::string filename =
997 par.output_name + "_" + std::to_string(cycle) + ".vtu";
998 data_out.write_vtu_in_parallel(par.output_directory + "/" + filename,
999 mpi_communicator);
1000 return filename;
1001}
1002
1003
1004
1005// template <int dim, int spacedim>
1006// std::string
1007// CoupledElasticityProblem<dim, spacedim>::output_particles() const
1008// {
1009// Particles::DataOut<spacedim> particles_out;
1010// particles_out.build_patches(inclusions.inclusions_as_particles);
1011// const std::string filename =
1012// par.output_name + "_particles_" + std::to_string(cycle) + ".vtu";
1013// particles_out.write_vtu_in_parallel(par.output_directory + "/" +
1014// filename,
1015// mpi_communicator);
1016// return filename;
1017// }
1018
1019
1020template <int dim, int spacedim>
1021void
1022CoupledElasticityProblem<dim, spacedim>::output_results() const
1023{
1024 TimerOutput::Scope t(computing_timer, "Postprocessing: Output results");
1025 static std::vector<std::pair<double, std::string>> cycles_and_solutions;
1026 static std::vector<std::pair<double, std::string>> cycles_and_particles;
1027
1028 if (cycles_and_solutions.size() == cycle)
1029 {
1030 cycles_and_solutions.push_back({(double)cycle, output_solution()});
1031 std::ofstream pvd_solutions(par.output_directory + "/" + par.output_name +
1032 ".pvd");
1033 DataOutBase::write_pvd_record(pvd_solutions, cycles_and_solutions);
1034
1035 if (cycle == 0)
1036 {
1037 const std::string particles_filename =
1038 par.output_name + "_particles.vtu";
1039
1040 inclusions.output_particles(par.output_directory + "/" +
1041 particles_filename);
1042 cycles_and_particles.push_back({(double)cycle, particles_filename});
1043
1044 std::ofstream pvd_particles(par.output_directory + "/" +
1045 par.output_name + "_particles.pvd");
1046 DataOutBase::write_pvd_record(pvd_particles, cycles_and_particles);
1047 }
1048 }
1049}
1050
1051template <int dim, int spacedim>
1052void
1053CoupledElasticityProblem<dim, spacedim>::print_parameters() const
1054{
1055# ifdef USE_PETSC_LA
1056 pcout << "Running CoupledElasticityProblem<"
1057 << Utilities::dim_string(dim, spacedim) << "> using PETSc."
1058 << std::endl;
1059# else
1060 pcout << "Running CoupledElasticityProblem<"
1061 << Utilities::dim_string(dim, spacedim) << "> using Trilinos."
1062 << std::endl;
1063# endif
1064 par.prm.print_parameters(par.output_directory + "/" + "used_parameters_" +
1065 std::to_string(dim) + std::to_string(spacedim) +
1066 ".prm",
1067 ParameterHandler::Short);
1068}
1069
1070template <int dim, int spacedim>
1071void
1072CoupledElasticityProblem<dim, spacedim>::check_boundary_ids()
1073{
1074 for (const auto id : par.dirichlet_ids)
1075 {
1076 for (const auto Nid : par.neumann_ids)
1077 if (id == Nid)
1078 AssertThrow(false,
1079 ExcNotImplemented("incoherent boundary conditions."));
1080 for (const auto noid : par.normal_flux_ids)
1081 if (id == noid)
1082 AssertThrow(false,
1083 ExcNotImplemented("incoherent boundary conditions."));
1084 }
1085}
1086/*
1087template <int dim, int spacedim>
1088void
1089CoupledElasticityProblem<dim, spacedim>::compute_face_stress(bool
1090openfilefirsttime) const
1091{
1092 TimerOutput::Scope t(computing_timer,
1093 "Postprocessing: Computing face stresses");
1094
1095
1096 const std::string full_filename(par.output_directory +
1097"/full_face_stress.csv"); std::ofstream full_face_stress_file; const
1098std::string filename(par.output_directory + "/face_stress.csv"); std::ofstream
1099face_stress_file;
1100
1101 if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
1102 {
1103 full_face_stress_file.open(full_filename);
1104 full_face_stress_file
1105 << "cellx celly cellz area normalx normaly normalz "
1106 //<< "face_stress_x face_stress_y face_stress_z"
1107 << "face_stress_11 face_stress_12 face_stress_13 face_stress_21
1108face_stress_22 face_stress_23 face_stress_31 face_stress_32 face_stress_33"
1109 << std::endl;
1110 face_stress_file.open(filename);
1111 face_stress_file
1112 << "cellx celly cellz area normalx normaly normalz "
1113 << "face_stress_x face_stress_y face_stress_z"
1114 << std::endl;
1115 }
1116 else
1117 {
1118 face_stress_file.open(filename, std::ios_base::app);
1119 full_face_stress_file.open(full_filename, std::ios_base::app);
1120 }
1121
1122 // unsigned int number_active_boundary_faces = 0;
1123 // for (const auto &cell : dh.active_cell_iterators())
1124 // if (cell->is_locally_owned())
1125 // for (unsigned int f = 0; f <
1126GeometryInfo<spacedim>::faces_per_cell;
1127 // ++f)
1128 // if (cell->face(f)->at_boundary())
1129 // number_active_boundary_faces ++;
1130
1131 // Utilities::MPI::sum(count, mpi_communicator);
1132
1133 // Vector<double> sigma(number_active_boundary_faces);
1134 // std::vector<Tensor<1, spacedim>> face_stresses_vector;
1135 // std::vector<Tensor<1, spacedim>> face_normals_vector;
1136 // std::vector<Point<spacedim>> baricenters_vector;
1137
1138 std::vector<std::vector<double>> cells_normals_full_stresses;
1139 std::vector<std::vector<double>> cells_normals_stresses;
1140
1141
1142 FEFaceValues<spacedim> fe_face_values(*fe,
1143 *face_quadrature_formula,
1144 update_values | update_gradients |
1145 update_JxW_values |
1146 update_quadrature_points |
1147 update_normal_vectors);
1148 const FEValuesExtractors::Vector displacement(0);
1149
1150 Tensor<2, spacedim> identity;
1151 for (unsigned int ix = 0; ix < spacedim; ++ix)
1152 identity[ix][ix] = 1;
1153
1154 std::vector<Tensor<2, spacedim>> displacement_gradient(
1155 face_quadrature_formula->size());
1156 std::vector<double> displacement_divergence(
1157 face_quadrature_formula->size());
1158 std::vector<Tensor<1, spacedim>> displacement_values(
1159 face_quadrature_formula->size());
1160
1161 for (const auto &cell : dh.active_cell_iterators())
1162 {
1163 if (cell->is_locally_owned())
1164 {
1165 for (unsigned int f = 0; f < GeometryInfo<spacedim>::faces_per_cell;
1166 ++f)
1167 {
1168 Tensor<1, spacedim, double> face_stress;
1169 Tensor<2, spacedim, double> full_face_stress;
1170 Tensor<1, spacedim, double> face_normal;
1171 Point<spacedim> baricenter(0.0, 0.0, 0.0);
1172 face_stress = 0.0;
1173 full_face_stress = 0.0;
1174 face_normal = 0.0;
1175 double cell_area = 0.0;
1176
1177 if (cell->face(f)->at_boundary())
1178 {
1179 fe_face_values.reinit(cell, f);
1180
1181 fe_face_values[displacement].get_function_gradients(
1182 locally_relevant_solution.block(0), displacement_gradient);
1183 fe_face_values[displacement].get_function_divergences(
1184 locally_relevant_solution.block(0), displacement_divergence);
1185
1186 for (unsigned int q = 0; q < fe_face_values.n_quadrature_points;
1187 ++q)
1188 {
1189 face_stress +=
1190 (2 * par.Lame_mu * displacement_gradient[q] +
1191 par.Lame_lambda *
1192 displacement_divergence[q] * identity
1193 ) *
1194 fe_face_values.JxW(q)* fe_face_values.normal_vector(q);
1195
1196 full_face_stress +=
1197 (2 * par.Lame_mu * displacement_gradient[q] +
1198 par.Lame_lambda *
1199 displacement_divergence[q] * identity
1200 ) *
1201 fe_face_values.JxW(q);
1202
1203 face_normal += fe_face_values.normal_vector(q);
1204 baricenter += fe_face_values.quadrature_point(q);
1205 cell_area += fe_face_values.JxW(q);
1206 }
1207
1208 face_normal /= fe_face_values.get_normal_vectors().size();
1209 baricenter /= fe_face_values.get_quadrature_points().size();
1210 //print on file
1211 // face_stress_file << baricenter << " " << cell_area << " " <<
1212face_normal << " " << face_stress << std::endl; std::vector<double>
1213temp(5*spacedim+1); baricenter.unroll(temp.begin(), temp.begin()+spacedim);
1214 temp[spacedim] = cell_area;
1215 face_normal.unroll(temp.begin()+spacedim+1,
1216temp.begin()+2*spacedim+1); full_face_stress.unroll(temp.begin()+2*spacedim+1,
1217temp.begin()+(3+2)*spacedim+1); cells_normals_full_stresses.push_back(temp);
1218
1219 std::vector<double> tump(3*spacedim+1);
1220 baricenter.unroll(tump.begin(), tump.begin()+spacedim);
1221 tump[spacedim] = cell_area;
1222 face_normal.unroll(tump.begin()+spacedim+1,
1223tump.begin()+2*spacedim+1); face_stress.unroll(tump.begin()+2*spacedim+1,
1224tump.begin()+3*spacedim+1); cells_normals_stresses.push_back(tump);
1225 }
1226 // else
1227 // face_stress = 0.0;
1228 // face_stresses_vector.push_back(face_stress);
1229 // face_normals_vector.push_back(face_normal);
1230 // baricenters_vector.push_back(baricenter);
1231 //sigma.add(cell->active_cell_index(), face_stress);
1232 }
1233 }
1234 }
1235 //print on file
1236
1237 for (unsigned int proc_id = 0; proc_id <
1238Utilities::MPI::n_mpi_processes(mpi_communicator); ++proc_id)
1239 {
1240 if (Utilities::MPI::this_mpi_process(mpi_communicator) == proc_id)
1241 {
1242 for (const auto & cell_data : cells_normals_stresses)
1243 {
1244 for (const auto & data_i : cell_data)
1245 {
1246 face_stress_file << data_i << " ";
1247 }
1248 face_stress_file << std::endl;
1249 }
1250 for (const auto & full_cell_data : cells_normals_full_stresses)
1251 {
1252 for (const auto & full_data_i : full_cell_data)
1253 {
1254 full_face_stress_file << full_data_i << " ";
1255 }
1256 full_face_stress_file << std::endl;
1257 }
1258 }
1259 MPI_Barrier(mpi_communicator);
1260 }
1261 face_stress_file.close();
1262 full_face_stress_file.close();
1263 // sigma_n.compress(VectorOperation::add);
1264
1265 return;
1266}
1267*/
1268template <int dim, int spacedim>
1269void
1270CoupledElasticityProblem<dim, spacedim>::compute_internal_and_boundary_stress(
1271 bool openfilefirsttime) const
1272{
1273 TimerOutput::Scope t(
1274 computing_timer,
1275 "Postprocessing: Computing internal and boundary stresses");
1276
1277 std::map<types::boundary_id, Tensor<1, spacedim>> boundary_stress;
1278 Tensor<2, spacedim> internal_stress;
1279 Tensor<1, spacedim> average_displacement;
1280 std::vector<double> u_dot_n(spacedim * spacedim);
1281
1282 auto all_ids = tria.get_boundary_ids();
1283 std::vector<double> perimeter;
1284 for (auto id : all_ids)
1285 // for (const auto id : par.dirichlet_ids)
1286 {
1287 // boundary_stress[id] = Tensor<1, spacedim>();
1288 boundary_stress[id] = 0.0;
1289 perimeter.push_back(0.0);
1290 }
1291 double internal_area = 0.;
1292 FEValues<spacedim> fe_values(*fe,
1293 *quadrature,
1294 update_values | update_gradients |
1295 update_quadrature_points | update_JxW_values);
1296 FEFaceValues<spacedim> fe_face_values(*fe,
1297 *face_quadrature_formula,
1298 update_values | update_gradients |
1299 update_JxW_values |
1300 update_quadrature_points |
1301 update_normal_vectors);
1302 const FEValuesExtractors::Vector displacement(0);
1303
1304 const unsigned int dofs_per_cell = fe->n_dofs_per_cell();
1305 const unsigned int n_q_points = quadrature->size();
1306 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1307 Tensor<2, spacedim> grad_phi_u;
1308 double div_phi_u;
1309 Tensor<2, spacedim> identity;
1310 for (unsigned int ix = 0; ix < spacedim; ++ix)
1311 identity[ix][ix] = 1;
1312 // std::vector<std::vector<Tensor<1,spacedim>>>
1313 // solution_gradient(face_quadrature_formula->size(),
1314 // std::vector<Tensor<1,spacedim> >(spacedim+1));
1315 std::vector<Tensor<2, spacedim>> displacement_gradient(
1316 face_quadrature_formula->size());
1317 std::vector<double> displacement_divergence(face_quadrature_formula->size());
1318 std::vector<Tensor<1, spacedim>> displacement_values(
1319 face_quadrature_formula->size());
1320
1321 for (const auto &cell : dh.active_cell_iterators())
1322 {
1323 if (cell->is_locally_owned())
1324 {
1325 if constexpr (spacedim == 2)
1326 {
1327 cell->get_dof_indices(local_dof_indices);
1328 fe_values.reinit(cell);
1329 for (unsigned int q = 0; q < n_q_points; ++q)
1330 {
1331 internal_area += fe_values.JxW(q);
1332 for (unsigned int k = 0; k < dofs_per_cell; ++k)
1333 {
1334 grad_phi_u =
1335 fe_values[displacement].symmetric_gradient(k, q);
1336 div_phi_u = fe_values[displacement].divergence(k, q);
1337 internal_stress +=
1338 (2 * par.Lame_mu * grad_phi_u +
1339 par.Lame_lambda * div_phi_u * identity) *
1340 locally_relevant_solution.block(
1341 0)[local_dof_indices[k]] *
1342 fe_values.JxW(q);
1343 average_displacement +=
1344 fe_values[displacement].value(k, q) *
1345 locally_relevant_solution.block(
1346 0)[local_dof_indices[k]] *
1347 fe_values.JxW(q);
1348 }
1349 }
1350 }
1351
1352 for (unsigned int f = 0; f < GeometryInfo<spacedim>::faces_per_cell;
1353 ++f)
1354 // for (const auto &f : cell->face_iterators())
1355 // for (const auto f : GeometryInfo<spacedim>::face_indices())
1356 if (cell->face(f)->at_boundary())
1357 {
1358 auto boundary_index = cell->face(f)->boundary_id();
1359 fe_face_values.reinit(cell, f);
1360
1361 fe_face_values[displacement].get_function_gradients(
1362 locally_relevant_solution.block(0), displacement_gradient);
1363 fe_face_values[displacement].get_function_values(
1364 locally_relevant_solution.block(0), displacement_values);
1365 fe_face_values[displacement].get_function_divergences(
1366 locally_relevant_solution.block(0), displacement_divergence);
1367
1368 for (unsigned int q = 0; q < fe_face_values.n_quadrature_points;
1369 ++q)
1370 {
1371 perimeter[boundary_index] += fe_face_values.JxW(q);
1372
1373 boundary_stress[boundary_index] +=
1374 (2 * par.Lame_mu * displacement_gradient[q] +
1375 par.Lame_lambda * displacement_divergence[q] *
1376 identity) *
1377 fe_face_values.JxW(q) * fe_face_values.normal_vector(q);
1378 u_dot_n[boundary_index] +=
1379 (displacement_values[q] *
1380 fe_face_values.normal_vector(q)) *
1381 fe_face_values.JxW(q);
1382 }
1383 }
1384 }
1385 }
1386
1387 if constexpr (spacedim == 2)
1388 {
1389 internal_stress = Utilities::MPI::sum(internal_stress, mpi_communicator);
1390 average_displacement =
1391 Utilities::MPI::sum(average_displacement, mpi_communicator);
1392 internal_area = Utilities::MPI::sum(internal_area, mpi_communicator);
1393
1394 internal_stress /= internal_area;
1395 average_displacement /= internal_area;
1396 }
1397 for (auto id : all_ids)
1398 {
1399 boundary_stress[id] =
1400 Utilities::MPI::sum(boundary_stress[id], mpi_communicator);
1401 perimeter[id] = Utilities::MPI::sum(perimeter[id], mpi_communicator);
1402 boundary_stress[id] /= perimeter[id];
1403 }
1404
1405 if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
1406 {
1407 const std::string filename(par.output_directory + "/forces.txt");
1408 std::ofstream forces_file;
1409 if (openfilefirsttime)
1410 {
1411 forces_file.open(filename);
1412 if constexpr (spacedim == 2)
1413 {
1414 forces_file
1415 << "cycle area perimeter meanInternalStressxx meanInternalStressxy meanInternalStressyx meanInternalStressyy avg_u_x avg_u_y";
1416 for (auto id : all_ids)
1417 forces_file << " boundaryStressX_" << id << " boundaryStressY_"
1418 << id << " uDotN_" << id;
1419 forces_file << std::endl;
1420 }
1421 else
1422 {
1423 forces_file << "cycle perimeter";
1424 for (auto id : all_ids)
1425 forces_file << " sigmanX_" << id << " sigmanY_"
1426 << " sigmanZ_" << id << " uDotN_" << id;
1427 forces_file << std::endl;
1428 }
1429 }
1430 else
1431 forces_file.open(filename, std::ios_base::app);
1432
1433 if constexpr (spacedim == 2)
1434 {
1435 forces_file << cycle << " " << internal_area << " ";
1436 for (auto id : all_ids)
1437 forces_file << perimeter[id] << " ";
1438 forces_file << internal_stress << " " << average_displacement << " ";
1439 for (auto id : all_ids)
1440 forces_file << boundary_stress[id] << " " << u_dot_n[id] << " ";
1441 forces_file << std::endl;
1442 }
1443 else // spacedim = 3
1444 {
1445 forces_file << cycle << " ";
1446 for (auto id : all_ids)
1447 forces_file << perimeter[id] << " ";
1448 for (auto id : all_ids)
1449 forces_file << boundary_stress[id] << " " << u_dot_n[id] << " ";
1450 forces_file << std::endl;
1451 }
1452 forces_file.close();
1453 }
1454
1455 return;
1456}
1457
1458template <int dim, int spacedim>
1459// TrilinosWrappers::MPI::Vector
1460void
1461CoupledElasticityProblem<dim, spacedim>::compute_coupling_pressure() /*const*/
1462{
1463 TimerOutput::Scope t(computing_timer, "Postprocessing: Computing Pressure");
1464 if (inclusions.n_inclusions() > 0 &&
1465 inclusions.get_offset_coefficients() == 1 &&
1466 inclusions.get_n_coefficients() >= 2)
1467 {
1468 const auto locally_owned_vessels =
1469 Utilities::MPI::create_evenly_distributed_partitioning(
1470 mpi_communicator, inclusions.get_n_vessels());
1471 // mpi_communicator,
1472 // std::min(Utilities::MPI::this_n_processes(mpi_communicator),
1473 // inclusions.get_n_vessels()), inclusions.get_n_vessels());
1474 const auto locally_owned_inclusions =
1475 Utilities::MPI::create_evenly_distributed_partitioning(
1476 mpi_communicator, inclusions.n_inclusions());
1477
1478 // TrilinosWrappers::MPI::Vector pressure(locally_owned_vessels,
1479 // mpi_communicator);
1480 coupling_pressure.reinit(locally_owned_vessels, mpi_communicator);
1481 auto &pressure = coupling_pressure;
1482 pressure = 0;
1483 TrilinosWrappers::MPI::Vector pressure_at_inc(locally_owned_inclusions,
1484 mpi_communicator);
1485 // coupling_pressure_at_inclusions.reinit(locally_owned_inclusions,mpi_communicator);
1486 // auto & pressure_at_inc = coupling_pressure_at_inclusions;
1487 pressure_at_inc = 0;
1488
1489 const auto &lambda_to_pressure = locally_relevant_solution.block(1);
1490
1491 // TODO: set the weight in a smarter way
1492 // std::vector<double> weights(inclusions.n_inclusions(),
1493 // inclusions.get_h3D1D());
1494 TrilinosWrappers::MPI::Vector inclusions_to_divide_by(
1495 locally_owned_vessels, mpi_communicator);
1496 inclusions_to_divide_by = 0;
1497
1498 const auto used_number_modes = inclusions.get_n_coefficients();
1499
1500 const auto local_lambda = lambda_to_pressure.locally_owned_elements();
1501 if constexpr (spacedim == 3)
1502 {
1503 unsigned int previous_inclusion_number =
1504 numbers::invalid_unsigned_int;
1505 auto tensorR = inclusions.get_rotation(0);
1506 for (const auto &ll : local_lambda)
1507 {
1508 const unsigned inclusion_number = (unsigned int)floor(
1509 ll / (inclusions.get_n_coefficients() * spacedim));
1510
1511 auto lii = ll - inclusion_number *
1512 inclusions.get_n_coefficients() * spacedim;
1513 const unsigned mode_number = (unsigned int)floor(lii / spacedim);
1514 const unsigned coor_number = lii % spacedim;
1515
1516 if (previous_inclusion_number != inclusion_number)
1517 tensorR = inclusions.get_rotation(inclusion_number);
1518
1519 if (mode_number == 0 || mode_number == 1)
1520 {
1521 AssertIndexRange(inclusion_number, inclusions.n_inclusions());
1522 pressure[inclusions.get_vesselID(inclusion_number)] +=
1523 lambda_to_pressure[ll] * tensorR[coor_number][mode_number] /
1524 used_number_modes;
1525 // * weights[inclusion_number];
1526 // inclusions_to_divide_by[inclusions.get_vesselID(inclusion_number)]
1527 // += 1;
1528 pressure_at_inc[inclusion_number] +=
1529 lambda_to_pressure[ll] * tensorR[coor_number][mode_number] /
1530 used_number_modes;
1531 }
1532 previous_inclusion_number = inclusion_number;
1533 }
1534 pressure.compress(VectorOperation::add);
1535 pressure_at_inc.compress(VectorOperation::add);
1536 }
1537 else // spacedim = 2
1538 {
1539 for (auto ll : local_lambda)
1540 {
1541 const unsigned inclusion_number = (unsigned int)floor(
1542 ll / (inclusions.get_n_coefficients() * spacedim));
1543
1544 auto lii = ll - inclusion_number *
1545 (inclusions.get_n_coefficients() * spacedim);
1546 if (lii == 0 || lii == 3)
1547 {
1548 AssertIndexRange(inclusion_number, inclusions.n_inclusions());
1549 pressure[inclusions.get_vesselID(inclusion_number)] +=
1550 lambda_to_pressure[ll] / used_number_modes;
1551 }
1552 }
1553 pressure_at_inc = pressure;
1554 pressure.print(std::cout);
1555 local_lambda.print(std::cout);
1556 }
1557 Utilities::MPI::gather(mpi_communicator, pressure_at_inc, 0);
1558 coupling_pressure_at_inclusions = pressure_at_inc;
1559
1560 output_coupling_pressure(cycle == 1 ? true : false);
1561 }
1562}
1563
1564template <int dim, int spacedim>
1565// TrilinosWrappers::MPI::Vector
1566void
1567CoupledElasticityProblem<dim, spacedim>::output_coupling_pressure(
1568 bool openfilefirsttime) const
1569{
1570 TimerOutput::Scope t(computing_timer, "Postprocessing: output Pressure");
1571 if (inclusions.n_inclusions() > 0 &&
1572 inclusions.get_offset_coefficients() == 1 &&
1573 inclusions.get_n_coefficients() >= 2)
1574 {
1575 const auto &pressure = coupling_pressure;
1576 // print .txt only sequential
1577 if (Utilities::MPI::n_mpi_processes(mpi_communicator) == 1)
1578 {
1579 const std::string filename(par.output_directory +
1580 "/externalPressure.txt");
1581 std::ofstream pressure_file;
1582 if (openfilefirsttime)
1583 pressure_file.open(filename);
1584 else
1585 pressure_file.open(filename, std::ios_base::app);
1586 // pressure_file << cycle << " ";
1587 for (unsigned int in = 0; in < pressure.size(); ++in)
1588 pressure_file << in << " " << pressure[in] << std::endl;
1589 // pressure.print(pressure_file);
1590 pressure_file.close();
1591 }
1592 else
1593 // print .h5
1594 if (par.initial_time == par.final_time)
1595 {
1596# ifdef DEAL_II_WITH_HDF5
1597 const std::string FILE_NAME(par.output_directory +
1598 "/externalPressure.h5");
1599
1600 auto accessMode = HDF5::File::FileAccessMode::create;
1601 if (!openfilefirsttime)
1602 accessMode = HDF5::File::FileAccessMode::open;
1603
1604 HDF5::File file_h5(FILE_NAME, accessMode, mpi_communicator);
1605 const std::string DATASET_NAME("externalPressure_" +
1606 std::to_string(cycle));
1607
1608 HDF5::DataSet dataset =
1609 file_h5.create_dataset<double>(DATASET_NAME,
1610 {inclusions.get_n_vessels()});
1611
1612 std::vector<double> data_to_write;
1613 // std::vector<hsize_t> coordinates;
1614 data_to_write.reserve(pressure.locally_owned_size());
1615 // coordinates.reserve(pressure.locally_owned_size());
1616 const auto locally_owned_vessels =
1617 Utilities::MPI::create_evenly_distributed_partitioning(
1618 mpi_communicator, inclusions.get_n_vessels());
1619
1620 for (const auto &el : locally_owned_vessels)
1621 {
1622 // coordinates.emplace_back(el);
1623 data_to_write.emplace_back(pressure[el]);
1624 }
1625 if (pressure.locally_owned_size() > 0)
1626 {
1627 hsize_t prefix = 0;
1628 hsize_t los = pressure.locally_owned_size();
1629 int ierr = MPI_Exscan(&los,
1630 &prefix,
1631 1,
1632 MPI_UNSIGNED_LONG_LONG,
1633 MPI_SUM,
1634 mpi_communicator);
1635 AssertThrowMPI(ierr);
1636
1637 std::vector<hsize_t> offset = {prefix, 1};
1638 std::vector<hsize_t> count = {pressure.locally_owned_size(), 1};
1639 // data.write_selection(data_to_write, coordinates);
1640 dataset.write_hyperslab(data_to_write, offset, count);
1641 }
1642 else
1643 dataset.write_none<int>();
1644# else
1645
1646 AssertThrow(false, ExcNeedsHDF5());
1647# endif
1648 }
1649 else
1650 {
1651 pcout
1652 << "output_pressure file for time dependent simulation not implemented"
1653 << std::endl;
1654 }
1655 // // coupling_pressure = pressure;
1656 // return pressure;
1657 }
1658 // else
1659 // {
1660 // pcout
1661 // << "inclusions parameters ('Start index of Fourier coefficients' or
1662 // 'Number of fourier coefficients') not compatible with the computation
1663 // of the pressure as intended, pressure.hdf5 not generated"
1664 // << std::endl;
1665 // }
1666 // TrilinosWrappers::MPI::Vector temp;
1667 // return temp;
1668}
1669
1670template <int dim, int spacedim>
1671std::vector<std::vector<double>>
1672CoupledElasticityProblem<dim, spacedim>::split_pressure_over_inclusions(
1673 std::vector<int> number_of_cells_per_vessel,
1674 Vector<double> /* full_press */) const
1675{
1676 Assert(number_of_cells_per_vessel.size() == inclusions.get_n_vessels(),
1677 ExcInternalError());
1678
1679 std::vector<std::vector<double>> split_pressure;
1680 unsigned starting_inclusion = 0;
1681
1682 for (unsigned int vessel = 0; vessel < number_of_cells_per_vessel.size();
1683 ++vessel)
1684 {
1685 auto N2 = number_of_cells_per_vessel[vessel];
1686 auto N1 = inclusions.get_inclusions_in_vessel(vessel);
1687
1688 std::vector<double> new_vector;
1689
1690 // const std::vector<double> pressure_of_inc_in_vessel(pressure);//
1691 // (*pressure[starting_inclusion], *pressure[starting_inclusion+N1]);
1692 const Vector<double> pressure_of_inc_in_vessel(
1693 coupling_pressure_at_inclusions);
1694
1695 new_vector.push_back(pressure_of_inc_in_vessel[starting_inclusion]);
1696
1697 for (auto cell = 1; cell < N2 - 1; ++cell)
1698 {
1699 auto X = cell / (N2 - 1) * (N1 - 1);
1700 auto j = floor(X);
1701 auto w = X - j;
1702
1703 new_vector.push_back(
1704 (1 - w) * pressure_of_inc_in_vessel[starting_inclusion + j] +
1705 (w)*pressure_of_inc_in_vessel[starting_inclusion + j + 1]);
1706 }
1707 new_vector.push_back(
1708 pressure_of_inc_in_vessel[starting_inclusion + N1 - 1]);
1709 starting_inclusion += N1;
1710 split_pressure.push_back(new_vector);
1711 }
1712
1713 return split_pressure;
1714}
1715
1716
1717template <int dim, int spacedim>
1718void
1719CoupledElasticityProblem<dim, spacedim>::run()
1720{
1721 if (par.initial_time == par.final_time) // time stationary
1722 {
1723 print_parameters();
1724 make_grid();
1725 setup_fe();
1726 check_boundary_ids();
1727 {
1728 TimerOutput::Scope t(computing_timer, "Setup inclusion");
1729 inclusions.setup_inclusions_particles(tria);
1730 }
1731 // setup_dofs(); // called inside refine_and_transfer
1732 for (cycle = 0; cycle < par.n_refinement_cycles; ++cycle)
1733 {
1734 setup_dofs();
1735
1736 assemble_elasticity_system();
1737
1738 assemble_coupling();
1739 solve();
1740 if (par.output_results)
1741 output_results();
1742
1743 {
1744 if constexpr (spacedim == 2)
1745 {
1746 FunctionParser<spacedim> weight(par.weight_expression);
1747 par.convergence_table.error_from_exact(
1748 dh,
1749 locally_relevant_solution.block(0),
1750 par.exact_solution,
1751 &weight);
1752 }
1753 else
1754 par.convergence_table.error_from_exact(
1755 dh, locally_relevant_solution.block(0), par.bc);
1756 par.convergence_table.output_table(pcout.get_stream());
1757 }
1758 if constexpr (spacedim == 2)
1759 {
1760 // output_pressure(cycle == 0 ? true : false);
1761 compute_coupling_pressure();
1762 output_coupling_pressure(cycle == 0 ? true : false);
1763 }
1764
1765 if (cycle != par.n_refinement_cycles - 1)
1766 refine_and_transfer();
1767 }
1768 // output_pressure(true);
1769 compute_coupling_pressure();
1770 output_coupling_pressure(true);
1771
1772 if (par.domain_type == "generate")
1773 compute_internal_and_boundary_stress(true);
1774 }
1775 else // Time dependent simulation
1776 {
1777 // TODO: add refinement as the first cycle,
1778 pcout << "time dependent simulation, refinement not implemented"
1779 << std::endl;
1780 print_parameters();
1781 make_grid();
1782 setup_fe();
1783 check_boundary_ids();
1784 cycle = 0;
1785 {
1786 TimerOutput::Scope t(computing_timer, "Setup inclusion");
1787 inclusions.setup_inclusions_particles(tria);
1788 }
1789 setup_dofs();
1790 assemble_elasticity_system();
1791 for (current_time = par.initial_time; current_time < par.final_time;
1792 current_time += par.dt, ++cycle)
1793 {
1794 pcout << "Time: " << current_time << std::endl;
1795 // assemble_elasticity_system();
1796 inclusions.inclusions_rhs.set_time(current_time);
1797 par.Neumann_bc.set_time(current_time);
1798 assemble_coupling();
1799 solve();
1800
1801 if (par.output_results)
1802 output_results();
1803 // output_pressure(cycle == 0 ? true : false);
1804 compute_coupling_pressure();
1805 output_coupling_pressure(cycle == 0 ? true : false);
1806
1807 if (par.domain_type == "generate")
1808 compute_internal_and_boundary_stress(cycle == 0 ? true : false);
1809 }
1810 }
1811}
1812
1813template <int dim, int spacedim>
1814void
1815CoupledElasticityProblem<dim, spacedim>::run_timestep0()
1816{
1817 print_parameters();
1818 make_grid();
1819 setup_fe();
1820 check_boundary_ids();
1821 {
1822 TimerOutput::Scope t(computing_timer, "Setup inclusion");
1823 inclusions.setup_inclusions_particles(tria);
1824 }
1825 cycle = 0;
1826 setup_dofs();
1827 // assemble_elasticity_system();
1828}
1829
1830template <int dim, int spacedim>
1831void
1832CoupledElasticityProblem<dim, spacedim>::run_timestep()
1833{
1834 if (cycle == 0) // at first timestep we refine
1835 {
1836 if (par.refinement_strategy == "inclusions")
1837 {
1838 refine_and_transfer_around_inclusions();
1839 std::cout << "refining around inclusions" << std::endl;
1840
1841 assemble_elasticity_system(); // questo mi serve perche sto raffinando
1842 assemble_coupling();
1843 solve();
1844 }
1845 else
1846 {
1847 for (unsigned int ref_cycle = 0; ref_cycle < par.n_refinement_cycles;
1848 ++ref_cycle)
1849 {
1850 assemble_elasticity_system(); // questo mi serve perche sto
1851 // raffinando
1852 assemble_coupling();
1853 solve();
1854 if (ref_cycle != par.n_refinement_cycles - 1)
1855 refine_and_transfer();
1856 }
1857 }
1858 }
1859 else
1860 {
1861 reassemble_coupling_rhs();
1862 solve();
1863 }
1864
1865
1866 if (par.output_results)
1867 output_results();
1868
1869 coupling_pressure.clear();
1870 coupling_pressure_at_inclusions.clear();
1871 // coupling_pressure = output_pressure(cycle == 0 ? true : false);
1872 // compute_coupling_pressure();
1873 // output_coupling_pressure(cycle == 0 ? true : false);
1874
1875 compute_internal_and_boundary_stress(cycle == 0 ? true : false);
1876 cycle++;
1877}
1878
1879template <int dim, int spacedim>
1880void
1881CoupledElasticityProblem<dim, spacedim>::update_inclusions_data(
1882 std::vector<double> new_data)
1883{
1884 inclusions.update_inclusions_data(new_data);
1885}
1886
1887template <int dim, int spacedim>
1888void
1889CoupledElasticityProblem<dim, spacedim>::update_inclusions_data(
1890 std::vector<double> new_data,
1891 std::vector<int> cells_per_vessel)
1892{
1893 Assert(cells_per_vessel.size() == inclusions.get_n_vessels(),
1894 ExcInternalError());
1895 std::vector<std::vector<double>> full_vector;
1896 unsigned int starting_point = 0;
1897 for (auto &N1 : cells_per_vessel)
1898 {
1899 std::vector<double> vessel_vector;
1900 for (auto j = 0; j < N1; ++j)
1901 {
1902 AssertIndexRange(starting_point + j, new_data.size());
1903 vessel_vector.push_back(new_data[starting_point + j]);
1904 }
1905 starting_point += N1;
1906
1907 full_vector.push_back(vessel_vector);
1908 }
1909 inclusions.update_inclusions_data(full_vector);
1910}
1911
1912// Template instantiations
1913template class CoupledElasticityProblemParameters<2>;
1914template class CoupledElasticityProblemParameters<2, 3>;
1915template class CoupledElasticityProblemParameters<3>;
1916
1917template class CoupledElasticityProblem<2>;
1918template class CoupledElasticityProblem<2, 3>; // dim != spacedim
1919template class CoupledElasticityProblem<3>;
1920#endif
DataOut::write_vtu_in_parallel
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
DataOut::type_dof_data
type_dof_data
DataOut::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
DataOut::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
FE_Q
GridIn::attach_triangulation
void attach_triangulation(Triangulation< dim, spacedim > &tria)
GridIn::read_msh
void read_msh(std::istream &in)
GridIn::read
void read(std::istream &in, Format format=Default)
HDF5::DataSet::write_hyperslab
void write_hyperslab(const Container &data, const std::vector< hsize_t > &offset, const std::vector< hsize_t > &count)
HDF5::DataSet::write_none
void write_none()
KellyErrorEstimator::estimate
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
Triangulation::locally_owned_subdomain
virtual types::subdomain_id locally_owned_subdomain() const
Triangulation::n_active_cells
unsigned int n_active_cells() const
Triangulation::refine_global
void refine_global(const unsigned int times=1)
Triangulation::prepare_coarsening_and_refinement
bool prepare_coarsening_and_refinement()
Triangulation::execute_coarsening_and_refinement
virtual void execute_coarsening_and_refinement()
Triangulation::get_boundary_ids
virtual std::vector< types::boundary_id > get_boundary_ids() const
coupled_elasticity.h
read_grid_and_cad_files
void read_grid_and_cad_files(const std::string &grid_file_name, const std::string &ids_and_cad_file_names, Triangulation< dim, spacedim > &tria)
Definition elasticity.cc:46
grid_tools.h
Triangulation::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Assert
#define Assert(cond, exc)
AssertThrowMPI
#define AssertThrowMPI(error_code)
AssertIndexRange
#define AssertIndexRange(index, range)
AssertThrow
#define AssertThrow(cond, exc)
linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
transpose_operator
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
inverse_operator
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
VectorTools::compute_nonzero_normal_flux_constraints
void compute_nonzero_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, const std::map< types::boundary_id, const Function< spacedim, double > * > &function_map, AffineConstraints< double > &constraints, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const bool use_manifold_for_normal=true)
update_values
update_values
update_normal_vectors
update_normal_vectors
update_JxW_values
update_JxW_values
update_gradients
update_gradients
update_quadrature_points
update_quadrature_points
Triangulation::set_manifold
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
complete_index_set
IndexSet complete_index_set(const IndexSet::size_type N)
InitializeLibrary::MPI
@ MPI
data
std::vector< index_type > data
DataOutBase::write_pvd_record
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > &times_and_names)
floor
Expression floor(const Expression &x)
DoFTools::extract_locally_relevant_dofs
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
DoFTools::extract_constant_modes
std::vector< std::vector< bool > > extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask={})
GridGenerator::generate_from_name_and_arguments
void generate_from_name_and_arguments(Triangulation< dim, spacedim > &tria, const std::string &grid_generator_function_name, const std::string &grid_generator_function_arguments)
GridGenerator::cheese
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
GridTools::find_cells_adjacent_to_vertex
std::vector< typename MeshType::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType &container, const unsigned int vertex_index)
LAPACKSupport::A
constexpr char A
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
OpenCASCADE::count_elements
std::tuple< unsigned int, unsigned int, unsigned int > count_elements(const TopoDS_Shape &shape)
OpenCASCADE::read_STEP
TopoDS_Shape read_STEP(const std::string &filename, const double scale_factor=1e-3)
OpenCASCADE::read_IGES
TopoDS_Shape read_IGES(const std::string &filename, const double scale_factor=1e-3)
w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Utilities::MPI::create_evenly_distributed_partitioning
IndexSet create_evenly_distributed_partitioning(const MPI_Comm comm, const types::global_dof_index total_size)
Utilities::MPI::gather
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
Utilities::dim_string
std::string dim_string(const int dim, const int spacedim)
VectorTools::interpolate
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::L1_norm)
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
parallel::distributed::SolutionTransfer
::SolutionTransfer< dim, VectorType, spacedim > SolutionTransfer
types::manifold_id
unsigned int manifold_id
types::material_id
unsigned int material_id
types::boundary_id
unsigned int boundary_id
identity
std_cxx20::type_identity< T > identity