--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2000 - 2020 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+
+// A slightly modified version of step-8 that tests the postprocessor
+// discussed in the documentation of DataPostprocessorTensor
+//
+// Unlike the _01 test, this one outputs in VTU format, which
+// uncovered a bug in the implementation of DataPostprocessorTensor,
+// namely that the components of the output field were not actually
+// described as one tensor, but dim*dim independent scalars.
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step8
+{
+ template <int dim>
+ class ElasticProblem
+ {
+ public:
+ ElasticProblem();
+ ~ElasticProblem();
+ void
+ run();
+
+ private:
+ void
+ setup_system();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ refine_grid();
+ void
+ output_results(const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+
+ FESystem<dim> fe;
+
+ AffineConstraints<double> hanging_node_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+ };
+
+
+
+ template <int dim>
+ void
+ right_hand_side(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<1, dim>> & values)
+ {
+ Assert(values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+ Assert(dim >= 2, ExcNotImplemented());
+
+ Point<dim> point_1, point_2;
+ point_1(0) = 0.5;
+ point_2(0) = -0.5;
+
+ for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
+ {
+ if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
+ ((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
+ values[point_n][0] = 1.0;
+ else
+ values[point_n][0] = 0.0;
+
+ if (points[point_n].norm_square() < 0.2 * 0.2)
+ values[point_n][1] = 1.0;
+ else
+ values[point_n][1] = 0.0;
+ }
+ }
+
+
+
+ template <int dim>
+ ElasticProblem<dim>::ElasticProblem()
+ : dof_handler(triangulation)
+ , fe(FE_Q<dim>(1), dim)
+ {}
+
+
+
+ template <int dim>
+ ElasticProblem<dim>::~ElasticProblem()
+ {
+ dof_handler.clear();
+ }
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::setup_system()
+ {
+ dof_handler.distribute_dofs(fe);
+ hanging_node_constraints.clear();
+ DoFTools::make_hanging_node_constraints(dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close();
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler,
+ dsp,
+ hanging_node_constraints,
+ /*keep_constrained_dofs = */ true);
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+ }
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::assemble_system()
+ {
+ QGauss<dim> quadrature_formula(2);
+
+ FEValues<dim> fe_values(fe,
+ quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+ std::vector<double> lambda_values(n_q_points);
+ std::vector<double> mu_values(n_q_points);
+
+ Functions::ConstantFunction<dim> lambda(1.), mu(1.);
+
+ std::vector<Tensor<1, dim>> rhs_values(n_q_points);
+
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ fe_values.reinit(cell);
+
+ lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
+ mu.value_list(fe_values.get_quadrature_points(), mu_values);
+ right_hand_side(fe_values.get_quadrature_points(), rhs_values);
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ const unsigned int component_j =
+ fe.system_to_component_index(j).first;
+
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ cell_matrix(i, j) +=
+ ((fe_values.shape_grad(i, q_point)[component_i] *
+ fe_values.shape_grad(j, q_point)[component_j] *
+ lambda_values[q_point]) +
+ (fe_values.shape_grad(i, q_point)[component_j] *
+ fe_values.shape_grad(j, q_point)[component_i] *
+ mu_values[q_point]) +
+ ((component_i == component_j) ?
+ (fe_values.shape_grad(i, q_point) *
+ fe_values.shape_grad(j, q_point) *
+ mu_values[q_point]) :
+ 0)) *
+ fe_values.JxW(q_point);
+ }
+ }
+ }
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ cell_rhs(i) += fe_values.shape_value(i, q_point) *
+ rhs_values[q_point][component_i] *
+ fe_values.JxW(q_point);
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i, j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ hanging_node_constraints.condense(system_matrix);
+ hanging_node_constraints.condense(system_rhs);
+
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 0,
+ Functions::ZeroFunction<dim>(dim),
+ boundary_values);
+ MatrixTools::apply_boundary_values(boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+ }
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::solve()
+ {
+ SolverControl solver_control(1000, 1e-12);
+ SolverCG<> cg(solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ cg.solve(system_matrix, solution, system_rhs, preconditioner);
+
+ hanging_node_constraints.distribute(solution);
+ }
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::refine_grid()
+ {
+ Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::estimate(
+ dof_handler,
+ QGauss<dim - 1>(2),
+ std::map<types::boundary_id, const Function<dim> *>(),
+ solution,
+ estimated_error_per_cell);
+
+ GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+ estimated_error_per_cell,
+ 0.3,
+ 0.03);
+
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+
+
+ template <int dim>
+ class StrainPostprocessor : public DataPostprocessorTensor<dim>
+ {
+ public:
+ StrainPostprocessor()
+ : DataPostprocessorTensor<dim>("strain", update_gradients)
+ {}
+
+ virtual void
+ evaluate_vector_field(
+ const DataPostprocessorInputs::Vector<dim> &input_data,
+ std::vector<Vector<double>> & computed_quantities) const
+ {
+ AssertDimension(input_data.solution_gradients.size(),
+ computed_quantities.size());
+
+ for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
+ {
+ AssertDimension(computed_quantities[p].size(),
+ (Tensor<2, dim>::n_independent_components));
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ computed_quantities[p]
+ [Tensor<2, dim>::component_to_unrolled_index(
+ TableIndices<2>(d, e))] =
+ (input_data.solution_gradients[p][d][e] +
+ input_data.solution_gradients[p][e][d]) /
+ 2;
+ }
+ }
+ };
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::output_results(const unsigned int cycle) const
+ {
+ std::string filename = "solution-";
+ filename += ('0' + cycle);
+ Assert(cycle < 10, ExcInternalError());
+
+ StrainPostprocessor<dim> grad_u;
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+
+
+
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation(
+ dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_out.add_data_vector(solution,
+ std::vector<std::string>(dim, "displacement"),
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ data_out.add_data_vector(solution, grad_u);
+ data_out.build_patches();
+ data_out.write_vtu(deallog.get_file_stream());
+ }
+
+
+
+ template <int dim>
+ void
+ ElasticProblem<dim>::run()
+ {
+ for (unsigned int cycle = 0; cycle < 1; ++cycle)
+ {
+ deallog << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(2);
+ }
+ else
+ refine_grid();
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells() << std::endl;
+
+ setup_system();
+
+ deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system();
+ solve();
+ output_results(cycle);
+ }
+ }
+} // namespace Step8
+
+
+int
+main()
+{
+ initlog();
+
+ Step8::ElasticProblem<2> elastic_problem_2d;
+ elastic_problem_2d.run();
+}
projects / dealii.git / commitdiff