OGS
SmallDeformationLocalAssemblerFracture-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
23
24namespace ProcessLib
25{
26namespace LIE
27{
28namespace SmallDeformation
29{
30template <typename ShapeFunction, int DisplacementDim>
33 MeshLib::Element const& e,
34 std::size_t const n_variables,
35 std::size_t const /*local_matrix_size*/,
36 std::vector<unsigned> const& dofIndex_to_localIndex,
37 NumLib::GenericIntegrationMethod const& integration_method,
38 bool const is_axially_symmetric,
41 n_variables * ShapeFunction::NPOINTS * DisplacementDim,
42 dofIndex_to_localIndex),
43 _process_data(process_data),
44 _integration_method(integration_method),
45 _shape_matrices(
46 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
47 DisplacementDim>(e, is_axially_symmetric,
48 _integration_method)),
49 _element(e)
50{
51 assert(_element.getDimension() == DisplacementDim - 1);
52
53 unsigned const n_integration_points =
55
56 _ip_data.reserve(n_integration_points);
57 _secondary_data.N.resize(n_integration_points);
58
59 auto mat_id = (*_process_data.mesh_prop_materialIDs)[e.getID()];
60 auto frac_id = _process_data.map_materialID_to_fractureID[mat_id];
61 _fracture_property = &_process_data.fracture_properties[frac_id];
62 for (auto fid : process_data.vec_ele_connected_fractureIDs[e.getID()])
63 {
64 _fracID_to_local.insert({fid, _fracture_props.size()});
65 _fracture_props.push_back(&_process_data.fracture_properties[fid]);
66 }
67
69 ranges::views::transform(
70 [&](auto const jid)
71 { return &_process_data.junction_properties[jid]; }) |
72 ranges::to<std::vector>;
73
75 x_position.setElementID(_element.getID());
76 for (unsigned ip = 0; ip < n_integration_points; ip++)
77 {
78 x_position.setIntegrationPoint(ip);
79
80 _ip_data.emplace_back(*_process_data.fracture_model);
81 auto const& sm = _shape_matrices[ip];
82 auto& ip_data = _ip_data[ip];
83 ip_data.integration_weight =
85 sm.integralMeasure * sm.detJ;
86 ip_data.h_matrices.setZero(DisplacementDim,
87 ShapeFunction::NPOINTS * DisplacementDim);
88
89 computeHMatrix<DisplacementDim, ShapeFunction::NPOINTS,
91 HMatrixType>(sm.N, ip_data.h_matrices);
92
93 // Initialize current time step values
94 ip_data.w.setZero(DisplacementDim);
95 ip_data.sigma.setZero(DisplacementDim);
96
97 // Previous time step values are not initialized and are set later.
98 ip_data.sigma_prev.resize(DisplacementDim);
99 ip_data.w_prev.resize(DisplacementDim);
100
101 ip_data.C.resize(DisplacementDim, DisplacementDim);
102
103 ip_data.aperture0 = _fracture_property->aperture0(0, x_position)[0];
104 ip_data.aperture_prev = ip_data.aperture0;
105
106 _secondary_data.N[ip] = sm.N;
107 }
108}
109
110template <typename ShapeFunction, int DisplacementDim>
112 assembleWithJacobian(double const t, double const /*dt*/,
113 Eigen::VectorXd const& local_u,
114 Eigen::VectorXd& local_b, Eigen::MatrixXd& local_J)
115{
116 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
117 auto const n_fractures = _fracture_props.size();
118 auto const n_junctions = _junction_props.size();
119 auto const n_enrich_var = n_fractures + n_junctions;
120
121 //--------------------------------------------------------------------------------------
122 // prepare sub vectors, matrices for regular displacement (u) and
123 // displacement jumps (g)
124 //
125 // example with two fractures with one intersection:
126 // b = |b(g1)|
127 // |b(g2)|
128 //
129 // J = |J(g1,g1) J(g1,g2)|
130 // |J(g2,g1) J(g2,g2)|
131 //--------------------------------------------------------------------------------------
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 std::vector<BlockVectorType> vec_local_b_g;
139 for (unsigned i = 0; i < n_enrich_var; i++)
140 {
141 vec_local_b_g.push_back(
142 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
143 }
144 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
145 for (unsigned i = 0; i < n_enrich_var; i++)
146 {
147 for (unsigned j = 0; j < n_enrich_var; j++)
148 {
149 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
150 N_DOF_PER_VAR * i, N_DOF_PER_VAR * j);
151 vec_local_J_gg[i].push_back(sub_gg);
152 }
153 }
154
155 std::vector<Eigen::VectorXd> vec_nodal_g;
156 for (unsigned i = 0; i < n_enrich_var; i++)
157 {
158 auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
159 N_DOF_PER_VAR * i);
160 vec_nodal_g.push_back(sub);
161 }
162
163 //------------------------------------------------
164 // integration
165 //------------------------------------------------
166 // the index of a normal (normal to a fracture plane) component
167 // in a displacement vector
168 int const index_normal = DisplacementDim - 1;
169 auto const& R = _fracture_property->R;
170
171 unsigned const n_integration_points =
172 _integration_method.getNumberOfPoints();
173
175 x_position.setElementID(_element.getID());
176
177 for (unsigned ip = 0; ip < n_integration_points; ip++)
178 {
179 x_position.setIntegrationPoint(ip);
180
181 auto& ip_data = _ip_data[ip];
182 auto const& integration_weight = ip_data.integration_weight;
183 auto const& H = ip_data.h_matrices;
184 auto& mat = ip_data.fracture_material;
185 auto& sigma = ip_data.sigma;
186 auto const& sigma_prev = ip_data.sigma_prev;
187 auto& w = ip_data.w;
188 auto const& w_prev = ip_data.w_prev;
189 auto& C = ip_data.C;
190 auto& state = *ip_data.material_state_variables;
191 auto const& N = _secondary_data.N[ip];
192
193 auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
194 std::vector<double> const levelsets(duGlobalEnrichments(
195 _fracture_property->fracture_id, _fracture_props, _junction_props,
196 _fracID_to_local, ip_physical_coords));
197
198 // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
199 // * [u]_i)
200 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
201 nodal_gap.setZero();
202 for (unsigned i = 0; i < n_enrich_var; i++)
203 {
204 nodal_gap += levelsets[i] * vec_nodal_g[i];
205 }
206
207 // displacement jumps
208 w.noalias() = R * H * nodal_gap;
209
210 // total aperture
211 ip_data.aperture = ip_data.aperture0 + w[index_normal];
212
213 // local C, local stress
214 mat.computeConstitutiveRelation(
215 t, x_position, ip_data.aperture0,
216 Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
217 // Replace with
218 // initial
219 // stress values
220 w_prev, w, sigma_prev, sigma, C, state);
221
222 // r_[u] += H^T*Stress
223 for (unsigned i = 0; i < n_enrich_var; i++)
224 {
225 vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
226 R.transpose() * sigma *
227 integration_weight;
228 }
229
230 // J_[u][u] += H^T*C*H
231 for (unsigned i = 0; i < n_enrich_var; i++)
232 {
233 for (unsigned j = 0; j < n_enrich_var; j++)
234 {
235 // J_[u][u] += (levelset * B)^T * C * (levelset * B)
236 vec_local_J_gg[i][j].noalias() +=
237 (levelsets[i] * H.transpose() * R.transpose()) * C *
238 (levelsets[j] * R * H) * integration_weight;
239 }
240 }
241 }
242}
243
244template <typename ShapeFunction, int DisplacementDim>
247 Eigen::VectorXd const& local_u)
248{
249 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
250 auto const n_fractures = _fracture_props.size();
251 auto const n_junctions = _junction_props.size();
252 auto const n_enrich_var = n_fractures + n_junctions;
253
254 auto const& R = _fracture_property->R;
255
256 // the index of a normal (normal to a fracture plane) component
257 // in a displacement vector
258 int const index_normal = DisplacementDim - 1;
259
260 unsigned const n_integration_points =
261 _integration_method.getNumberOfPoints();
262
264 auto const e_id = _element.getID();
265 x_position.setElementID(e_id);
266
267 std::vector<Eigen::VectorXd> vec_nodal_g;
268 for (unsigned i = 0; i < n_enrich_var; i++)
269 {
270 auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
271 N_DOF_PER_VAR * i);
272 vec_nodal_g.push_back(sub);
273 }
274
275 for (unsigned ip = 0; ip < n_integration_points; ip++)
276 {
277 x_position.setIntegrationPoint(ip);
278
279 auto& ip_data = _ip_data[ip];
280 auto const& H = ip_data.h_matrices;
281 auto& mat = ip_data.fracture_material;
282 auto& sigma = ip_data.sigma;
283 auto const& sigma_prev = ip_data.sigma_prev;
284 auto& w = ip_data.w;
285 auto const& w_prev = ip_data.w_prev;
286 auto& C = ip_data.C;
287 auto& state = *ip_data.material_state_variables;
288 auto& b_m = ip_data.aperture;
289 auto const& N = _secondary_data.N[ip];
290
291 auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
292 std::vector<double> const levelsets(duGlobalEnrichments(
293 _fracture_property->fracture_id, _fracture_props, _junction_props,
294 _fracID_to_local, ip_physical_coords));
295
296 // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
297 // * [u]_i)
298 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
299 nodal_gap.setZero();
300 for (unsigned i = 0; i < n_enrich_var; i++)
301 {
302 nodal_gap += levelsets[i] * vec_nodal_g[i];
303 }
304
305 // displacement jumps in local coordinates
306 w.noalias() = R * H * nodal_gap;
307
308 // aperture
309 b_m = ip_data.aperture0 + w[index_normal];
310 if (b_m < 0.0)
311 {
312 OGS_FATAL(
313 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
314 "be non-negative.",
315 _element.getID(), ip, b_m);
316 }
317
318 // local C, local stress
319 mat.computeConstitutiveRelation(
320 t, x_position, ip_data.aperture0,
321 Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
322 // Replace with
323 // initial
324 // stress values
325 w_prev, w, sigma_prev, sigma, C, state);
326 }
327
328 double ele_b = 0;
329 ForceVectorType ele_sigma = ForceVectorType::Zero(DisplacementDim);
330 ForceVectorType ele_w = ForceVectorType::Zero(DisplacementDim);
331
332 for (unsigned ip = 0; ip < n_integration_points; ip++)
333 {
334 auto const& ip_data = _ip_data[ip];
335 ele_b += ip_data.aperture;
336 ele_w += ip_data.w;
337 ele_sigma += ip_data.sigma;
338 }
339 ele_b /= n_integration_points;
340 ele_w /= n_integration_points;
341 ele_sigma /= n_integration_points;
342 (*_process_data.mesh_prop_b)[_element.getID()] = ele_b;
343
344 Eigen::Map<GlobalDimVectorType>(
345 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
346 ele_sigma;
347
348 Eigen::Map<GlobalDimVectorType>(
349 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;
350}
351
352template <typename ShapeFunction, int DisplacementDim>
353std::vector<double> const&
356 const double /*t*/,
357 std::vector<GlobalVector*> const& /*x*/,
358 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
359 std::vector<double>& cache) const
360{
361 unsigned const n_integration_points = _ip_data.size();
362 cache.clear();
363 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
364 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
365 cache, DisplacementDim, n_integration_points);
366
367 for (unsigned ip = 0; ip < n_integration_points; ip++)
368 {
369 cache_matrix.col(ip).noalias() = _ip_data[ip].sigma;
370 }
371
372 return cache;
373}
374
375template <typename ShapeFunction, int DisplacementDim>
376std::vector<double> const&
379 const double /*t*/,
380 std::vector<GlobalVector*> const& /*x*/,
381 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
382 std::vector<double>& cache) const
383{
385 _ip_data, &IntegrationPointDataType::aperture, cache);
386}
387
388} // namespace SmallDeformation
389} // namespace LIE
390} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
double getWeight() const
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
std::size_t getID() const
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)
std::vector< double > const & getIntPtFractureAperture(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void computeSecondaryVariableConcreteWithVector(const double t, Eigen::VectorXd const &local_u) override
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
void assembleWithJacobian(double const t, double const dt, Eigen::VectorXd const &local_u, Eigen::VectorXd &local_b, Eigen::MatrixXd &local_J) override
SmallDeformationLocalAssemblerFracture(SmallDeformationLocalAssemblerFracture const &)=delete
std::vector< ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > _shape_matrices
std::vector< double > const & getIntPtFractureStress(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
std::vector< double > duGlobalEnrichments(std::size_t this_frac_id, 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)
Eigen::Vector3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
Definition Utils.h:24
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
ParameterLib::Parameter< double > const & aperture0
Initial aperture.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N