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);
56 local_Jac_data, displacement_size, displacement_size);
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++)
73 auto const& w = _ip_data[ip].integration_weight;
74 auto const& N = _ip_data[ip].N;
75 auto const& dNdx = _ip_data[ip].dNdx;
83 ShapeFunction::NPOINTS,
85 dNdx, N, x_coord, _is_axially_symmetric);
87 auto& eps = _ip_data[ip].eps;
88 eps.noalias() = B * u;
89 double const k = _process_data.residual_stiffness(t, x_position)[0];
90 double const ls = _process_data.crack_length_scale(t, x_position)[0];
91 double const d_ip = N.dot(d);
92 double const degradation =
93 _process_data.degradation_derivative->degradation(d_ip, k, ls);
94 _ip_data[ip].updateConstitutiveRelation(
95 t, x_position, dt, u, degradation,
96 _process_data.energy_split_model);
98 auto& sigma = _ip_data[ip].sigma;
99 auto const& D = _ip_data[ip].D;
101 typename ShapeMatricesType::template MatrixType<DisplacementDim,
103 N_u = ShapeMatricesType::template MatrixType<
104 DisplacementDim, displacement_size>::Zero(DisplacementDim,
107 for (
int i = 0; i < DisplacementDim; ++i)
109 N_u.template block<1, displacement_size / DisplacementDim>(
110 i, i * displacement_size / DisplacementDim)
114 auto const rho_sr = _process_data.solid_density(t, x_position)[0];
115 auto const& b = _process_data.specific_body_force;
117 local_rhs.noalias() -=
118 (B.transpose() * sigma - N_u.transpose() * rho_sr * b -
119 local_pressure * N_u.transpose() * dNdx * d) *
121 local_Jac.noalias() += B.transpose() * D * B * w;
128 Eigen::VectorXd
const& local_x,
129 std::vector<double>& local_b_data,
130 std::vector<double>& local_Jac_data)
132 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
134 local_x.template segment<displacement_size>(displacement_index);
137 local_Jac_data, phasefield_size, phasefield_size);
139 local_b_data, phasefield_size);
144 auto local_pressure = 0.0;
145 if (_process_data.static_pressurized_crack)
147 local_pressure = _process_data.unity_pressure;
149 else if (_process_data.propagating_pressurized_crack)
151 local_pressure = _process_data.pressure;
154 double const k = _process_data.residual_stiffness(t, x_position)[0];
155 double const ls = _process_data.crack_length_scale(t, x_position)[0];
156 double const gc = _process_data.crack_resistance(t, x_position)[0];
158 int const n_integration_points = _integration_method.getNumberOfPoints();
159 for (
int ip = 0; ip < n_integration_points; ip++)
161 auto const& w = _ip_data[ip].integration_weight;
162 auto const& N = _ip_data[ip].N;
163 auto const& dNdx = _ip_data[ip].dNdx;
165 double const d_ip = N.dot(d);
166 double const degradation =
167 _process_data.degradation_derivative->degradation(d_ip, k, ls);
168 double const degradation_df1 =
169 _process_data.degradation_derivative->degradation_df1(d_ip, ls);
170 double const degradation_df2 =
171 _process_data.degradation_derivative->degradation_df2(d_ip, ls);
177 ShapeFunction::NPOINTS,
179 dNdx, N, x_coord, _is_axially_symmetric);
181 auto& eps = _ip_data[ip].eps;
182 eps.noalias() = B * u;
183 _ip_data[ip].updateConstitutiveRelation(
184 t, x_position, dt, u, degradation,
185 _process_data.energy_split_model);
187 auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
189 auto& ip_data = _ip_data[ip];
190 ip_data.strain_energy_tensile = strain_energy_tensile;
192 typename ShapeMatricesType::template MatrixType<DisplacementDim,
194 N_u = ShapeMatricesType::template MatrixType<
195 DisplacementDim, displacement_size>::Zero(DisplacementDim,
198 for (
int i = 0; i < DisplacementDim; ++i)
200 N_u.template block<1, displacement_size / DisplacementDim>(
201 i, i * displacement_size / DisplacementDim)
205 local_Jac.noalias() +=
206 (N.transpose() * N * degradation_df2 * strain_energy_tensile) * w;
208 local_rhs.noalias() -=
209 (N.transpose() * degradation_df1 * strain_energy_tensile -
210 local_pressure * dNdx.transpose() * N_u * u) *
213 switch (_process_data.phasefield_model)
217 auto const local_Jac_AT1 =
218 (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
219 local_Jac.noalias() += local_Jac_AT1;
221 local_rhs.noalias() -=
222 gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
227 auto const local_Jac_AT2 =
229 (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
232 local_Jac.noalias() += local_Jac_AT2;
234 local_rhs.noalias() -=
235 local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
240 auto const local_Jac_COHESIVE =
241 (2.0 / std::numbers::pi * gc *
242 (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
246 local_Jac.noalias() += local_Jac_COHESIVE;
248 local_rhs.noalias() -= local_Jac_COHESIVE * d;
258 std::size_t mesh_item_id,
259 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
260 std::vector<GlobalVector*>
const& x,
double const ,
261 double& crack_volume)
263 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
264 indices_of_processes.reserve(dof_tables.size());
265 std::transform(dof_tables.begin(), dof_tables.end(),
266 std::back_inserter(indices_of_processes),
267 [&](
auto const dof_table)
268 { return NumLib::getIndices(mesh_item_id, *dof_table); });
271 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
273 auto const d = Eigen::Map<PhaseFieldVector const>(
274 &local_coupled_xs[phasefield_index], phasefield_size);
275 auto const u = Eigen::Map<DeformationVector const>(
276 &local_coupled_xs[displacement_index], displacement_size);
281 int const n_integration_points = _integration_method.getNumberOfPoints();
282 for (
int ip = 0; ip < n_integration_points; ip++)
284 auto const& w = _ip_data[ip].integration_weight;
285 auto const& N = _ip_data[ip].N;
286 auto const& dNdx = _ip_data[ip].dNdx;
288 typename ShapeMatricesType::template MatrixType<DisplacementDim,
290 N_u = ShapeMatricesType::template MatrixType<
291 DisplacementDim, displacement_size>::Zero(DisplacementDim,
294 for (
int i = 0; i < DisplacementDim; ++i)
296 N_u.template block<1, displacement_size / DisplacementDim>(
297 i, i * displacement_size / DisplacementDim)
301 crack_volume += (N_u * u).dot(dNdx * d) * w;
307 std::size_t mesh_item_id,
308 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
309 std::vector<GlobalVector*>
const& x,
double const t,
double& elastic_energy,
310 double& surface_energy,
double& pressure_work)
312 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
313 indices_of_processes.reserve(dof_tables.size());
314 std::transform(dof_tables.begin(), dof_tables.end(),
315 std::back_inserter(indices_of_processes),
316 [&](
auto const dof_table)
317 { return NumLib::getIndices(mesh_item_id, *dof_table); });
319 auto const local_coupled_xs =
321 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
323 auto const d = Eigen::Map<PhaseFieldVector const>(
324 &local_coupled_xs[phasefield_index], phasefield_size);
325 auto const u = Eigen::Map<DeformationVector const>(
326 &local_coupled_xs[displacement_index], displacement_size);
331 double element_elastic_energy = 0.0;
332 double element_surface_energy = 0.0;
333 double element_pressure_work = 0.0;
335 double const gc = _process_data.crack_resistance(t, x_position)[0];
336 double const ls = _process_data.crack_length_scale(t, x_position)[0];
337 int const n_integration_points = _integration_method.getNumberOfPoints();
338 for (
int ip = 0; ip < n_integration_points; ip++)
340 auto const& w = _ip_data[ip].integration_weight;
341 auto const& N = _ip_data[ip].N;
342 auto const& dNdx = _ip_data[ip].dNdx;
343 auto pressure_ip = _process_data.pressure;
345 typename ShapeMatricesType::template MatrixType<DisplacementDim,
347 N_u = ShapeMatricesType::template MatrixType<
348 DisplacementDim, displacement_size>::Zero(DisplacementDim,
351 for (
int i = 0; i < DisplacementDim; ++i)
353 N_u.template block<1, displacement_size / DisplacementDim>(
354 i, i * displacement_size / DisplacementDim)
358 element_elastic_energy += _ip_data[ip].elastic_energy * w;
360 double const d_ip = N.dot(d);
361 switch (_process_data.phasefield_model)
365 element_surface_energy +=
367 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
373 element_surface_energy += 0.5 * gc *
374 ((1 - d_ip) * (1 - d_ip) / ls +
375 (dNdx * d).dot((dNdx * d)) * ls) *
381 element_surface_energy +=
382 gc / std::numbers::pi *
383 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
389 if (_process_data.pressurized_crack)
391 element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
396 int const n_all_nodes = indices_of_processes[1].size();
397 int const n_regular_nodes = std::count_if(
398 begin(indices_of_processes[1]), end(indices_of_processes[1]),
400 if (n_all_nodes != n_regular_nodes)
402 element_elastic_energy *=
403 static_cast<double>(n_regular_nodes) / n_all_nodes;
404 element_surface_energy *=
405 static_cast<double>(n_regular_nodes) / n_all_nodes;
406 element_pressure_work *=
407 static_cast<double>(n_regular_nodes) / n_all_nodes;
410 elastic_energy += element_elastic_energy;
411 surface_energy += element_surface_energy;
412 pressure_work += element_pressure_work;