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>;
74 for (
unsigned ip = 0; ip < n_integration_points; ip++)
86 ip_data.integration_weight =
88 sm.integralMeasure * sm.detJ;
89 ip_data.h_matrices.setZero(DisplacementDim,
90 ShapeFunction::NPOINTS * DisplacementDim);
97 ip_data.w.setZero(DisplacementDim);
98 ip_data.sigma.setZero(DisplacementDim);
101 ip_data.sigma_prev.resize(DisplacementDim);
102 ip_data.w_prev.resize(DisplacementDim);
104 ip_data.C.resize(DisplacementDim, DisplacementDim);
107 ip_data.aperture_prev = ip_data.aperture0;
116 Eigen::VectorXd
const& local_u,
117 Eigen::VectorXd& local_b, Eigen::MatrixXd& local_J)
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;
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>;
141 std::vector<BlockVectorType> vec_local_b_g;
142 for (
unsigned i = 0; i < n_enrich_var; i++)
144 vec_local_b_g.push_back(
145 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
147 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
148 for (
unsigned i = 0; i < n_enrich_var; i++)
150 for (
unsigned j = 0; j < n_enrich_var; j++)
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);
158 std::vector<Eigen::VectorXd> vec_nodal_g;
159 for (
unsigned i = 0; i < n_enrich_var; i++)
161 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
163 vec_nodal_g.push_back(sub);
171 int const index_normal = DisplacementDim - 1;
172 auto const& R = _fracture_property->R;
174 unsigned const n_integration_points =
175 _integration_method.getNumberOfPoints();
180 for (
unsigned ip = 0; ip < n_integration_points; ip++)
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;
189 auto const& w_prev = ip_data.w_prev;
191 auto& state = *ip_data.material_state_variables;
192 auto const& N = _secondary_data.N[ip];
196 _fracture_property->fracture_id, _fracture_props, _junction_props,
197 _fracID_to_local, ip_physical_coords));
201 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
203 for (
unsigned i = 0; i < n_enrich_var; i++)
205 nodal_gap += levelsets[i] * vec_nodal_g[i];
209 w.noalias() = R * H * nodal_gap;
212 ip_data.aperture = ip_data.aperture0 + w[index_normal];
215 mat.computeConstitutiveRelation(
216 t, x_position, ip_data.aperture0,
217 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
221 w_prev, w, sigma_prev, sigma, C, state);
224 for (
unsigned i = 0; i < n_enrich_var; i++)
226 vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
227 R.transpose() * sigma *
232 for (
unsigned i = 0; i < n_enrich_var; i++)
234 for (
unsigned j = 0; j < n_enrich_var; j++)
237 vec_local_J_gg[i][j].noalias() +=
238 (levelsets[i] * H.transpose() * R.transpose()) * C *
239 (levelsets[j] * R * H) * integration_weight;
248 Eigen::VectorXd
const& local_u)
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;
255 auto const& R = _fracture_property->R;
259 int const index_normal = DisplacementDim - 1;
261 unsigned const n_integration_points =
262 _integration_method.getNumberOfPoints();
265 auto const e_id = _element.getID();
268 std::vector<Eigen::VectorXd> vec_nodal_g;
269 for (
unsigned i = 0; i < n_enrich_var; i++)
271 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
273 vec_nodal_g.push_back(sub);
276 for (
unsigned ip = 0; ip < n_integration_points; ip++)
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;
284 auto const& w_prev = ip_data.w_prev;
286 auto& state = *ip_data.material_state_variables;
287 auto& b_m = ip_data.aperture;
288 auto const& N = _secondary_data.N[ip];
292 _fracture_property->fracture_id, _fracture_props, _junction_props,
293 _fracID_to_local, ip_physical_coords));
297 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
299 for (
unsigned i = 0; i < n_enrich_var; i++)
301 nodal_gap += levelsets[i] * vec_nodal_g[i];
305 w.noalias() = R * H * nodal_gap;
308 b_m = ip_data.aperture0 + w[index_normal];
312 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
314 _element.getID(), ip, b_m);
318 mat.computeConstitutiveRelation(
319 t, x_position, ip_data.aperture0,
320 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
324 w_prev, w, sigma_prev, sigma, C, state);
331 for (
unsigned ip = 0; ip < n_integration_points; ip++)
333 auto const& ip_data = _ip_data[ip];
334 ele_b += ip_data.aperture;
336 ele_sigma += ip_data.sigma;
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;
343 Eigen::Map<GlobalDimVectorType>(
344 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
347 Eigen::Map<GlobalDimVectorType>(
348 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;