From f1050cd3a9744f54b685e0ba54ce0e19eed22d03 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 16 Oct 2017 08:37:15 +0200 Subject: [PATCH] Add Kelvin notation to the physics module. --- doc/doxygen/headers/physics.h | 30 + .../changes/major/20171016Jean-PaulPelteret | 4 + include/deal.II/physics/notation.h | 1463 +++++++++++++++++ tests/physics/notation-kelvin_01.cc | 310 ++++ tests/physics/notation-kelvin_01.output | 2 + tests/physics/notation-kelvin_02.cc | 532 ++++++ tests/physics/notation-kelvin_02.output | 2 + tests/physics/notation-kelvin_03.cc | 295 ++++ tests/physics/notation-kelvin_03.output | 119 ++ 9 files changed, 2757 insertions(+) create mode 100644 doc/news/changes/major/20171016Jean-PaulPelteret create mode 100644 include/deal.II/physics/notation.h create mode 100644 tests/physics/notation-kelvin_01.cc create mode 100644 tests/physics/notation-kelvin_01.output create mode 100644 tests/physics/notation-kelvin_02.cc create mode 100644 tests/physics/notation-kelvin_02.output create mode 100644 tests/physics/notation-kelvin_03.cc create mode 100644 tests/physics/notation-kelvin_03.output diff --git a/doc/doxygen/headers/physics.h b/doc/doxygen/headers/physics.h index eda1175e7f..95577792f1 100644 --- a/doc/doxygen/headers/physics.h +++ b/doc/doxygen/headers/physics.h @@ -28,6 +28,36 @@ namespace Physics { + /** + * Notations that reduce the order of tensors, effectively storing them in some + * sort of consistent compressed storage pattern. An example is storing the + * 6 independent components of $3\times 3$ symmetric tensors of rank 2 as a + * vector with 6 components, and then representing the 36 independent elements + * of symmetric $3\times 3 \times 3\times 3$ tensors of rank 4 (which when + * applied to a symmetric rank-2 tensor yields another symmetric rank-2 tensor) + * as a $6 \times 6$ matrix. + * + * Although this method of representing tensors is most regularly associated with + * the efficient storage of the fourth-order elasticity tensor, with its + * generalization it has wider applicability. This representation is also common + * in the physics, material science and FEM literature. + * + * There are several variations of tensor notation, each a slightly different + * structure. The primary difference between the various forms of tensor notation + * is the weighting prescribed to the various elements of the compressed tensors. + * This wikipedia article has + * some further general insights on this topic. + * + * @ingroup physics + * + * @author Jean-Paul Pelteret, 2017 + */ + namespace Notation + { + + } + /** * A collection of operations to assist in the transformation of tensor * quantities from the reference to spatial configuration, and vice versa. diff --git a/doc/news/changes/major/20171016Jean-PaulPelteret b/doc/news/changes/major/20171016Jean-PaulPelteret new file mode 100644 index 0000000000..1683280470 --- /dev/null +++ b/doc/news/changes/major/20171016Jean-PaulPelteret @@ -0,0 +1,4 @@ +New: Using the Physics::Notation::Kelvin class, it is possible to store and +retrieve tensors and symmetric tensors in or from a compressed format. +
+(Jean-Paul Pelteret, 2017/10/16) diff --git a/include/deal.II/physics/notation.h b/include/deal.II/physics/notation.h new file mode 100644 index 0000000000..120fc9050a --- /dev/null +++ b/include/deal.II/physics/notation.h @@ -0,0 +1,1463 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_physics_notation_h +#define dealii_physics_notation_h + +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace Physics +{ + namespace Notation + { + /** + * @brief A namespace with functions that assist in the conversion of + * vectors and tensors to and from a compressed format using Kelvin notation + * and weighting. + * + * Both Kelvin and Voigt notation adopt the same indexing convention. + * With specific reference to the spatial dimension 3 case, for + * a rank-2 symmetric tensor $\mathbf{S}$ we enumerate its tensor + * components + * @f[ + * \mathbf{S} := \left[ \begin{array}{ccc} + * S_{00} & S_{01} & S_{02} \\ + * S_{10} = S_{01} & S_{11} & S_{12} \\ + * S_{20} = S_{02} & S_{21} = S_{12} & S_{22} + * \end{array}\right] + * \quad \Rightarrow \quad + * \left[ \begin{array}{ccc} + * n = 0 & n = 5 & n = 4 \\ + * sym & n = 1 & n = 3 \\ + * sym & sym & n = 2} + * \end{array}\right] , + * @f] + * where $n$ denotes the Kelvin index for the tensor component, + * while for a general rank-2 tensor $\mathbf{T}$ + * @f[ + * \mathbf{T} := \left[ \begin{array}{ccc} + * T_{00} & T_{01} & T_{02} \\ + * T_{10} & T_{11} & T_{12} \\ + * T_{20} & T_{21} & T_{22} + * \end{array}\right] + * \quad \Rightarrow \quad + * \left[ \begin{array}{ccc} + * n = 0 & n = 5 & n = 4 \\ + * n = 6 & n = 1 & n = 3 \\ + * n = 7 & n = 8 & n = 2} + * \end{array}\right] , + * @f] + * and for a rank-1 tensor $\mathbf{v}$ + * @f[ + * \mathbf{v} := \left[ \begin{array}{ccc} + * v_{0} & v_{1} & v_{2} + * \end{array}\right]^{T} + * \quad \Rightarrow \quad + * \left[ \begin{array}{ccc} + * n = 0 & n = 1 & n = 2 + * \end{array}\right]^{T} . + * @f] + * To summarize, the relationship between tensor and Kelvin indices for both + * the three-dimensional case and the analogously discerned two-dimensional + * case outlined in the following table: + * + * + * + * + * + * + * + * + * + *
Dimension 2 Dimension 3
+ * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
Tensor index pairsKelvin index
000
111
122
213
+ * 		+ * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
Tensor index pairsKelvin index
000
111
222
123
024
015
106
207
218
+ *
+ * + * To illustrate the purpose of this notation, consider the rank-2 symmetric + * tensors $\mathbf{S}$ and $\mathbf{E}$ that are related to one another by + * $\mathbf{S} = \cal{C} : \mathbf{E}$, where the operator $\cal{C}$ is a fourth-order + * symmetric tensor. + * As opposed to the commonly used Voigt notation, Kelvin (or Mandel) notation + * keeps the same definition of the inner product $\mathbf{S} : \mathbf{E}$ + * when both $\mathbf{S}$ and $\mathbf{E}$ are symmetric. In general, the inner product + * of all symmetric and general tensors remain the same regardless of the notation + * with which it is represented. + * + * To achieve these two properties, namely that + * @f[ + * \mathbf{S} = \cal{C} : \mathbf{E} + * \quad \Rightarrow \quad + * \tilde{\mathbf{S}} = \tilde{\cal{C}} \tilde{\mathbf{E}} + * @f] + * and + * @f[ + * \mathbf{S} : \mathbf{E} + * \, \equiv \, + * \tilde{\mathbf{S}} \cdot \tilde{\mathbf{E}} , + * @f] + * it holds that the Kelvin-condensed equivalents of the previously defined + * symmetric tensors, indicated by the $\tilde{\left(\bullet\right)}$, must be + * defined as + * @f[ + * \tilde{\mathbf{S}} + * = \left[ \begin{array}{c} + * S_{00} \\ S_{11} \\ S_{22} \\ \sqrt{2} S_{12} \\ \sqrt{2} S_{02} \\ \sqrt{2} S_{01} + * \end{array}\right] + * \quad \text{and} \quad + * \tilde{\mathbf{E}} + * = \left[ \begin{array}{c} + * E_{00} \\ E_{11} \\ E_{22} \\ \sqrt{2} E_{12} \\ \sqrt{2} E_{02} \\ \sqrt{2} E_{01} + * \end{array}\right] . + * @f] + * The corresponding and consistent condensed fourth-order symmetric tensor is + * @f[ + * \tilde{\cal{C}} + * = \left[ \begin{array}{cccccc} + * \tilde{\cal{C}}_{00} & \tilde{\cal{C}}_{01} & \tilde{\cal{C}}_{02} & \tilde{\cal{C}}_{03} & \tilde{\cal{C}}_{04} & \tilde{\cal{C}}_{05} \\ + * \tilde{\cal{C}}_{10} & \tilde{\cal{C}}_{11} & \tilde{\cal{C}}_{12} & \tilde{\cal{C}}_{13} & \tilde{\cal{C}}_{14} & \tilde{\cal{C}}_{15} \\ + * \tilde{\cal{C}}_{20} & \tilde{\cal{C}}_{21} & \tilde{\cal{C}}_{22} & \tilde{\cal{C}}_{23} & \tilde{\cal{C}}_{24} & \tilde{\cal{C}}_{25} \\ + * \tilde{\cal{C}}_{30} & \tilde{\cal{C}}_{31} & \tilde{\cal{C}}_{32} & \tilde{\cal{C}}_{33} & \tilde{\cal{C}}_{34} & \tilde{\cal{C}}_{35} \\ + * \tilde{\cal{C}}_{40} & \tilde{\cal{C}}_{41} & \tilde{\cal{C}}_{42} & \tilde{\cal{C}}_{43} & \tilde{\cal{C}}_{44} & \tilde{\cal{C}}_{45} \\ + * \tilde{\cal{C}}_{50} & \tilde{\cal{C}}_{51} & \tilde{\cal{C}}_{52} & \tilde{\cal{C}}_{53} & \tilde{\cal{C}}_{54} & \tilde{\cal{C}}_{55} + * \end{array}\right] + * \equiv + * \left[ \begin{array}{cccccc} + * {\cal{C}}_{0000} & {\cal{C}}_{0011} & {\cal{C}}_{0022} & \sqrt{2} {\cal{C}}_{0012} & \sqrt{2} {\cal{C}}_{0002} & \sqrt{2} {\cal{C}}_{0001} \\ + * {\cal{C}}_{1100} & {\cal{C}}_{1111} & {\cal{C}}_{1122} & \sqrt{2} {\cal{C}}_{1112} & \sqrt{2} {\cal{C}}_{1102} & \sqrt{2} {\cal{C}}_{1101} \\ + * {\cal{C}}_{2200} & {\cal{C}}_{2211} & {\cal{C}}_{2222} & \sqrt{2} {\cal{C}}_{2212} & \sqrt{2} {\cal{C}}_{2202} & \sqrt{2} {\cal{C}}_{2201} \\ + * \sqrt{2} {\cal{C}}_{1200} & \sqrt{2} {\cal{C}}_{1211} & \sqrt{2} {\cal{C}}_{1222} & 2 {\cal{C}}_{1212} & 2 {\cal{C}}_{1202} & 2 {\cal{C}}_{1201} \\ + * \sqrt{2} {\cal{C}}_{0200} & \sqrt{2} {\cal{C}}_{0211} & \sqrt{2} {\cal{C}}_{0222} & 2 {\cal{C}}_{0212} & 2 {\cal{C}}_{0202} & 2 {\cal{C}}_{0201} \\ + * \sqrt{2} {\cal{C}}_{0100} & \sqrt{2} {\cal{C}}_{0111} & \sqrt{2} {\cal{C}}_{0122} & 2 {\cal{C}}_{0112} & 2 {\cal{C}}_{0102} & 2 {\cal{C}}_{0101} + * \end{array}\right] . + * @f] + * The mapping from the two Kelvin indices of the FullMatrix $\tilde{\cal{C}}$ to the + * rank-4 SymmetricTensor $\cal{C}$ can be inferred using the table shown above. + * + * An important observation is that both the left-hand side tensor $\tilde{\mathbf{S}}$ + * and right-hand side tensor $\tilde{\mathbf{E}}$ have the same form; this is a property + * that is not present in Voigt notation. + * The various factors introduced into $\tilde{\mathbf{S}}$, $\tilde{\mathbf{E}}$ + * and $\tilde{\cal{C}}$ account for the symmetry of the tensors. The Kelvin + * description of their non-symmetric counterparts include no such factors. + * + * Some useful references that show how this notation works include, amongst others, + * @code{.bib} + * @Article{Nagel2016, + * author = {Nagel, T. and G{\"o}rke, U-J. and Moerman, K. and Kolditz, O.}, + * title = {On advantages of the Kelvin mapping in finite element implementations of deformation processes}, + * journal = {Environmental Earth Sciences}, + * year = {2016}, + * volume = {75}, + * number = {11}, + * pages = {937} + * } + * @endcode + * and + * @code{.bib} + * @Article{Dellinger1998, + * author = {Dellinger, J. and Vasicek, D. and Sondergeld, C.}, + * title = {Kelvin notation for stabilizing elastic-constant inversion}, + * journal = {Revue de l'Institut Fran{\c{c}}ais du P{\'e}trole}, + * year = {1998}, + * volume = {53}, + * number = {5}, + * pages = {709--719}, + * url = {http://sepwww.stanford.edu/oldsep/joe/Reprints/8IWSA.pdf}, + * } + * @endcode + * as well as the online reference found on + * this wikipedia page + * and the unit tests. + * + * @author Jean-Paul Pelteret, 2017 + */ + namespace Kelvin + { + + /** + * Input matrix has incorrect number of rows. + */ + DeclException3 (ExcNotationExcFullMatrixToTensorRowSize2, int, int, int, + << "The number of rows in the input matrix is " << arg1 + << ", but needs to be either " << arg2 + << " or " << arg3 << "."); + + + /** + * Input matrix has incorrect number of rows. + */ + DeclException4 (ExcNotationExcFullMatrixToTensorRowSize3, int, int, int, int, + << "The number of rows in the input matrix is " << arg1 + << ", but needs to be either " << arg2 << "," << arg3 + << ", or " << arg4 << "."); + + + /** + * Input matrix has incorrect number of columns. + */ + DeclException3 (ExcNotationExcFullMatrixToTensorColSize2, int, int, int, + << "The number of columns in the input matrix is " << arg1 + << ", but needs to be either " << arg2 + << " or " << arg3 << "."); + + + /** + * Input matrix has incorrect number of columns. + */ + DeclException4 (ExcNotationExcFullMatrixToTensorColSize3, int, int, int, int, + << "The number of columns in the input matrix is " << arg1 + << ", but needs to be either " << arg2 << "," << arg3 + << ", or " << arg4 << "."); + + + /** + * @name Forward operation: Tensor notation to Kelvin notation + */ +//@{ + + /** + * Convert a scalar value to its compressed vector equivalent. + */ + template + Vector + to_vector (const Number &s); + + + /** + * Convert a rank-0 tensor to its compressed vector equivalent. + */ + template + Vector + to_vector (const Tensor<0,dim,Number> &s); + + + /** + * Convert a rank-1 tensor to its compressed vector equivalent. + */ + template + Vector + to_vector (const Tensor<1,dim,Number> &v); + + + /** + * Convert a rank-2 tensor to its compressed vector equivalent. + */ + template + Vector + to_vector (const Tensor<2,dim,Number> &t); + + + /** + * Convert a rank-2 symmetric tensor to its compressed vector equivalent. + */ + template + Vector + to_vector (const SymmetricTensor<2,dim,Number> &st); + + + /** + * Convert a scalar value to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const Number &s); + + + /** + * Convert a rank-0 tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const Tensor<0,dim,Number> &s); + + + /** + * Convert a rank-1 tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const Tensor<1,dim,Number> &v); + + + /** + * Convert a rank-2 tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const Tensor<2,dim,Number> &t); + + + /** + * Convert a rank-2 symmetric tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const SymmetricTensor<2,dim,Number> &st); + + /** + * Convert a rank-3 tensor to its compressed matrix equivalent. + * + * The template arguments @p SubTensor1 and @p SubTensor2 determine how + * the unrolling occurs, in particular how the elements of the rank-3 + * tensor are to be interpreted. + * + * So, for example, with the following two conversions + * @code + * Tensor<3,dim> r3_tnsr; // All elements filled differently + * Tensor<3,dim> r3_symm_tnsr; // Some elements filled symmetrically + * + * const FullMatrix mtrx_1 + * = Physics::Notation::to_matrix,Tensor<1,dim> >(r3_tnsr); + * const FullMatrix mtrx_2 + * = Physics::Notation::to_matrix,SymmetricTensor<2,dim> >(r3_symm_tnsr); + * @endcode + * the matrix @p mtrx_1 will have $dim \times dim$ rows and $dim$ columns + * (i.e. size Tensor<2,dim>::n_independent_components $\times$ Tensor<1,dim>::n_independent_components), + * while those of the matrix @p mtrx_2 will have $dim$ rows and + * $(dim \times dim + dim)/2$ columns + * (i.e. size Tensor<1,dim>::n_independent_components $\times$ SymmetricTensor<2,dim>::n_independent_components), + * as it is assumed that the entries corresponding to the alternation of the + * second and third indices are equal. + * That is to say that r3_symm_tnsr[i][j][k] == r3_symm_tnsr[i][k][j]. + */ + template, + typename SubTensor2 = Tensor<1,dim>, + typename Number> + FullMatrix + to_matrix (const Tensor<3,dim,Number> &t); + + + /** + * Convert a rank-4 tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const Tensor<4,dim,Number> &t); + + + /** + * Convert a rank-4 symmetric tensor to its compressed matrix equivalent. + */ + template + FullMatrix + to_matrix (const SymmetricTensor<4,dim,Number> &st); + +//@} + + /** + * @name Reverse operation: Kelvin notation to tensor notation + */ +//@{ + + /** + * Convert a compressed vector to its equivalent scalar value. + */ + template + void + to_tensor (Number &s, + const Vector &vec); + + + /** + * Convert a compressed vector to its equivalent rank-0 tensor. + */ + template + void + to_tensor (Tensor<0,dim,Number> &s, + const Vector &vec); + + + /** + * Convert a compressed vector to its equivalent rank-1 tensor. + */ + template + void + to_tensor (Tensor<1,dim,Number> &v, + const Vector &vec); + + + /** + * Convert a compressed vector to its equivalent rank-2 tensor. + */ + template + void + to_tensor (Tensor<2,dim,Number> &t, + const Vector &vec); + + + /** + * Convert a compressed vector to its equivalent rank-2 symmetric tensor. + */ + template + void + to_tensor (SymmetricTensor<2,dim,Number> &st, + const Vector &vec); + + + /** + * Convert a compressed matrix to its equivalent scalar value. + */ + template + void + to_tensor (Number &s, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-0 tensor. + */ + template + void + to_tensor (Tensor<0,dim,Number> &s, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-1 tensor. + */ + template + void + to_tensor (Tensor<1,dim,Number> &v, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-2 tensor. + */ + template + void + to_tensor (Tensor<2,dim,Number> &t, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-2 symmetric tensor. + */ + template + void + to_tensor (SymmetricTensor<2,dim,Number> &st, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-3 tensor. + * + * @note Based on the size of the matrix @p mtrx, some of the + * components of @p t may be interpreted as having symmetric + * counterparts. This is the reverse of the operation explained + * in the documentation of the counterpart to_matrix() + * function. + */ + template + void + to_tensor (Tensor<3,dim,Number> &t, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-4 tensor. + */ + template + void + to_tensor (Tensor<4,dim,Number> &t, + const FullMatrix &mtrx); + + + /** + * Convert a compressed matrix to its equivalent rank-4 symmetric tensor. + */ + template + void + to_tensor (SymmetricTensor<4,dim,Number> &st, + const FullMatrix &mtrx); + + + /** + * A generic helper function that will convert a compressed vector + * to its equivalent @p TensorType. + */ + template + TensorType + to_tensor (const Vector &vec); + + + /** + * A generic helper function that will convert a compressed matrix + * to its equivalent @p TensorType. + */ + template + TensorType + to_tensor (const FullMatrix &vec); +//@} + + }; + + } +} + + +#ifndef DOXYGEN + + +// ------------------------- inline functions ------------------------ + + +namespace Physics +{ + namespace Notation + { + namespace Kelvin + { + + namespace internal + { + + /** + * Return the tensor indices + * associated with a condensed component index. The + * @p symmetric flag indicates whether or not the + * @p component_n index is associated with a tensor that + * has symmetric entries. + */ + template + std::pair + indices_from_component (const unsigned int component_n, + const bool symmetric); + + + template + std::pair + indices_from_component (const unsigned int component_n, + const bool) + { + AssertThrow(false, ExcNotImplemented()); + return std::pair (); + } + + + template<> + inline + std::pair + indices_from_component<1> (const unsigned int component_n, + const bool) + { + Assert(component_n < 1, ExcIndexRange(component_n,0,1)); + + return std::make_pair(0u,0u); + } + + + template<> + inline + std::pair + indices_from_component<2> (const unsigned int component_n, + const bool symmetric) + { + if (symmetric == true) + { + Assert((component_n < SymmetricTensor<2,2>::n_independent_components), + ExcIndexRange(component_n,0,SymmetricTensor<2,2>::n_independent_components)); + } + else + { + Assert((component_n < Tensor<2,2>::n_independent_components), + ExcIndexRange(component_n,0,Tensor<2,2>::n_independent_components)); + } + + static const unsigned int indices[4][2] = + { + {0,0}, {1,1}, + {0,1}, {1,0} + }; + return std::make_pair(indices[component_n][0], indices[component_n][1]); + } + + + template<> + inline + std::pair + indices_from_component<3> (const unsigned int component_n, + const bool symmetric) + { + if (symmetric == true) + { + Assert((component_n < SymmetricTensor<2,3>::n_independent_components), + ExcIndexRange(component_n,0,SymmetricTensor<2,3>::n_independent_components)); + } + else + { + Assert((component_n < Tensor<2,3>::n_independent_components), + ExcIndexRange(component_n,0,Tensor<2,3>::n_independent_components)); + } + + static const unsigned int indices[9][2] = + { + {0,0}, {1,1}, {2,2}, + {1,2}, {0,2}, {0,1}, + {1,0}, {2,0}, {2,1} + }; + return std::make_pair(indices[component_n][0], indices[component_n][1]); + } + + + /** + * Return the scaling factor to be applied to the + * entry in the condensed vector. + */ + template + double + vector_component_factor (const unsigned int component_i, + const bool symmetric) + { + if (symmetric == false) + return 1.0; + + if (component_i < dim) + return 1.0; + else + return numbers::SQRT2; + } + + + /** + * Return the scaling factor to be applied to the + * entry in the condensed matrix. + */ + template + double + matrix_component_factor (const unsigned int component_i, + const unsigned int component_j, + const bool symmetric) + { + if (symmetric == false) + return 1.0; + + // This case check returns equivalent of this result: + // internal::vector_component_factor(component_i,symmetric)*internal::vector_component_factor(component_j,symmetric); + if (component_i < dim && component_j < dim) + return 1.0; + else if (component_i >= dim && component_j >= dim) + return 2.0; + else // ((component_i >= dim && component_j < dim) || (component_i < dim && component_j >= dim)) + return numbers::SQRT2; + } + + } + + + template + Vector + to_vector (const Number &s) + { + Vector out (1); + const unsigned int n_rows = out.size(); + for (unsigned int r=0; r + Vector + to_vector (const Tensor<0,dim,Number> &s) + { + return to_vector(s.operator const Number &()); + } + + + template + Vector + to_vector (const Tensor<1,dim,Number> &v) + { + Vector out (v.n_independent_components); + const unsigned int n_rows = out.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + const unsigned int i = indices.first; + out(r) = v[i]; + } + return out; + } + + + template + Vector + to_vector (const Tensor<2,dim,Number> &t) + { + Vector out (t.n_independent_components); + const unsigned int n_rows = out.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + Assert(indices.second < dim, ExcInternalError()); + const unsigned int i = indices.first; + const unsigned int j = indices.second; + out(r) = t[i][j]; + } + return out; + } + + + template + Vector + to_vector (const SymmetricTensor<2,dim,Number> &st) + { + Vector out (st.n_independent_components); + const unsigned int n_rows = out.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,true); + Assert(indices.first < dim, ExcInternalError()); + Assert(indices.second < dim, ExcInternalError()); + Assert(indices.second >= indices.first, ExcInternalError()); + const unsigned int i = indices.first; + const unsigned int j = indices.second; + + const double factor = internal::vector_component_factor(r,true); + + out(r) = factor*st[i][j]; + } + return out; + } + + + template + FullMatrix + to_matrix (const Number &s) + { + FullMatrix out (1,1); + out(0,0) = s; + return out; + } + + + template + FullMatrix + to_matrix (const Tensor<0,dim,Number> &s) + { + return to_matrix(s.operator const Number &()); + } + + + template + FullMatrix + to_matrix (const Tensor<1,dim,Number> &v) + { + FullMatrix out (v.n_independent_components,1); + const unsigned int n_rows = out.m(); + const unsigned int n_cols = out.n(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + const unsigned int i = indices.first; + + for (unsigned int c=0; c + FullMatrix + to_matrix (const Tensor<2,dim,Number> &t) + { + FullMatrix out (dim,dim); + const unsigned int n_rows = out.m(); + const unsigned int n_cols = out.n(); + for (unsigned int r=0; r indices_i = internal::indices_from_component(r,false); + Assert(indices_i.first < dim, ExcInternalError()); + Assert(indices_i.second < dim, ExcInternalError()); + const unsigned int i = indices_i.first; + + for (unsigned int c=0; c indices_j = internal::indices_from_component(c,false); + Assert(indices_j.first < dim, ExcInternalError()); + Assert(indices_j.second < dim, ExcInternalError()); + const unsigned int j = indices_j.second; + + out(r,c) = t[i][j]; + } + } + return out; + } + + + template + FullMatrix + to_matrix (const SymmetricTensor<2,dim,Number> &st) + { + return to_matrix(Tensor<2,dim,Number>(st)); + } + + + namespace internal + { + namespace + { + template + struct is_rank_2_symmetric_tensor : std::false_type + {}; + + template + struct is_rank_2_symmetric_tensor > + : std::true_type + {}; + } + } + + + template + FullMatrix + to_matrix (const Tensor<3,dim,Number> &t) + { + static_assert( + (SubTensor1::dimension == dim && SubTensor2::dimension == dim), + "Sub-tensor spatial dimension is different from those of the input tensor."); + + static_assert( + (SubTensor1::rank == 2 && SubTensor2::rank == 1) || + (SubTensor1::rank == 1 && SubTensor2::rank == 2), + "Cannot build a rank 3 tensor from the given combination of sub-tensors."); + + FullMatrix out (SubTensor1::n_independent_components,SubTensor2::n_independent_components); + const unsigned int n_rows = out.m(); + const unsigned int n_cols = out.n(); + + if (SubTensor1::rank == 2 && SubTensor2::rank == 1) + { + const bool subtensor_is_rank_2_symmetric_tensor = internal::is_rank_2_symmetric_tensor::value; + + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,subtensor_is_rank_2_symmetric_tensor); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + if (subtensor_is_rank_2_symmetric_tensor) + { + Assert(indices_ij.second >= indices_ij.first, ExcInternalError()); + } + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + const double factor = internal::vector_component_factor(r,subtensor_is_rank_2_symmetric_tensor); + + for (unsigned int c=0; c indices_k = internal::indices_from_component(c,false); + Assert(indices_k.first < dim, ExcInternalError()); + const unsigned int k = indices_k.first; + + if (subtensor_is_rank_2_symmetric_tensor) + out(r,c) = factor*t[i][j][k]; + else + out(r,c) = t[i][j][k]; + } + } + } + else if (SubTensor1::rank == 1 && SubTensor2::rank == 2) + { + const bool subtensor_is_rank_2_symmetric_tensor = internal::is_rank_2_symmetric_tensor::value; + + for (unsigned int r=0; r indices_k = internal::indices_from_component(r,false); + Assert(indices_k.first < dim, ExcInternalError()); + const unsigned int k = indices_k.first; + + for (unsigned int c=0; c indices_ij = internal::indices_from_component(c,subtensor_is_rank_2_symmetric_tensor); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + if (subtensor_is_rank_2_symmetric_tensor) + { + Assert(indices_ij.second >= indices_ij.first, ExcInternalError()); + } + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + if (subtensor_is_rank_2_symmetric_tensor) + { + const double factor = internal::vector_component_factor(c,subtensor_is_rank_2_symmetric_tensor); + out(r,c) = factor*t[k][i][j]; + } + else + out(r,c) = t[k][i][j]; + } + } + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + + return out; + } + + + template + FullMatrix + to_matrix (const Tensor<4,dim,Number> &t) + { + FullMatrix out (Tensor<2,dim,Number>::n_independent_components, + Tensor<2,dim,Number>::n_independent_components); + const unsigned int n_rows = out.m(); + const unsigned int n_cols = out.n(); + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,false); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + for (unsigned int c=0; c indices_kl = internal::indices_from_component(c,false); + Assert(indices_kl.first < dim, ExcInternalError()); + Assert(indices_kl.second < dim, ExcInternalError()); + const unsigned int k = indices_kl.first; + const unsigned int l = indices_kl.second; + + out(r,c) = t[i][j][k][l]; + } + } + return out; + } + + + template + FullMatrix + to_matrix (const SymmetricTensor<4,dim,Number> &st) + { + FullMatrix out (SymmetricTensor<2,dim,Number>::n_independent_components, + SymmetricTensor<2,dim,Number>::n_independent_components); + const unsigned int n_rows = out.m(); + const unsigned int n_cols = out.n(); + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,true); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + Assert(indices_ij.second >= indices_ij.first, ExcInternalError()); + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + for (unsigned int c=0; c indices_kl = internal::indices_from_component(c,true); + Assert(indices_kl.first < dim, ExcInternalError()); + Assert(indices_kl.second < dim, ExcInternalError()); + Assert(indices_kl.second >= indices_kl.first, ExcInternalError()); + const unsigned int k = indices_kl.first; + const unsigned int l = indices_kl.second; + + const double factor = internal::matrix_component_factor(r,c,true); + + out(r,c) = factor*st[i][j][k][l]; + } + } + return out; + } + + + template + void + to_tensor (Number &s, + const Vector &vec) + { + Assert(vec.size() == 1, ExcDimensionMismatch(vec.size(), 1)); + s = vec(0); + } + + + template + void + to_tensor (Tensor<0,dim,Number> &s, + const Vector &vec) + { + return to_tensor(s.operator Number &(), vec); + } + + + template + void + to_tensor (Tensor<1,dim,Number> &v, + const Vector &vec) + { + Assert(vec.size() == v.n_independent_components, + ExcDimensionMismatch(vec.size(), v.n_independent_components)); + const unsigned int n_rows = vec.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + const unsigned int i = indices.first; + v[i] = vec(r); + } + } + + + template + void + to_tensor (Tensor<2,dim,Number> &t, + const Vector &vec) + { + Assert(vec.size() == t.n_independent_components, + ExcDimensionMismatch(vec.size(), t.n_independent_components)); + const unsigned int n_rows = vec.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + Assert(indices.second < dim, ExcInternalError()); + const unsigned int i = indices.first; + const unsigned int j = indices.second; + t[i][j] = vec(r); + } + } + + + template + void + to_tensor (SymmetricTensor<2,dim,Number> &st, + const Vector &vec) + { + Assert(vec.size() == st.n_independent_components, + ExcDimensionMismatch(vec.size(), st.n_independent_components)); + const unsigned int n_rows = vec.size(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,true); + Assert(indices.first < dim, ExcInternalError()); + Assert(indices.second < dim, ExcInternalError()); + Assert(indices.second >= indices.first, ExcInternalError()); + const unsigned int i = indices.first; + const unsigned int j = indices.second; + + const double inv_factor = 1.0/internal::vector_component_factor(r,true); + + st[i][j] = inv_factor*vec(r); + } + } + + + template + void + to_tensor (Number &s, + const FullMatrix &mtrx) + { + Assert(mtrx.m() == 1, ExcDimensionMismatch(mtrx.m(), 1)); + Assert(mtrx.n() == 1, ExcDimensionMismatch(mtrx.n(), 1)); + Assert(mtrx.n_elements() == 1, ExcDimensionMismatch(mtrx.n_elements(), 1)); + s = mtrx(0,0); + } + + + template + void + to_tensor (Tensor<0,dim,Number> &s, + const FullMatrix &mtrx) + { + return to_tensor(s.operator Number &(), mtrx); + } + + + template + void + to_tensor (Tensor<1,dim,Number> &v, + const FullMatrix &mtrx) + { + Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim)); + Assert(mtrx.n() == 1, ExcDimensionMismatch(mtrx.n(), 1)); + Assert(mtrx.n_elements() == v.n_independent_components, + ExcDimensionMismatch(mtrx.n_elements(), v.n_independent_components)); + + const unsigned int n_rows = mtrx.m(); + const unsigned int n_cols = mtrx.n(); + for (unsigned int r=0; r indices = internal::indices_from_component(r,false); + Assert(indices.first < dim, ExcInternalError()); + Assert(indices.second == 0, ExcInternalError()); + const unsigned int i = indices.first; + + for (unsigned int c=0; c + void + to_tensor (Tensor<2,dim,Number> &t, + const FullMatrix &mtrx) + { + Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim)); + Assert(mtrx.n() == dim, ExcDimensionMismatch(mtrx.n(), dim)); + Assert(mtrx.n_elements() == t.n_independent_components, + ExcDimensionMismatch(mtrx.n_elements(), t.n_independent_components)); + + const unsigned int n_rows = mtrx.m(); + const unsigned int n_cols = mtrx.n(); + for (unsigned int r=0; r indices_i = internal::indices_from_component(r,false); + Assert(indices_i.first < dim, ExcInternalError()); + Assert(indices_i.second < dim, ExcInternalError()); + const unsigned int i = indices_i.first; + + for (unsigned int c=0; c indices_j = internal::indices_from_component(c,false); + Assert(indices_j.first < dim, ExcInternalError()); + Assert(indices_j.second < dim, ExcInternalError()); + const unsigned int j = indices_j.second; + + t[i][j] = mtrx(r,c); + } + } + + } + + + template + void + to_tensor (SymmetricTensor<2,dim,Number> &st, + const FullMatrix &mtrx) + { + // Its impossible to fit the (dim^2 + dim)/2 entries into a square matrix + // We therefore assume that its been converted to a standard tensor format + // using to_matrix (SymmetricTensor<2,dim,Number>) at some point... + Assert(mtrx.m() == dim, ExcDimensionMismatch(mtrx.m(), dim)); + Assert(mtrx.n() == dim, ExcDimensionMismatch(mtrx.n(), dim)); + Assert((mtrx.n_elements() == Tensor<2,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.n_elements(), Tensor<2,dim,Number>::n_independent_components)); + + Tensor<2,dim,Number> tmp; + to_tensor(tmp,mtrx); + st = symmetrize(tmp); + Assert((Tensor<2,dim,Number>(st) - tmp).norm() < 1e-12, + ExcMessage("The entries stored inside the matrix were not symmetric")); + } + + + template + void + to_tensor (Tensor<3,dim,Number> &t, + const FullMatrix &mtrx) + { + Assert((mtrx.m() == Tensor<1,dim,Number>::n_independent_components) || + (mtrx.m() == Tensor<2,dim,Number>::n_independent_components) || + (mtrx.m() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcNotationExcFullMatrixToTensorColSize3( + mtrx.m(), + Tensor<1,dim,Number>::n_independent_components, + Tensor<2,dim,Number>::n_independent_components, + SymmetricTensor<2,dim,Number>::n_independent_components)); + Assert((mtrx.n() == Tensor<1,dim,Number>::n_independent_components) || + (mtrx.n() == Tensor<2,dim,Number>::n_independent_components) || + (mtrx.n() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcNotationExcFullMatrixToTensorColSize3( + mtrx.n(), + Tensor<1,dim,Number>::n_independent_components, + Tensor<2,dim,Number>::n_independent_components, + SymmetricTensor<2,dim,Number>::n_independent_components)); + + const unsigned int n_rows = mtrx.m(); + const unsigned int n_cols = mtrx.n(); + if (mtrx.n() == Tensor<1,dim,Number>::n_independent_components) + { + Assert((mtrx.m() == Tensor<2,dim,Number>::n_independent_components) || + (mtrx.m() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcNotationExcFullMatrixToTensorRowSize2( + mtrx.m(), + Tensor<2,dim,Number>::n_independent_components, + SymmetricTensor<2,dim,Number>::n_independent_components + )); + + const bool subtensor_is_rank_2_symmetric_tensor + = (mtrx.m() == SymmetricTensor<2,dim,Number>::n_independent_components); + + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,subtensor_is_rank_2_symmetric_tensor); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + if (subtensor_is_rank_2_symmetric_tensor) + { + Assert(indices_ij.second >= indices_ij.first, ExcInternalError()); + } + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + const double inv_factor = 1.0/internal::vector_component_factor(r,subtensor_is_rank_2_symmetric_tensor); + + for (unsigned int c=0; c indices_k = internal::indices_from_component(c,false); + Assert(indices_k.first < dim, ExcInternalError()); + const unsigned int k = indices_k.first; + + if (subtensor_is_rank_2_symmetric_tensor) + { + t[i][j][k] = inv_factor*mtrx(r,c); + t[j][i][k] = t[i][j][k]; + } + else + t[i][j][k] = mtrx(r,c); + } + } + } + else + { + Assert((mtrx.m() == Tensor<1,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.m(), Tensor<1,dim,Number>::n_independent_components)); + Assert((mtrx.n() == Tensor<2,dim,Number>::n_independent_components) || + (mtrx.n() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcNotationExcFullMatrixToTensorColSize2( + mtrx.n(), + Tensor<2,dim,Number>::n_independent_components, + SymmetricTensor<2,dim,Number>::n_independent_components + )); + + const bool subtensor_is_rank_2_symmetric_tensor + = (mtrx.n() == SymmetricTensor<2,dim,Number>::n_independent_components); + + for (unsigned int r=0; r indices_k = internal::indices_from_component(r,false); + Assert(indices_k.first < dim, ExcInternalError()); + const unsigned int k = indices_k.first; + + for (unsigned int c=0; c indices_ij = internal::indices_from_component(c,subtensor_is_rank_2_symmetric_tensor); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + if (subtensor_is_rank_2_symmetric_tensor) + { + Assert(indices_ij.second >= indices_ij.first, ExcInternalError()); + } + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + if (subtensor_is_rank_2_symmetric_tensor) + { + const double inv_factor = 1.0/internal::vector_component_factor(c,subtensor_is_rank_2_symmetric_tensor); + t[k][i][j] = inv_factor*mtrx(r,c); + t[k][j][i] = t[k][i][j]; + } + else + t[k][i][j] = mtrx(r,c); + } + } + } + } + + + template + void + to_tensor (Tensor<4,dim,Number> &t, + const FullMatrix &mtrx) + { + Assert((mtrx.m() == Tensor<2,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.m(), Tensor<2,dim,Number>::n_independent_components)); + Assert((mtrx.n() == Tensor<2,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.n(), Tensor<2,dim,Number>::n_independent_components)); + Assert(mtrx.n_elements() == t.n_independent_components, + ExcDimensionMismatch(mtrx.n_elements(), t.n_independent_components)); + + const unsigned int n_rows = mtrx.m(); + const unsigned int n_cols = mtrx.n(); + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,false); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + for (unsigned int c=0; c indices_kl = internal::indices_from_component(c,false); + Assert(indices_kl.first < dim, ExcInternalError()); + Assert(indices_kl.second < dim, ExcInternalError()); + const unsigned int k = indices_kl.first; + const unsigned int l = indices_kl.second; + + t[i][j][k][l] = mtrx(r,c); + } + } + } + + + template + void + to_tensor (SymmetricTensor<4,dim,Number> &st, + const FullMatrix &mtrx) + { + Assert((mtrx.m() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.m(), SymmetricTensor<2,dim,Number>::n_independent_components)); + Assert((mtrx.n() == SymmetricTensor<2,dim,Number>::n_independent_components), + ExcDimensionMismatch(mtrx.n(), SymmetricTensor<2,dim,Number>::n_independent_components)); + Assert(mtrx.n_elements() == st.n_independent_components, + ExcDimensionMismatch(mtrx.n_elements(), st.n_independent_components)); + + const unsigned int n_rows = mtrx.m(); + const unsigned int n_cols = mtrx.n(); + for (unsigned int r=0; r indices_ij = internal::indices_from_component(r,false); + Assert(indices_ij.first < dim, ExcInternalError()); + Assert(indices_ij.second < dim, ExcInternalError()); + const unsigned int i = indices_ij.first; + const unsigned int j = indices_ij.second; + + for (unsigned int c=0; c indices_kl = internal::indices_from_component(c,false); + Assert(indices_kl.first < dim, ExcInternalError()); + Assert(indices_kl.second < dim, ExcInternalError()); + const unsigned int k = indices_kl.first; + const unsigned int l = indices_kl.second; + + const double inv_factor = 1.0/internal::matrix_component_factor(r,c,true); + + st[i][j][k][l] = inv_factor*mtrx(r,c); + } + } + } + + + template + inline TensorType + to_tensor (const Vector &vec) + { + TensorType out; + to_tensor(out, vec); + return out; + } + + + template + inline TensorType + to_tensor (const FullMatrix &mtrx) + { + TensorType out; + to_tensor(out, mtrx); + return out; + } + + } + } +} + + +#endif // DOXYGEN + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/physics/notation-kelvin_01.cc b/tests/physics/notation-kelvin_01.cc new file mode 100644 index 0000000000..c658cd90c4 --- /dev/null +++ b/tests/physics/notation-kelvin_01.cc @@ -0,0 +1,310 @@ +// --------------------------------------------------------------------- +// +// 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 that rolling to, and unrolling from, Kelvin notation. + +#include "../tests.h" + +#include +#include +#include +#include + +#include + + +using namespace dealii; +using namespace dealii::Physics; + +template +void initialize (Tensor<1,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (Tensor<2,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (SymmetricTensor<2,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (Tensor<3,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (Tensor<3,dim,Number> &x, + const bool left_components_are_symmetric) +{ + Tensor<1,dim,Number> v; + initialize(v); + SymmetricTensor<2,dim,Number> st; + initialize(st); + const Tensor<2,dim,Number> t (st); + if (left_components_are_symmetric == true) + x = outer_product(t,v); + else + x = outer_product(v,t); +} + +template +void initialize (Tensor<4,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (SymmetricTensor<4,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void +test_scalars () +{ + const double A = 5; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto vA_conv = Notation::Kelvin::to_tensor(vA); + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert(std::abs(vA_conv - A) < 1e-12, ExcMessage("Different result for vector conversion")); + Assert(std::abs(mA_conv - A) < 1e-12, ExcMessage("Different result for matrix conversion")); +} + +template +void +test_rank_0_tensors () +{ + const Tensor<0,dim,double> A = 5; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto vA_conv = Notation::Kelvin::to_tensor(vA); + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert(std::abs(vA_conv - A) < 1e-12, ExcMessage("Different result for vector conversion")); + Assert(std::abs(mA_conv - A) < 1e-12, ExcMessage("Different result for matrix conversion")); +} + +template +void +test_rank_1_tensors () +{ + Tensor<1,dim,double> A; + initialize(A); + + const Vector vA = Notation::Kelvin::to_vector(A); + + typedef typename std::decay::type InpType; + const auto vA_conv = Notation::Kelvin::to_tensor(vA); + + Assert((vA_conv - A).norm() < 1e-12, ExcMessage("Different result for vector conversion")); +} + +template +void +test_rank_2_tensors () +{ + // Non-symmetric tensor + { + Tensor<2,dim,double> A; + initialize(A); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto vA_conv = Notation::Kelvin::to_tensor(vA); + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((vA_conv - A).norm() < 1e-12, ExcMessage("Different result for vector conversion")); + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } + + // Symmetric tensor + { + SymmetricTensor<2,dim,double> A; + initialize(A); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto vA_conv = Notation::Kelvin::to_tensor(vA); + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((vA_conv - A).norm() < 1e-12, ExcMessage("Different result for vector conversion")); + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } +} + +template +void +test_rank_3_tensors () +{ + // Non-symmetric tensor: Version 1 + { + Tensor<3,dim,double> A; + initialize(A); + + const FullMatrix mA = Notation::Kelvin::to_matrix< dim, Tensor<1,dim>, Tensor<2,dim> >(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } + + // Non-symmetric tensor: Version 2 + { + Tensor<3,dim,double> A; + initialize(A); + + const FullMatrix mA = Notation::Kelvin::to_matrix< dim, Tensor<2,dim>, Tensor<1,dim> >(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } + + // Symmetric tensor: Version 1 + { + Tensor<3,dim,double> A; + initialize(A,true); // Specialised constructor + + const FullMatrix mA = Notation::Kelvin::to_matrix< dim, SymmetricTensor<2,dim>, Tensor<1,dim> >(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } + + // Symmetric tensor: Version 2 + { + Tensor<3,dim,double> A; + initialize(A,false); // Specialised constructor + + const FullMatrix mA = Notation::Kelvin::to_matrix< dim, Tensor<1,dim>, SymmetricTensor<2,dim> >(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } +} + +template +void +test_rank_4_tensors () +{ + // Non-symmetric tensor + { + Tensor<4,dim,double> A; + initialize(A); + + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } + + // Symmetric tensor + { + SymmetricTensor<4,dim,double> A; + initialize(A); + + const FullMatrix mA = Notation::Kelvin::to_matrix(A); + + typedef typename std::decay::type InpType; + const auto mA_conv = Notation::Kelvin::to_tensor(mA); + + Assert((mA_conv - A).norm() < 1e-12, ExcMessage("Different result for matrix conversion")); + } +} + +template +void +test_tensors () +{ + test_scalars(); + test_rank_0_tensors(); + test_rank_1_tensors(); + test_rank_2_tensors(); + test_rank_3_tensors(); + test_rank_4_tensors(); +} + +int main () +{ + initlog(); + + test_tensors<2> (); + test_tensors<3> (); + + deallog << "OK" << std::endl; +} diff --git a/tests/physics/notation-kelvin_01.output b/tests/physics/notation-kelvin_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/notation-kelvin_01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/physics/notation-kelvin_02.cc b/tests/physics/notation-kelvin_02.cc new file mode 100644 index 0000000000..350e551686 --- /dev/null +++ b/tests/physics/notation-kelvin_02.cc @@ -0,0 +1,532 @@ +// --------------------------------------------------------------------- +// +// 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 that Kelvin notation works as expected. The following properties +// should hold: +// - B = C:A --> kelvin(B) = kelvin(C)*kelvin(A) +// - A:A == kelvin(A)*kelvin(A) + +#include "../tests.h" + +#include +#include +#include +#include + +#include + + +using namespace dealii; +using namespace dealii::Physics; + +template +void initialize (Tensor<1,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (Tensor<2,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (SymmetricTensor<2,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (Tensor<3,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (Tensor<3,dim,Number> &x, + const bool left_components_are_symmetric) +{ + Tensor<1,dim,Number> v; + initialize(v); + SymmetricTensor<2,dim,Number> st; + initialize(st); + const Tensor<2,dim,Number> t (st); + if (left_components_are_symmetric == true) + x = outer_product(t,v); + else + x = outer_product(v,t); +} + +template +void initialize (Tensor<4,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void initialize (SymmetricTensor<4,dim,Number> &x) +{ + unsigned int c=1; + for (unsigned int i=0; i +void +test_scalars () +{ + const double A = 5; + const double C = 12; + const double B = C*A; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Scalar" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert(std::abs(A_conv - A) < 1e-12, ExcMessage("Different result for input A")); + Assert(std::abs(C_conv - C) < 1e-12, ExcMessage("Different result for input C")); + Assert(std::abs(B_conv - B) < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - std::abs(A)) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - std::abs(C)) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - std::abs(B)) < 1e-12, ExcMessage("Different norm for output B")); +} + +template +void +test_rank_0_tensors () +{ + const Tensor<0,dim,double> A = 5; + const Tensor<0,dim,double> C = 12; + const Tensor<0,dim,double> B = C*A; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 0" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); +} + +template +void +test_rank_1_2_tensors () +{ + // Non-symmetric tensor + { + Tensor<1,dim,double> A; + Tensor<2,dim,double> C; + initialize(A); + initialize(C); + const Tensor<1,dim,double> B = C*A; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 1 (non-symm)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Symmetric tensor + { + Tensor<1,dim,double> A; + SymmetricTensor<2,dim,double> C; + initialize(A); + initialize(C); + const Tensor<1,dim,double> B = C*A; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 1 (symm)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } +} + +template +void +test_rank_2_4_tensors () +{ + // Non-symmetric tensor + { + Tensor<2,dim,double> A; + Tensor<4,dim,double> C; + initialize(A); + initialize(C); + const Tensor<2,dim,double> B = double_contract<2,0,3,1>(C,A); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 2 (non-symm)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Symmetric tensor + { + SymmetricTensor<2,dim,double> A; + SymmetricTensor<4,dim,double> C; + initialize(A); + initialize(C); + const SymmetricTensor<2,dim,double> B = C*A; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 2 (symm)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Non-symmetric tensor from symmetric tensor + { + SymmetricTensor<2,dim,double> A_symm; + SymmetricTensor<4,dim,double> C_symm; + initialize(A_symm); + initialize(C_symm); + + Tensor<2,dim,double> A (A_symm); + Tensor<4,dim,double> C (C_symm); + const Tensor<2,dim,double> B = double_contract<2,0,3,1>(C,A); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.m()); + mC.vmult(vB, vA); + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 2 (non-symm from symm)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } +} + +template +void +test_rank_3_tensors () +{ + // Non-symmetric tensor: Version 1 + { + Tensor<2,dim,double> A; + Tensor<3,dim,double> C; + initialize(A); + initialize(C); + const Tensor<1,dim,double> B = double_contract<0,0,1,1>(C,A); // This implies that a Tvmult is necessary + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + Vector vB (mC.n()); // Note result size + mC.Tvmult(vB, vA); // Note transpose vmult + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 3 (non-symm 1)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Non-symmetric tensor: Version 2 + { + Tensor<2,dim,double> A; + Tensor<3,dim,double> C; + initialize(A); + initialize(C); + const Tensor<1,dim,double> B = double_contract<1,0,2,1>(C,A); // This implies that a standard vmult is necessary + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix,Tensor<2,dim>>(C); // Define subtensor representation Tensor<1,dim> \otimes Tensor<2,dim> + Vector vB (mC.m()); // Note result size + mC.vmult(vB, vA); // Note transpose vmult + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 3 (non-symm 2)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Symmetric tensor: Version 1 + { + SymmetricTensor<2,dim,double> A; + Tensor<3,dim,double> C; + initialize(A); + initialize(C,true); // Specialised constructor + const Tensor<2,dim,double> A_ns (A); + const Tensor<1,dim,double> B = double_contract<0,0,1,1>(C,A_ns); // This implies that a Tvmult is necessary + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix, Tensor<1,dim>>(C); // Define subtensor representation SymmetricTensor<2,dim> \otimes Tensor<1,dim> + Vector vB (mC.n()); // Note result size + mC.Tvmult(vB, vA); // Note transpose vmult + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 3 (symm 1)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } + + // Symmetric tensor: Version 2 + { + SymmetricTensor<2,dim,double> A; + Tensor<3,dim,double> C; + initialize(A); + initialize(C,false); // Specialised constructor + const Tensor<2,dim,double> A_ns (A); + const Tensor<1,dim,double> B = double_contract<1,0,2,1>(C,A_ns); // This implies that a standard vmult is necessary + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix,SymmetricTensor<2,dim>>(C); // Define subtensor representation Tensor<1,dim> \otimes SymmetricTensor<2,dim> + Vector vB (mC.m()); // Note result size + mC.vmult(vB, vA); // Note transpose vmult + + typedef typename std::decay::type InpVecType; + typedef typename std::decay::type ResVecType; + typedef typename std::decay::type InpMatType; + const auto A_conv = Notation::Kelvin::to_tensor(vA); + const auto B_conv = Notation::Kelvin::to_tensor(vB); + const auto C_conv = Notation::Kelvin::to_tensor(mC); + + std::cout << "Rank 3 (symm 2)" << std::endl; + std::cout << "A: " << A << " A_conv: " << A_conv << std::endl; + std::cout << "B: " << B << " B_conv: " << B_conv << std::endl; + std::cout << "C: " << C << " C_conv: " << C_conv << std::endl; + + Assert((A_conv - A).norm() < 1e-12, ExcMessage("Different result for input A")); + Assert((C_conv - C).norm() < 1e-12, ExcMessage("Different result for input C")); + Assert((B_conv - B).norm() < 1e-12, ExcMessage("Different result for output B")); + + Assert(std::abs(vA.l2_norm() - A.norm()) < 1e-12, ExcMessage("Different norm for input A")); + Assert(std::abs(mC.frobenius_norm() - C.norm()) < 1e-12, ExcMessage("Different norm for input C")); + Assert(std::abs(vB.l2_norm() - B.norm()) < 1e-12, ExcMessage("Different norm for output B")); + } +} + +template +void +test_tensors () +{ + test_scalars(); + test_rank_0_tensors(); + test_rank_1_2_tensors(); + test_rank_2_4_tensors(); + test_rank_3_tensors(); +} + +int main () +{ + initlog(); + + test_tensors<2> (); + test_tensors<3> (); + + deallog << "OK" << std::endl; +} diff --git a/tests/physics/notation-kelvin_02.output b/tests/physics/notation-kelvin_02.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/physics/notation-kelvin_02.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/physics/notation-kelvin_03.cc b/tests/physics/notation-kelvin_03.cc new file mode 100644 index 0000000000..295eb41f63 --- /dev/null +++ b/tests/physics/notation-kelvin_03.cc @@ -0,0 +1,295 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// Check the scaling factors used during conversion between +// Tensor and Kelvin notation + +#include "../tests.h" + +#include +#include +#include +#include + +#include + + +using namespace dealii; +using namespace dealii::Physics; + +template +void initialize (Tensor<0,dim,Number> &x) +{ + x = 1.0; +} + +template +void initialize (Tensor<1,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (Tensor<2,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (SymmetricTensor<2,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (Tensor<3,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (Tensor<4,dim,Number> &x) +{ + for (unsigned int i=0; i +void initialize (SymmetricTensor<4,dim,Number> &x) +{ + for (unsigned int i=0; i +void +test_scalars () +{ + const double A = 1.0; + const double C = 1.0; + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + + deallog.push("Scalar"); + { + deallog << "Vector" << std::endl; + vA.print(deallog.get_file_stream()); + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); +} + +template +void +test_tensors_012 () +{ + static_assert((rank==0 || rank==1 || rank==2), "Must be rank 0, 1 or 2"); + + Tensor A; + Tensor C; + initialize(A); + initialize(C); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + + const std::string name = std::string("Rank-") + Utilities::int_to_string(rank) + std::string (" tensor"); + deallog.push(name); + { + deallog << "Vector" << std::endl; + vA.print(deallog.get_file_stream()); + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); +} + +template +void +test_tensor_3 () +{ + // Non-symmetric tensor: Version 1 + { + Tensor<3,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix< dim, Tensor<1,dim>, Tensor<2,dim> >(C); + + deallog.push("Rank-3 tensor: NSymm 1"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); + } + + // Non-symmetric tensor: Version 2 + { + Tensor<3,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix< dim, Tensor<2,dim>, Tensor<1,dim> >(C); + + deallog.push("Rank-3 tensor: NSymm 2"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); + } + + // Symmetric tensor: Version 1 + { + Tensor<3,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix< dim, Tensor<1,dim>, SymmetricTensor<2,dim> >(C); + + deallog.push("Rank-3 tensor: Symm 1"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); + } + + // Symmetric tensor: Version 2 + { + Tensor<3,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix< dim, SymmetricTensor<2,dim>, Tensor<1,dim> >(C); + + deallog.push("Rank-3 tensor: Symm 2"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); + } +} + +template +void +test_tensor_4 () +{ + Tensor<4,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + + deallog.push("Rank-4 tensor"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); +} + +template +void +test_symmetric_tensor_2 () +{ + SymmetricTensor<2,dim,double> A; + SymmetricTensor<2,dim,double> C; + initialize(A); + initialize(C); + + const Vector vA = Notation::Kelvin::to_vector(A); + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + + deallog.push("Rank-2 symm tensor"); + { + deallog << "Vector" << std::endl; + vA.print(deallog.get_file_stream()); + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); +} + +template +void +test_symmetric_tensor_4 () +{ + SymmetricTensor<4,dim,double> C; + initialize(C); + + const FullMatrix mC = Notation::Kelvin::to_matrix(C); + + deallog.push("Rank-4 symm tensor"); + { + deallog << "Matrix" << std::endl; + mC.print_formatted(deallog.get_file_stream()); + } + deallog.pop(); +} + +template +void +test_tensors () +{ + deallog.push(Utilities::int_to_string(dim)); + { + test_scalars(); + test_tensors_012<0,dim>(); + test_tensors_012<1,dim>(); + test_tensors_012<2,dim>(); + test_tensor_3(); + test_tensor_4(); + test_symmetric_tensor_2(); + test_symmetric_tensor_4(); + + deallog << "OK" << std::endl; + } + deallog.pop(); +} + +int main () +{ + initlog(); + + test_tensors<2> (); + test_tensors<3> (); + + deallog << "OK" << std::endl; +} diff --git a/tests/physics/notation-kelvin_03.output b/tests/physics/notation-kelvin_03.output new file mode 100644 index 0000000000..4f1863c621 --- /dev/null +++ b/tests/physics/notation-kelvin_03.output @@ -0,0 +1,119 @@ + +DEAL:2:Scalar::Vector +1.000e+00 +DEAL:2:Scalar::Matrix +1.000e+00 +DEAL:2:Rank-0 tensor::Vector +1.000e+00 +DEAL:2:Rank-0 tensor::Matrix +1.000e+00 +DEAL:2:Rank-1 tensor::Vector +1.000e+00 1.000e+00 +DEAL:2:Rank-1 tensor::Matrix +1.000e+00 +1.000e+00 +DEAL:2:Rank-2 tensor::Vector +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:2:Rank-2 tensor::Matrix +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +DEAL:2:Rank-3 tensor: NSymm 1::Matrix +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:2:Rank-3 tensor: NSymm 2::Matrix +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +DEAL:2:Rank-3 tensor: Symm 1::Matrix +1.000e+00 1.000e+00 1.414e+00 +1.000e+00 1.000e+00 1.414e+00 +DEAL:2:Rank-3 tensor: Symm 2::Matrix +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +1.414e+00 1.414e+00 +DEAL:2:Rank-4 tensor::Matrix +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:2:Rank-2 symm tensor::Vector +1.000e+00 1.000e+00 1.414e+00 +DEAL:2:Rank-2 symm tensor::Matrix +1.000e+00 1.000e+00 +1.000e+00 1.000e+00 +DEAL:2:Rank-4 symm tensor::Matrix +1.000e+00 1.000e+00 1.414e+00 +1.000e+00 1.000e+00 1.414e+00 +1.414e+00 1.414e+00 2.000e+00 +DEAL:2::OK +DEAL:3:Scalar::Vector +1.000e+00 +DEAL:3:Scalar::Matrix +1.000e+00 +DEAL:3:Rank-0 tensor::Vector +1.000e+00 +DEAL:3:Rank-0 tensor::Matrix +1.000e+00 +DEAL:3:Rank-1 tensor::Vector +1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-1 tensor::Matrix +1.000e+00 +1.000e+00 +1.000e+00 +DEAL:3:Rank-2 tensor::Vector +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-2 tensor::Matrix +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-3 tensor: NSymm 1::Matrix +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-3 tensor: NSymm 2::Matrix +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-3 tensor: Symm 1::Matrix +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +DEAL:3:Rank-3 tensor: Symm 2::Matrix +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.414e+00 1.414e+00 1.414e+00 +1.414e+00 1.414e+00 1.414e+00 +1.414e+00 1.414e+00 1.414e+00 +DEAL:3:Rank-4 tensor::Matrix +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-2 symm tensor::Vector +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +DEAL:3:Rank-2 symm tensor::Matrix +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +1.000e+00 1.000e+00 1.000e+00 +DEAL:3:Rank-4 symm tensor::Matrix +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +1.000e+00 1.000e+00 1.000e+00 1.414e+00 1.414e+00 1.414e+00 +1.414e+00 1.414e+00 1.414e+00 2.000e+00 2.000e+00 2.000e+00 +1.414e+00 1.414e+00 1.414e+00 2.000e+00 2.000e+00 2.000e+00 +1.414e+00 1.414e+00 1.414e+00 2.000e+00 2.000e+00 2.000e+00 +DEAL:3::OK +DEAL::OK -- 2.39.5