OGS
PhaseFieldFEM-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <numbers>
14
17
18namespace ProcessLib
19{
20namespace PhaseField
21{
22template <typename ShapeFunction, int DisplacementDim>
25 double const t, double const dt, Eigen::VectorXd const& local_x,
26 Eigen::VectorXd const& /*local_x_prev*/, int const process_id,
27 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
28{
29 // For the equations with phase field.
30 if (process_id == phase_process_id)
31 {
32 assembleWithJacobianPhaseFieldEquations(t, dt, local_x, local_b_data,
33 local_Jac_data);
34 return;
35 }
36
37 // For the equations with deformation
38 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
39 local_Jac_data);
40}
41
42template <typename ShapeFunction, int DisplacementDim>
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)
47{
48 using DeformationMatrix =
49 typename ShapeMatricesType::template MatrixType<displacement_size,
50 displacement_size>;
51 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
52 auto const u =
53 local_x.template segment<displacement_size>(displacement_index);
54
56 local_Jac_data, displacement_size, displacement_size);
57
59 local_b_data, displacement_size);
60
62 x_position.setElementID(_element.getID());
63
64 auto local_pressure = 0.0;
65 if (_process_data.pressurized_crack)
66 {
67 local_pressure = _process_data.unity_pressure;
68 }
69
70 int const n_integration_points = _integration_method.getNumberOfPoints();
71 for (int ip = 0; ip < n_integration_points; ip++)
72 {
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;
76
77 auto const x_coord =
79 _element, N);
80
81 auto const& B =
82 LinearBMatrix::computeBMatrix<DisplacementDim,
83 ShapeFunction::NPOINTS,
85 dNdx, N, x_coord, _is_axially_symmetric);
86
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);
97
98 auto& sigma = _ip_data[ip].sigma;
99 auto const& D = _ip_data[ip].D;
100
101 typename ShapeMatricesType::template MatrixType<DisplacementDim,
102 displacement_size>
103 N_u = ShapeMatricesType::template MatrixType<
104 DisplacementDim, displacement_size>::Zero(DisplacementDim,
105 displacement_size);
106
107 for (int i = 0; i < DisplacementDim; ++i)
108 {
109 N_u.template block<1, displacement_size / DisplacementDim>(
110 i, i * displacement_size / DisplacementDim)
111 .noalias() = N;
112 }
113
114 auto const rho_sr = _process_data.solid_density(t, x_position)[0];
115 auto const& b = _process_data.specific_body_force;
116
117 local_rhs.noalias() -=
118 (B.transpose() * sigma - N_u.transpose() * rho_sr * b -
119 local_pressure * N_u.transpose() * dNdx * d) *
120 w;
121 local_Jac.noalias() += B.transpose() * D * B * w;
122 }
123}
124
125template <typename ShapeFunction, int DisplacementDim>
127 assembleWithJacobianPhaseFieldEquations(double const t, double const dt,
128 Eigen::VectorXd const& local_x,
129 std::vector<double>& local_b_data,
130 std::vector<double>& local_Jac_data)
131{
132 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
133 auto const u =
134 local_x.template segment<displacement_size>(displacement_index);
135
137 local_Jac_data, phasefield_size, phasefield_size);
139 local_b_data, phasefield_size);
140
142 x_position.setElementID(_element.getID());
143
144 auto local_pressure = 0.0;
145 if (_process_data.static_pressurized_crack)
146 {
147 local_pressure = _process_data.unity_pressure;
148 }
149 else if (_process_data.propagating_pressurized_crack)
150 {
151 local_pressure = _process_data.pressure;
152 }
153
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];
157
158 int const n_integration_points = _integration_method.getNumberOfPoints();
159 for (int ip = 0; ip < n_integration_points; ip++)
160 {
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;
164
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);
172 auto const x_coord =
174 _element, N);
175 auto const& B =
176 LinearBMatrix::computeBMatrix<DisplacementDim,
177 ShapeFunction::NPOINTS,
179 dNdx, N, x_coord, _is_axially_symmetric);
180
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);
186
187 auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
188
189 auto& ip_data = _ip_data[ip];
190 ip_data.strain_energy_tensile = strain_energy_tensile;
191
192 typename ShapeMatricesType::template MatrixType<DisplacementDim,
193 displacement_size>
194 N_u = ShapeMatricesType::template MatrixType<
195 DisplacementDim, displacement_size>::Zero(DisplacementDim,
196 displacement_size);
197
198 for (int i = 0; i < DisplacementDim; ++i)
199 {
200 N_u.template block<1, displacement_size / DisplacementDim>(
201 i, i * displacement_size / DisplacementDim)
202 .noalias() = N;
203 }
204
205 local_Jac.noalias() +=
206 (N.transpose() * N * degradation_df2 * strain_energy_tensile) * w;
207
208 local_rhs.noalias() -=
209 (N.transpose() * degradation_df1 * strain_energy_tensile -
210 local_pressure * dNdx.transpose() * N_u * u) *
211 w;
212
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);
218 }
219}
220
221template <typename ShapeFunction, int DisplacementDim>
224 std::size_t mesh_item_id,
225 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
226 std::vector<GlobalVector*> const& x, double const /*t*/,
227 double& crack_volume)
228{
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); });
235
236 auto local_coupled_xs = getCoupledLocalSolutions(x, indices_of_processes);
237 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
238
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);
243
245 x_position.setElementID(_element.getID());
246
247 int const n_integration_points = _integration_method.getNumberOfPoints();
248 for (int ip = 0; ip < n_integration_points; ip++)
249 {
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;
253
254 typename ShapeMatricesType::template MatrixType<DisplacementDim,
255 displacement_size>
256 N_u = ShapeMatricesType::template MatrixType<
257 DisplacementDim, displacement_size>::Zero(DisplacementDim,
258 displacement_size);
259
260 for (int i = 0; i < DisplacementDim; ++i)
261 {
262 N_u.template block<1, displacement_size / DisplacementDim>(
263 i, i * displacement_size / DisplacementDim)
264 .noalias() = N;
265 }
266
267 crack_volume += (N_u * u).dot(dNdx * d) * w;
268 }
269}
270
271template <typename ShapeFunction, int DisplacementDim>
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)
277{
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); });
284
285 auto const local_coupled_xs =
286 getCoupledLocalSolutions(x, indices_of_processes);
287 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
288
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);
293
295 x_position.setElementID(_element.getID());
296
297 double element_elastic_energy = 0.0;
298 double element_surface_energy = 0.0;
299 double element_pressure_work = 0.0;
300
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++)
305 {
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;
310
311 typename ShapeMatricesType::template MatrixType<DisplacementDim,
312 displacement_size>
313 N_u = ShapeMatricesType::template MatrixType<
314 DisplacementDim, displacement_size>::Zero(DisplacementDim,
315 displacement_size);
316
317 for (int i = 0; i < DisplacementDim; ++i)
318 {
319 N_u.template block<1, displacement_size / DisplacementDim>(
320 i, i * displacement_size / DisplacementDim)
321 .noalias() = N;
322 }
323
324 element_elastic_energy += _ip_data[ip].elastic_energy * w;
325
326 double const d_ip = N.dot(d);
327 switch (_process_data.phasefield_model)
328 {
330 {
331 element_surface_energy +=
332 gc * 0.375 *
333 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
334
335 break;
336 }
338 {
339 element_surface_energy += 0.5 * gc *
340 ((1 - d_ip) * (1 - d_ip) / ls +
341 (dNdx * d).dot((dNdx * d)) * ls) *
342 w;
343 break;
344 }
346 {
347 element_surface_energy +=
348 gc / std::numbers::pi *
349 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
350 w;
351 break;
352 }
353 }
354
355 if (_process_data.pressurized_crack)
356 {
357 element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
358 }
359 }
360
361#ifdef USE_PETSC
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]),
365 [](GlobalIndexType const& index) { return index >= 0; });
366 if (n_all_nodes != n_regular_nodes)
367 {
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;
374 }
375#endif // USE_PETSC
376 elastic_energy += element_elastic_energy;
377 surface_energy += element_surface_energy;
378 pressure_work += element_pressure_work;
379}
380
381template <typename ShapeFunctionDisplacement, int DisplacementDim>
382std::size_t PhaseFieldLocalAssembler<
383 ShapeFunctionDisplacement,
384 DisplacementDim>::setIPDataInitialConditions(std::string_view const name,
385 double const* values,
386 int const integration_order)
387{
388 if (integration_order !=
389 static_cast<int>(_integration_method.getIntegrationOrder()))
390 {
391 OGS_FATAL(
392 "Setting integration point initial conditions; The integration "
393 "order of the local assembler for element {:d} is different from "
394 "the integration order in the initial condition.",
395 _element.getID());
396 }
397
398 if (name == "sigma")
399 {
400 if (_process_data.initial_stress.value != nullptr)
401 {
402 OGS_FATAL(
403 "Setting initial conditions for stress from integration "
404 "point data and from a parameter '{:s}' is not possible "
405 "simultaneously.",
406 _process_data.initial_stress.value->name);
407 }
408
410 values, _ip_data, &IpData::sigma);
411 }
412
413 return 0;
414}
415
416template <typename ShapeFunctionDisplacement, int DisplacementDim>
417std::vector<double> PhaseFieldLocalAssembler<ShapeFunctionDisplacement,
418 DisplacementDim>::getSigma() const
419{
421 _ip_data, &IpData::sigma);
422}
423
424template <typename ShapeFunctionDisplacement, int DisplacementDim>
425std::vector<double> PhaseFieldLocalAssembler<
426 ShapeFunctionDisplacement, DisplacementDim>::getEpsilon() const
427{
428 auto const kelvin_vector_size =
430 unsigned const n_integration_points =
431 _integration_method.getNumberOfPoints();
432
433 std::vector<double> ip_epsilon_values;
434 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
435 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
436 ip_epsilon_values, n_integration_points, kelvin_vector_size);
437
438 for (unsigned ip = 0; ip < n_integration_points; ++ip)
439 {
440 auto const& eps = _ip_data[ip].eps;
441 cache_mat.row(ip) =
443 }
444
445 return ip_epsilon_values;
446}
447
448} // namespace PhaseField
449} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void setElementID(std::size_t element_id)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianPhaseFieldEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void computeCrackIntegral(std::size_t mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double &crack_volume) override
void computeEnergy(std::size_t mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double &elastic_energy, double &surface_energy, double &pressure_work) override
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)