From cff7df1fef999b9e51c01ce51df293e4770b5b23 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 6 Aug 2020 17:27:35 -0600 Subject: [PATCH] Add a new test. --- .../data_out_postprocessor_tensor_02.cc | 423 ++++++++++++++++++ .../data_out_postprocessor_tensor_02.output | 40 ++ 2 files changed, 463 insertions(+) create mode 100644 tests/data_out/data_out_postprocessor_tensor_02.cc create mode 100644 tests/data_out/data_out_postprocessor_tensor_02.output diff --git a/tests/data_out/data_out_postprocessor_tensor_02.cc b/tests/data_out/data_out_postprocessor_tensor_02.cc new file mode 100644 index 0000000000..163157ad38 --- /dev/null +++ b/tests/data_out/data_out_postprocessor_tensor_02.cc @@ -0,0 +1,423 @@ +/* --------------------------------------------------------------------- + * + * 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 +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + +namespace Step8 +{ + template + 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 triangulation; + DoFHandler dof_handler; + + FESystem fe; + + AffineConstraints hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + + + template + void + right_hand_side(const std::vector> &points, + std::vector> & values) + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + Assert(dim >= 2, ExcNotImplemented()); + + Point 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 + ElasticProblem::ElasticProblem() + : dof_handler(triangulation) + , fe(FE_Q(1), dim) + {} + + + + template + ElasticProblem::~ElasticProblem() + { + dof_handler.clear(); + } + + + + template + void + ElasticProblem::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 + void + ElasticProblem::assemble_system() + { + QGauss quadrature_formula(2); + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + std::vector lambda_values(n_q_points); + std::vector mu_values(n_q_points); + + Functions::ConstantFunction lambda(1.), mu(1.); + + std::vector> rhs_values(n_q_points); + + typename DoFHandler::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 boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(dim), + boundary_values); + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + system_rhs); + } + + + + template + void + ElasticProblem::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 + void + ElasticProblem::refine_grid() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, + QGauss(2), + std::map *>(), + 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 + class StrainPostprocessor : public DataPostprocessorTensor + { + public: + StrainPostprocessor() + : DataPostprocessorTensor("strain", update_gradients) + {} + + virtual void + evaluate_vector_field( + const DataPostprocessorInputs::Vector &input_data, + std::vector> & 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 + void + ElasticProblem::output_results(const unsigned int cycle) const + { + std::string filename = "solution-"; + filename += ('0' + cycle); + Assert(cycle < 10, ExcInternalError()); + + StrainPostprocessor grad_u; + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + + + + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_out.add_data_vector(solution, + std::vector(dim, "displacement"), + DataOut::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 + void + ElasticProblem::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(); +} diff --git a/tests/data_out/data_out_postprocessor_tensor_02.output b/tests/data_out/data_out_postprocessor_tensor_02.output new file mode 100644 index 0000000000..738d772f10 --- /dev/null +++ b/tests/data_out/data_out_postprocessor_tensor_02.output @@ -0,0 +1,40 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 16 +DEAL:: Number of degrees of freedom: 50 +DEAL:cg::Starting value 0.279006 +DEAL:cg::Convergence step 11 value 3.55007e-13 + + + + + + + +AQAAAAADAAAAAwAAdQAAAA==eNqFkdENACEIQ7tZV2MzV7uED66QoiTG2tQXUCAOcmWJzl18mEzlZga/3jhhsjdO80YViwuXXdeM7HchPl4cLQ7tOJIJmt7HjOlTfD7e4Zg+Lpz8d/GDd07LjBnsO+iZQzuGZIJLD1zus2vLkUyU/gCxEj+h + + + + + +AQAAAAABAAAAAQAAaQAAAA==eNoNwwVSAkAAAMADpaQ7pEGluyT//yt2ZzaEECJ+GPXTmAnjJk2Z9suMWfPmLFi0bMmKVevWbNj025ZtO/bs2nfgyKFjf/zz14lT585cuHTtyo1b9+48ePTsyYv/3rx69+HLp2+r3gfh + + +AQAAAEAAAABAAAAAMQAAAA==eNpjYWBg4ABiHiAWAGIRIJYAYhkgVgBiFSDWAGIdIDYAYhMgtgBiGyB2AGIAM0ACIQ== + + +AQAAABAAAAAQAAAACwAAAA==eNrj5EQFAATYAJE= + + + + +AQAAAAADAAAAAwAAegAAAA==eNpjYCAMWLsbbVTfeNsQqybya6eNCuNSG2LNaWI3tw3MsrbBZw6yGhl1ddttRba2+NyDrJeJ4ZX1fNm51vjU41KDyxxkNyC7DZc5yGpwAVxqkM3EFQ7IanCFP6nhgMscZIArHJD1IscpLj/iShu4zKEkrSIDAI4LYGU= + +AQAAAAAJAAAACQAAIgIAAA==eNqdVr8vREEQfufH0ZAoRCSu4golpbenQiQqhYR/gUgkGi0V0ZHoUGklChHebCIhBFEqqCUS0YlCI+Ybs42bzW1um+/lvW+/mdmZnXlZlrbat9ddwOrHg6v/XqpVP2bkfUCL07f+Jt+Alk6KLayFry0XcKhUrf3/PnFUzvv2Lv9sMVo64Cxld/IeaOmkxtVZG5a9wP3utuI/Z6PjrDa7mMv+gBZnd+DPB6ClA583b3uFA4xxvk9KwgFOv87VxTVf/il6ljfkPdDSAWfnekU4QEsnJRdYLdl7HvCwUqmL/+T9gG7Pr8bwDLR0wFlbLROegZZOiq2U1ZIt0s11v+gALR1wDiv7cm4BLZ8nx7tkL9A6Q3Cmnk6FA2w9W3BGLvL74xGJHRjJaf75NSo+Ay0d+Pw4WJa9wBinUVxcq/7l+Un2Ai0dcLjWSWueYjWfci/YBqktYltk+OzY11x9ziM5dRyz19i9pYOa57Mr9AyLGIdz4DQXyKk3cuo4l6Q5NX0Gh2uCtDbI0kmNi2u00FotuFYpa3LxnSkCWjrwme+e0zvoYhy+wxd6ly/4LvtmbGlPoYCWDmqee5PXHuWtMwSHe5zXXue51+VGzyy4V3rtmaYOONxzvfZeUwc1z73baw+PchrdC54pnmeJ15li6oDDM4l0Npk68JlnG+mMM88QHJ6RpLOSeFY2lQud2RTQ0oHPPPtJ/wGinEZxpdhKWb+GgHOA + + + + -- 2.39.5