From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 6 Aug 2020 23:27:35 +0000 (-0600)
Subject: Add a new test.
X-Git-Tag: v9.3.0-rc1~1208^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cff7df1fef999b9e51c01ce51df293e4770b5b23;p=dealii.git

Add a new test.
---

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 <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();
+}
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
+<?xml version="1.0" ?> 
+<!-- 
+# vtk DataFile Version 3.0
+#This file was generated 
+-->
+<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="64" NumberOfCells="16" >
+  <Points>
+    <DataArray type="Float32" NumberOfComponents="3" format="binary">
+AQAAAAADAAAAAwAAdQAAAA==eNqFkdENACEIQ7tZV2MzV7uED66QoiTG2tQXUCAOcmWJzl18mEzlZga/3jhhsjdO80YViwuXXdeM7HchPl4cLQ7tOJIJmt7HjOlTfD7e4Zg+Lpz8d/GDd07LjBnsO+iZQzuGZIJLD1zus2vLkUyU/gCxEj+h
+    </DataArray>
+  </Points>
+
+  <Cells>
+    <DataArray type="Int32" Name="connectivity" format="binary">
+AQAAAAABAAAAAQAAaQAAAA==eNoNwwVSAkAAAMADpaQ7pEGluyT//yt2ZzaEECJ+GPXTmAnjJk2Z9suMWfPmLFi0bMmKVevWbNj025ZtO/bs2nfgyKFjf/zz14lT585cuHTtyo1b9+48ePTsyYv/3rx69+HLp2+r3gfh
+    </DataArray>
+    <DataArray type="Int32" Name="offsets" format="binary">
+AQAAAEAAAABAAAAAMQAAAA==eNpjYWBg4ABiHiAWAGIRIJYAYhkgVgBiFSDWAGIdIDYAYhMgtgBiGyB2AGIAM0ACIQ==
+    </DataArray>
+    <DataArray type="UInt8" Name="types" format="binary">
+AQAAABAAAAAQAAAACwAAAA==eNrj5EQFAATYAJE=
+    </DataArray>
+  </Cells>
+  <PointData Scalars="scalars">
+    <DataArray type="Float32" Name="displacement" NumberOfComponents="3" format="binary">
+AQAAAAADAAAAAwAAegAAAA==eNpjYCAMWLsbbVTfeNsQqybya6eNCuNSG2LNaWI3tw3MsrbBZw6yGhl1ddttRba2+NyDrJeJ4ZX1fNm51vjU41KDyxxkNyC7DZc5yGpwAVxqkM3EFQ7IanCFP6nhgMscZIArHJD1IscpLj/iShu4zKEkrSIDAI4LYGU=    </DataArray>
+    <DataArray type="Float32" Name="strain_xx__strain_xy__strain_yx__strain_yy" NumberOfComponents="9" format="binary">
+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    </DataArray>
+  </PointData>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>