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->degradationDf1(d_ip, k, ls);
170 double const degradation_df2 =
171 _process_data.degradation_derivative->degradationDf2(d_ip, k, 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 calculateCrackLocalJacobianAndResidual<
214 decltype(dNdx),
decltype(N),
decltype(w),
decltype(d),
215 decltype(local_Jac),
decltype(local_rhs)>(
216 dNdx, N, w, d, local_Jac, local_rhs, gc, ls,
217 _process_data.phasefield_model);
224 std::size_t mesh_item_id,
225 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
226 std::vector<GlobalVector*>
const& x,
double const ,
227 double& crack_volume)
229 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
230 indices_of_processes.reserve(dof_tables.size());
231 std::transform(dof_tables.begin(), dof_tables.end(),
232 std::back_inserter(indices_of_processes),
233 [&](
auto const dof_table)
234 { return NumLib::getIndices(mesh_item_id, *dof_table); });
237 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
239 auto const d = Eigen::Map<PhaseFieldVector const>(
240 &local_coupled_xs[phasefield_index], phasefield_size);
241 auto const u = Eigen::Map<DeformationVector const>(
242 &local_coupled_xs[displacement_index], displacement_size);
247 int const n_integration_points = _integration_method.getNumberOfPoints();
248 for (
int ip = 0; ip < n_integration_points; ip++)
250 auto const& w = _ip_data[ip].integration_weight;
251 auto const& N = _ip_data[ip].N;
252 auto const& dNdx = _ip_data[ip].dNdx;
254 typename ShapeMatricesType::template MatrixType<DisplacementDim,
256 N_u = ShapeMatricesType::template MatrixType<
257 DisplacementDim, displacement_size>::Zero(DisplacementDim,
260 for (
int i = 0; i < DisplacementDim; ++i)
262 N_u.template block<1, displacement_size / DisplacementDim>(
263 i, i * displacement_size / DisplacementDim)
267 crack_volume += (N_u * u).dot(dNdx * d) * w;
273 std::size_t mesh_item_id,
274 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
275 std::vector<GlobalVector*>
const& x,
double const t,
double& elastic_energy,
276 double& surface_energy,
double& pressure_work)
278 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
279 indices_of_processes.reserve(dof_tables.size());
280 std::transform(dof_tables.begin(), dof_tables.end(),
281 std::back_inserter(indices_of_processes),
282 [&](
auto const dof_table)
283 { return NumLib::getIndices(mesh_item_id, *dof_table); });
285 auto const local_coupled_xs =
287 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
289 auto const d = Eigen::Map<PhaseFieldVector const>(
290 &local_coupled_xs[phasefield_index], phasefield_size);
291 auto const u = Eigen::Map<DeformationVector const>(
292 &local_coupled_xs[displacement_index], displacement_size);
297 double element_elastic_energy = 0.0;
298 double element_surface_energy = 0.0;
299 double element_pressure_work = 0.0;
301 double const gc = _process_data.crack_resistance(t, x_position)[0];
302 double const ls = _process_data.crack_length_scale(t, x_position)[0];
303 int const n_integration_points = _integration_method.getNumberOfPoints();
304 for (
int ip = 0; ip < n_integration_points; ip++)
306 auto const& w = _ip_data[ip].integration_weight;
307 auto const& N = _ip_data[ip].N;
308 auto const& dNdx = _ip_data[ip].dNdx;
309 auto pressure_ip = _process_data.pressure;
311 typename ShapeMatricesType::template MatrixType<DisplacementDim,
313 N_u = ShapeMatricesType::template MatrixType<
314 DisplacementDim, displacement_size>::Zero(DisplacementDim,
317 for (
int i = 0; i < DisplacementDim; ++i)
319 N_u.template block<1, displacement_size / DisplacementDim>(
320 i, i * displacement_size / DisplacementDim)
324 element_elastic_energy += _ip_data[ip].elastic_energy * w;
326 double const d_ip = N.dot(d);
327 switch (_process_data.phasefield_model)
331 element_surface_energy +=
333 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
339 element_surface_energy += 0.5 * gc *
340 ((1 - d_ip) * (1 - d_ip) / ls +
341 (dNdx * d).dot((dNdx * d)) * ls) *
347 element_surface_energy +=
348 gc / std::numbers::pi *
349 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
355 if (_process_data.pressurized_crack)
357 element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
362 int const n_all_nodes = indices_of_processes[1].size();
363 int const n_regular_nodes = std::count_if(
364 begin(indices_of_processes[1]), end(indices_of_processes[1]),
366 if (n_all_nodes != n_regular_nodes)
368 element_elastic_energy *=
369 static_cast<double>(n_regular_nodes) / n_all_nodes;
370 element_surface_energy *=
371 static_cast<double>(n_regular_nodes) / n_all_nodes;
372 element_pressure_work *=
373 static_cast<double>(n_regular_nodes) / n_all_nodes;
376 elastic_energy += element_elastic_energy;
377 surface_energy += element_surface_energy;
378 pressure_work += element_pressure_work;