OGS
LinearBMatrix.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <cmath>
14#include <optional>
15
18
19namespace ProcessLib
20{
21namespace LinearBMatrix
22{
23namespace detail
24{
25template <int NPOINTS, typename DNDX_Type, typename BMatrixType>
26void fillBMatrix2DCartesianPart(DNDX_Type const& dNdx, BMatrixType& B)
27{
28 for (int i = 0; i < NPOINTS; ++i)
29 {
30 B(1, NPOINTS + i) = dNdx(1, i);
31 B(3, i) = dNdx(1, i) / std::sqrt(2);
32 B(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
33 B(0, i) = dNdx(0, i);
34 }
35}
36} // namespace detail
37
39template <int DisplacementDim,
40 int NPOINTS,
41 typename BMatrixType,
42 typename N_Type,
43 typename DNDX_Type>
44BMatrixType computeBMatrix(DNDX_Type const& dNdx,
45 N_Type const& N,
46 const double radius,
47 const bool is_axially_symmetric)
48{
49 static_assert(0 < DisplacementDim && DisplacementDim <= 3,
50 "LinearBMatrix::computeBMatrix: DisplacementDim must be in "
51 "range [1,3].");
52
53 BMatrixType B = BMatrixType::Zero(
55 NPOINTS * DisplacementDim);
56
57 switch (DisplacementDim)
58 {
59 case 3:
60 for (int i = 0; i < NPOINTS; ++i)
61 {
62 B(2, 2 * NPOINTS + i) = dNdx(2, i);
63 B(4, NPOINTS + i) = dNdx(2, i) / std::sqrt(2);
64 B(4, 2 * NPOINTS + i) = dNdx(1, i) / std::sqrt(2);
65 B(5, i) = dNdx(2, i) / std::sqrt(2);
66 B(5, 2 * NPOINTS + i) = dNdx(0, i) / std::sqrt(2);
67 }
69 break;
70 case 2:
72 if (is_axially_symmetric)
73 {
74 for (int i = 0; i < NPOINTS; ++i)
75 {
76 B(2, i) = N[i] / radius;
77 }
78 }
79 break;
80 default:
81 break;
82 }
83
84 return B;
85}
86
87template <int DisplacementDim, int NPOINTS, typename ShapeFunction,
88 typename BBarMatrixType, typename ShapeMatricesType, typename IpData>
90 std::vector<IpData, Eigen::aligned_allocator<IpData>> const& ip_data,
91 MeshLib::Element const& element,
92 NumLib::GenericIntegrationMethod const& integration_method,
93 const bool is_axially_symmetric)
94{
95 unsigned const n_integration_points =
96 integration_method.getNumberOfPoints();
97
98 BBarMatrixType B_bar = BBarMatrixType::Zero(3, ShapeFunction::NPOINTS);
99
100 // Compute the element volume
101 double volume = 0.0;
102 for (unsigned ip = 0; ip < n_integration_points; ip++)
103 {
104 auto const& w = ip_data[ip].integration_weight;
105 volume += w;
106 }
107
108 for (int i = 0; i < NPOINTS; i++)
109 {
110 B_bar.col(i).noalias() +=
111 NumLib::averageGradShapeFunction<DisplacementDim, ShapeFunction,
112 ShapeMatricesType, IpData>(
113 i, element, integration_method, ip_data, is_axially_symmetric);
114 }
115
116 return B_bar / volume;
117}
118
119namespace detail
120{
121template <int DisplacementDim, int NPOINTS, typename BBarMatrixType,
122 typename BMatrixType>
123void applyBbar(BBarMatrixType const& B_bar, BMatrixType& B,
124 const bool is_axially_symmetric)
125{
126 if constexpr (DisplacementDim == 3)
127 {
128 for (int i = 0; i < NPOINTS; ++i)
129 {
130 auto const B_bar_i = B_bar.col(i);
131
132 // The following loop is based on the following facts:
133 // B1 (dN/dx) is the 1st element of the 1st column of B,
134 // B2 (dN/dy) is the 2nd element of the 2nd column of B,
135 // B3 (dN/dz) is the 3rd element of the 3rd column of B.
136 for (int k = 0; k < 3; k++)
137 {
138 // k-th column of B matrix at node i
139 auto B_i_k = B.col(k * NPOINTS + i);
140
141 B_i_k.template segment<3>(0) -=
142 Eigen::Vector3d::Constant((B_i_k[k] - B_bar_i[k]) / 3.0);
143 }
144 }
145
146 return;
147 }
148
149 // 2D or axisymmetry
150 for (int i = 0; i < NPOINTS; ++i)
151 {
152 auto B_i_0 = B.col(i);
153 auto B_i_1 = B.col(NPOINTS + i);
154
155 auto const B_bar_i = B_bar.col(i);
156
157 if (is_axially_symmetric)
158 {
159 double const b0_dil_pertubation =
160 (B_i_0[0] - B_bar_i[0] + B_i_0[2] - B_bar_i[2]);
161 B_i_0.template segment<3>(0) -=
162 Eigen::Vector3d::Constant((b0_dil_pertubation) / 3.);
163 B_i_1.template segment<3>(0) -=
164 Eigen::Vector3d::Constant((B_i_1[1] - B_bar_i[1]) / 3.);
165 continue;
166 }
167
168 // Plane strain
169 B_i_0.template segment<2>(0) -=
170 Eigen::Vector2d::Constant((B_i_0[0] - B_bar_i[0]) / 2.);
171 B_i_1.template segment<2>(0) -=
172 Eigen::Vector2d::Constant((B_i_1[1] - B_bar_i[1]) / 2.);
173 }
174}
175} // namespace detail
176
178template <int DisplacementDim, int NPOINTS, typename BBarMatrixType,
179 typename BMatrixType, typename N_Type, typename DNDX_Type>
181 DNDX_Type const& dNdx,
182 N_Type const& N,
183 std::optional<BBarMatrixType> const& B_dil_bar,
184 const double radius,
185 const bool is_axially_symmetric)
186{
187 auto B = computeBMatrix<DisplacementDim, NPOINTS, BMatrixType, N_Type,
188 DNDX_Type>(dNdx, N, radius, is_axially_symmetric);
189
190 if (!B_dil_bar)
191 {
192 return B;
193 }
194
196 *B_dil_bar, B, is_axially_symmetric);
197
198 return B;
199}
200
201} // namespace LinearBMatrix
202} // namespace ProcessLib
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Vector3d averageGradShapeFunction(int const local_node_id, MeshLib::Element const &element, NumLib::GenericIntegrationMethod const &integration_method, std::vector< IpData, Eigen::aligned_allocator< IpData > > const &ip_data, const bool is_axially_symmetric)
void fillBMatrix2DCartesianPart(DNDX_Type const &dNdx, BMatrixType &B)
void applyBbar(BBarMatrixType const &B_bar, BMatrixType &B, const bool is_axially_symmetric)
BBarMatrixType computeDilatationalBbar(std::vector< IpData, Eigen::aligned_allocator< IpData > > const &ip_data, MeshLib::Element const &element, NumLib::GenericIntegrationMethod const &integration_method, const bool is_axially_symmetric)
BMatrixType computeBMatrixPossiblyWithBbar(DNDX_Type const &dNdx, N_Type const &N, std::optional< BBarMatrixType > const &B_dil_bar, const double radius, const bool is_axially_symmetric)
Fills a B matrix, or a B bar matrix if required.
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.