/**
* Return the result of the push forward transformation on a rank-4
- * contravariant tensor, i.e. (in index notation)
+ * contravariant tensor, i.e. (in index notation):
* @f[
* \left[ \chi\left(\bullet\right)^{\sharp} \right]_{ijkl}
* \dealcoloneq F_{iI} F_{jJ}
/**
* Return the result of the push forward transformation on a rank-4
- * contravariant symmetric tensor, i.e. (in index notation)
+ * contravariant symmetric tensor, i.e. (in index notation):
* @f[
* \left[ \chi\left(\bullet\right)^{\sharp} \right]_{ijkl}
* \dealcoloneq F_{iI} F_{jJ}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant tensor, i.e. (in index notation)
+ * contravariant tensor, i.e. (in index notation):
* @f[
* \left[ \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL}
* \dealcoloneq F^{-1}_{Ii} F^{-1}_{Jj}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant symmetric tensor, i.e. (in index notation)
+ * contravariant symmetric tensor, i.e. (in index notation):
* @f[
* \left[ \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL}
* \dealcoloneq F^{-1}_{Ii} F^{-1}_{Jj}
/**
* Return the result of the push forward transformation on a rank-4
- * covariant tensor, i.e. (in index notation)
+ * covariant tensor, i.e. (in index notation):
* @f[
* \left[ \chi\left(\bullet\right)^{\flat} \right]_{ijkl}
* \dealcoloneq F^{-T}_{iI} F^{-T}_{jJ}
/**
* Return the result of the push forward transformation on a rank-4
- * covariant symmetric tensor, i.e. (in index notation)
+ * covariant symmetric tensor, i.e. (in index notation):
* @f[
* \left[ \chi\left(\bullet\right)^{\flat} \right]_{ijkl}
* \dealcoloneq F^{-T}_{iI} F^{-T}_{jJ}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant tensor, i.e. (in index notation)
+ * contravariant tensor, i.e. (in index notation):
* @f[
* \left[ \chi^{-1}\left(\bullet\right)^{\flat} \right]_{IJKL}
* \dealcoloneq F^{T}_{Ii} F^{T}_{Jj}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant symmetric tensor, i.e. (in index notation)
+ * contravariant symmetric tensor, i.e. (in index notation):
* @f[
* \left[ \chi^{-1}\left(\bullet\right)^{\flat} \right]_{IJKL}
* \dealcoloneq F^{T}_{Ii} F^{T}_{Jj}
/**
* Return the result of the push forward transformation on a rank-4
- * contravariant tensor, i.e. (in index notation)
+ * contravariant tensor, i.e. (in index notation):
* @f[
* \textrm{det} \mathbf{F}^{-1} \; \left[
* \chi\left(\bullet\right)^{\sharp} \right]_{ijkl}
/**
* Return the result of the push forward transformation on a rank-4
- * contravariant symmetric tensor, i.e. (in index notation)
+ * contravariant symmetric tensor, i.e. (in index notation):
* @f[
* \textrm{det} \mathbf{F}^{-1} \; \left[
* \chi\left(\bullet\right)^{\sharp} \right]_{ijkl}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant tensor, i.e. (in index notation)
+ * contravariant tensor, i.e. (in index notation):
* @f[
* \textrm{det} \mathbf{F} \; \left[
* \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL}
/**
* Return the result of the pull back transformation on a rank-4
- * contravariant symmetric tensor, i.e. (in index notation)
+ * contravariant symmetric tensor, i.e. (in index notation):
* @f[
* \textrm{det} \mathbf{F} \; \left[
* \chi^{-1}\left(\bullet\right)^{\sharp} \right]_{IJKL}
const Tensor<2, dim, Number> &F);
//@}
- } // namespace Transformations
-} // namespace Physics
-
-
-#ifndef DOXYGEN
-
-// ------------------------- inline functions ------------------------
+ /**
+ * @name Basis transformations
+ */
+ //@{
-namespace internal
-{
- namespace Physics
- {
+ /**
+ * Return a vector with a changed basis, i.e.
+ * @f[
+ * \mathbf{V}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{V}
+ * @f]
+ *
+ * @param[in] V The vector to be transformed $\mathbf{V}$
+ * @param[in] B The transformation matrix $\mathbf{B}$
+ * @return $\mathbf{V}^{\prime}$
+ */
template <int dim, typename Number>
- inline Tensor<1, dim, Number>
- transformation_contraction(const Tensor<1, dim, Number> &V,
- const Tensor<2, dim, Number> &F)
- {
- return contract<1, 0>(F, V);
- }
-
-
+ Tensor<1, dim, Number>
+ basis_transformation(const Tensor<1, dim, Number> &V,
+ const Tensor<2, dim, Number> &B);
+ /**
+ * Return a rank-2 tensor with a changed basis, i.e.
+ * @f[
+ * \mathbf{T}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{T} \cdot
+ * \mathbf{B}^{T}
+ * @f]
+ *
+ * @param[in] T The tensor to be transformed $\mathbf{T}$
+ * @param[in] B The transformation matrix $\mathbf{B}$
+ * @return $\mathbf{T}^{\prime}$
+ */
template <int dim, typename Number>
- inline Tensor<2, dim, Number>
- transformation_contraction(const Tensor<2, dim, Number> &T,
- const Tensor<2, dim, Number> &F)
- {
- return contract<1, 0>(F, contract<1, 1>(T, F));
- }
-
-
+ Tensor<2, dim, Number>
+ basis_transformation(const Tensor<2, dim, Number> &T,
+ const Tensor<2, dim, Number> &B);
+ /**
+ * Return a symmetric rank-2 tensor with a changed basis, i.e.
+ * @f[
+ * \mathbf{T}^{\prime} \dealcoloneq \mathbf{B} \cdot \mathbf{T} \cdot
+ * \mathbf{B}^{T}
+ * @f]
+ *
+ * @param[in] T The tensor to be transformed $\mathbf{T}$
+ * @param[in] B The transformation matrix $\mathbf{B}$
+ * @return $\mathbf{T}^{\prime}$
+ */
template <int dim, typename Number>
- inline dealii::SymmetricTensor<2, dim, Number>
- transformation_contraction(const dealii::SymmetricTensor<2, dim, Number> &T,
- const Tensor<2, dim, Number> & F)
- {
- Tensor<2, dim, Number> tmp_1;
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int J = 0; J < dim; ++J)
- // Loop over I but complex.h defines a macro I, so use I_ instead
- for (unsigned int I_ = 0; I_ < dim; ++I_)
- tmp_1[i][J] += F[i][I_] * T[I_][J];
-
- dealii::SymmetricTensor<2, dim, Number> out;
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int j = i; j < dim; ++j)
- for (unsigned int J = 0; J < dim; ++J)
- out[i][j] += F[j][J] * tmp_1[i][J];
-
- return out;
- }
-
-
+ SymmetricTensor<2, dim, Number>
+ basis_transformation(const SymmetricTensor<2, dim, Number> &T,
+ const Tensor<2, dim, Number> & B);
+ /**
+ * Return a rank-4 tensor with a changed basis, i.e. (in index notation):
+ * @f[
+ * H_{ijkl}^{\prime} \dealcoloneq B_{iI} B_{jJ} H_{IJKL} B_{kK} B_{lL}
+ * @f]
+ *
+ * @param[in] H The tensor to be transformed $\mathbf{T}$
+ * @param[in] B The transformation matrix $\mathbf{B}$
+ * @return $\mathbf{H}^{\prime}$
+ */
template <int dim, typename Number>
- inline Tensor<4, dim, Number>
- transformation_contraction(const Tensor<4, dim, Number> &H,
- const Tensor<2, dim, Number> &F)
- {
- // This contraction order and indexing might look a bit dubious, so a
- // quick explanation as to what's going on is probably in order:
- //
- // When the contract() function operates on the inner indices, the
- // result has the inner index and outer index transposed, i.e.
- // contract<2,1>(H,F) implies
- // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
- // rather than T_{IJkL} (the desired result).
- // So, in effect, contraction of the 3rd (inner) index with F as the
- // second argument results in its transposition with respect to its
- // adjacent neighbor. This is due to the position of the argument F,
- // leading to the free index being on the right hand side of the result.
- // However, given that we can do two transformations from the LHS of H
- // and two from the right we can undo the otherwise erroneous
- // swapping of the outer indices upon application of the second
- // sets of contractions.
- //
- // Note: Its significantly quicker (in 3d) to push forward
- // each index individually
- return contract<1, 1>(
- F, contract<1, 1>(F, contract<2, 1>(contract<2, 1>(H, F), F)));
- }
+ Tensor<4, dim, Number>
+ basis_transformation(const Tensor<4, dim, Number> &H,
+ const Tensor<2, dim, Number> &B);
+ /**
+ * Return a symmetric rank-4 tensor with a changed basis, i.e. (in index
+ * notation):
+ * @f[
+ * H_{ijkl}^{\prime} \dealcoloneq B_{iI} B_{jJ} H_{IJKL} B_{kK} B_{lL}
+ * @f]
+ *
+ * @param[in] H The tensor to be transformed $\mathbf{T}$
+ * @param[in] B The transformation matrix $\mathbf{B}$
+ * @return $\mathbf{H}^{\prime}$
+ */
+ template <int dim, typename Number>
+ SymmetricTensor<4, dim, Number>
+ basis_transformation(const SymmetricTensor<4, dim, Number> &H,
+ const Tensor<2, dim, Number> & B);
+ //@}
- template <int dim, typename Number>
- inline dealii::SymmetricTensor<4, dim, Number>
- transformation_contraction(const dealii::SymmetricTensor<4, dim, Number> &H,
- const Tensor<2, dim, Number> & F)
- {
- // The first and last transformation operations respectively
- // break and recover the symmetry properties of the tensors.
- // We also want to perform a minimal number of operations here
- // and avoid some complications related to the transposition of
- // tensor indices when contracting inner indices using the contract()
- // function. (For an explanation of the contraction operations,
- // please see the note in the equivalent function for standard
- // Tensors.) So what we'll do here is manually perform the first
- // and last contractions that break/recover the tensor symmetries
- // on the inner indices, and use the contract() function only on
- // the outer indices.
- //
- // Note: Its significantly quicker (in 3d) to push forward
- // each index individually
-
- // Push forward (inner) index 1
- Tensor<4, dim, Number> tmp;
- // Loop over I but complex.h defines a macro I, so use I_ instead
- for (unsigned int I_ = 0; I_ < dim; ++I_)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int K = 0; K < dim; ++K)
- for (unsigned int L = 0; L < dim; ++L)
- for (unsigned int J = 0; J < dim; ++J)
- tmp[I_][j][K][L] += F[j][J] * H[I_][J][K][L];
+ } // namespace Transformations
+} // namespace Physics
- // Push forward (outer) indices 0 and 3
- tmp = contract<1, 0>(F, contract<3, 1>(tmp, F));
- // Push forward (inner) index 2
- dealii::SymmetricTensor<4, dim, Number> out;
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int j = i; j < dim; ++j)
- for (unsigned int k = 0; k < dim; ++k)
- for (unsigned int l = k; l < dim; ++l)
- for (unsigned int K = 0; K < dim; ++K)
- out[i][j][k][l] += F[k][K] * tmp[i][j][K][l];
- return out;
- }
- } // namespace Physics
-} // namespace internal
+#ifndef DOXYGEN
const Tensor<1, dim, Number> &V,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(V, F);
+ return Physics::Transformations::basis_transformation(V, F);
}
const Tensor<2, dim, Number> &T,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(T, F);
+ return Physics::Transformations::basis_transformation(T, F);
}
const SymmetricTensor<2, dim, Number> &T,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(T, F);
+ return Physics::Transformations::basis_transformation(T, F);
}
const Tensor<4, dim, Number> &H,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(H, F);
+ return Physics::Transformations::basis_transformation(H, F);
}
const SymmetricTensor<4, dim, Number> &H,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(H, F);
+ return Physics::Transformations::basis_transformation(H, F);
}
const Tensor<1, dim, Number> &v,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(v, invert(F));
+ return Physics::Transformations::basis_transformation(v, invert(F));
}
const Tensor<2, dim, Number> &t,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(t, invert(F));
+ return Physics::Transformations::basis_transformation(t, invert(F));
}
const SymmetricTensor<2, dim, Number> &t,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(t, invert(F));
+ return Physics::Transformations::basis_transformation(t, invert(F));
}
const Tensor<4, dim, Number> &h,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(h, invert(F));
+ return Physics::Transformations::basis_transformation(h, invert(F));
}
const SymmetricTensor<4, dim, Number> &h,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(h, invert(F));
+ return Physics::Transformations::basis_transformation(h, invert(F));
}
const Tensor<1, dim, Number> &V,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(V, transpose(invert(F)));
+ return Physics::Transformations::basis_transformation(V,
+ transpose(invert(F)));
}
const Tensor<2, dim, Number> &T,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(T, transpose(invert(F)));
+ return Physics::Transformations::basis_transformation(T,
+ transpose(invert(F)));
}
const SymmetricTensor<2, dim, Number> &T,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(T, transpose(invert(F)));
+ return Physics::Transformations::basis_transformation(T,
+ transpose(invert(F)));
}
const Tensor<4, dim, Number> &H,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(H, transpose(invert(F)));
+ return Physics::Transformations::basis_transformation(H,
+ transpose(invert(F)));
}
const SymmetricTensor<4, dim, Number> &H,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(H, transpose(invert(F)));
+ return Physics::Transformations::basis_transformation(H,
+ transpose(invert(F)));
}
Physics::Transformations::Covariant::pull_back(const Tensor<1, dim, Number> &v,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(v, transpose(F));
+ return Physics::Transformations::basis_transformation(v, transpose(F));
}
Physics::Transformations::Covariant::pull_back(const Tensor<2, dim, Number> &t,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(t, transpose(F));
+ return Physics::Transformations::basis_transformation(t, transpose(F));
}
const SymmetricTensor<2, dim, Number> &t,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(t, transpose(F));
+ return Physics::Transformations::basis_transformation(t, transpose(F));
}
Physics::Transformations::Covariant::pull_back(const Tensor<4, dim, Number> &h,
const Tensor<2, dim, Number> &F)
{
- return internal::Physics::transformation_contraction(h, transpose(F));
+ return Physics::Transformations::basis_transformation(h, transpose(F));
}
const SymmetricTensor<4, dim, Number> &h,
const Tensor<2, dim, Number> & F)
{
- return internal::Physics::transformation_contraction(h, transpose(F));
+ return Physics::Transformations::basis_transformation(h, transpose(F));
}
return cofactor(F) * N;
}
+
+template <int dim, typename Number>
+inline Tensor<1, dim, Number>
+Physics::Transformations::basis_transformation(const Tensor<1, dim, Number> &V,
+ const Tensor<2, dim, Number> &B)
+{
+ return contract<1, 0>(B, V);
+}
+
+
+
+template <int dim, typename Number>
+inline Tensor<2, dim, Number>
+Physics::Transformations::basis_transformation(const Tensor<2, dim, Number> &T,
+ const Tensor<2, dim, Number> &B)
+{
+ return contract<1, 0>(B, contract<1, 1>(T, B));
+}
+
+
+
+template <int dim, typename Number>
+inline SymmetricTensor<2, dim, Number>
+Physics::Transformations::basis_transformation(
+ const SymmetricTensor<2, dim, Number> &T,
+ const Tensor<2, dim, Number> & B)
+{
+ Tensor<2, dim, Number> tmp_1;
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int J = 0; J < dim; ++J)
+ // Loop over I but complex.h defines a macro I, so use I_ instead
+ for (unsigned int I_ = 0; I_ < dim; ++I_)
+ tmp_1[i][J] += B[i][I_] * T[I_][J];
+
+ SymmetricTensor<2, dim, Number> out;
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = i; j < dim; ++j)
+ for (unsigned int J = 0; J < dim; ++J)
+ out[i][j] += B[j][J] * tmp_1[i][J];
+
+ return out;
+}
+
+
+
+template <int dim, typename Number>
+inline Tensor<4, dim, Number>
+Physics::Transformations::basis_transformation(const Tensor<4, dim, Number> &H,
+ const Tensor<2, dim, Number> &B)
+{
+ // This contraction order and indexing might look a bit dubious, so a
+ // quick explanation as to what's going on is probably in order:
+ //
+ // When the contract() function operates on the inner indices, the
+ // result has the inner index and outer index transposed, i.e.
+ // contract<2,1>(H,F) implies
+ // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
+ // rather than T_{IJkL} (the desired result).
+ // So, in effect, contraction of the 3rd (inner) index with F as the
+ // second argument results in its transposition with respect to its
+ // adjacent neighbor. This is due to the position of the argument F,
+ // leading to the free index being on the right hand side of the result.
+ // However, given that we can do two transformations from the LHS of H
+ // and two from the right we can undo the otherwise erroneous
+ // swapping of the outer indices upon application of the second
+ // sets of contractions.
+ //
+ // Note: Its significantly quicker (in 3d) to push forward
+ // each index individually
+ return contract<1, 1>(
+ B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B)));
+}
+
+
+
+template <int dim, typename Number>
+inline SymmetricTensor<4, dim, Number>
+Physics::Transformations::basis_transformation(
+ const SymmetricTensor<4, dim, Number> &H,
+ const Tensor<2, dim, Number> & B)
+{
+ // The first and last transformation operations respectively
+ // break and recover the symmetry properties of the tensors.
+ // We also want to perform a minimal number of operations here
+ // and avoid some complications related to the transposition of
+ // tensor indices when contracting inner indices using the contract()
+ // function. (For an explanation of the contraction operations,
+ // please see the note in the equivalent function for standard
+ // Tensors.) So what we'll do here is manually perform the first
+ // and last contractions that break/recover the tensor symmetries
+ // on the inner indices, and use the contract() function only on
+ // the outer indices.
+ //
+ // Note: Its significantly quicker (in 3d) to push forward
+ // each index individually
+
+ // Push forward (inner) index 1
+ Tensor<4, dim, Number> tmp;
+ // Loop over I but complex.h defines a macro I, so use I_ instead
+ for (unsigned int I_ = 0; I_ < dim; ++I_)
+ for (unsigned int j = 0; j < dim; ++j)
+ for (unsigned int K = 0; K < dim; ++K)
+ for (unsigned int L = 0; L < dim; ++L)
+ for (unsigned int J = 0; J < dim; ++J)
+ tmp[I_][j][K][L] += B[j][J] * H[I_][J][K][L];
+
+ // Push forward (outer) indices 0 and 3
+ tmp = contract<1, 0>(B, contract<3, 1>(tmp, B));
+
+ // Push forward (inner) index 2
+ SymmetricTensor<4, dim, Number> out;
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = i; j < dim; ++j)
+ for (unsigned int k = 0; k < dim; ++k)
+ for (unsigned int l = k; l < dim; ++l)
+ for (unsigned int K = 0; K < dim; ++K)
+ out[i][j][k][l] += B[k][K] * tmp[i][j][K][l];
+
+ return out;
+}
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff