OGS
SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Core>
14#include <range/v3/range/conversion.hpp>
15#include <range/v3/view/transform.hpp>
16#include <vector>
17
22#include "MathLib/Point3d.h"
23#include "MeshLib/Node.h"
31#include "SecondaryData.h"
34
35namespace ProcessLib
36{
37namespace LIE
38{
39namespace SmallDeformation
40{
41template <typename ShapeFunction,
42
43 int DisplacementDim>
44SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
45
46 DisplacementDim>::
47 SmallDeformationLocalAssemblerMatrixNearFracture(
48 MeshLib::Element const& e,
49 std::size_t const n_variables,
50 std::size_t const /*local_matrix_size*/,
51 std::vector<unsigned> const& dofIndex_to_localIndex,
52 NumLib::GenericIntegrationMethod const& integration_method,
53 bool const is_axially_symmetric,
56 n_variables * ShapeFunction::NPOINTS * DisplacementDim,
57 dofIndex_to_localIndex),
58 _process_data(process_data),
59 _integration_method(integration_method),
60 _element(e),
61 _is_axially_symmetric(is_axially_symmetric)
62{
63 std::vector<ShapeMatrices, Eigen::aligned_allocator<
65 shape_matrices =
67 DisplacementDim>(e, is_axially_symmetric,
69
70 unsigned const n_integration_points =
72
73 _ip_data.reserve(n_integration_points);
74 _secondary_data.N.resize(n_integration_points);
75
77 _process_data.solid_materials, _process_data.material_ids, e.getID());
78
79 for (unsigned ip = 0; ip < n_integration_points; ip++)
80 {
81 _ip_data.emplace_back(solid_material);
82 auto& ip_data = _ip_data[ip];
83 auto const& sm = shape_matrices[ip];
84 ip_data.N = sm.N;
85 ip_data.dNdx = sm.dNdx;
86 ip_data.integration_weight =
88 sm.integralMeasure * sm.detJ;
89
90 // Initialize current time step values
91 static const int kelvin_vector_size =
93 ip_data._sigma.setZero(kelvin_vector_size);
94 ip_data._eps.setZero(kelvin_vector_size);
95
96 // Previous time step values are not initialized and are set later.
97 ip_data._sigma_prev.resize(kelvin_vector_size);
98 ip_data._eps_prev.resize(kelvin_vector_size);
99
100 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
101
102 _secondary_data.N[ip] = sm.N;
103 }
104
105 for (auto fid : process_data._vec_ele_connected_fractureIDs[e.getID()])
106 {
107 _fracID_to_local.insert({fid, _fracture_props.size()});
108 _fracture_props.push_back(&_process_data.fracture_properties[fid]);
109 }
110
112 ranges::views::transform(
113 [&](auto const jid)
114 { return &_process_data.junction_properties[jid]; }) |
115 ranges::to<std::vector>;
116}
117
118template <typename ShapeFunction, int DisplacementDim>
120 ShapeFunction,
121 DisplacementDim>::assembleWithJacobian(double const t, double const dt,
122 Eigen::VectorXd const& local_u,
123 Eigen::VectorXd& local_b,
124 Eigen::MatrixXd& local_J)
125{
126 assert(_element.getDimension() == DisplacementDim);
127
128 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
129 auto const n_fractures = _fracture_props.size();
130 auto const n_junctions = _junction_props.size();
131 auto const n_enrich_var = n_fractures + n_junctions;
132
133 using BlockVectorType =
134 typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
135 using BlockMatrixType =
136 Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
137
138 //--------------------------------------------------------------------------------------
139 // prepare sub vectors, matrices for regular displacement (u) and
140 // displacement jumps (g)
141 //
142 // example with two fractures with one intersection:
143 // |b(u)|
144 // b = |b(g1)|
145 // |b(g2)|
146 // |b(j1)|
147 //
148 // |J(u,u) J(u,g1) J(u,g2) J(u,j1) |
149 // J = |J(g1,u) J(g1,g1) J(g1,g2) J(g1,j1)|
150 // |J(g2,u) J(g2,g1) J(g2,g2) J(g2,j1)|
151 // |J(j1,u) J(j1,g1) J(j1,g2) J(j1,j1)|
152 //--------------------------------------------------------------------------------------
153 auto local_b_u = local_b.segment<N_DOF_PER_VAR>(0);
154 std::vector<BlockVectorType> vec_local_b_g;
155 for (unsigned i = 0; i < n_enrich_var; i++)
156 {
157 vec_local_b_g.push_back(
158 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * (i + 1)));
159 }
160
161 auto local_J_uu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(0, 0);
162 std::vector<BlockMatrixType> vec_local_J_ug;
163 std::vector<BlockMatrixType> vec_local_J_gu;
164 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
165 for (unsigned i = 0; i < n_enrich_var; i++)
166 {
167 auto sub_ug = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
168 0, N_DOF_PER_VAR * (i + 1));
169 vec_local_J_ug.push_back(sub_ug);
170
171 auto sub_gu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
172 N_DOF_PER_VAR * (i + 1), 0);
173 vec_local_J_gu.push_back(sub_gu);
174
175 for (unsigned j = 0; j < n_enrich_var; j++)
176 {
177 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
178 N_DOF_PER_VAR * (i + 1), N_DOF_PER_VAR * (j + 1));
179 vec_local_J_gg[i].push_back(sub_gg);
180 }
181 }
182
183 auto const nodal_u = local_u.segment<N_DOF_PER_VAR>(0);
184 std::vector<BlockVectorType> vec_nodal_g;
185 for (unsigned i = 0; i < n_enrich_var; i++)
186 {
187 auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
188 N_DOF_PER_VAR * (i + 1));
189 vec_nodal_g.push_back(sub);
190 }
191
192 //------------------------------------------------
193 // integration
194 //------------------------------------------------
195 unsigned const n_integration_points =
196 _integration_method.getNumberOfPoints();
197
198 MPL::VariableArray variables;
199 MPL::VariableArray variables_prev;
201 x_position.setElementID(_element.getID());
202
203 for (unsigned ip = 0; ip < n_integration_points; ip++)
204 {
205 x_position.setIntegrationPoint(ip);
206
207 auto& ip_data = _ip_data[ip];
208 auto const& w = _ip_data[ip].integration_weight;
209
210 auto const& N = ip_data.N;
211 auto const& dNdx = ip_data.dNdx;
212
213 // levelset functions
214 Eigen::Vector3d const ip_physical_coords(
215 computePhysicalCoordinates(_element, N).data());
216 std::vector<double> const levelsets(
217 uGlobalEnrichments(_fracture_props, _junction_props,
218 _fracID_to_local, ip_physical_coords));
219
220 // u = u^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x) *
221 // [u]_i)
222 NodalDisplacementVectorType nodal_total_u = nodal_u;
223 for (unsigned i = 0; i < n_enrich_var; i++)
224 {
225 nodal_total_u += levelsets[i] * vec_nodal_g[i];
226 }
227
228 auto const x_coord =
229 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
230 _element, N);
231 auto const B =
232 LinearBMatrix::computeBMatrix<DisplacementDim,
233 ShapeFunction::NPOINTS,
235 dNdx, N, x_coord, _is_axially_symmetric);
236
237 // strain, stress
238 auto const& eps_prev = ip_data._eps_prev;
239 auto const& sigma_prev = ip_data._sigma_prev;
240
241 auto& eps = ip_data._eps;
242 auto& sigma = ip_data._sigma;
243 auto& state = ip_data._material_state_variables;
244
245 eps.noalias() = B * nodal_total_u;
246
247 variables.mechanical_strain
249 eps);
250
251 variables_prev.stress
253 sigma_prev);
254 variables_prev.mechanical_strain
256 eps_prev);
257 variables_prev.temperature = _process_data._reference_temperature;
258
259 auto&& solution = _ip_data[ip]._solid_material.integrateStress(
260 variables_prev, variables, t, x_position, dt, *state);
261
262 if (!solution)
263 {
264 OGS_FATAL("Computation of local constitutive relation failed.");
265 }
266
268 std::tie(sigma, state, C) = std::move(*solution);
269
270 // r_u = B^T * Sigma = B^T * C * B * (u+phi*[u])
271 // r_[u] = (phi*B)^T * Sigma = (phi*B)^T * C * B * (u+phi*[u])
272 local_b_u.noalias() -= B.transpose() * sigma * w;
273 for (unsigned i = 0; i < n_enrich_var; i++)
274 {
275 vec_local_b_g[i].noalias() -=
276 levelsets[i] * B.transpose() * sigma * w;
277 }
278
279 // J_uu += B^T * C * B
280 local_J_uu.noalias() += B.transpose() * C * B * w;
281
282 for (unsigned i = 0; i < n_enrich_var; i++)
283 {
284 // J_u[u] += B^T * C * (levelset * B)
285 vec_local_J_ug[i].noalias() +=
286 B.transpose() * C * (levelsets[i] * B) * w;
287
288 // J_[u]u += (levelset * B)^T * C * B
289 vec_local_J_gu[i].noalias() +=
290 (levelsets[i] * B.transpose()) * C * B * w;
291
292 for (unsigned j = 0; j < n_enrich_var; j++)
293 {
294 // J_[u][u] += (levelset * B)^T * C * (levelset * B)
295 vec_local_J_gg[i][j].noalias() +=
296 (levelsets[i] * B.transpose()) * C * (levelsets[j] * B) * w;
297 }
298 }
299 }
300}
301
302template <typename ShapeFunction,
303
304 int DisplacementDim>
306
307 DisplacementDim>::
308 computeSecondaryVariableConcreteWithVector(
309 double const /*t*/, Eigen::VectorXd const& /*local_x*/)
310{
311 // Compute average value per element
312 const int n = DisplacementDim == 2 ? 4 : 6;
313 Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
314 Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
315
316 unsigned const n_integration_points =
317 _integration_method.getNumberOfPoints();
318 for (unsigned ip = 0; ip < n_integration_points; ip++)
319 {
320 auto const& ip_data = _ip_data[ip];
321
322 ele_stress += ip_data._sigma;
323 ele_strain += ip_data._eps;
324 }
325 ele_stress /= n_integration_points;
326 ele_strain /= n_integration_points;
327
328 (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
329 (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
330 (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
331 (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
332 if (DisplacementDim == 3)
333 {
334 (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
335 (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
336 }
337
338 (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
339 (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
340 (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
341 (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
342 if (DisplacementDim == 3)
343 {
344 (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
345 (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
346 }
347}
348
349} // namespace SmallDeformation
350} // namespace LIE
351} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
Definition of the Node class.
Definition of the Point3d class.
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
Definition: VariableType.h:192
double getWeight() const
Definition: WeightedPoint.h:80
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
std::vector< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim > > > _ip_data
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:57
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Eigen::Vector3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
Definition: Utils.h:24
std::vector< double > uGlobalEnrichments(std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
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.
Definition: LinearBMatrix.h:42
Coordinates mapping matrices at particular location.
Definition: ShapeMatrices.h:45
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:27