]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added functions to compute rotation matrices in 2d and 3d.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Mon, 9 Jan 2017 17:06:18 +0000 (18:06 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 11 Jan 2017 16:11:14 +0000 (17:11 +0100)
doc/news/changes/minor/20170109Jean-PaulPelteret [new file with mode: 0644]
include/deal.II/physics/transformations.h
tests/physics/step-18-rotation_matrix.cc [new file with mode: 0644]
tests/physics/step-18-rotation_matrix.with_petsc=true.output [new file with mode: 0644]
tests/physics/step-18.cc [new file with mode: 0644]
tests/physics/step-18.with_petsc=true.output [new file with mode: 0644]
tests/physics/transformations-rotations_01.cc [new file with mode: 0644]
tests/physics/transformations-rotations_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170109Jean-PaulPelteret b/doc/news/changes/minor/20170109Jean-PaulPelteret
new file mode 100644 (file)
index 0000000..9bddb7b
--- /dev/null
@@ -0,0 +1,4 @@
+New: Functions to compute the 2d and 3d rotation matrix have been implemented
+and can be found in the Physics::Transformations::Rotations namespace.
+<br>
+(Jean-Paul Pelteret, 2017/01/09)
index b77d6d3cc64350a1246bb74e44fce909bef4b04a..e7e6ebff06b3734fcf7c1b48d3ead42e5e687daf 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2016 by the deal.II authors
+// Copyright (C) 2016 - 2017 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -16,6 +16,7 @@
 #ifndef dealii__transformations_h
 #define dealii__transformations_h
 
+#include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/symmetric_tensor.h>
 
@@ -28,6 +29,75 @@ namespace Physics
   namespace Transformations
   {
 
+    /**
+     * Transformation functions and tensors that are defined in terms of
+     * rotation angles and axes of rotation.
+     *
+     * @author Jean-Paul Pelteret, 2017
+    */
+    namespace Rotations
+    {
+      /**
+       * @name Rotation matrices
+       */
+//@{
+
+      /**
+      * Returns the rotation matrix for 2-d Euclidean space, namely
+      * @f[
+      *  \mathbf{R} := \left[ \begin{array}{cc}
+      *  cos(\theta) & sin(\theta) \\
+      *  -sin(\theta) & cos(\theta)
+      * \end{array}\right]
+      * @f]
+      * where $\theta$ is the rotation angle given in radians.
+      * In particular, this describes the counter-clockwise rotation of a vector
+      * relative to a
+      * <a href="http://mathworld.wolfram.com/RotationMatrix.html">fixed set of right-handed axes</a>.
+      *
+      * @param[in] angle The rotation angle (about the z-axis) in radians
+      */
+      template<typename Number>
+      Tensor<2,2,Number>
+      rotation_matrix_2d (const Number &angle);
+
+
+      /**
+      * Returns the rotation matrix for 3-d Euclidean space.
+      * Most concisely stated using the Rodrigues' rotation formula, this
+      * function returns the equivalent of
+      * @f[
+      *  \mathbf{R} := cos(\theta)\mathbf{I} + sin(\theta)\mathbf{W}
+      *              + (1-cos(\theta))\mathbf{u}\otimes\mathbf{u}
+      * @f]
+      * where $\mathbf{u}$ is the axial vector (an axial vector) and $\theta$
+      * is the rotation angle given in radians, $\mathbf{I}$ is the identity
+      * tensor and $\mathbf{W}$ is the skew symmetric tensor of $\mathbf{u}$.
+      *
+      * @dealiiWriggersA{374,9.194}
+      * This presents Rodrigues' rotation formula, but the implementation used
+      * in this function is described in this
+      * <a href="https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle">wikipedia link</a>.
+      * In particular, this describes the counter-clockwise rotation of a vector
+      * <a href="http://mathworld.wolfram.com/RotationMatrix.html">in a plane with its normal</a>.
+      * defined by the @p axis of rotation.
+      * An alternative implementation is discussed at
+      * <a href="https://www.gamedev.net/resources/_/technical/math-and-physics/do-we-really-need-quaternions-r1199">this link</a>,
+      * but is inconsistent (sign-wise) with the Rodrigues' rotation formula as
+      * it describes the rotation of a coordinate system.
+      *
+      * @param[in] axis  A unit vector that defines the axis of rotation
+      * @param[in] angle The rotation angle in radians
+      */
+      template<typename Number>
+      Tensor<2,3,Number>
+      rotation_matrix_3d (const Point<3,Number> &axis,
+                          const Number          &angle);
+
+//@}
+
+    }
+
     /**
      * Transformation of tensors that are defined in terms of a set of
      * contravariant bases. Rank-1 and rank-2 contravariant tensors
@@ -789,6 +859,55 @@ namespace internal
 
 
 
+template<typename Number>
+Tensor<2,2,Number>
+Physics::Transformations::Rotations::rotation_matrix_2d (const Number &angle)
+{
+  const double rotation[2][2]
+  = {{
+      std::cos(angle) , -std::sin(angle)
+    },
+    {
+      std::sin(angle) , std::cos(angle)
+    }
+  };
+  return Tensor<2,2> (rotation);
+}
+
+
+
+template<typename Number>
+Tensor<2,3,Number>
+Physics::Transformations::Rotations::rotation_matrix_3d (const Point<3,Number> &axis,
+                                                         const Number          &angle)
+{
+  Assert(std::abs(axis.norm() - 1.0) < 1e-9,
+         ExcMessage("The supplied axial vector is not a unit vector."));
+  const Number c = std::cos(angle);
+  const Number s = std::sin(angle);
+  const Number t = 1.-c;
+  const double rotation[3][3]
+  = {{
+      t *axis[0] *axis[0] + c,
+      t *axis[0] *axis[1] - s *axis[2],
+      t *axis[0] *axis[2] + s *axis[1]
+    },
+    {
+      t *axis[0] *axis[1] + s *axis[2],
+      t *axis[1] *axis[1] + c,
+      t *axis[1] *axis[2] - s *axis[0]
+    },
+    {
+      t *axis[0] *axis[2] - s *axis[1],
+      t *axis[1] *axis[2] + s *axis[0],
+      t *axis[2] *axis[2] + c
+    }
+  };
+  return Tensor<2,3,Number>(rotation);
+}
+
+
+
 template <int dim, typename Number>
 inline
 Tensor<1,dim,Number>
diff --git a/tests/physics/step-18-rotation_matrix.cc b/tests/physics/step-18-rotation_matrix.cc
new file mode 100644 (file)
index 0000000..27a58b0
--- /dev/null
@@ -0,0 +1,812 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2017 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 at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// This is a version of step-18 that employs the rotation matrices defined
+// in the Physics::Transformations::Angular namespace.
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/physics/transformations.h>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+namespace Step18
+{
+  using namespace dealii;
+  template <int dim>
+  struct PointHistory
+  {
+    SymmetricTensor<2,dim> old_stress;
+  };
+  template <int dim>
+  SymmetricTensor<4,dim>
+  get_stress_strain_tensor (const double lambda, const double mu)
+  {
+    SymmetricTensor<4,dim> tmp;
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
+        for (unsigned int k=0; k<dim; ++k)
+          for (unsigned int l=0; l<dim; ++l)
+            tmp[i][j][k][l] = (((i==k) && (j==l) ? mu : 0.0) +
+                               ((i==l) && (j==k) ? mu : 0.0) +
+                               ((i==j) && (k==l) ? lambda : 0.0));
+    return tmp;
+  }
+  template <int dim>
+  inline
+  SymmetricTensor<2,dim>
+  get_strain (const FEValues<dim> &fe_values,
+              const unsigned int   shape_func,
+              const unsigned int   q_point)
+  {
+    SymmetricTensor<2,dim> tmp;
+    for (unsigned int i=0; i<dim; ++i)
+      tmp[i][i] = fe_values.shape_grad_component (shape_func,q_point,i)[i];
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        tmp[i][j]
+          = (fe_values.shape_grad_component (shape_func,q_point,i)[j] +
+             fe_values.shape_grad_component (shape_func,q_point,j)[i]) / 2;
+    return tmp;
+  }
+  template <int dim>
+  inline
+  SymmetricTensor<2,dim>
+  get_strain (const std::vector<Tensor<1,dim> > &grad)
+  {
+    Assert (grad.size() == dim, ExcInternalError());
+    SymmetricTensor<2,dim> strain;
+    for (unsigned int i=0; i<dim; ++i)
+      strain[i][i] = grad[i][i];
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
+    return strain;
+  }
+  Tensor<2,2>
+  get_rotation_matrix (const std::vector<Tensor<1,2> > &grad_u)
+  {
+    const double curl = (grad_u[1][0] - grad_u[0][1]);
+    // Note: Here the negative angle suggests that we're computing the rotation
+    // of the coordinate system around a fixed point
+    const double angle = -std::atan (curl);
+    return Physics::Transformations::Rotations::rotation_matrix_2d(angle);
+  }
+  Tensor<2,3>
+  get_rotation_matrix (const std::vector<Tensor<1,3> > &grad_u)
+  {
+    const Point<3> curl (grad_u[2][1] - grad_u[1][2],
+                         grad_u[0][2] - grad_u[2][0],
+                         grad_u[1][0] - grad_u[0][1]);
+    const double tan_angle = std::sqrt(curl*curl);
+    // Note: Here the negative angle suggests that we're computing the rotation
+    // of the coordinate system around a fixed point
+    const double angle = -std::atan (tan_angle);
+    const Point<3> axis = curl/tan_angle;
+    return Physics::Transformations::Rotations::rotation_matrix_3d(axis,angle);
+  }
+  template <int dim>
+  class TopLevel
+  {
+  public:
+    TopLevel ();
+    ~TopLevel ();
+    void run ();
+  private:
+    void create_coarse_grid ();
+    void setup_system ();
+    void assemble_system ();
+    void solve_timestep ();
+    unsigned int solve_linear_problem ();
+    void output_results () const;
+    void do_initial_timestep ();
+    void do_timestep ();
+    void refine_initial_grid ();
+    void move_mesh ();
+    void setup_quadrature_point_history ();
+    void update_quadrature_point_history ();
+    parallel::shared::Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    ConstraintMatrix     hanging_node_constraints;
+    const QGauss<dim>          quadrature_formula;
+    std::vector<PointHistory<dim> > quadrature_point_history;
+    PETScWrappers::MPI::SparseMatrix system_matrix;
+    PETScWrappers::MPI::Vector       system_rhs;
+    Vector<double>                   incremental_displacement;
+    double       present_time;
+    double       present_timestep;
+    double       end_time;
+    unsigned int timestep_no;
+    MPI_Comm mpi_communicator;
+    const unsigned int n_mpi_processes;
+    const unsigned int this_mpi_process;
+    ConditionalOStream pcout;
+    std::vector<types::global_dof_index> local_dofs_per_process;
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+    unsigned int         n_local_cells;
+    static const SymmetricTensor<4,dim> stress_strain_tensor;
+    int monitored_vertex_first_dof;
+  };
+  template <int dim>
+  class BodyForce :  public Function<dim>
+  {
+  public:
+    BodyForce ();
+    virtual
+    void
+    vector_value (const Point<dim> &p,
+                  Vector<double>   &values) const;
+    virtual
+    void
+    vector_value_list (const std::vector<Point<dim> > &points,
+                       std::vector<Vector<double> >   &value_list) const;
+  };
+  template <int dim>
+  BodyForce<dim>::BodyForce ()
+    :
+    Function<dim> (dim)
+  {}
+  template <int dim>
+  inline
+  void
+  BodyForce<dim>::vector_value (const Point<dim> &/*p*/,
+                                Vector<double>   &values) const
+  {
+    Assert (values.size() == dim,
+            ExcDimensionMismatch (values.size(), dim));
+    const double g   = 9.81;
+    const double rho = 7700;
+    values = 0;
+    values(dim-1) = -rho * g;
+  }
+  template <int dim>
+  void
+  BodyForce<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<Vector<double> >   &value_list) const
+  {
+    const unsigned int n_points = points.size();
+    Assert (value_list.size() == n_points,
+            ExcDimensionMismatch (value_list.size(), n_points));
+    for (unsigned int p=0; p<n_points; ++p)
+      BodyForce<dim>::vector_value (points[p],
+                                    value_list[p]);
+  }
+  template <int dim>
+  class IncrementalBoundaryValues :  public Function<dim>
+  {
+  public:
+    IncrementalBoundaryValues (const double present_time,
+                               const double present_timestep);
+    virtual
+    void
+    vector_value (const Point<dim> &p,
+                  Vector<double>   &values) const;
+    virtual
+    void
+    vector_value_list (const std::vector<Point<dim> > &points,
+                       std::vector<Vector<double> >   &value_list) const;
+  private:
+    const double velocity;
+    const double present_time;
+    const double present_timestep;
+  };
+  template <int dim>
+  IncrementalBoundaryValues<dim>::
+  IncrementalBoundaryValues (const double present_time,
+                             const double present_timestep)
+    :
+    Function<dim> (dim),
+    velocity (.1),
+    present_time (present_time),
+    present_timestep (present_timestep)
+  {}
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::
+  vector_value (const Point<dim> &/*p*/,
+                Vector<double>   &values) const
+  {
+    Assert (values.size() == dim,
+            ExcDimensionMismatch (values.size(), dim));
+    values = 0;
+    values(2) = -present_timestep * velocity;
+  }
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::
+  vector_value_list (const std::vector<Point<dim> > &points,
+                     std::vector<Vector<double> >   &value_list) const
+  {
+    const unsigned int n_points = points.size();
+    Assert (value_list.size() == n_points,
+            ExcDimensionMismatch (value_list.size(), n_points));
+    for (unsigned int p=0; p<n_points; ++p)
+      IncrementalBoundaryValues<dim>::vector_value (points[p],
+                                                    value_list[p]);
+  }
+  template <int dim>
+  const SymmetricTensor<4,dim>
+  TopLevel<dim>::stress_strain_tensor
+    = get_stress_strain_tensor<dim> (/*lambda = */ 9.695e10,
+                                                   /*mu     = */ 7.617e10);
+  template <int dim>
+  TopLevel<dim>::TopLevel ()
+    :
+    triangulation(MPI_COMM_WORLD),
+    fe (FE_Q<dim>(1), dim),
+    dof_handler (triangulation),
+    quadrature_formula (2),
+    mpi_communicator (MPI_COMM_WORLD),
+    n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
+    this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
+    pcout (std::cout, this_mpi_process == 0),
+    monitored_vertex_first_dof(0)
+  {}
+  template <int dim>
+  TopLevel<dim>::~TopLevel ()
+  {
+    dof_handler.clear ();
+  }
+  template <int dim>
+  void TopLevel<dim>::run ()
+  {
+    present_time = 0;
+    present_timestep = 1;
+    end_time = 10;
+    timestep_no = 0;
+    do_initial_timestep ();
+    while (present_time < end_time)
+      do_timestep ();
+  }
+  template <int dim>
+  void TopLevel<dim>::create_coarse_grid ()
+  {
+    const double inner_radius = 0.8,
+                 outer_radius = 1;
+    GridGenerator::cylinder_shell (triangulation,
+                                   3, inner_radius, outer_radius);
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell=triangulation.begin_active();
+         cell!=triangulation.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+        if (cell->face(f)->at_boundary())
+          {
+            const Point<dim> face_center = cell->face(f)->center();
+            if (face_center[2] == 0)
+              cell->face(f)->set_boundary_id (0);
+            else if (face_center[2] == 3)
+              cell->face(f)->set_boundary_id (1);
+            else if (std::sqrt(face_center[0]*face_center[0] +
+                               face_center[1]*face_center[1])
+                     <
+                     (inner_radius + outer_radius) / 2)
+              cell->face(f)->set_boundary_id (2);
+            else
+              cell->face(f)->set_boundary_id (3);
+          }
+    static const CylindricalManifold<dim> cylindrical_manifold (2);
+    triangulation.set_all_manifold_ids(0);
+    triangulation.set_manifold (0, cylindrical_manifold);
+    // triangulation.refine_global (1);
+    setup_quadrature_point_history ();
+  }
+  template <int dim>
+  void TopLevel<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs);
+    n_local_cells
+      = GridTools::count_cells_with_subdomain_association (triangulation,
+                                                           triangulation.locally_owned_subdomain ());
+    local_dofs_per_process = dof_handler.n_locally_owned_dofs_per_processor();
+    hanging_node_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                             hanging_node_constraints);
+    hanging_node_constraints.close ();
+    DynamicSparsityPattern sparsity_pattern (locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
+                                     hanging_node_constraints, /*keep constrained dofs*/ false);
+    SparsityTools::distribute_sparsity_pattern (sparsity_pattern,
+                                                local_dofs_per_process,
+                                                mpi_communicator,
+                                                locally_relevant_dofs);
+    system_matrix.reinit (locally_owned_dofs,
+                          locally_owned_dofs,
+                          sparsity_pattern,
+                          mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs,mpi_communicator);
+    incremental_displacement.reinit (dof_handler.n_dofs());
+  }
+  template <int dim>
+  void TopLevel<dim>::assemble_system ()
+  {
+    system_rhs = 0;
+    system_matrix = 0;
+    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);
+    BodyForce<dim>      body_force;
+    std::vector<Vector<double> > body_force_values (n_q_points,
+                                                    Vector<double>(dim));
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0;
+          cell_rhs = 0;
+          fe_values.reinit (cell);
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              for (unsigned int q_point=0; q_point<n_q_points;
+                   ++q_point)
+                {
+                  const SymmetricTensor<2,dim>
+                  eps_phi_i = get_strain (fe_values, i, q_point),
+                  eps_phi_j = get_strain (fe_values, j, q_point);
+                  cell_matrix(i,j)
+                  += (eps_phi_i * stress_strain_tensor * eps_phi_j
+                      *
+                      fe_values.JxW (q_point));
+                }
+          const PointHistory<dim> *local_quadrature_points_data
+            = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+          body_force.vector_value_list (fe_values.get_quadrature_points(),
+                                        body_force_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 q_point=0; q_point<n_q_points; ++q_point)
+                {
+                  const SymmetricTensor<2,dim> &old_stress
+                    = local_quadrature_points_data[q_point].old_stress;
+                  cell_rhs(i) += (body_force_values[q_point](component_i) *
+                                  fe_values.shape_value (i,q_point)
+                                  -
+                                  old_stress *
+                                  get_strain (fe_values,i,q_point))
+                                 *
+                                 fe_values.JxW (q_point);
+                }
+            }
+          cell->get_dof_indices (local_dof_indices);
+          hanging_node_constraints
+          .distribute_local_to_global (cell_matrix, cell_rhs,
+                                       local_dof_indices,
+                                       system_matrix, system_rhs);
+        }
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+    FEValuesExtractors::Scalar z_component (dim-1);
+    std::map<types::global_dof_index,double> boundary_values;
+    VectorTools::
+    interpolate_boundary_values (dof_handler,
+                                 0,
+                                 ZeroFunction<dim> (dim),
+                                 boundary_values);
+    VectorTools::
+    interpolate_boundary_values (dof_handler,
+                                 1,
+                                 IncrementalBoundaryValues<dim>(present_time,
+                                                                present_timestep),
+                                 boundary_values,
+                                 fe.component_mask(z_component));
+    PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                        system_matrix, tmp,
+                                        system_rhs, false);
+    incremental_displacement = tmp;
+  }
+  template <int dim>
+  void TopLevel<dim>::solve_timestep ()
+  {
+    pcout << "    Assembling system..." << std::flush;
+    assemble_system ();
+    pcout << " norm of rhs is " << system_rhs.l2_norm()
+          << std::endl;
+    const unsigned int n_iterations = solve_linear_problem ();
+    pcout << "    Solver converged in " << n_iterations
+          << " iterations." << std::endl;
+    pcout << "    Updating quadrature point data..." << std::flush;
+    update_quadrature_point_history ();
+    pcout << std::endl;
+  }
+  template <int dim>
+  unsigned int TopLevel<dim>::solve_linear_problem ()
+  {
+    PETScWrappers::MPI::Vector
+    distributed_incremental_displacement (locally_owned_dofs,mpi_communicator);
+    distributed_incremental_displacement = incremental_displacement;
+    SolverControl           solver_control (dof_handler.n_dofs(),
+                                            1e-16*system_rhs.l2_norm(),
+                                            false,false);
+    PETScWrappers::SolverCG cg (solver_control,
+                                mpi_communicator);
+    PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
+    cg.solve (system_matrix, distributed_incremental_displacement, system_rhs,
+              preconditioner);
+    incremental_displacement = distributed_incremental_displacement;
+    hanging_node_constraints.distribute (incremental_displacement);
+    return solver_control.last_step();
+  }
+  template <int dim>
+  void TopLevel<dim>::output_results () const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    std::vector<std::string> solution_names;
+    switch (dim)
+      {
+      case 1:
+        solution_names.push_back ("delta_x");
+        break;
+      case 2:
+        solution_names.push_back ("delta_x");
+        solution_names.push_back ("delta_y");
+        break;
+      case 3:
+        solution_names.push_back ("delta_x");
+        solution_names.push_back ("delta_y");
+        solution_names.push_back ("delta_z");
+        break;
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+    data_out.add_data_vector (incremental_displacement,
+                              solution_names);
+    Vector<double> norm_of_stress (triangulation.n_active_cells());
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = triangulation.begin_active(),
+      endc = triangulation.end();
+      for (; cell!=endc; ++cell)
+        if (cell->is_locally_owned())
+          {
+            SymmetricTensor<2,dim> accumulated_stress;
+            for (unsigned int q=0;
+                 q<quadrature_formula.size();
+                 ++q)
+              accumulated_stress +=
+                reinterpret_cast<PointHistory<dim>*>(cell->user_pointer())[q]
+                .old_stress;
+            norm_of_stress(cell->active_cell_index())
+              = (accumulated_stress /
+                 quadrature_formula.size()).norm();
+          }
+        else
+          norm_of_stress(cell->active_cell_index()) = -1e+20;
+    }
+    data_out.add_data_vector (norm_of_stress, "norm_of_stress");
+    std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
+    GridTools::get_subdomain_association (triangulation, partition_int);
+    const Vector<double> partitioning(partition_int.begin(),
+                                      partition_int.end());
+    data_out.add_data_vector (partitioning, "partitioning");
+    data_out.build_patches ();
+    std::string filename = "solution-" + Utilities::int_to_string(timestep_no,4)
+                           + "." + Utilities::int_to_string(this_mpi_process,3)
+                           + ".vtu";
+    AssertThrow (n_mpi_processes < 1000, ExcNotImplemented());
+    std::ofstream output (filename.c_str());
+    data_out.write_vtu (output);
+    if (this_mpi_process==0)
+      {
+        std::vector<std::string> filenames;
+        for (unsigned int i=0; i<n_mpi_processes; ++i)
+          filenames.push_back ("solution-" + Utilities::int_to_string(timestep_no,4)
+                               + "." + Utilities::int_to_string(i,3)
+                               + ".vtu");
+        const std::string
+        visit_master_filename = ("solution-" +
+                                 Utilities::int_to_string(timestep_no,4) +
+                                 ".visit");
+        std::ofstream visit_master (visit_master_filename.c_str());
+        DataOutBase::write_visit_record (visit_master, filenames);
+        const std::string
+        pvtu_master_filename = ("solution-" +
+                                Utilities::int_to_string(timestep_no,4) +
+                                ".pvtu");
+        std::ofstream pvtu_master (pvtu_master_filename.c_str());
+        data_out.write_pvtu_record (pvtu_master, filenames);
+        static std::vector<std::pair<double,std::string> > times_and_names;
+        times_and_names.push_back (std::pair<double,std::string> (present_time, pvtu_master_filename));
+        std::ofstream pvd_output ("solution.pvd");
+        DataOutBase::write_pvd_record (pvd_output, times_and_names);
+      }
+  }
+  template <int dim>
+  void TopLevel<dim>::do_initial_timestep ()
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    pcout << "Timestep " << timestep_no << " at time " << present_time
+          << std::endl;
+    for (unsigned int cycle=0; cycle<2; ++cycle)
+      {
+        pcout << "  Cycle " << cycle << ':' << std::endl;
+        if (cycle == 0)
+          create_coarse_grid ();
+        else
+          refine_initial_grid ();
+        pcout << "    Number of active cells:       "
+              << triangulation.n_active_cells()
+              << " (by partition:";
+        for (unsigned int p=0; p<n_mpi_processes; ++p)
+          pcout << (p==0 ? ' ' : '+')
+                << (GridTools::
+                    count_cells_with_subdomain_association (triangulation,p));
+        pcout << ")" << std::endl;
+        setup_system ();
+        pcout << "    Number of degrees of freedom: "
+              << dof_handler.n_dofs()
+              << " (by partition:";
+        for (unsigned int p=0; p<n_mpi_processes; ++p)
+          pcout << (p==0 ? ' ' : '+')
+                << (DoFTools::
+                    count_dofs_with_subdomain_association (dof_handler,p));
+        pcout << ")" << std::endl;
+        // Get point at which to output displacement
+        // (outer radius of displaced surface)
+        {
+          const Point<dim> soln_pt (1.0, 0.0, 3.0);
+          typename DoFHandler<dim>::active_cell_iterator cell =
+            dof_handler.begin_active(), endc = dof_handler.end();
+          for (; cell != endc; ++cell)
+            for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+              if (cell->vertex(v).distance(soln_pt) < 1e-6)
+                {
+                  monitored_vertex_first_dof = cell->vertex_dof_index(v,0);
+                }
+        }
+        solve_timestep ();
+      }
+    move_mesh ();
+    // output_results ();
+    pcout << std::endl;
+  }
+  template <int dim>
+  void TopLevel<dim>::do_timestep ()
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    pcout << "Timestep " << timestep_no << " at time " << present_time
+          << std::endl;
+    if (present_time > end_time)
+      {
+        present_timestep -= (present_time - end_time);
+        present_time = end_time;
+      }
+    solve_timestep ();
+    move_mesh ();
+    // output_results ();
+    // Output displacement at edge of displaced surface
+    {
+      static Tensor<1,dim> soln;
+      for (unsigned int d=0; d<dim; ++d)
+        soln[d] += incremental_displacement(monitored_vertex_first_dof+d);
+      deallog
+          << "Timestep " << timestep_no
+          << ": " << soln << std::endl;
+    }
+    pcout << std::endl;
+  }
+  template <int dim>
+  void TopLevel<dim>::refine_initial_grid ()
+  {
+    Vector<float> error_per_cell (triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                        QGauss<dim-1>(2),
+                                        typename FunctionMap<dim>::type(),
+                                        incremental_displacement,
+                                        error_per_cell,
+                                        ComponentMask(),
+                                        0,
+                                        MultithreadInfo::n_threads(),
+                                        this_mpi_process);
+    const unsigned int n_local_cells = triangulation.n_locally_owned_active_cells ();
+    PETScWrappers::MPI::Vector
+    distributed_error_per_cell (mpi_communicator,
+                                triangulation.n_active_cells(),
+                                n_local_cells);
+    for (unsigned int i=0; i<error_per_cell.size(); ++i)
+      if (error_per_cell(i) != 0)
+        distributed_error_per_cell(i) = error_per_cell(i);
+    distributed_error_per_cell.compress (VectorOperation::insert);
+    error_per_cell = distributed_error_per_cell;
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                     error_per_cell,
+                                                     0.25, 0.03);
+    triangulation.execute_coarsening_and_refinement ();
+    setup_quadrature_point_history ();
+  }
+  template <int dim>
+  void TopLevel<dim>::move_mesh ()
+  {
+    pcout << "    Moving mesh..." << std::endl;
+    std::vector<bool> vertex_touched (triangulation.n_vertices(),
+                                      false);
+    for (typename DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active ();
+         cell != dof_handler.end(); ++cell)
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+        if (vertex_touched[cell->vertex_index(v)] == false)
+          {
+            vertex_touched[cell->vertex_index(v)] = true;
+            Point<dim> vertex_displacement;
+            for (unsigned int d=0; d<dim; ++d)
+              vertex_displacement[d]
+                = incremental_displacement(cell->vertex_dof_index(v,d));
+            cell->vertex(v) += vertex_displacement;
+          }
+  }
+  template <int dim>
+  void TopLevel<dim>::setup_quadrature_point_history ()
+  {
+    unsigned int our_cells = 0;
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell)
+      if (cell->is_locally_owned())
+        ++our_cells;
+    triangulation.clear_user_data();
+    {
+      std::vector<PointHistory<dim> > tmp;
+      tmp.swap (quadrature_point_history);
+    }
+    quadrature_point_history.resize (our_cells *
+                                     quadrature_formula.size());
+    unsigned int history_index = 0;
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell)
+      if (cell->is_locally_owned())
+        {
+          cell->set_user_pointer (&quadrature_point_history[history_index]);
+          history_index += quadrature_formula.size();
+        }
+    Assert (history_index == quadrature_point_history.size(),
+            ExcInternalError());
+  }
+  template <int dim>
+  void TopLevel<dim>::update_quadrature_point_history ()
+  {
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values | update_gradients);
+    std::vector<std::vector<Tensor<1,dim> > >
+    displacement_increment_grads (quadrature_formula.size(),
+                                  std::vector<Tensor<1,dim> >(dim));
+    for (typename DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell->is_locally_owned())
+        {
+          PointHistory<dim> *local_quadrature_points_history
+            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+          Assert (local_quadrature_points_history >=
+                  &quadrature_point_history.front(),
+                  ExcInternalError());
+          Assert (local_quadrature_points_history <
+                  &quadrature_point_history.back(),
+                  ExcInternalError());
+          fe_values.reinit (cell);
+          fe_values.get_function_gradients (incremental_displacement,
+                                            displacement_increment_grads);
+          for (unsigned int q=0; q<quadrature_formula.size(); ++q)
+            {
+              const SymmetricTensor<2,dim> new_stress
+                = (local_quadrature_points_history[q].old_stress
+                   +
+                   (stress_strain_tensor *
+                    get_strain (displacement_increment_grads[q])));
+              const Tensor<2,dim> rotation
+                = get_rotation_matrix (displacement_increment_grads[q]);
+              const SymmetricTensor<2,dim> rotated_new_stress
+                = symmetrize(transpose(rotation) *
+                             static_cast<Tensor<2,dim> >(new_stress) *
+                             rotation);
+              local_quadrature_points_history[q].old_stress
+                = rotated_new_stress;
+            }
+        }
+  }
+}
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step18;
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      TopLevel<3> elastic_problem;
+      elastic_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/physics/step-18-rotation_matrix.with_petsc=true.output b/tests/physics/step-18-rotation_matrix.with_petsc=true.output
new file mode 100644 (file)
index 0000000..91219e3
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Timestep 2: 0.00982 -1.35e-06 -0.100
+DEAL::Timestep 3: 0.0201 -3.32e-06 -0.200
+DEAL::Timestep 4: 0.0308 -5.19e-06 -0.300
+DEAL::Timestep 5: 0.0424 -5.68e-06 -0.400
+DEAL::Timestep 6: 0.0554 7.87e-07 -0.500
+DEAL::Timestep 7: 0.0720 3.84e-05 -0.600
+DEAL::Timestep 8: 0.0979 0.000188 -0.700
+DEAL::Timestep 9: 0.139 0.000618 -0.800
+DEAL::Timestep 10: 0.178 0.00126 -0.900
diff --git a/tests/physics/step-18.cc b/tests/physics/step-18.cc
new file mode 100644 (file)
index 0000000..8f24216
--- /dev/null
@@ -0,0 +1,839 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2017 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 at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// This is a copy of step-18 (git rev ecafd3f), but with a corrected rotation
+// matrix, to use as a base-line for results produced via different approaches
+// to be compared to.
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+namespace Step18
+{
+  using namespace dealii;
+  template <int dim>
+  struct PointHistory
+  {
+    SymmetricTensor<2,dim> old_stress;
+  };
+  template <int dim>
+  SymmetricTensor<4,dim>
+  get_stress_strain_tensor (const double lambda, const double mu)
+  {
+    SymmetricTensor<4,dim> tmp;
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=0; j<dim; ++j)
+        for (unsigned int k=0; k<dim; ++k)
+          for (unsigned int l=0; l<dim; ++l)
+            tmp[i][j][k][l] = (((i==k) && (j==l) ? mu : 0.0) +
+                               ((i==l) && (j==k) ? mu : 0.0) +
+                               ((i==j) && (k==l) ? lambda : 0.0));
+    return tmp;
+  }
+  template <int dim>
+  inline
+  SymmetricTensor<2,dim>
+  get_strain (const FEValues<dim> &fe_values,
+              const unsigned int   shape_func,
+              const unsigned int   q_point)
+  {
+    SymmetricTensor<2,dim> tmp;
+    for (unsigned int i=0; i<dim; ++i)
+      tmp[i][i] = fe_values.shape_grad_component (shape_func,q_point,i)[i];
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        tmp[i][j]
+          = (fe_values.shape_grad_component (shape_func,q_point,i)[j] +
+             fe_values.shape_grad_component (shape_func,q_point,j)[i]) / 2;
+    return tmp;
+  }
+  template <int dim>
+  inline
+  SymmetricTensor<2,dim>
+  get_strain (const std::vector<Tensor<1,dim> > &grad)
+  {
+    Assert (grad.size() == dim, ExcInternalError());
+    SymmetricTensor<2,dim> strain;
+    for (unsigned int i=0; i<dim; ++i)
+      strain[i][i] = grad[i][i];
+    for (unsigned int i=0; i<dim; ++i)
+      for (unsigned int j=i+1; j<dim; ++j)
+        strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
+    return strain;
+  }
+  Tensor<2,2>
+  get_rotation_matrix (const std::vector<Tensor<1,2> > &grad_u)
+  {
+    const double curl = (grad_u[1][0] - grad_u[0][1]);
+    const double angle = std::atan (curl);
+    const double t[2][2] = {{ cos(angle), sin(angle) },
+      {-sin(angle), cos(angle) }
+    };
+    return Tensor<2,2>(t);
+  }
+  Tensor<2,3>
+  get_rotation_matrix (const std::vector<Tensor<1,3> > &grad_u)
+  {
+    const Point<3> curl (grad_u[2][1] - grad_u[1][2],
+                         grad_u[0][2] - grad_u[2][0],
+                         grad_u[1][0] - grad_u[0][1]);
+    const double tan_angle = std::sqrt(curl*curl);
+    const double angle = std::atan (tan_angle);
+    if (angle < 1e-9)
+      {
+        static const double rotation[3][3]
+        = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
+        static const Tensor<2,3> rot(rotation);
+        return rot;
+      }
+    const double c = std::cos(angle);
+    const double s = std::sin(angle);
+    const double t = 1-c;
+    const Point<3> axis = curl/tan_angle;
+    // Rotation matrix fixed: See issue #468
+    const double rotation[3][3]
+    = {{
+        t *axis[0] *axis[0]+c,
+        t *axis[0] *axis[1]+s *axis[2],
+        t *axis[0] *axis[2]-s *axis[1]
+      },
+      {
+        t *axis[0] *axis[1]-s *axis[2],
+        t *axis[1] *axis[1]+c,
+        t *axis[1] *axis[2]+s *axis[0]
+      },
+      {
+        t *axis[0] *axis[2]+s *axis[1],
+        t *axis[1] *axis[2]-s *axis[0],
+        t *axis[2] *axis[2]+c
+      }
+    };
+    return Tensor<2,3>(rotation);
+  }
+  template <int dim>
+  class TopLevel
+  {
+  public:
+    TopLevel ();
+    ~TopLevel ();
+    void run ();
+  private:
+    void create_coarse_grid ();
+    void setup_system ();
+    void assemble_system ();
+    void solve_timestep ();
+    unsigned int solve_linear_problem ();
+    void output_results () const;
+    void do_initial_timestep ();
+    void do_timestep ();
+    void refine_initial_grid ();
+    void move_mesh ();
+    void setup_quadrature_point_history ();
+    void update_quadrature_point_history ();
+    parallel::shared::Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    ConstraintMatrix     hanging_node_constraints;
+    const QGauss<dim>          quadrature_formula;
+    std::vector<PointHistory<dim> > quadrature_point_history;
+    PETScWrappers::MPI::SparseMatrix system_matrix;
+    PETScWrappers::MPI::Vector       system_rhs;
+    Vector<double>                   incremental_displacement;
+    double       present_time;
+    double       present_timestep;
+    double       end_time;
+    unsigned int timestep_no;
+    MPI_Comm mpi_communicator;
+    const unsigned int n_mpi_processes;
+    const unsigned int this_mpi_process;
+    ConditionalOStream pcout;
+    std::vector<types::global_dof_index> local_dofs_per_process;
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+    unsigned int         n_local_cells;
+    static const SymmetricTensor<4,dim> stress_strain_tensor;
+    int monitored_vertex_first_dof;
+  };
+  template <int dim>
+  class BodyForce :  public Function<dim>
+  {
+  public:
+    BodyForce ();
+    virtual
+    void
+    vector_value (const Point<dim> &p,
+                  Vector<double>   &values) const;
+    virtual
+    void
+    vector_value_list (const std::vector<Point<dim> > &points,
+                       std::vector<Vector<double> >   &value_list) const;
+  };
+  template <int dim>
+  BodyForce<dim>::BodyForce ()
+    :
+    Function<dim> (dim)
+  {}
+  template <int dim>
+  inline
+  void
+  BodyForce<dim>::vector_value (const Point<dim> &/*p*/,
+                                Vector<double>   &values) const
+  {
+    Assert (values.size() == dim,
+            ExcDimensionMismatch (values.size(), dim));
+    const double g   = 9.81;
+    const double rho = 7700;
+    values = 0;
+    values(dim-1) = -rho * g;
+  }
+  template <int dim>
+  void
+  BodyForce<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<Vector<double> >   &value_list) const
+  {
+    const unsigned int n_points = points.size();
+    Assert (value_list.size() == n_points,
+            ExcDimensionMismatch (value_list.size(), n_points));
+    for (unsigned int p=0; p<n_points; ++p)
+      BodyForce<dim>::vector_value (points[p],
+                                    value_list[p]);
+  }
+  template <int dim>
+  class IncrementalBoundaryValues :  public Function<dim>
+  {
+  public:
+    IncrementalBoundaryValues (const double present_time,
+                               const double present_timestep);
+    virtual
+    void
+    vector_value (const Point<dim> &p,
+                  Vector<double>   &values) const;
+    virtual
+    void
+    vector_value_list (const std::vector<Point<dim> > &points,
+                       std::vector<Vector<double> >   &value_list) const;
+  private:
+    const double velocity;
+    const double present_time;
+    const double present_timestep;
+  };
+  template <int dim>
+  IncrementalBoundaryValues<dim>::
+  IncrementalBoundaryValues (const double present_time,
+                             const double present_timestep)
+    :
+    Function<dim> (dim),
+    velocity (.1),
+    present_time (present_time),
+    present_timestep (present_timestep)
+  {}
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::
+  vector_value (const Point<dim> &/*p*/,
+                Vector<double>   &values) const
+  {
+    Assert (values.size() == dim,
+            ExcDimensionMismatch (values.size(), dim));
+    values = 0;
+    values(2) = -present_timestep * velocity;
+  }
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::
+  vector_value_list (const std::vector<Point<dim> > &points,
+                     std::vector<Vector<double> >   &value_list) const
+  {
+    const unsigned int n_points = points.size();
+    Assert (value_list.size() == n_points,
+            ExcDimensionMismatch (value_list.size(), n_points));
+    for (unsigned int p=0; p<n_points; ++p)
+      IncrementalBoundaryValues<dim>::vector_value (points[p],
+                                                    value_list[p]);
+  }
+  template <int dim>
+  const SymmetricTensor<4,dim>
+  TopLevel<dim>::stress_strain_tensor
+    = get_stress_strain_tensor<dim> (/*lambda = */ 9.695e10,
+                                                   /*mu     = */ 7.617e10);
+  template <int dim>
+  TopLevel<dim>::TopLevel ()
+    :
+    triangulation(MPI_COMM_WORLD),
+    fe (FE_Q<dim>(1), dim),
+    dof_handler (triangulation),
+    quadrature_formula (2),
+    mpi_communicator (MPI_COMM_WORLD),
+    n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
+    this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
+    pcout (std::cout, this_mpi_process == 0),
+    monitored_vertex_first_dof(0)
+  {}
+  template <int dim>
+  TopLevel<dim>::~TopLevel ()
+  {
+    dof_handler.clear ();
+  }
+  template <int dim>
+  void TopLevel<dim>::run ()
+  {
+    present_time = 0;
+    present_timestep = 1;
+    end_time = 10;
+    timestep_no = 0;
+    do_initial_timestep ();
+    while (present_time < end_time)
+      do_timestep ();
+  }
+  template <int dim>
+  void TopLevel<dim>::create_coarse_grid ()
+  {
+    const double inner_radius = 0.8,
+                 outer_radius = 1;
+    GridGenerator::cylinder_shell (triangulation,
+                                   3, inner_radius, outer_radius);
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell=triangulation.begin_active();
+         cell!=triangulation.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+        if (cell->face(f)->at_boundary())
+          {
+            const Point<dim> face_center = cell->face(f)->center();
+            if (face_center[2] == 0)
+              cell->face(f)->set_boundary_id (0);
+            else if (face_center[2] == 3)
+              cell->face(f)->set_boundary_id (1);
+            else if (std::sqrt(face_center[0]*face_center[0] +
+                               face_center[1]*face_center[1])
+                     <
+                     (inner_radius + outer_radius) / 2)
+              cell->face(f)->set_boundary_id (2);
+            else
+              cell->face(f)->set_boundary_id (3);
+          }
+    static const CylindricalManifold<dim> cylindrical_manifold (2);
+    triangulation.set_all_manifold_ids(0);
+    triangulation.set_manifold (0, cylindrical_manifold);
+    // triangulation.refine_global (1);
+    setup_quadrature_point_history ();
+  }
+  template <int dim>
+  void TopLevel<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs);
+    n_local_cells
+      = GridTools::count_cells_with_subdomain_association (triangulation,
+                                                           triangulation.locally_owned_subdomain ());
+    local_dofs_per_process = dof_handler.n_locally_owned_dofs_per_processor();
+    hanging_node_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                             hanging_node_constraints);
+    hanging_node_constraints.close ();
+    DynamicSparsityPattern sparsity_pattern (locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
+                                     hanging_node_constraints, /*keep constrained dofs*/ false);
+    SparsityTools::distribute_sparsity_pattern (sparsity_pattern,
+                                                local_dofs_per_process,
+                                                mpi_communicator,
+                                                locally_relevant_dofs);
+    system_matrix.reinit (locally_owned_dofs,
+                          locally_owned_dofs,
+                          sparsity_pattern,
+                          mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs,mpi_communicator);
+    incremental_displacement.reinit (dof_handler.n_dofs());
+  }
+  template <int dim>
+  void TopLevel<dim>::assemble_system ()
+  {
+    system_rhs = 0;
+    system_matrix = 0;
+    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);
+    BodyForce<dim>      body_force;
+    std::vector<Vector<double> > body_force_values (n_q_points,
+                                                    Vector<double>(dim));
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0;
+          cell_rhs = 0;
+          fe_values.reinit (cell);
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              for (unsigned int q_point=0; q_point<n_q_points;
+                   ++q_point)
+                {
+                  const SymmetricTensor<2,dim>
+                  eps_phi_i = get_strain (fe_values, i, q_point),
+                  eps_phi_j = get_strain (fe_values, j, q_point);
+                  cell_matrix(i,j)
+                  += (eps_phi_i * stress_strain_tensor * eps_phi_j
+                      *
+                      fe_values.JxW (q_point));
+                }
+          const PointHistory<dim> *local_quadrature_points_data
+            = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
+          body_force.vector_value_list (fe_values.get_quadrature_points(),
+                                        body_force_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 q_point=0; q_point<n_q_points; ++q_point)
+                {
+                  const SymmetricTensor<2,dim> &old_stress
+                    = local_quadrature_points_data[q_point].old_stress;
+                  cell_rhs(i) += (body_force_values[q_point](component_i) *
+                                  fe_values.shape_value (i,q_point)
+                                  -
+                                  old_stress *
+                                  get_strain (fe_values,i,q_point))
+                                 *
+                                 fe_values.JxW (q_point);
+                }
+            }
+          cell->get_dof_indices (local_dof_indices);
+          hanging_node_constraints
+          .distribute_local_to_global (cell_matrix, cell_rhs,
+                                       local_dof_indices,
+                                       system_matrix, system_rhs);
+        }
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+    FEValuesExtractors::Scalar z_component (dim-1);
+    std::map<types::global_dof_index,double> boundary_values;
+    VectorTools::
+    interpolate_boundary_values (dof_handler,
+                                 0,
+                                 ZeroFunction<dim> (dim),
+                                 boundary_values);
+    VectorTools::
+    interpolate_boundary_values (dof_handler,
+                                 1,
+                                 IncrementalBoundaryValues<dim>(present_time,
+                                                                present_timestep),
+                                 boundary_values,
+                                 fe.component_mask(z_component));
+    PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                        system_matrix, tmp,
+                                        system_rhs, false);
+    incremental_displacement = tmp;
+  }
+  template <int dim>
+  void TopLevel<dim>::solve_timestep ()
+  {
+    pcout << "    Assembling system..." << std::flush;
+    assemble_system ();
+    pcout << " norm of rhs is " << system_rhs.l2_norm()
+          << std::endl;
+    const unsigned int n_iterations = solve_linear_problem ();
+    pcout << "    Solver converged in " << n_iterations
+          << " iterations." << std::endl;
+    pcout << "    Updating quadrature point data..." << std::flush;
+    update_quadrature_point_history ();
+    pcout << std::endl;
+  }
+  template <int dim>
+  unsigned int TopLevel<dim>::solve_linear_problem ()
+  {
+    PETScWrappers::MPI::Vector
+    distributed_incremental_displacement (locally_owned_dofs,mpi_communicator);
+    distributed_incremental_displacement = incremental_displacement;
+    SolverControl           solver_control (dof_handler.n_dofs(),
+                                            1e-16*system_rhs.l2_norm(),
+                                            false,false);
+    PETScWrappers::SolverCG cg (solver_control,
+                                mpi_communicator);
+    PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
+    cg.solve (system_matrix, distributed_incremental_displacement, system_rhs,
+              preconditioner);
+    incremental_displacement = distributed_incremental_displacement;
+    hanging_node_constraints.distribute (incremental_displacement);
+    return solver_control.last_step();
+  }
+  template <int dim>
+  void TopLevel<dim>::output_results () const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (dof_handler);
+    std::vector<std::string> solution_names;
+    switch (dim)
+      {
+      case 1:
+        solution_names.push_back ("delta_x");
+        break;
+      case 2:
+        solution_names.push_back ("delta_x");
+        solution_names.push_back ("delta_y");
+        break;
+      case 3:
+        solution_names.push_back ("delta_x");
+        solution_names.push_back ("delta_y");
+        solution_names.push_back ("delta_z");
+        break;
+      default:
+        Assert (false, ExcNotImplemented());
+      }
+    data_out.add_data_vector (incremental_displacement,
+                              solution_names);
+    Vector<double> norm_of_stress (triangulation.n_active_cells());
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = triangulation.begin_active(),
+      endc = triangulation.end();
+      for (; cell!=endc; ++cell)
+        if (cell->is_locally_owned())
+          {
+            SymmetricTensor<2,dim> accumulated_stress;
+            for (unsigned int q=0;
+                 q<quadrature_formula.size();
+                 ++q)
+              accumulated_stress +=
+                reinterpret_cast<PointHistory<dim>*>(cell->user_pointer())[q]
+                .old_stress;
+            norm_of_stress(cell->active_cell_index())
+              = (accumulated_stress /
+                 quadrature_formula.size()).norm();
+          }
+        else
+          norm_of_stress(cell->active_cell_index()) = -1e+20;
+    }
+    data_out.add_data_vector (norm_of_stress, "norm_of_stress");
+    std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
+    GridTools::get_subdomain_association (triangulation, partition_int);
+    const Vector<double> partitioning(partition_int.begin(),
+                                      partition_int.end());
+    data_out.add_data_vector (partitioning, "partitioning");
+    data_out.build_patches ();
+    std::string filename = "solution-" + Utilities::int_to_string(timestep_no,4)
+                           + "." + Utilities::int_to_string(this_mpi_process,3)
+                           + ".vtu";
+    AssertThrow (n_mpi_processes < 1000, ExcNotImplemented());
+    std::ofstream output (filename.c_str());
+    data_out.write_vtu (output);
+    if (this_mpi_process==0)
+      {
+        std::vector<std::string> filenames;
+        for (unsigned int i=0; i<n_mpi_processes; ++i)
+          filenames.push_back ("solution-" + Utilities::int_to_string(timestep_no,4)
+                               + "." + Utilities::int_to_string(i,3)
+                               + ".vtu");
+        const std::string
+        visit_master_filename = ("solution-" +
+                                 Utilities::int_to_string(timestep_no,4) +
+                                 ".visit");
+        std::ofstream visit_master (visit_master_filename.c_str());
+        DataOutBase::write_visit_record (visit_master, filenames);
+        const std::string
+        pvtu_master_filename = ("solution-" +
+                                Utilities::int_to_string(timestep_no,4) +
+                                ".pvtu");
+        std::ofstream pvtu_master (pvtu_master_filename.c_str());
+        data_out.write_pvtu_record (pvtu_master, filenames);
+        static std::vector<std::pair<double,std::string> > times_and_names;
+        times_and_names.push_back (std::pair<double,std::string> (present_time, pvtu_master_filename));
+        std::ofstream pvd_output ("solution.pvd");
+        DataOutBase::write_pvd_record (pvd_output, times_and_names);
+      }
+  }
+  template <int dim>
+  void TopLevel<dim>::do_initial_timestep ()
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    pcout << "Timestep " << timestep_no << " at time " << present_time
+          << std::endl;
+    for (unsigned int cycle=0; cycle<2; ++cycle)
+      {
+        pcout << "  Cycle " << cycle << ':' << std::endl;
+        if (cycle == 0)
+          create_coarse_grid ();
+        else
+          refine_initial_grid ();
+        pcout << "    Number of active cells:       "
+              << triangulation.n_active_cells()
+              << " (by partition:";
+        for (unsigned int p=0; p<n_mpi_processes; ++p)
+          pcout << (p==0 ? ' ' : '+')
+                << (GridTools::
+                    count_cells_with_subdomain_association (triangulation,p));
+        pcout << ")" << std::endl;
+        setup_system ();
+        pcout << "    Number of degrees of freedom: "
+              << dof_handler.n_dofs()
+              << " (by partition:";
+        for (unsigned int p=0; p<n_mpi_processes; ++p)
+          pcout << (p==0 ? ' ' : '+')
+                << (DoFTools::
+                    count_dofs_with_subdomain_association (dof_handler,p));
+        pcout << ")" << std::endl;
+        // Get point at which to output displacement
+        // (outer radius of displaced surface)
+        {
+          const Point<dim> soln_pt (1.0, 0.0, 3.0);
+          typename DoFHandler<dim>::active_cell_iterator cell =
+            dof_handler.begin_active(), endc = dof_handler.end();
+          for (; cell != endc; ++cell)
+            for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+              if (cell->vertex(v).distance(soln_pt) < 1e-6)
+                {
+                  monitored_vertex_first_dof = cell->vertex_dof_index(v,0);
+                }
+        }
+        solve_timestep ();
+      }
+    move_mesh ();
+    // output_results ();
+    pcout << std::endl;
+  }
+  template <int dim>
+  void TopLevel<dim>::do_timestep ()
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    pcout << "Timestep " << timestep_no << " at time " << present_time
+          << std::endl;
+    if (present_time > end_time)
+      {
+        present_timestep -= (present_time - end_time);
+        present_time = end_time;
+      }
+    solve_timestep ();
+    move_mesh ();
+    // output_results ();
+    // Output displacement at edge of displaced surface
+    {
+      static Tensor<1,dim> soln;
+      for (unsigned int d=0; d<dim; ++d)
+        soln[d] += incremental_displacement(monitored_vertex_first_dof+d);
+      deallog
+          << "Timestep " << timestep_no
+          << ": " << soln << std::endl;
+    }
+    pcout << std::endl;
+  }
+  template <int dim>
+  void TopLevel<dim>::refine_initial_grid ()
+  {
+    Vector<float> error_per_cell (triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                        QGauss<dim-1>(2),
+                                        typename FunctionMap<dim>::type(),
+                                        incremental_displacement,
+                                        error_per_cell,
+                                        ComponentMask(),
+                                        0,
+                                        MultithreadInfo::n_threads(),
+                                        this_mpi_process);
+    const unsigned int n_local_cells = triangulation.n_locally_owned_active_cells ();
+    PETScWrappers::MPI::Vector
+    distributed_error_per_cell (mpi_communicator,
+                                triangulation.n_active_cells(),
+                                n_local_cells);
+    for (unsigned int i=0; i<error_per_cell.size(); ++i)
+      if (error_per_cell(i) != 0)
+        distributed_error_per_cell(i) = error_per_cell(i);
+    distributed_error_per_cell.compress (VectorOperation::insert);
+    error_per_cell = distributed_error_per_cell;
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                     error_per_cell,
+                                                     0.25, 0.03);
+    triangulation.execute_coarsening_and_refinement ();
+    setup_quadrature_point_history ();
+  }
+  template <int dim>
+  void TopLevel<dim>::move_mesh ()
+  {
+    pcout << "    Moving mesh..." << std::endl;
+    std::vector<bool> vertex_touched (triangulation.n_vertices(),
+                                      false);
+    for (typename DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active ();
+         cell != dof_handler.end(); ++cell)
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+        if (vertex_touched[cell->vertex_index(v)] == false)
+          {
+            vertex_touched[cell->vertex_index(v)] = true;
+            Point<dim> vertex_displacement;
+            for (unsigned int d=0; d<dim; ++d)
+              vertex_displacement[d]
+                = incremental_displacement(cell->vertex_dof_index(v,d));
+            cell->vertex(v) += vertex_displacement;
+          }
+  }
+  template <int dim>
+  void TopLevel<dim>::setup_quadrature_point_history ()
+  {
+    unsigned int our_cells = 0;
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell)
+      if (cell->is_locally_owned())
+        ++our_cells;
+    triangulation.clear_user_data();
+    {
+      std::vector<PointHistory<dim> > tmp;
+      tmp.swap (quadrature_point_history);
+    }
+    quadrature_point_history.resize (our_cells *
+                                     quadrature_formula.size());
+    unsigned int history_index = 0;
+    for (typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+         cell != triangulation.end(); ++cell)
+      if (cell->is_locally_owned())
+        {
+          cell->set_user_pointer (&quadrature_point_history[history_index]);
+          history_index += quadrature_formula.size();
+        }
+    Assert (history_index == quadrature_point_history.size(),
+            ExcInternalError());
+  }
+  template <int dim>
+  void TopLevel<dim>::update_quadrature_point_history ()
+  {
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values | update_gradients);
+    std::vector<std::vector<Tensor<1,dim> > >
+    displacement_increment_grads (quadrature_formula.size(),
+                                  std::vector<Tensor<1,dim> >(dim));
+    for (typename DoFHandler<dim>::active_cell_iterator
+         cell = dof_handler.begin_active();
+         cell != dof_handler.end(); ++cell)
+      if (cell->is_locally_owned())
+        {
+          PointHistory<dim> *local_quadrature_points_history
+            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+          Assert (local_quadrature_points_history >=
+                  &quadrature_point_history.front(),
+                  ExcInternalError());
+          Assert (local_quadrature_points_history <
+                  &quadrature_point_history.back(),
+                  ExcInternalError());
+          fe_values.reinit (cell);
+          fe_values.get_function_gradients (incremental_displacement,
+                                            displacement_increment_grads);
+          for (unsigned int q=0; q<quadrature_formula.size(); ++q)
+            {
+              const SymmetricTensor<2,dim> new_stress
+                = (local_quadrature_points_history[q].old_stress
+                   +
+                   (stress_strain_tensor *
+                    get_strain (displacement_increment_grads[q])));
+              const Tensor<2,dim> rotation
+                = get_rotation_matrix (displacement_increment_grads[q]);
+              const SymmetricTensor<2,dim> rotated_new_stress
+                = symmetrize(transpose(rotation) *
+                             static_cast<Tensor<2,dim> >(new_stress) *
+                             rotation);
+              local_quadrature_points_history[q].old_stress
+                = rotated_new_stress;
+            }
+        }
+  }
+}
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step18;
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      TopLevel<3> elastic_problem;
+      elastic_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/physics/step-18.with_petsc=true.output b/tests/physics/step-18.with_petsc=true.output
new file mode 100644 (file)
index 0000000..91219e3
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Timestep 2: 0.00982 -1.35e-06 -0.100
+DEAL::Timestep 3: 0.0201 -3.32e-06 -0.200
+DEAL::Timestep 4: 0.0308 -5.19e-06 -0.300
+DEAL::Timestep 5: 0.0424 -5.68e-06 -0.400
+DEAL::Timestep 6: 0.0554 7.87e-07 -0.500
+DEAL::Timestep 7: 0.0720 3.84e-05 -0.600
+DEAL::Timestep 8: 0.0979 0.000188 -0.700
+DEAL::Timestep 9: 0.139 0.000618 -0.800
+DEAL::Timestep 10: 0.178 0.00126 -0.900
diff --git a/tests/physics/transformations-rotations_01.cc b/tests/physics/transformations-rotations_01.cc
new file mode 100644 (file)
index 0000000..3114853
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test rotation matrix definitions
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/physics/transformations.h>
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+using namespace dealii::Physics;
+
+
+void
+test_rotation_matrix_3d_z_axis (const double angle)
+{
+  Tensor<2,3> R_z = unit_symmetric_tensor<3>();
+  const Tensor<2,2> R_2d = Transformations::Rotations::rotation_matrix_2d(angle);
+  for (unsigned int i=0; i<2; ++i)
+    for (unsigned int j=0; j<2; ++j)
+      R_z[i][j] = R_2d[i][j];
+
+  Assert(std::abs(determinant(R_z) - 1.0) < 1e-9,
+         ExcMessage("Rodrigues rotation matrix determinant is not unity"));
+  const Tensor<2,3> R = Transformations::Rotations::rotation_matrix_3d(Point<3>({0,0,1}),angle);
+  Assert(std::abs(determinant(R) - 1.0) < 1e-9,
+         ExcMessage("Rotation matrix determinant is not unity"));
+
+  Assert((transpose(R)*R - unit_symmetric_tensor<3>()).norm() < 1e-9,
+         ExcMessage("Matrix is not a rotation matrix"));
+  Assert((R - R_z).norm() < 1e-12,
+         ExcMessage("Incorrect computation of R in 3d"));
+}
+
+void
+test_rotation_matrix_3d (const Point<3> &axis,
+                         const double    angle)
+{
+  // http://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
+  // http://en.wikipedia.org/wiki/Rotation_matrix
+  // NOTE: Angle in radians
+  const Tensor<1,3> u = axis / axis.norm(); // Ensure unit vector
+  const Tensor<2,3> u_dyad_u = outer_product (u, u);
+  const double u_skew_array [3][3] =
+  {
+    {0.0, -u[2], u[1]},
+    {u[2], 0.0, -u[0]},
+    {-u[1], u[0], 0.0}
+  };
+
+  const Tensor<2,3> R_rodrigues
+    = u_dyad_u
+      + std::cos(angle)*(static_cast < Tensor<2,3> > (unit_symmetric_tensor<3>()) - u_dyad_u)
+      + std::sin(angle)*Tensor<2,3> (u_skew_array);
+
+
+  Assert(std::abs(determinant(R_rodrigues) - 1.0) < 1e-9,
+         ExcMessage("Rodrigues rotation matrix determinant is not unity"));
+  const Tensor<2,3> R = Transformations::Rotations::rotation_matrix_3d(axis,angle);
+  Assert(std::abs(determinant(R) - 1.0) < 1e-9,
+         ExcMessage("Rotation matrix determinant is not unity"));
+
+  Assert((transpose(R)*R - unit_symmetric_tensor<3>()).norm() < 1e-9,
+         ExcMessage("Matrix is not a rotation matrix"));
+  Assert((R - R_rodrigues).norm() < 1e-12,
+         ExcMessage("Incorrect computation of R in 3d"));
+}
+
+Point<3>
+normalise(const Point<3> &p)
+{
+  Assert(p.norm() > 0.0, ExcMessage("Point vector has zero norm"));
+  return p/p.norm();
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  const double deg_to_rad = M_PI/180.0;
+
+  // 2-d
+  {
+    const double angle = 90.0*deg_to_rad;
+    const Tensor<1,2> in ({1,0});
+    const Tensor<2,2> R = Transformations::Rotations::rotation_matrix_2d(angle);
+    Assert((transpose(R)*R - unit_symmetric_tensor<2>()).norm() < 1e-9, ExcMessage("Matrix is not a rotation matrix"));
+    Assert(std::abs(determinant(R) - 1.0) < 1e-9, ExcMessage("Rotation matrix determinant is not unity"));
+    const Tensor<1,2> out = R*in;
+    Assert((out - Tensor<1,2>({0,1})).norm() < 1e-12,
+           ExcMessage("Incorrect computation of 90 degree R in 2d"));
+  }
+  {
+    const double angle = 135.0*deg_to_rad;
+    const Tensor<1,2> in ({1,0});
+    const Tensor<2,2> R = Transformations::Rotations::rotation_matrix_2d(angle);
+    Assert((transpose(R)*R - unit_symmetric_tensor<2>()).norm() < 1e-9, ExcMessage("Matrix is not a rotation matrix"));
+    Assert(std::abs(determinant(R) - 1.0) < 1e-9, ExcMessage("Rotation matrix determinant is not unity"));
+    const Tensor<1,2> out = R*in;
+    Assert((out - Tensor<1,2>({-1.0/std::sqrt(2.0),1.0/std::sqrt(2.0)})).norm() < 1e-12,
+           ExcMessage("Incorrect computation of 135 degree R in 2d"));
+  }
+  {
+    const double angle = 240.0*deg_to_rad;
+    const Tensor<1,2> in ({1,0});
+    const Tensor<2,2> R = Transformations::Rotations::rotation_matrix_2d(angle);
+    Assert((transpose(R)*R - unit_symmetric_tensor<2>()).norm() < 1e-9, ExcMessage("Matrix is not a rotation matrix"));
+    Assert(std::abs(determinant(R) - 1.0) < 1e-9, ExcMessage("Rotation matrix determinant is not unity"));
+    const Tensor<1,2> out = R*in;
+    Assert((out - Tensor<1,2>({-0.5,-std::sqrt(3.0)/2.0})).norm() < 1e-12,
+           ExcMessage("Incorrect computation of 240 degree R in 2d"));
+  }
+
+  // 3-d
+  test_rotation_matrix_3d_z_axis (90.0*deg_to_rad);
+  test_rotation_matrix_3d_z_axis (45.0*deg_to_rad);
+  test_rotation_matrix_3d_z_axis (60.0*deg_to_rad);
+
+  test_rotation_matrix_3d (normalise(Point<3>({1,1,1})), 90.0*deg_to_rad);
+  test_rotation_matrix_3d (normalise(Point<3>({0,2,1})), 45.0*deg_to_rad);
+  test_rotation_matrix_3d (normalise(Point<3>({-1,3,2})), 60.0*deg_to_rad);
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/physics/transformations-rotations_01.output b/tests/physics/transformations-rotations_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.