Loading [MathJax]/extensions/tex2jax.js
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
74 for (unsigned ip = 0; ip < n_integration_points; ip++)
75 {
76 _ip_data.emplace_back(*_process_data.fracture_model);
77 auto const& sm = _shape_matrices[ip];
78 auto& ip_data = _ip_data[ip];
79
80 ParameterLib::SpatialPosition const x_position{
81 std::nullopt, _element.getID(),
84 _element, sm.N))};
85
86 ip_data.integration_weight =
88 sm.integralMeasure * sm.detJ;
89 ip_data.h_matrices.setZero(DisplacementDim,
90 ShapeFunction::NPOINTS * DisplacementDim);
91
92 computeHMatrix<DisplacementDim, ShapeFunction::NPOINTS,
94 HMatrixType>(sm.N, ip_data.h_matrices);
95
96 // Initialize current time step values
97 ip_data.w.setZero(DisplacementDim);
98 ip_data.sigma.setZero(DisplacementDim);
99
100 // Previous time step values are not initialized and are set later.
101 ip_data.sigma_prev.resize(DisplacementDim);
102 ip_data.w_prev.resize(DisplacementDim);
103
104 ip_data.C.resize(DisplacementDim, DisplacementDim);
105
106 ip_data.aperture0 = _fracture_property->aperture0(0, x_position)[0];
107 ip_data.aperture_prev = ip_data.aperture0;
108
109 _secondary_data.N[ip] = sm.N;
110 }
111}
112
113template <typename ShapeFunction, int DisplacementDim>
115 assembleWithJacobian(double const t, double const /*dt*/,
116 Eigen::VectorXd const& local_u,
117 Eigen::VectorXd& local_b, Eigen::MatrixXd& local_J)
118{
119 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
120 auto const n_fractures = _fracture_props.size();
121 auto const n_junctions = _junction_props.size();
122 auto const n_enrich_var = n_fractures + n_junctions;
123
124 //--------------------------------------------------------------------------------------
125 // prepare sub vectors, matrices for regular displacement (u) and
126 // displacement jumps (g)
127 //
128 // example with two fractures with one intersection:
129 // b = |b(g1)|
130 // |b(g2)|
131 //
132 // J = |J(g1,g1) J(g1,g2)|
133 // |J(g2,g1) J(g2,g2)|
134 //--------------------------------------------------------------------------------------
135
136 using BlockVectorType =
137 typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
138 using BlockMatrixType =
139 Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
140
141 std::vector<BlockVectorType> vec_local_b_g;
142 for (unsigned i = 0; i < n_enrich_var; i++)
143 {
144 vec_local_b_g.push_back(
145 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
146 }
147 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
148 for (unsigned i = 0; i < n_enrich_var; i++)
149 {
150 for (unsigned j = 0; j < n_enrich_var; j++)
151 {
152 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
153 N_DOF_PER_VAR * i, N_DOF_PER_VAR * j);
154 vec_local_J_gg[i].push_back(sub_gg);
155 }
156 }
157
158 std::vector<Eigen::VectorXd> vec_nodal_g;
159 for (unsigned i = 0; i < n_enrich_var; i++)
160 {
161 auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
162 N_DOF_PER_VAR * i);
163 vec_nodal_g.push_back(sub);
164 }
165
166 //------------------------------------------------
167 // integration
168 //------------------------------------------------
169 // the index of a normal (normal to a fracture plane) component
170 // in a displacement vector
171 int const index_normal = DisplacementDim - 1;
172 auto const& R = _fracture_property->R;
173
174 unsigned const n_integration_points =
175 _integration_method.getNumberOfPoints();
176
178 x_position.setElementID(_element.getID());
179
180 for (unsigned ip = 0; ip < n_integration_points; ip++)
181 {
182 auto& ip_data = _ip_data[ip];
183 auto const& integration_weight = ip_data.integration_weight;
184 auto const& H = ip_data.h_matrices;
185 auto& mat = ip_data.fracture_material;
186 auto& sigma = ip_data.sigma;
187 auto const& sigma_prev = ip_data.sigma_prev;
188 auto& w = ip_data.w;
189 auto const& w_prev = ip_data.w_prev;
190 auto& C = ip_data.C;
191 auto& state = *ip_data.material_state_variables;
192 auto const& N = _secondary_data.N[ip];
193
194 auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
195 std::vector<double> const levelsets(duGlobalEnrichments(
196 _fracture_property->fracture_id, _fracture_props, _junction_props,
197 _fracID_to_local, ip_physical_coords));
198
199 // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
200 // * [u]_i)
201 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
202 nodal_gap.setZero();
203 for (unsigned i = 0; i < n_enrich_var; i++)
204 {
205 nodal_gap += levelsets[i] * vec_nodal_g[i];
206 }
207
208 // displacement jumps
209 w.noalias() = R * H * nodal_gap;
210
211 // total aperture
212 ip_data.aperture = ip_data.aperture0 + w[index_normal];
213
214 // local C, local stress
215 mat.computeConstitutiveRelation(
216 t, x_position, ip_data.aperture0,
217 Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
218 // Replace with
219 // initial
220 // stress values
221 w_prev, w, sigma_prev, sigma, C, state);
222
223 // r_[u] += H^T*Stress
224 for (unsigned i = 0; i < n_enrich_var; i++)
225 {
226 vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
227 R.transpose() * sigma *
228 integration_weight;
229 }
230
231 // J_[u][u] += H^T*C*H
232 for (unsigned i = 0; i < n_enrich_var; i++)
233 {
234 for (unsigned j = 0; j < n_enrich_var; j++)
235 {
236 // J_[u][u] += (levelset * B)^T * C * (levelset * B)
237 vec_local_J_gg[i][j].noalias() +=
238 (levelsets[i] * H.transpose() * R.transpose()) * C *
239 (levelsets[j] * R * H) * integration_weight;
240 }
241 }
242 }
243}
244
245template <typename ShapeFunction, int DisplacementDim>
248 Eigen::VectorXd const& local_u)
249{
250 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
251 auto const n_fractures = _fracture_props.size();
252 auto const n_junctions = _junction_props.size();
253 auto const n_enrich_var = n_fractures + n_junctions;
254
255 auto const& R = _fracture_property->R;
256
257 // the index of a normal (normal to a fracture plane) component
258 // in a displacement vector
259 int const index_normal = DisplacementDim - 1;
260
261 unsigned const n_integration_points =
262 _integration_method.getNumberOfPoints();
263
265 auto const e_id = _element.getID();
266 x_position.setElementID(e_id);
267
268 std::vector<Eigen::VectorXd> vec_nodal_g;
269 for (unsigned i = 0; i < n_enrich_var; i++)
270 {
271 auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
272 N_DOF_PER_VAR * i);
273 vec_nodal_g.push_back(sub);
274 }
275
276 for (unsigned ip = 0; ip < n_integration_points; ip++)
277 {
278 auto& ip_data = _ip_data[ip];
279 auto const& H = ip_data.h_matrices;
280 auto& mat = ip_data.fracture_material;
281 auto& sigma = ip_data.sigma;
282 auto const& sigma_prev = ip_data.sigma_prev;
283 auto& w = ip_data.w;
284 auto const& w_prev = ip_data.w_prev;
285 auto& C = ip_data.C;
286 auto& state = *ip_data.material_state_variables;
287 auto& b_m = ip_data.aperture;
288 auto const& N = _secondary_data.N[ip];
289
290 auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
291 std::vector<double> const levelsets(duGlobalEnrichments(
292 _fracture_property->fracture_id, _fracture_props, _junction_props,
293 _fracID_to_local, ip_physical_coords));
294
295 // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
296 // * [u]_i)
297 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
298 nodal_gap.setZero();
299 for (unsigned i = 0; i < n_enrich_var; i++)
300 {
301 nodal_gap += levelsets[i] * vec_nodal_g[i];
302 }
303
304 // displacement jumps in local coordinates
305 w.noalias() = R * H * nodal_gap;
306
307 // aperture
308 b_m = ip_data.aperture0 + w[index_normal];
309 if (b_m < 0.0)
310 {
311 OGS_FATAL(
312 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
313 "be non-negative.",
314 _element.getID(), ip, b_m);
315 }
316
317 // local C, local stress
318 mat.computeConstitutiveRelation(
319 t, x_position, ip_data.aperture0,
320 Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
321 // Replace with
322 // initial
323 // stress values
324 w_prev, w, sigma_prev, sigma, C, state);
325 }
326
327 double ele_b = 0;
328 ForceVectorType ele_sigma = ForceVectorType::Zero(DisplacementDim);
329 ForceVectorType ele_w = ForceVectorType::Zero(DisplacementDim);
330
331 for (unsigned ip = 0; ip < n_integration_points; ip++)
332 {
333 auto const& ip_data = _ip_data[ip];
334 ele_b += ip_data.aperture;
335 ele_w += ip_data.w;
336 ele_sigma += ip_data.sigma;
337 }
338 ele_b /= n_integration_points;
339 ele_w /= n_integration_points;
340 ele_sigma /= n_integration_points;
341 (*_process_data.mesh_prop_b)[_element.getID()] = ele_b;
342
343 Eigen::Map<GlobalDimVectorType>(
344 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
345 ele_sigma;
346
347 Eigen::Map<GlobalDimVectorType>(
348 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;
349}
350
351template <typename ShapeFunction, int DisplacementDim>
352std::vector<double> const&
355 const double /*t*/,
356 std::vector<GlobalVector*> const& /*x*/,
357 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
358 std::vector<double>& cache) const
359{
360 unsigned const n_integration_points = _ip_data.size();
361 cache.clear();
362 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
363 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
364 cache, DisplacementDim, n_integration_points);
365
366 for (unsigned ip = 0; ip < n_integration_points; ip++)
367 {
368 cache_matrix.col(ip).noalias() = _ip_data[ip].sigma;
369 }
370
371 return cache;
372}
373
374template <typename ShapeFunction, int DisplacementDim>
375std::vector<double> const&
378 const double /*t*/,
379 std::vector<GlobalVector*> const& /*x*/,
380 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
381 std::vector<double>& cache) const
382{
384 _ip_data, &IntegrationPointDataType::aperture, cache);
385}
386
387} // namespace SmallDeformation
388} // namespace LIE
389} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
constexpr 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:91
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
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::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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