From: Jean-Paul Pelteret Date: Mon, 9 Jan 2017 17:06:18 +0000 (+0100) Subject: Added functions to compute rotation matrices in 2d and 3d. X-Git-Tag: v8.5.0-rc1~268^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5b61da2ba47cabdaeb70e14bbf16321a3ca9f79e;p=dealii.git Added functions to compute rotation matrices in 2d and 3d. --- diff --git a/doc/news/changes/minor/20170109Jean-PaulPelteret b/doc/news/changes/minor/20170109Jean-PaulPelteret new file mode 100644 index 0000000000..9bddb7b187 --- /dev/null +++ b/doc/news/changes/minor/20170109Jean-PaulPelteret @@ -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. +
+(Jean-Paul Pelteret, 2017/01/09) diff --git a/include/deal.II/physics/transformations.h b/include/deal.II/physics/transformations.h index b77d6d3cc6..e7e6ebff06 100644 --- a/include/deal.II/physics/transformations.h +++ b/include/deal.II/physics/transformations.h @@ -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 #include #include @@ -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 + * fixed set of right-handed axes. + * + * @param[in] angle The rotation angle (about the z-axis) in radians + */ + template + 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 + * wikipedia link. + * In particular, this describes the counter-clockwise rotation of a vector + * in a plane with its normal. + * defined by the @p axis of rotation. + * An alternative implementation is discussed at + * this link, + * 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 + 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 +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 +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 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 index 0000000000..27a58b0e90 --- /dev/null +++ b/tests/physics/step-18-rotation_matrix.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace Step18 +{ + using namespace dealii; + template + struct PointHistory + { + SymmetricTensor<2,dim> old_stress; + }; + template + SymmetricTensor<4,dim> + get_stress_strain_tensor (const double lambda, const double mu) + { + SymmetricTensor<4,dim> tmp; + for (unsigned int i=0; i + inline + SymmetricTensor<2,dim> + get_strain (const FEValues &fe_values, + const unsigned int shape_func, + const unsigned int q_point) + { + SymmetricTensor<2,dim> tmp; + for (unsigned int i=0; i + inline + SymmetricTensor<2,dim> + get_strain (const std::vector > &grad) + { + Assert (grad.size() == dim, ExcInternalError()); + SymmetricTensor<2,dim> strain; + for (unsigned int i=0; i + get_rotation_matrix (const std::vector > &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 > &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 + 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 triangulation; + FESystem fe; + DoFHandler dof_handler; + ConstraintMatrix hanging_node_constraints; + const QGauss quadrature_formula; + std::vector > quadrature_point_history; + PETScWrappers::MPI::SparseMatrix system_matrix; + PETScWrappers::MPI::Vector system_rhs; + Vector 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 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 + class BodyForce : public Function + { + public: + BodyForce (); + virtual + void + vector_value (const Point &p, + Vector &values) const; + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + }; + template + BodyForce::BodyForce () + : + Function (dim) + {} + template + inline + void + BodyForce::vector_value (const Point &/*p*/, + Vector &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 + void + BodyForce::vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], + value_list[p]); + } + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double present_timestep); + virtual + void + vector_value (const Point &p, + Vector &values) const; + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + const double velocity; + const double present_time; + const double present_timestep; + }; + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double present_timestep) + : + Function (dim), + velocity (.1), + present_time (present_time), + present_timestep (present_timestep) + {} + template + void + IncrementalBoundaryValues:: + vector_value (const Point &/*p*/, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + values = 0; + values(2) = -present_timestep * velocity; + } + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], + value_list[p]); + } + template + const SymmetricTensor<4,dim> + TopLevel::stress_strain_tensor + = get_stress_strain_tensor (/*lambda = */ 9.695e10, + /*mu = */ 7.617e10); + template + TopLevel::TopLevel () + : + triangulation(MPI_COMM_WORLD), + fe (FE_Q(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 + TopLevel::~TopLevel () + { + dof_handler.clear (); + } + template + void TopLevel::run () + { + present_time = 0; + present_timestep = 1; + end_time = 10; + timestep_no = 0; + do_initial_timestep (); + while (present_time < end_time) + do_timestep (); + } + template + void TopLevel::create_coarse_grid () + { + const double inner_radius = 0.8, + outer_radius = 1; + GridGenerator::cylinder_shell (triangulation, + 3, inner_radius, outer_radius); + for (typename Triangulation::active_cell_iterator + cell=triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + const Point 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 cylindrical_manifold (2); + triangulation.set_all_manifold_ids(0); + triangulation.set_manifold (0, cylindrical_manifold); + // triangulation.refine_global (1); + setup_quadrature_point_history (); + } + template + void TopLevel::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 + void TopLevel::assemble_system () + { + system_rhs = 0; + system_matrix = 0; + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + BodyForce body_force; + std::vector > body_force_values (n_q_points, + Vector(dim)); + typename DoFHandler::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 + 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 *local_quadrature_points_data + = reinterpret_cast*>(cell->user_pointer()); + body_force.vector_value_list (fe_values.get_quadrature_points(), + body_force_values); + for (unsigned int i=0; i &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 boundary_values; + VectorTools:: + interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (dim), + boundary_values); + VectorTools:: + interpolate_boundary_values (dof_handler, + 1, + IncrementalBoundaryValues(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 + void TopLevel::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 + unsigned int TopLevel::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 + void TopLevel::output_results () const + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + std::vector 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 norm_of_stress (triangulation.n_active_cells()); + { + typename Triangulation::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*>(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 partition_int (triangulation.n_active_cells()); + GridTools::get_subdomain_association (triangulation, partition_int); + const Vector 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 filenames; + for (unsigned int i=0; i > times_and_names; + times_and_names.push_back (std::pair (present_time, pvtu_master_filename)); + std::ofstream pvd_output ("solution.pvd"); + DataOutBase::write_pvd_record (pvd_output, times_and_names); + } + } + template + void TopLevel::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 soln_pt (1.0, 0.0, 3.0); + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + for (; cell != endc; ++cell) + for (unsigned int v=0; v::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 + void TopLevel::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 + void TopLevel::refine_initial_grid () + { + Vector error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss(2), + typename FunctionMap::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 + void TopLevel::move_mesh () + { + pcout << " Moving mesh..." << std::endl; + std::vector vertex_touched (triangulation.n_vertices(), + false); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (); + cell != dof_handler.end(); ++cell) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + Point vertex_displacement; + for (unsigned int d=0; dvertex_dof_index(v,d)); + cell->vertex(v) += vertex_displacement; + } + } + template + void TopLevel::setup_quadrature_point_history () + { + unsigned int our_cells = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned()) + ++our_cells; + triangulation.clear_user_data(); + { + std::vector > tmp; + tmp.swap (quadrature_point_history); + } + quadrature_point_history.resize (our_cells * + quadrature_formula.size()); + unsigned int history_index = 0; + for (typename Triangulation::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 + void TopLevel::update_quadrature_point_history () + { + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients); + std::vector > > + displacement_increment_grads (quadrature_formula.size(), + std::vector >(dim)); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(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 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 >(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 index 0000000000..91219e30f8 --- /dev/null +++ b/tests/physics/step-18-rotation_matrix.with_petsc=true.output @@ -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 index 0000000000..8f24216ae4 --- /dev/null +++ b/tests/physics/step-18.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace Step18 +{ + using namespace dealii; + template + struct PointHistory + { + SymmetricTensor<2,dim> old_stress; + }; + template + SymmetricTensor<4,dim> + get_stress_strain_tensor (const double lambda, const double mu) + { + SymmetricTensor<4,dim> tmp; + for (unsigned int i=0; i + inline + SymmetricTensor<2,dim> + get_strain (const FEValues &fe_values, + const unsigned int shape_func, + const unsigned int q_point) + { + SymmetricTensor<2,dim> tmp; + for (unsigned int i=0; i + inline + SymmetricTensor<2,dim> + get_strain (const std::vector > &grad) + { + Assert (grad.size() == dim, ExcInternalError()); + SymmetricTensor<2,dim> strain; + for (unsigned int i=0; i + get_rotation_matrix (const std::vector > &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 > &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 + 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 triangulation; + FESystem fe; + DoFHandler dof_handler; + ConstraintMatrix hanging_node_constraints; + const QGauss quadrature_formula; + std::vector > quadrature_point_history; + PETScWrappers::MPI::SparseMatrix system_matrix; + PETScWrappers::MPI::Vector system_rhs; + Vector 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 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 + class BodyForce : public Function + { + public: + BodyForce (); + virtual + void + vector_value (const Point &p, + Vector &values) const; + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + }; + template + BodyForce::BodyForce () + : + Function (dim) + {} + template + inline + void + BodyForce::vector_value (const Point &/*p*/, + Vector &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 + void + BodyForce::vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], + value_list[p]); + } + template + class IncrementalBoundaryValues : public Function + { + public: + IncrementalBoundaryValues (const double present_time, + const double present_timestep); + virtual + void + vector_value (const Point &p, + Vector &values) const; + virtual + void + vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + private: + const double velocity; + const double present_time; + const double present_timestep; + }; + template + IncrementalBoundaryValues:: + IncrementalBoundaryValues (const double present_time, + const double present_timestep) + : + Function (dim), + velocity (.1), + present_time (present_time), + present_timestep (present_timestep) + {} + template + void + IncrementalBoundaryValues:: + vector_value (const Point &/*p*/, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + values = 0; + values(2) = -present_timestep * velocity; + } + template + void + IncrementalBoundaryValues:: + vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], + value_list[p]); + } + template + const SymmetricTensor<4,dim> + TopLevel::stress_strain_tensor + = get_stress_strain_tensor (/*lambda = */ 9.695e10, + /*mu = */ 7.617e10); + template + TopLevel::TopLevel () + : + triangulation(MPI_COMM_WORLD), + fe (FE_Q(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 + TopLevel::~TopLevel () + { + dof_handler.clear (); + } + template + void TopLevel::run () + { + present_time = 0; + present_timestep = 1; + end_time = 10; + timestep_no = 0; + do_initial_timestep (); + while (present_time < end_time) + do_timestep (); + } + template + void TopLevel::create_coarse_grid () + { + const double inner_radius = 0.8, + outer_radius = 1; + GridGenerator::cylinder_shell (triangulation, + 3, inner_radius, outer_radius); + for (typename Triangulation::active_cell_iterator + cell=triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + const Point 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 cylindrical_manifold (2); + triangulation.set_all_manifold_ids(0); + triangulation.set_manifold (0, cylindrical_manifold); + // triangulation.refine_global (1); + setup_quadrature_point_history (); + } + template + void TopLevel::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 + void TopLevel::assemble_system () + { + system_rhs = 0; + system_matrix = 0; + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + BodyForce body_force; + std::vector > body_force_values (n_q_points, + Vector(dim)); + typename DoFHandler::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 + 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 *local_quadrature_points_data + = reinterpret_cast*>(cell->user_pointer()); + body_force.vector_value_list (fe_values.get_quadrature_points(), + body_force_values); + for (unsigned int i=0; i &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 boundary_values; + VectorTools:: + interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (dim), + boundary_values); + VectorTools:: + interpolate_boundary_values (dof_handler, + 1, + IncrementalBoundaryValues(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 + void TopLevel::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 + unsigned int TopLevel::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 + void TopLevel::output_results () const + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + std::vector 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 norm_of_stress (triangulation.n_active_cells()); + { + typename Triangulation::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*>(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 partition_int (triangulation.n_active_cells()); + GridTools::get_subdomain_association (triangulation, partition_int); + const Vector 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 filenames; + for (unsigned int i=0; i > times_and_names; + times_and_names.push_back (std::pair (present_time, pvtu_master_filename)); + std::ofstream pvd_output ("solution.pvd"); + DataOutBase::write_pvd_record (pvd_output, times_and_names); + } + } + template + void TopLevel::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 soln_pt (1.0, 0.0, 3.0); + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + for (; cell != endc; ++cell) + for (unsigned int v=0; v::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 + void TopLevel::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 + void TopLevel::refine_initial_grid () + { + Vector error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss(2), + typename FunctionMap::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 + void TopLevel::move_mesh () + { + pcout << " Moving mesh..." << std::endl; + std::vector vertex_touched (triangulation.n_vertices(), + false); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (); + cell != dof_handler.end(); ++cell) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + Point vertex_displacement; + for (unsigned int d=0; dvertex_dof_index(v,d)); + cell->vertex(v) += vertex_displacement; + } + } + template + void TopLevel::setup_quadrature_point_history () + { + unsigned int our_cells = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned()) + ++our_cells; + triangulation.clear_user_data(); + { + std::vector > tmp; + tmp.swap (quadrature_point_history); + } + quadrature_point_history.resize (our_cells * + quadrature_formula.size()); + unsigned int history_index = 0; + for (typename Triangulation::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 + void TopLevel::update_quadrature_point_history () + { + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients); + std::vector > > + displacement_increment_grads (quadrature_formula.size(), + std::vector >(dim)); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + { + PointHistory *local_quadrature_points_history + = reinterpret_cast *>(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 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 >(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 index 0000000000..91219e30f8 --- /dev/null +++ b/tests/physics/step-18.with_petsc=true.output @@ -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 index 0000000000..3114853361 --- /dev/null +++ b/tests/physics/transformations-rotations_01.cc @@ -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 + +#include +#include +#include + +#include + +#include +#include + +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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/transformations-rotations_01.output @@ -0,0 +1,2 @@ + +DEAL::OK