45 double const t,
double const dt, Eigen::VectorXd
const& local_x,
46 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
48 using DeformationMatrix =
49 typename ShapeMatricesType::template MatrixType<displacement_size,
51 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
53 local_x.template segment<displacement_size>(displacement_index);
55 auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
56 local_Jac_data, displacement_size, displacement_size);
58 auto local_rhs = MathLib::createZeroedVector<DeformationVector>(
59 local_b_data, displacement_size);
64 auto local_pressure = 0.0;
65 if (_process_data.pressurized_crack)
67 local_pressure = _process_data.unity_pressure;
70 int const n_integration_points = _integration_method.getNumberOfPoints();
71 for (
int ip = 0; ip < n_integration_points; ip++)
74 auto const& w = _ip_data[ip].integration_weight;
75 auto const& N = _ip_data[ip].N;
76 auto const& dNdx = _ip_data[ip].dNdx;
79 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
84 ShapeFunction::NPOINTS,
86 dNdx, N, x_coord, _is_axially_symmetric);
88 auto& eps = _ip_data[ip].eps;
89 eps.noalias() = B * u;
90 double const k = _process_data.residual_stiffness(t, x_position)[0];
91 double const ls = _process_data.crack_length_scale(t, x_position)[0];
92 double const d_ip = N.dot(d);
93 double const degradation =
94 _process_data.degradation_derivative->degradation(d_ip, k, ls);
95 _ip_data[ip].updateConstitutiveRelation(
96 t, x_position, dt, u, degradation,
97 _process_data.energy_split_model);
99 auto& sigma = _ip_data[ip].sigma;
100 auto const& D = _ip_data[ip].D;
102 typename ShapeMatricesType::template MatrixType<DisplacementDim,
104 N_u = ShapeMatricesType::template MatrixType<
105 DisplacementDim, displacement_size>::Zero(DisplacementDim,
108 for (
int i = 0; i < DisplacementDim; ++i)
110 N_u.template block<1, displacement_size / DisplacementDim>(
111 i, i * displacement_size / DisplacementDim)
115 auto const rho_sr = _process_data.solid_density(t, x_position)[0];
116 auto const& b = _process_data.specific_body_force;
118 local_rhs.noalias() -=
119 (B.transpose() * sigma - N_u.transpose() * rho_sr * b -
120 local_pressure * N_u.transpose() * dNdx * d) *
122 local_Jac.noalias() += B.transpose() * D * B * w;
129 Eigen::VectorXd
const& local_x,
130 std::vector<double>& local_b_data,
131 std::vector<double>& local_Jac_data)
133 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
135 local_x.template segment<displacement_size>(displacement_index);
137 auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
138 local_Jac_data, phasefield_size, phasefield_size);
139 auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
140 local_b_data, phasefield_size);
145 auto local_pressure = 0.0;
146 if (_process_data.static_pressurized_crack)
148 local_pressure = _process_data.unity_pressure;
150 else if (_process_data.propagating_pressurized_crack)
152 local_pressure = _process_data.pressure;
155 double const k = _process_data.residual_stiffness(t, x_position)[0];
156 double const ls = _process_data.crack_length_scale(t, x_position)[0];
157 double const gc = _process_data.crack_resistance(t, x_position)[0];
159 int const n_integration_points = _integration_method.getNumberOfPoints();
160 for (
int ip = 0; ip < n_integration_points; ip++)
163 auto const& w = _ip_data[ip].integration_weight;
164 auto const& N = _ip_data[ip].N;
165 auto const& dNdx = _ip_data[ip].dNdx;
167 double const d_ip = N.dot(d);
168 double const degradation =
169 _process_data.degradation_derivative->degradation(d_ip, k, ls);
170 double const degradation_df1 =
171 _process_data.degradation_derivative->degradation_df1(d_ip, ls);
172 double const degradation_df2 =
173 _process_data.degradation_derivative->degradation_df2(d_ip, ls);
175 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
179 ShapeFunction::NPOINTS,
181 dNdx, N, x_coord, _is_axially_symmetric);
183 auto& eps = _ip_data[ip].eps;
184 eps.noalias() = B * u;
185 _ip_data[ip].updateConstitutiveRelation(
186 t, x_position, dt, u, degradation,
187 _process_data.energy_split_model);
189 auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
191 auto& ip_data = _ip_data[ip];
192 ip_data.strain_energy_tensile = strain_energy_tensile;
194 typename ShapeMatricesType::template MatrixType<DisplacementDim,
196 N_u = ShapeMatricesType::template MatrixType<
197 DisplacementDim, displacement_size>::Zero(DisplacementDim,
200 for (
int i = 0; i < DisplacementDim; ++i)
202 N_u.template block<1, displacement_size / DisplacementDim>(
203 i, i * displacement_size / DisplacementDim)
207 local_Jac.noalias() +=
208 (N.transpose() * N * degradation_df2 * strain_energy_tensile) * w;
210 local_rhs.noalias() -=
211 (N.transpose() * degradation_df1 * strain_energy_tensile -
212 local_pressure * dNdx.transpose() * N_u * u) *
215 switch (_process_data.phasefield_model)
219 auto const local_Jac_AT1 =
220 (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
221 local_Jac.noalias() += local_Jac_AT1;
223 local_rhs.noalias() -=
224 gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
229 auto const local_Jac_AT2 =
231 (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
234 local_Jac.noalias() += local_Jac_AT2;
236 local_rhs.noalias() -=
237 local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
242 auto const local_Jac_COHESIVE =
243 (2.0 / std::numbers::pi * gc *
244 (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
248 local_Jac.noalias() += local_Jac_COHESIVE;
250 local_rhs.noalias() -= local_Jac_COHESIVE * d;
260 std::size_t mesh_item_id,
261 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
262 std::vector<GlobalVector*>
const& x,
double const ,
263 double& crack_volume)
265 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
266 indices_of_processes.reserve(dof_tables.size());
267 std::transform(dof_tables.begin(), dof_tables.end(),
268 std::back_inserter(indices_of_processes),
269 [&](
auto const dof_table)
270 { return NumLib::getIndices(mesh_item_id, *dof_table); });
273 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
275 auto const d = Eigen::Map<PhaseFieldVector const>(
276 &local_coupled_xs[phasefield_index], phasefield_size);
277 auto const u = Eigen::Map<DeformationVector const>(
278 &local_coupled_xs[displacement_index], displacement_size);
283 int const n_integration_points = _integration_method.getNumberOfPoints();
284 for (
int ip = 0; ip < n_integration_points; ip++)
287 auto const& w = _ip_data[ip].integration_weight;
288 auto const& N = _ip_data[ip].N;
289 auto const& dNdx = _ip_data[ip].dNdx;
291 typename ShapeMatricesType::template MatrixType<DisplacementDim,
293 N_u = ShapeMatricesType::template MatrixType<
294 DisplacementDim, displacement_size>::Zero(DisplacementDim,
297 for (
int i = 0; i < DisplacementDim; ++i)
299 N_u.template block<1, displacement_size / DisplacementDim>(
300 i, i * displacement_size / DisplacementDim)
304 crack_volume += (N_u * u).dot(dNdx * d) * w;
310 std::size_t mesh_item_id,
311 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
312 std::vector<GlobalVector*>
const& x,
double const t,
double& elastic_energy,
313 double& surface_energy,
double& pressure_work)
315 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
316 indices_of_processes.reserve(dof_tables.size());
317 std::transform(dof_tables.begin(), dof_tables.end(),
318 std::back_inserter(indices_of_processes),
319 [&](
auto const dof_table)
320 { return NumLib::getIndices(mesh_item_id, *dof_table); });
322 auto const local_coupled_xs =
324 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
326 auto const d = Eigen::Map<PhaseFieldVector const>(
327 &local_coupled_xs[phasefield_index], phasefield_size);
328 auto const u = Eigen::Map<DeformationVector const>(
329 &local_coupled_xs[displacement_index], displacement_size);
334 double element_elastic_energy = 0.0;
335 double element_surface_energy = 0.0;
336 double element_pressure_work = 0.0;
338 double const gc = _process_data.crack_resistance(t, x_position)[0];
339 double const ls = _process_data.crack_length_scale(t, x_position)[0];
340 int const n_integration_points = _integration_method.getNumberOfPoints();
341 for (
int ip = 0; ip < n_integration_points; ip++)
344 auto const& w = _ip_data[ip].integration_weight;
345 auto const& N = _ip_data[ip].N;
346 auto const& dNdx = _ip_data[ip].dNdx;
347 auto pressure_ip = _process_data.pressure;
349 typename ShapeMatricesType::template MatrixType<DisplacementDim,
351 N_u = ShapeMatricesType::template MatrixType<
352 DisplacementDim, displacement_size>::Zero(DisplacementDim,
355 for (
int i = 0; i < DisplacementDim; ++i)
357 N_u.template block<1, displacement_size / DisplacementDim>(
358 i, i * displacement_size / DisplacementDim)
362 element_elastic_energy += _ip_data[ip].elastic_energy * w;
364 double const d_ip = N.dot(d);
365 switch (_process_data.phasefield_model)
369 element_surface_energy +=
371 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
377 element_surface_energy += 0.5 * gc *
378 ((1 - d_ip) * (1 - d_ip) / ls +
379 (dNdx * d).dot((dNdx * d)) * ls) *
385 element_surface_energy +=
386 gc / std::numbers::pi *
387 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
393 if (_process_data.pressurized_crack)
395 element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
400 int const n_all_nodes = indices_of_processes[1].size();
401 int const n_regular_nodes = std::count_if(
402 begin(indices_of_processes[1]), end(indices_of_processes[1]),
404 if (n_all_nodes != n_regular_nodes)
406 element_elastic_energy *=
407 static_cast<double>(n_regular_nodes) / n_all_nodes;
408 element_surface_energy *=
409 static_cast<double>(n_regular_nodes) / n_all_nodes;
410 element_pressure_work *=
411 static_cast<double>(n_regular_nodes) / n_all_nodes;
414 elastic_energy += element_elastic_energy;
415 surface_energy += element_surface_energy;
416 pressure_work += element_pressure_work;