From 7e618d28ae4bbae15c2d5e6e48a11013eace0df8 Mon Sep 17 00:00:00 2001 From: Mathias Mentler Date: Fri, 3 Mar 2017 09:58:59 +0100 Subject: [PATCH] Added support for vectorization in Physics::Elasticity::Kinematics --- doc/news/changes/minor/20170301MathiasMentler | 4 + .../deal.II/physics/elasticity/kinematics.h | 9 +- tests/physics/elasticity-kinematics_02.cc | 101 ++++++++++++++++++ ...sticity-kinematics_02.with_cxx11=on.output | 2 + 4 files changed, 112 insertions(+), 4 deletions(-) create mode 100644 doc/news/changes/minor/20170301MathiasMentler create mode 100644 tests/physics/elasticity-kinematics_02.cc create mode 100644 tests/physics/elasticity-kinematics_02.with_cxx11=on.output diff --git a/doc/news/changes/minor/20170301MathiasMentler b/doc/news/changes/minor/20170301MathiasMentler new file mode 100644 index 0000000000..80ca434df2 --- /dev/null +++ b/doc/news/changes/minor/20170301MathiasMentler @@ -0,0 +1,4 @@ +New: Reworked Physics::Elasticity::Kinematics +that it is usable with vectorization +
+(Mathias Mentler, 2017/03/01) diff --git a/include/deal.II/physics/elasticity/kinematics.h b/include/deal.II/physics/elasticity/kinematics.h index 3e45a198ec..bc7913b537 100644 --- a/include/deal.II/physics/elasticity/kinematics.h +++ b/include/deal.II/physics/elasticity/kinematics.h @@ -19,6 +19,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -299,7 +300,7 @@ inline SymmetricTensor<2, dim, Number> Physics::Elasticity::Kinematics::F_vol (const Tensor<2, dim, Number> &F) { - return Number(std::pow(determinant(F),1.0/dim))*static_cast< SymmetricTensor<2,dim,Number> >(unit_symmetric_tensor()); + return NumberType::value(std::pow(determinant(F),1.0/dim))*static_cast< SymmetricTensor<2,dim,Number> >(unit_symmetric_tensor()); } @@ -329,7 +330,7 @@ inline SymmetricTensor<2, dim, Number> Physics::Elasticity::Kinematics::E (const Tensor<2, dim, Number> &F) { - return Number(0.5)*(C(F) - static_cast >(StandardTensors::I)); + return NumberType::value(0.5)*(C(F) - static_cast >(StandardTensors::I)); } @@ -351,7 +352,7 @@ SymmetricTensor<2, dim, Number> Physics::Elasticity::Kinematics::e (const Tensor<2, dim, Number> &F) { const Tensor<2, dim, Number> F_inv = invert(F); - return Number(0.5)*symmetrize(static_cast >(StandardTensors::I) - transpose(F_inv)*F_inv); + return NumberType::value(0.5)*symmetrize(static_cast >(StandardTensors::I) - transpose(F_inv)*F_inv); } @@ -390,7 +391,7 @@ Physics::Elasticity::Kinematics::w ( // This could be implemented as w = l-d, but that would mean computing "l" // a second time. const Tensor<2,dim> grad_v = l(F,dF_dt); - return 0.5*(grad_v - transpose(grad_v)) ; + return NumberType::value*(grad_v - transpose(grad_v)) ; } #endif // DOXYGEN diff --git a/tests/physics/elasticity-kinematics_02.cc b/tests/physics/elasticity-kinematics_02.cc new file mode 100644 index 0000000000..bcccbb5fb0 --- /dev/null +++ b/tests/physics/elasticity-kinematics_02.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 vectorization capabilities of the dealii::Physics::Elasticity::kinematics quantities. + +#include "../tests.h" + +#include +#include +#include +#include +#include + +#include + +using namespace dealii; +using namespace dealii::Physics; +using namespace dealii::Physics::Elasticity; + +int main() +{ + std::ofstream logfile("output"); + deallog << std::setprecision(3); + deallog.attach(logfile); + + const int dim = 3; + + Tensor<2,dim,VectorizedArray > grad_u; + + // grad_u = [2 -1 0; -1 2 -1; 0 -1 2] + // which is a random s.p.d. matrix -> F will be s.p.d + // -> det(F) > 0 which has to be satisfied for + // F to make any sense (no negative volume). + grad_u[0][0] = 2.0; + grad_u[0][1] = -1.0; + grad_u[0][2] = 0.0; + grad_u[1][0] = -1.0; + grad_u[1][1] = 2.0; + grad_u[1][2] = -1.0; + grad_u[2][0] = 0.0; + grad_u[2][1] = -1.0; + grad_u[2][2] = 2.0; + + // Scale the gradients along the vectorization-index so that each grad_u is unique. + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; v++) + for (unsigned int i = 0; i < dim; i++) + for (unsigned int j = 0; j < dim; j++) + grad_u[i][j][v] *= (v+1); + + Tensor<2,dim,VectorizedArray > F_solution; + F_solution = grad_u; + for (unsigned int i = 0; i < dim; i++) + F_solution[i][i] = F_solution[i][i] + 1.0; + + Tensor<2,dim,VectorizedArray > F_test = Kinematics::F(grad_u); + + // You can't use .norm() on some difference-tensor of the two so we compare element-wise! + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; v++) + for (unsigned int i = 0; i < dim; i++) + for (unsigned int j = 0; j < dim; j++) + if ( F_solution[i][j][v] - F_test[i][j][v] != 0.0) + deallog << "Not OK" << std::endl; + + SymmetricTensor<2,dim,VectorizedArray > E_solution; + E_solution = 0.5*(symmetrize(transpose(F_solution)*F_solution) - static_cast > >(StandardTensors::I)); + + SymmetricTensor<2,dim,VectorizedArray > E_test = Kinematics::E(F_test); + + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; v++) + for (unsigned int i = 0; i < dim; i++) + for (unsigned int j = 0; j < dim; j++) + if ( E_test[i][j][v] - E_solution[i][j][v] != 0.0) + deallog << "Not OK" << std::endl; + + Tensor<2,dim,VectorizedArray > F_iso_solution; + F_iso_solution = std::pow(determinant(F_solution),-1.0/dim)*F_solution; + + Tensor<2,dim,VectorizedArray > F_iso_test; + F_iso_test = Kinematics::F_iso(F_test); + + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; v++) + for (unsigned int i = 0; i < dim; i++) + for (unsigned int j = 0; j < dim; j++) + if ( F_iso_test[i][j][v] - F_iso_solution[i][j][v] != 0.0) + deallog << "Not OK" << std::endl; + + deallog << "OK" << std::endl; + logfile.close(); +} \ No newline at end of file diff --git a/tests/physics/elasticity-kinematics_02.with_cxx11=on.output b/tests/physics/elasticity-kinematics_02.with_cxx11=on.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/elasticity-kinematics_02.with_cxx11=on.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5