34 std::size_t
const n_variables,
36 std::vector<unsigned>
const& dofIndex_to_localIndex,
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),
47 DisplacementDim>(e, is_axially_symmetric,
48 _integration_method)),
53 unsigned const n_integration_points =
56 _ip_data.reserve(n_integration_points);
60 auto frac_id =
_process_data.map_materialID_to_fractureID[mat_id];
69 ranges::views::transform(
72 ranges::to<std::vector>;
76 for (
unsigned ip = 0; ip < n_integration_points; ip++)
78 x_position.setIntegrationPoint(ip);
83 ip_data.integration_weight =
85 sm.integralMeasure * sm.detJ;
86 ip_data.h_matrices.setZero(DisplacementDim,
87 ShapeFunction::NPOINTS * DisplacementDim);
94 ip_data.w.setZero(DisplacementDim);
95 ip_data.sigma.setZero(DisplacementDim);
98 ip_data.sigma_prev.resize(DisplacementDim);
99 ip_data.w_prev.resize(DisplacementDim);
101 ip_data.C.resize(DisplacementDim, DisplacementDim);
104 ip_data.aperture_prev = ip_data.aperture0;
113 Eigen::VectorXd
const& local_u,
114 Eigen::VectorXd& local_b, Eigen::MatrixXd& local_J)
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;
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>;
138 std::vector<BlockVectorType> vec_local_b_g;
139 for (
unsigned i = 0; i < n_enrich_var; i++)
141 vec_local_b_g.push_back(
142 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
144 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
145 for (
unsigned i = 0; i < n_enrich_var; i++)
147 for (
unsigned j = 0; j < n_enrich_var; j++)
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);
155 std::vector<Eigen::VectorXd> vec_nodal_g;
156 for (
unsigned i = 0; i < n_enrich_var; i++)
158 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
160 vec_nodal_g.push_back(sub);
168 int const index_normal = DisplacementDim - 1;
169 auto const& R = _fracture_property->R;
171 unsigned const n_integration_points =
172 _integration_method.getNumberOfPoints();
177 for (
unsigned ip = 0; ip < n_integration_points; ip++)
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;
188 auto const& w_prev = ip_data.w_prev;
190 auto& state = *ip_data.material_state_variables;
191 auto const& N = _secondary_data.N[ip];
195 _fracture_property->fracture_id, _fracture_props, _junction_props,
196 _fracID_to_local, ip_physical_coords));
200 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
202 for (
unsigned i = 0; i < n_enrich_var; i++)
204 nodal_gap += levelsets[i] * vec_nodal_g[i];
208 w.noalias() = R * H * nodal_gap;
211 ip_data.aperture = ip_data.aperture0 + w[index_normal];
214 mat.computeConstitutiveRelation(
215 t, x_position, ip_data.aperture0,
216 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
220 w_prev, w, sigma_prev, sigma, C, state);
223 for (
unsigned i = 0; i < n_enrich_var; i++)
225 vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
226 R.transpose() * sigma *
231 for (
unsigned i = 0; i < n_enrich_var; i++)
233 for (
unsigned j = 0; j < n_enrich_var; j++)
236 vec_local_J_gg[i][j].noalias() +=
237 (levelsets[i] * H.transpose() * R.transpose()) * C *
238 (levelsets[j] * R * H) * integration_weight;
247 Eigen::VectorXd
const& local_u)
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;
254 auto const& R = _fracture_property->R;
258 int const index_normal = DisplacementDim - 1;
260 unsigned const n_integration_points =
261 _integration_method.getNumberOfPoints();
264 auto const e_id = _element.getID();
267 std::vector<Eigen::VectorXd> vec_nodal_g;
268 for (
unsigned i = 0; i < n_enrich_var; i++)
270 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
272 vec_nodal_g.push_back(sub);
275 for (
unsigned ip = 0; ip < n_integration_points; ip++)
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;
285 auto const& w_prev = ip_data.w_prev;
287 auto& state = *ip_data.material_state_variables;
288 auto& b_m = ip_data.aperture;
289 auto const& N = _secondary_data.N[ip];
293 _fracture_property->fracture_id, _fracture_props, _junction_props,
294 _fracID_to_local, ip_physical_coords));
298 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
300 for (
unsigned i = 0; i < n_enrich_var; i++)
302 nodal_gap += levelsets[i] * vec_nodal_g[i];
306 w.noalias() = R * H * nodal_gap;
309 b_m = ip_data.aperture0 + w[index_normal];
313 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
315 _element.getID(), ip, b_m);
319 mat.computeConstitutiveRelation(
320 t, x_position, ip_data.aperture0,
321 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
325 w_prev, w, sigma_prev, sigma, C, state);
332 for (
unsigned ip = 0; ip < n_integration_points; ip++)
334 auto const& ip_data = _ip_data[ip];
335 ele_b += ip_data.aperture;
337 ele_sigma += ip_data.sigma;
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;
344 Eigen::Map<GlobalDimVectorType>(
345 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
348 Eigen::Map<GlobalDimVectorType>(
349 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;