--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2007, 2008, 2010, 2011, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// a complete testcase by Denis Davydov for FEValuesExtractors::Tensor
+
+#include "../tests.h"
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <numerics/vector_tools.h>
+#include <numerics/matrix_tools.h>
+#include <numerics/data_out.h>
+#include <fe/fe_q.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+template <int dim>
+class MixedElastoPlasticity
+{
+public:
+ MixedElastoPlasticity(const unsigned int degree);
+ void run();
+private:
+ void make_grid_and_dofs();
+ void assemble_system();
+
+ const unsigned int degree;
+ const unsigned int n_stress_components; // components of stress
+ const unsigned int n_gamma_components; // scalar plastic multiplier
+
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+
+};
+
+template <int dim>
+MixedElastoPlasticity<dim>::MixedElastoPlasticity(const unsigned int degree):
+degree(degree),
+n_stress_components(dim*dim),
+n_gamma_components(1),
+fe( FE_Q<dim>(degree), n_stress_components,
+ FE_Q<dim>(degree), n_gamma_components),
+dof_handler(triangulation)
+{
+}
+
+template <int dim>
+void MixedElastoPlasticity<dim>::make_grid_and_dofs()
+{
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(0);
+
+ dof_handler.distribute_dofs(fe);
+
+ deallog << "Number of stress components: "
+ << n_stress_components
+ << "\tNumber of gamma components: "
+ << n_gamma_components << std::endl;
+
+ // stress -> 0 gamma -> 1
+ std::vector<unsigned int> block_component(n_stress_components + n_gamma_components, 1);
+ for (unsigned int ii = 0; ii < n_stress_components; ii++)
+ block_component[ii] = 0;
+
+ DoFRenumbering::component_wise(dof_handler);
+
+ // total number of dof per block component
+ std::vector<unsigned int> dofs_per_block(2);
+ DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, block_component);
+
+ const unsigned int n_stress_dof = dofs_per_block[0];
+ const unsigned int n_gamma_dof = dofs_per_block[1];
+
+ deallog << "Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << "Total number of cells: "
+ << triangulation.n_cells()
+ << std::endl
+ << "Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << " = (" << n_stress_dof << " + "
+ << n_gamma_dof << ")"
+ << std::endl;
+
+ // following step-22 use of simple compressed block sparsity pattern for efficiency
+ {
+ BlockCompressedSimpleSparsityPattern csp(2, 2);
+
+ csp.block(0, 0).reinit(n_stress_dof, n_stress_dof);
+ csp.block(1, 0).reinit(n_gamma_dof, n_stress_dof);
+ csp.block(0, 1).reinit(n_stress_dof, n_gamma_dof);
+ csp.block(1, 1).reinit(n_gamma_dof, n_gamma_dof);
+ csp.collect_sizes();
+
+ DoFTools::make_sparsity_pattern(dof_handler, csp);
+ sparsity_pattern.copy_from(csp);
+ }
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(2);
+ solution.block(0).reinit(n_stress_dof);
+ solution.block(1).reinit(n_gamma_dof);
+ solution.collect_sizes();
+
+ system_rhs.reinit(2);
+ system_rhs.block(0).reinit(n_stress_dof);
+ system_rhs.block(1).reinit(n_gamma_dof);
+ system_rhs.collect_sizes();
+}
+
+template <int dim>
+void MixedElastoPlasticity<dim>::assemble_system()
+{
+ QGauss<dim> quadrature_formula(1);
+
+ 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;
+ deallog << "dofs_per_cell: " << fe.dofs_per_cell << std::endl;
+ //return;
+ 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<unsigned int> local_dof_indices(dofs_per_cell);
+
+
+ const FEValuesExtractors::Tensor<2> stress_extr(0);//rank2
+ const FEValuesExtractors::Scalar gamma_extr(n_stress_components);
+
+ deallog << "fe.dofs_per_cell: " << fe.dofs_per_cell
+ << "\tquadrature_formula.size(): " << quadrature_formula.size()
+ << std::endl;
+
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+
+ std::vector<Tensor<1,dim>> div_values(n_q_points);
+ std::vector<Tensor<2,dim>> stress_values(n_q_points);
+
+ unsigned int cc = 0;
+ for (; cell!=endc; ++cell) //loop over all cells
+ {
+ deallog<<++cc<<" ";
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ fe_values.reinit (cell); //compute requested values for a given cell
+
+ fe_values[stress_extr].get_function_values(solution,stress_values);
+ fe_values[stress_extr].get_function_divergences(solution,div_values);
+
+ }
+}
+
+template <int dim>
+void MixedElastoPlasticity<dim>::run()
+{
+ make_grid_and_dofs();
+ assemble_system();
+ deallog << "here" << std::endl;
+
+}
+
+void check()
+{
+ static const unsigned int dim = 3;
+ {
+ for (unsigned int i = 0; i < dim; i++) {
+ for (unsigned int j = 0; j < dim; j++) {
+ TableIndices<2> indices(i,j);
+ unsigned int unrolled = Tensor<2,dim>::component_to_unrolled_index(indices);
+ deallog<<i<<" "<<j<<" -> "<<unrolled<<std::endl;
+ indices = Tensor<2,dim>::unrolled_to_component_indices(unrolled);
+ deallog<<unrolled<<" -> "<<indices[0]<<" "<<indices[1]<<std::endl;
+ }
+ }
+ }
+
+ MixedElastoPlasticity < 3 > elasto_plasticity(1);
+ elasto_plasticity.run();
+
+ deallog << "Analysis complete" << std::endl;
+
+}
+
+
+int main()
+{
+ std::ofstream logfile ("fe_values_view_tensor_01/output");
+ deallog << std::setprecision (3);
+
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+ deallog.threshold_double(1.e-7);
+
+ check ();
+}
--- /dev/null
+
+DEAL::FE=FESystem<2>[FE_Q<2>(1)-FE_RaviartThomas<2>(1)-FE_Nedelec<2>(1)]
+DEAL::component=0
+DEAL::0.807 4.36
+DEAL::0.807 5.29
+DEAL::1.31 4.07
+DEAL::1.31 5.59
+DEAL::component=1
+DEAL::18.2 352.
+DEAL::-18.2 337.
+DEAL::333. 180.
+DEAL::-295. 296.
+DEAL::component=2
+DEAL::43.2 -31.7
+DEAL::43.2 118.
+DEAL::125. -61.1
+DEAL::125. -40.8
+DEAL::component=3
+DEAL::2.58e-16 5.84
+DEAL::2.58e-16 5.84
+DEAL::3.58e-16 9.52
+DEAL::3.58e-16 9.52
+DEAL::component=4
+DEAL::5.84 -6.75
+DEAL::5.84 6.75
+DEAL::9.52 -11.0
+DEAL::9.52 11.0
+DEAL::FE=FESystem<3>[FE_Q<3>(1)-FE_RaviartThomas<3>(1)-FE_Nedelec<3>(1)]
+DEAL::component=0
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::2.37 4.73 9.46
+DEAL::component=1
+DEAL::8.18e-13 -275. -138.
+DEAL::-5.71e-13 -275. -138.
+DEAL::-367. -275. 4.71e+03
+DEAL::367. -275. 4.92e+03
+DEAL::-184. 4.57e+03 -138.
+DEAL::184. 4.79e+03 -138.
+DEAL::4.80e+03 4.57e+03 4.71e+03
+DEAL::-4.38e+03 4.79e+03 4.92e+03
+DEAL::component=2
+DEAL::-275. 1.30e-12 -138.
+DEAL::-275. -459. 4.80e+03
+DEAL::-275. -1.40e-12 -138.
+DEAL::-275. 459. 5.02e+03
+DEAL::4.67e+03 -91.8 -138.
+DEAL::4.67e+03 3.52e+03 4.80e+03
+DEAL::4.88e+03 91.8 -138.
+DEAL::4.88e+03 -3.09e+03 5.02e+03
+DEAL::component=3
+DEAL::-275. -138. 1.36e-12
+DEAL::-275. 4.90e+03 -367.
+DEAL::4.76e+03 -138. -184.
+DEAL::4.76e+03 4.90e+03 2.23e+03
+DEAL::-275. -138. -3.64e-12
+DEAL::-275. 5.11e+03 367.
+DEAL::4.97e+03 -138. 184.
+DEAL::4.97e+03 5.11e+03 -1.81e+03
+DEAL::component=4
+DEAL::-7.70e-16 5.60 22.4
+DEAL::-7.70e-16 5.60 22.4
+DEAL::-6.69e-16 5.60 22.4
+DEAL::-6.69e-16 5.60 22.4
+DEAL::1.58e-16 5.60 22.4
+DEAL::1.58e-16 5.60 22.4
+DEAL::9.52e-16 5.60 22.4
+DEAL::9.52e-16 5.60 22.4
+DEAL::component=5
+DEAL::5.60 5.78e-16 22.4
+DEAL::5.60 -1.39e-15 22.4
+DEAL::5.60 5.78e-16 22.4
+DEAL::5.60 -1.39e-15 22.4
+DEAL::5.60 -6.60e-16 22.4
+DEAL::5.60 3.53e-16 22.4
+DEAL::5.60 -6.60e-16 22.4
+DEAL::5.60 3.53e-16 22.4
+DEAL::component=6
+DEAL::5.60 11.2 -1.38e-15
+DEAL::5.60 11.2 1.96e-16
+DEAL::5.60 11.2 1.23e-17
+DEAL::5.60 11.2 1.50e-15
+DEAL::5.60 11.2 -1.38e-15
+DEAL::5.60 11.2 1.96e-16
+DEAL::5.60 11.2 1.23e-17
+DEAL::5.60 11.2 1.50e-15