OGS
PhaseFieldProcess.cpp
Go to the documentation of this file.
1
11#include "PhaseFieldProcess.h"
12
13#include <cassert>
14
16#include "PhaseFieldFEM.h"
17#include "ProcessLib/Process.h"
19
20namespace ProcessLib
21{
22namespace PhaseField
23{
24template <int DisplacementDim>
26 std::string name,
27 MeshLib::Mesh& mesh,
28 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
29 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
30 unsigned const integration_order,
31 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
32 process_variables,
34 SecondaryVariableCollection&& secondary_variables,
35 bool const use_monolithic_scheme)
36 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
37 integration_order, std::move(process_variables),
38 std::move(secondary_variables), use_monolithic_scheme),
39 _process_data(std::move(process_data))
40{
41 if (use_monolithic_scheme)
42 {
44 "Monolithic scheme is not implemented for the PhaseField process.");
45 }
46
47 _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
48 mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
49}
50
51template <int DisplacementDim>
53{
54 return false;
55}
56
57template <int DisplacementDim>
60 const int process_id) const
61{
62 // For the M process (deformation) in the staggered scheme.
63 if (process_id == 0)
64 {
65 auto const& l = *_local_to_global_index_map;
66 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
67 &l.getGhostIndices(), &this->_sparsity_pattern};
68 }
69
70 // For staggered scheme and phase field process.
71 auto const& l = *_local_to_global_index_map_single_component;
72 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
73 &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
74}
75
76template <int DisplacementDim>
79{
80 // For the M process (deformation) in the staggered scheme.
81 if (process_id == 0)
82 {
83 return *_local_to_global_index_map;
84 }
85
86 // For the equation of phasefield
87 return *_local_to_global_index_map_single_component;
88}
89
90template <int DisplacementDim>
92{
93 // For displacement equation.
94 const int mechanics_process_id = 0;
95 constructDofTableOfSpecifiedProcessStaggeredScheme(mechanics_process_id);
96
97 // TODO move the two data members somewhere else.
98 // for extrapolation of secondary variables of stress or strain
99 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
100 *_mesh_subset_all_nodes};
101 _local_to_global_index_map_single_component =
102 std::make_unique<NumLib::LocalToGlobalIndexMap>(
103 std::move(all_mesh_subsets_single_component),
104 // by location order is needed for output
106
107 assert(_local_to_global_index_map_single_component);
108
109 // For phase field equation.
110 _sparsity_pattern_with_single_component = NumLib::computeSparsityPattern(
111 *_local_to_global_index_map_single_component, _mesh);
112}
113
114template <int DisplacementDim>
116 NumLib::LocalToGlobalIndexMap const& dof_table,
117 MeshLib::Mesh const& mesh,
118 unsigned const integration_order)
119{
121 DisplacementDim, PhaseFieldLocalAssembler>(
122 mesh.getElements(), dof_table, _local_assemblers,
123 mesh.isAxiallySymmetric(), integration_order, _process_data);
124
125 _secondary_variables.addSecondaryVariable(
126 "sigma",
128 DisplacementDim>::RowsAtCompileTime,
129 getExtrapolator(), _local_assemblers,
130 &LocalAssemblerInterface::getIntPtSigma));
131
132 _secondary_variables.addSecondaryVariable(
133 "epsilon",
135 DisplacementDim>::RowsAtCompileTime,
136 getExtrapolator(), _local_assemblers,
137 &LocalAssemblerInterface::getIntPtEpsilon));
138
139 // Initialize local assemblers after all variables have been set.
141 &LocalAssemblerInterface::initialize, _local_assemblers,
142 *_local_to_global_index_map);
143}
144
145template <int DisplacementDim>
147{
148 // Staggered scheme:
149 // for the equations of deformation.
150 const int mechanical_process_id = 0;
151 initializeProcessBoundaryConditionsAndSourceTerms(
152 *_local_to_global_index_map, mechanical_process_id);
153 // for the phase field
154 const int phasefield_process_id = 1;
155 initializeProcessBoundaryConditionsAndSourceTerms(
156 *_local_to_global_index_map_single_component, phasefield_process_id);
157}
158
159template <int DisplacementDim>
161 const double t, double const dt, std::vector<GlobalVector*> const& x,
162 std::vector<GlobalVector*> const& xdot, int const process_id,
164{
165 DBUG("Assemble PhaseFieldProcess.");
166
167 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
168 dof_tables;
169
170 // For the staggered scheme
171 if (process_id == 1)
172 {
173 DBUG(
174 "Assemble the equations of phase field in "
175 "PhaseFieldProcess for the staggered scheme.");
176 }
177 else
178 {
179 DBUG(
180 "Assemble the equations of deformation in "
181 "PhaseFieldProcess for the staggered scheme.");
182 }
183 dof_tables.emplace_back(*_local_to_global_index_map_single_component);
184 dof_tables.emplace_back(*_local_to_global_index_map);
185
186 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
187
188 // Call global assembler for each local assembly item.
190 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
191 pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
192 b);
193}
194
195template <int DisplacementDim>
197 const double t, double const dt, std::vector<GlobalVector*> const& x,
198 std::vector<GlobalVector*> const& xdot, int const process_id,
200{
201 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
202 dof_tables;
203
204 // For the staggered scheme
205 if (process_id == 1)
206 {
207 DBUG(
208 "Assemble the Jacobian equations of phase field in "
209 "PhaseFieldProcess for the staggered scheme.");
210 }
211 else
212 {
213 DBUG(
214 "Assemble the Jacobian equations of deformation in "
215 "PhaseFieldProcess for the staggered scheme.");
216 }
217 dof_tables.emplace_back(*_local_to_global_index_map);
218 dof_tables.emplace_back(*_local_to_global_index_map_single_component);
219
220 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
221
222 // Call global assembler for each local assembly item.
225 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
226 process_id, M, K, b, Jac);
227
228 if (process_id == 0)
229 {
230 b.copyValues(*_nodal_forces);
231 std::transform(_nodal_forces->begin(), _nodal_forces->end(),
232 _nodal_forces->begin(), [](double val) { return -val; });
233 }
234}
235
236template <int DisplacementDim>
238 std::vector<GlobalVector*> const& x, double const t, double const dt,
239 const int process_id)
240{
241 DBUG("PreTimestep PhaseFieldProcess {:d}.", process_id);
242
243 _process_data.injected_volume = t;
244
245 _x_previous_timestep =
247
248 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
249
251 &LocalAssemblerInterface::preTimestep, _local_assemblers,
252 pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
253 dt);
254}
255
256template <int DisplacementDim>
258 std::vector<GlobalVector*> const& x, const double t,
259 const double /*delta_t*/, int const process_id)
260{
261 if (isPhaseFieldProcess(process_id))
262 {
263 DBUG("PostTimestep PhaseFieldProcess.");
264
265 _process_data.elastic_energy = 0.0;
266 _process_data.surface_energy = 0.0;
267 _process_data.pressure_work = 0.0;
268
269 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
270 dof_tables;
271
272 dof_tables.emplace_back(*_local_to_global_index_map);
273 dof_tables.emplace_back(*_local_to_global_index_map_single_component);
274
276 getProcessVariables(process_id)[0];
277
279 &LocalAssemblerInterface::computeEnergy, _local_assemblers,
280 pv.getActiveElementIDs(), dof_tables, *x[process_id], t,
281 _process_data.elastic_energy, _process_data.surface_energy,
282 _process_data.pressure_work, _coupled_solutions);
283
284 INFO("Elastic energy: {:g} Surface energy: {:g} Pressure work: {:g} ",
285 _process_data.elastic_energy, _process_data.surface_energy,
286 _process_data.pressure_work);
287 }
288}
289
290template <int DisplacementDim>
292 GlobalVector const& x, GlobalVector const& /*xdot*/, const double t,
293 double const /*dt*/, const int process_id)
294{
295 _process_data.crack_volume = 0.0;
296
297 if (!isPhaseFieldProcess(process_id))
298 {
299 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
300 dof_tables;
301
302 dof_tables.emplace_back(*_local_to_global_index_map);
303 dof_tables.emplace_back(*_local_to_global_index_map_single_component);
304
305 DBUG("PostNonLinearSolver crack volume computation.");
306
308 getProcessVariables(process_id)[0];
310 &LocalAssemblerInterface::computeCrackIntegral, _local_assemblers,
311 pv.getActiveElementIDs(), dof_tables, x, t,
312 _process_data.crack_volume, _coupled_solutions);
313
314 INFO("Integral of crack: {:g}", _process_data.crack_volume);
315
316 if (_process_data.hydro_crack)
317 {
318 _process_data.pressure_old = _process_data.pressure;
319 _process_data.pressure =
320 _process_data.injected_volume / _process_data.crack_volume;
321 _process_data.pressure_error =
322 std::abs(_process_data.pressure_old - _process_data.pressure) /
323 _process_data.pressure;
324 INFO("Internal pressure: {:g} and Pressure error: {:.4e}",
325 _process_data.pressure, _process_data.pressure_error);
326 auto& u = *_coupled_solutions->coupled_xs[0];
328 _process_data.pressure);
329 }
330 }
331 else
332 {
333 if (_process_data.hydro_crack)
334 {
335 auto& u = *_coupled_solutions->coupled_xs[0];
337 1 / _process_data.pressure);
338 }
339 }
340}
341
342template <int DisplacementDim>
344 GlobalVector& lower, GlobalVector& upper, int const /*process_id*/)
345{
346 lower.setZero();
347 MathLib::LinAlg::setLocalAccessibleVector(*_x_previous_timestep);
348 MathLib::LinAlg::copy(*_x_previous_timestep, upper);
349
350 GlobalIndexType const x_begin = _x_previous_timestep->getRangeBegin();
351 GlobalIndexType const x_end = _x_previous_timestep->getRangeEnd();
352
353 for (GlobalIndexType i = x_begin; i < x_end; i++)
354 {
355 if ((*_x_previous_timestep)[i] > _process_data.irreversible_threshold)
356 {
357 upper.set(i, 1.0);
358 }
359 }
360}
361
362template <int DisplacementDim>
364 int const process_id) const
365{
366 return process_id == 1;
367}
368
369template class PhaseFieldProcess<2>;
370template class PhaseFieldProcess<3>;
371
372} // namespace PhaseField
373} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void copyValues(std::vector< double > &u) const
Definition: EigenVector.h:106
void set(IndexType rowId, double v)
set entry
Definition: EigenVector.h:76
bool isAxiallySymmetric() const
Definition: Mesh.h:131
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
bool isPhaseFieldProcess(int const process_id) const
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void updateConstraints(GlobalVector &lower, GlobalVector &upper, int const process_id) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
MeshLib::PropertyVector< double > * _nodal_forces
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double delta_t, int const process_id) override
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
PhaseFieldProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, PhaseFieldProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void postNonLinearSolverConcreteProcess(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id) override
std::vector< std::size_t > const & getActiveElementIDs() const
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
void scale(PETScVector &x, PetscScalar const a)
Definition: LinAlg.cpp:44
static const double u
static const double t
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
void createLocalAssemblers(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, ExtraCtorArgs &&... extra_ctor_args)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)