From a157752e6fe84ec034f8a85c31c1ee060925429a Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 14 Apr 2015 13:39:28 +0200 Subject: [PATCH] Make most parts of SymmetricTensor dimension-independent --- include/deal.II/base/symmetric_tensor.h | 371 ++++++++---------------- 1 file changed, 123 insertions(+), 248 deletions(-) diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 685de47e38..e50518ea4e 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -958,7 +958,16 @@ SymmetricTensor::SymmetricTensor (const Tensor<2,dim,Number> &t break; default: - Assert (false, ExcNotImplemented()); + for (unsigned int d=0; d::operator = (const Number d) -// helper function to convert symmetric tensor -// to generic tensor -namespace internal -{ - template - inline - Tensor<2,1,Number> - conversion (const Tensor<1,1,Number> &data) - { - const Number t[1][1] = {{data[0]}}; - return Tensor<2,1,Number>(t); - } - - template - inline - Tensor<2,2,Number> - conversion (const Tensor<1,3,Number> &data) - { - const Number t[2][2] = {{data[0], data[2]}, - {data[2], data[1]} - }; - return Tensor<2,2,Number>(t); - } - - template - inline - Tensor<2,3,Number> - conversion (const Tensor<1,6,Number> &data) - { - const Number t[3][3] = {{data[0], data[3], data[4]}, - {data[3], data[1], data[5]}, - {data[4], data[5], data[2]} - }; - return Tensor<2,3,Number>(t); - } -} - - - template inline SymmetricTensor:: operator Tensor () const { Assert (rank == 2, ExcNotImplemented()); - return internal::conversion(data); + Number t[dim][dim]; + for (unsigned int d=0; d(t); } @@ -1204,19 +1183,12 @@ namespace internal return (data[0] * sdata[0] + data[1] * sdata[1] + 2*data[2] * sdata[2]); - case 3: - return (data[0] * sdata[0] + - data[1] * sdata[1] + - data[2] * sdata[2] + - 2*data[3] * sdata[3] + - 2*data[4] * sdata[4] + - 2*data[5] * sdata[5]); default: - Number sum = 0; + Number sum = Number(); for (unsigned int d=0; d::base_tensor_type &data, const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata) { - Number tmp [SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components]; - switch (dim) - { - case 1: - tmp[0] = data[0][0] * sdata[0]; - break; - case 2: - for (unsigned int i=0; i<3; ++i) - tmp[i] = (data[i][0] * sdata[0] + - data[i][1] * sdata[1] + - 2 * data[i][2] * sdata[2]); - break; - case 3: - for (unsigned int i=0; i<6; ++i) - tmp[i] = (data[i][0] * sdata[0] + - data[i][1] * sdata[1] + - data[i][2] * sdata[2] + - 2 * data[i][3] * sdata[3] + - 2 * data[i][4] * sdata[4] + - 2 * data[i][5] * sdata[5]); - break; - default: - Assert (false, ExcNotImplemented()); - } + const unsigned int data_dim = + SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components; + Number tmp [data_dim]; + for (unsigned int i=0; i(data[i], sdata); return SymmetricTensor<2,dim,Number>(tmp); } @@ -1265,28 +1218,12 @@ namespace internal const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata) { typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type tmp; - switch (dim) + for (unsigned int i=0; i::base_tensor_type &data, const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata) { + const unsigned int data_dim = + SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components; typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type tmp; - switch (dim) - { - case 1: - tmp[0][0] = data[0][0] * sdata[0][0]; - break; - case 2: - for (unsigned int i=0; i<3; ++i) - for (unsigned int j=0; j<3; ++j) - tmp[i][j] = (data[i][0] * sdata[0][j] + - data[i][1] * sdata[1][j] + - 2*data[i][2] * sdata[2][j]); - break; - case 3: - for (unsigned int i=0; i<6; ++i) - for (unsigned int j=0; j<6; ++j) - tmp[i][j] = (data[i][0] * sdata[0][j] + - data[i][1] * sdata[1][j] + - data[i][2] * sdata[2][j] + - 2*data[i][3] * sdata[3][j] + - 2*data[i][4] * sdata[4][j] + - 2*data[i][5] * sdata[5][j]); - break; - default: - Assert (false, ExcNotImplemented()); - } + for (unsigned int i=0; i &indices, typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data) { + // 1d is very simple and done first + if (dim == 1) + return data[0]; + + // first treat the main diagonal elements, which are stored consecutively + // at the beginning + if (indices[0] == indices[1]) + return data[indices[0]]; + + // the rest is messier and requires a few switches. switch (dim) { - case 1: - return data[0]; - case 2: - // first treat the main diagonal - // elements, which are stored - // consecutively at the beginning - if (indices[0] == indices[1]) - return data[indices[0]]; - - // the rest is messier and requires a few - // switches. at least for the 2x2 case it - // is reasonably simple + // at least for the 2x2 case it is reasonably simple Assert (((indices[0]==1) && (indices[1]==0)) || ((indices[0]==0) && (indices[1]==1)), ExcInternalError()); return data[2]; - case 3: - // first treat the main diagonal - // elements, which are stored - // consecutively at the beginning - if (indices[0] == indices[1]) - return data[indices[0]]; - - // the rest is messier and requires a few - // switches, but simpler if we just sort - // our indices - { - TableIndices<2> sorted_indices (indices); - sorted_indices.sort (); - - if ((sorted_indices[0]==0) && (sorted_indices[1]==1)) - return data[3]; - else if ((sorted_indices[0]==0) && (sorted_indices[1]==2)) - return data[4]; - else if ((sorted_indices[0]==1) && (sorted_indices[1]==2)) - return data[5]; - else - Assert (false, ExcInternalError()); - } + default: + // to do the rest, sort our indices before comparing + { + TableIndices<2> sorted_indices (indices); + sorted_indices.sort (); + + for (unsigned int d=0, c=0; d &indices, const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data) { + // 1d is very simple and done first + if (dim == 1) + return data[0]; + + // first treat the main diagonal elements, which are stored consecutively + // at the beginning + if (indices[0] == indices[1]) + return data[indices[0]]; + + // the rest is messier and requires a few switches. switch (dim) { - case 1: - return data[0]; - case 2: - // first treat the main diagonal - // elements, which are stored - // consecutively at the beginning - if (indices[0] == indices[1]) - return data[indices[0]]; - - // the rest is messier and requires a few - // switches. at least for the 2x2 case it - // is reasonably simple + // at least for the 2x2 case it is reasonably simple Assert (((indices[0]==1) && (indices[1]==0)) || ((indices[0]==0) && (indices[1]==1)), ExcInternalError()); return data[2]; - case 3: - // first treat the main diagonal - // elements, which are stored - // consecutively at the beginning - if (indices[0] == indices[1]) - return data[indices[0]]; - - // the rest is messier and requires a few - // switches, but simpler if we just sort - // our indices - { - TableIndices<2> sorted_indices (indices); - sorted_indices.sort (); - - if ((sorted_indices[0]==0) && (sorted_indices[1]==1)) - return data[3]; - else if ((sorted_indices[0]==0) && (sorted_indices[1]==2)) - return data[4]; - else if ((sorted_indices[0]==1) && (sorted_indices[1]==2)) - return data[5]; - else - Assert (false, ExcInternalError()); - } + default: + // to do the rest, sort our indices before comparing + { + TableIndices<2> sorted_indices (indices); + sorted_indices.sort (); + + for (unsigned int d=0, c=0; d sorted_indices (indices); + sorted_indices.sort (); + + for (unsigned int d=0, c=0; d(0,0); + if (i (i,i); + + for (unsigned int d=0, c=0; d(d,e); } } @@ -2606,63 +2521,23 @@ outer_product (const SymmetricTensor<2,dim,Number> &t1, /** * Return the symmetrized version of a full rank-2 tensor, i.e. * (t+transpose(t))/2, as a symmetric rank-2 tensor. This is the version for - * dim==1. + * general dimensions. * * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template -inline -SymmetricTensor<2,1,Number> -symmetrize (const Tensor<2,1,Number> &t) -{ - const Number array[1] - = { t[0][0] }; - return SymmetricTensor<2,1,Number>(array); -} - - - -/** - * Return the symmetrized version of a full rank-2 tensor, i.e. - * (t+transpose(t))/2, as a symmetric rank-2 tensor. This is the version for - * dim==2. - * - * @relates SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -template +template inline -SymmetricTensor<2,2,Number> -symmetrize (const Tensor<2,2,Number> &t) -{ - const Number array[3] - = { t[0][0], t[1][1], (t[0][1] + t[1][0])/2 }; - return SymmetricTensor<2,2,Number>(array); -} - - - -/** - * Return the symmetrized version of a full rank-2 tensor, i.e. - * (t+transpose(t))/2, as a symmetric rank-2 tensor. This is the version for - * dim==3. - * - * @relates SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -template -inline -SymmetricTensor<2,3,Number> -symmetrize (const Tensor<2,3,Number> &t) -{ - const Number array[6] - = { t[0][0], t[1][1], t[2][2], - (t[0][1] + t[1][0])/2, - (t[0][2] + t[2][0])/2, - (t[1][2] + t[2][1])/2 - }; - return SymmetricTensor<2,3,Number>(array); +SymmetricTensor<2,dim,Number> +symmetrize (const Tensor<2,dim,Number> &t) +{ + Number array[(dim*dim+dim)/2]; + for (unsigned int d=0; d(array); } -- 2.39.5