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->degradation_df1(d_ip, ls);
170 double const degradation_df2 =
171 _process_data.degradation_derivative->degradation_df2(d_ip, 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 switch (_process_data.phasefield_model)
214 {
216 {
217 auto const local_Jac_AT1 =
218 (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
219 local_Jac.noalias() += local_Jac_AT1;
220
221 local_rhs.noalias() -=
222 gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
223 break;
224 }
226 {
227 auto const local_Jac_AT2 =
228 (gc *
229 (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
230 w)
231 .eval();
232 local_Jac.noalias() += local_Jac_AT2;
233
234 local_rhs.noalias() -=
235 local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
236 break;
237 }
239 {
240 auto const local_Jac_COHESIVE =
241 (2.0 / std::numbers::pi * gc *
242 (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
243 w)
244 .eval();
245
246 local_Jac.noalias() += local_Jac_COHESIVE;
247
248 local_rhs.noalias() -= local_Jac_COHESIVE * d;
249 break;
250 }
251 }
252 }
253}
254
255template <typename ShapeFunction, int DisplacementDim>
258 std::size_t mesh_item_id,
259 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
260 std::vector<GlobalVector*> const& x, double const /*t*/,
261 double& crack_volume)
262{
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); });
269
270 auto local_coupled_xs = getCoupledLocalSolutions(x, indices_of_processes);
271 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
272
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);
277
279 x_position.setElementID(_element.getID());
280
281 int const n_integration_points = _integration_method.getNumberOfPoints();
282 for (int ip = 0; ip < n_integration_points; ip++)
283 {
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;
287
288 typename ShapeMatricesType::template MatrixType<DisplacementDim,
289 displacement_size>
290 N_u = ShapeMatricesType::template MatrixType<
291 DisplacementDim, displacement_size>::Zero(DisplacementDim,
292 displacement_size);
293
294 for (int i = 0; i < DisplacementDim; ++i)
295 {
296 N_u.template block<1, displacement_size / DisplacementDim>(
297 i, i * displacement_size / DisplacementDim)
298 .noalias() = N;
299 }
300
301 crack_volume += (N_u * u).dot(dNdx * d) * w;
302 }
303}
304
305template <typename ShapeFunction, int DisplacementDim>
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)
311{
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); });
318
319 auto const local_coupled_xs =
320 getCoupledLocalSolutions(x, indices_of_processes);
321 assert(local_coupled_xs.size() == displacement_size + phasefield_size);
322
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);
327
329 x_position.setElementID(_element.getID());
330
331 double element_elastic_energy = 0.0;
332 double element_surface_energy = 0.0;
333 double element_pressure_work = 0.0;
334
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++)
339 {
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;
344
345 typename ShapeMatricesType::template MatrixType<DisplacementDim,
346 displacement_size>
347 N_u = ShapeMatricesType::template MatrixType<
348 DisplacementDim, displacement_size>::Zero(DisplacementDim,
349 displacement_size);
350
351 for (int i = 0; i < DisplacementDim; ++i)
352 {
353 N_u.template block<1, displacement_size / DisplacementDim>(
354 i, i * displacement_size / DisplacementDim)
355 .noalias() = N;
356 }
357
358 element_elastic_energy += _ip_data[ip].elastic_energy * w;
359
360 double const d_ip = N.dot(d);
361 switch (_process_data.phasefield_model)
362 {
364 {
365 element_surface_energy +=
366 gc * 0.375 *
367 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
368
369 break;
370 }
372 {
373 element_surface_energy += 0.5 * gc *
374 ((1 - d_ip) * (1 - d_ip) / ls +
375 (dNdx * d).dot((dNdx * d)) * ls) *
376 w;
377 break;
378 }
380 {
381 element_surface_energy +=
382 gc / std::numbers::pi *
383 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
384 w;
385 break;
386 }
387 }
388
389 if (_process_data.pressurized_crack)
390 {
391 element_pressure_work += pressure_ip * (N_u * u).dot(dNdx * d) * w;
392 }
393 }
394
395#ifdef USE_PETSC
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]),
399 [](GlobalIndexType const& index) { return index >= 0; });
400 if (n_all_nodes != n_regular_nodes)
401 {
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;
408 }
409#endif // USE_PETSC
410 elastic_energy += element_elastic_energy;
411 surface_energy += element_surface_energy;
412 pressure_work += element_pressure_work;
413}
414
415template <typename ShapeFunctionDisplacement, int DisplacementDim>
416std::size_t PhaseFieldLocalAssembler<
417 ShapeFunctionDisplacement,
418 DisplacementDim>::setIPDataInitialConditions(std::string_view const name,
419 double const* values,
420 int const integration_order)
421{
422 if (integration_order !=
423 static_cast<int>(_integration_method.getIntegrationOrder()))
424 {
425 OGS_FATAL(
426 "Setting integration point initial conditions; The integration "
427 "order of the local assembler for element {:d} is different from "
428 "the integration order in the initial condition.",
429 _element.getID());
430 }
431
432 if (name == "sigma")
433 {
434 if (_process_data.initial_stress.value != nullptr)
435 {
436 OGS_FATAL(
437 "Setting initial conditions for stress from integration "
438 "point data and from a parameter '{:s}' is not possible "
439 "simultaneously.",
440 _process_data.initial_stress.value->name);
441 }
442
444 values, _ip_data, &IpData::sigma);
445 }
446
447 return 0;
448}
449
450template <typename ShapeFunctionDisplacement, int DisplacementDim>
451std::vector<double> PhaseFieldLocalAssembler<ShapeFunctionDisplacement,
452 DisplacementDim>::getSigma() const
453{
455 _ip_data, &IpData::sigma);
456}
457
458template <typename ShapeFunctionDisplacement, int DisplacementDim>
459std::vector<double> PhaseFieldLocalAssembler<
460 ShapeFunctionDisplacement, DisplacementDim>::getEpsilon() const
461{
462 auto const kelvin_vector_size =
464 unsigned const n_integration_points =
465 _integration_method.getNumberOfPoints();
466
467 std::vector<double> ip_epsilon_values;
468 auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
469 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
470 ip_epsilon_values, n_integration_points, kelvin_vector_size);
471
472 for (unsigned ip = 0; ip < n_integration_points; ++ip)
473 {
474 auto const& eps = _ip_data[ip].eps;
475 cache_mat.row(ip) =
477 }
478
479 return ip_epsilon_values;
480}
481
482} // namespace PhaseField
483} // 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)