[PL] Implement BMatrices and Kelvin vec/mat types.

ProcessLib/Deformation/BMatrixPolicy.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_FEM_BMATRIXPOLICY_H_
+#define PROCESSLIB_FEM_BMATRIXPOLICY_H_
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+namespace ProcessLib
+{
+/// Kelvin vector dimensions for given displacement dimension.
+template <int DisplacementDim>
+struct KelvinVectorDimensions;
+
+template <>
+struct KelvinVectorDimensions<2>
+{
+    static int const value = 4;
+};
+
+template <>
+struct KelvinVectorDimensions<3>
+{
+    static int const value = 6;
+};
+
+//
+// Kelvin vector and matrix templates for given displacement dimension.
+//
+
+/// \TODO Maybe better to move the KelvinVector/MatrixType into
+/// MaterialLib/SolidModels/KelvinVector.h.
+
+/// Kelvin vector type for given displacement dimension.
+/// \note The Eigen vector is always a fixed size vector in contrast to the
+/// BMatrixPolicyType::StressVectorType.
+template <int DisplacementDim>
+using KelvinVectorType =
+    Eigen::Matrix<double, KelvinVectorDimensions<DisplacementDim>::value, 1,
+                  Eigen::ColMajor>;
+
+/// Kelvin matrix type for given displacement dimension.
+template <int DisplacementDim>
+using KelvinMatrixType =
+    Eigen::Matrix<double, KelvinVectorDimensions<DisplacementDim>::value,
+                  KelvinVectorDimensions<DisplacementDim>::value,
+                  Eigen::RowMajor>;
+
+/// An implementation of B-Matrix policy using same matrix and vector types
+/// (fixed size or dynamic) as in the ShapeMatrixPolicyType.
+template <typename ShapeFunction, unsigned DisplacementDim>
+class BMatrixPolicyType
+{
+private:
+    /// Reusing the ShapeMatrixPolicy vector type.
+    template <int N>
+    using VectorType =
+        typename ShapeMatrixPolicyType<ShapeFunction,
+                                       DisplacementDim>::template VectorType<N>;
+
+    /// Reusing the ShapeMatrixPolicy matrix type.
+    template <int N, int M>
+    using MatrixType = typename ShapeMatrixPolicyType<
+        ShapeFunction, DisplacementDim>::template MatrixType<N, M>;
+
+    // Dimensions of specific b-matrix for n-points and displacement dimension.
+    static int const _number_of_dof = ShapeFunction::NPOINTS * DisplacementDim;
+    static int const _kelvin_vector_size =
+        KelvinVectorDimensions<DisplacementDim>::value;
+
+public:
+    using StiffnessMatrixType = MatrixType<_number_of_dof, _number_of_dof>;
+
+    /// Rhs residual
+    using NodalForceVectorType = VectorType<_number_of_dof>;
+
+    /// Sigma
+    using StressVectorType = VectorType<_kelvin_vector_size>;
+
+    /// C
+    using ModulusMatrixType =
+        MatrixType<_kelvin_vector_size, _kelvin_vector_size>;
+
+    using BMatrixType = MatrixType<_kelvin_vector_size, _number_of_dof>;
+};
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_FEM_BMATRIXPOLICY_H_

ProcessLib/Deformation/LinearBMatrix.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_FEM_DEFORMATION_LINEARBMATRIX_H_
+#define PROCESSLIB_FEM_DEFORMATION_LINEARBMATRIX_H_
+
+#include <cmath>
+
+namespace ProcessLib
+{
+namespace LinearBMatrix
+{
+/// Fills a B-matrix based on given shape function dN/dx values.
+template <int DisplacementDim, int NPOINTS, typename DNDX_Type,
+          typename BMatrixType>
+void computeBMatrix(DNDX_Type const& dNdx, BMatrixType& b_matrix)
+{
+    static_assert(0 < DisplacementDim && DisplacementDim <= 3,
+                  "LinearBMatrix::computeBMatrix: DisplacementDim must be in "
+                  "range [1,3].");
+
+    b_matrix.setZero();
+
+    switch (DisplacementDim)
+    {
+        case 3:
+            for (int i = 0; i < NPOINTS; ++i)
+            {
+                b_matrix(2, 2 * NPOINTS + i) = dNdx(2, i);
+                b_matrix(4, NPOINTS + i) = dNdx(2, i) / std::sqrt(2);
+                b_matrix(4, 2 * NPOINTS + i) = dNdx(1, i) / std::sqrt(2);
+                b_matrix(5, i) = dNdx(2, i) / std::sqrt(2);
+                b_matrix(5, 2 * NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
+            }
+        // no break for fallthrough.
+        case 2:
+            for (int i = 0; i < NPOINTS; ++i)
+            {
+                b_matrix(1, NPOINTS + i) = dNdx(1, i);
+                b_matrix(3, i) = dNdx(1, i) / std::sqrt(2);
+                b_matrix(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
+                b_matrix(0, i) = dNdx(0, i);
+            }
+        default:
+            break;
+    }
+}
+
+}  // namespace LinearBMatrix
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_FEM_DEFORMATION_LINEARBMATRIX_H_