OGS
ThermoMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
20#include "ThermoMechanicsFEM.h"
21
22namespace ProcessLib
23{
24namespace ThermoMechanics
25{
26template <int DisplacementDim>
28 std::string name, MeshLib::Mesh& mesh,
29 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
30 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
31 unsigned const integration_order,
32 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
33 process_variables,
35 SecondaryVariableCollection&& secondary_variables,
36 bool const use_monolithic_scheme)
37 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
38 integration_order, std::move(process_variables),
39 std::move(secondary_variables), use_monolithic_scheme),
40 _process_data(std::move(process_data))
41{
42 _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
43 mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
44
45 _heat_flux = MeshLib::getOrCreateMeshProperty<double>(
46 mesh, "HeatFlowRate", MeshLib::MeshItemType::Node, 1);
47
48 _integration_point_writer.emplace_back(
49 std::make_unique<IntegrationPointWriter>(
50 "sigma_ip",
51 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
52 integration_order, _local_assemblers,
54
55 _integration_point_writer.emplace_back(
56 std::make_unique<IntegrationPointWriter>(
57 "epsilon_ip",
58 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
59 integration_order, _local_assemblers,
61
62 _integration_point_writer.emplace_back(
63 std::make_unique<IntegrationPointWriter>(
64 "epsilon_m_ip",
65 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
66 integration_order, _local_assemblers,
68}
69
70template <int DisplacementDim>
72{
73 return false;
74}
75
76template <int DisplacementDim>
79 const int process_id) const
80{
81 // For the monolithic scheme or the M process (deformation) in the staggered
82 // scheme.
83 if (_use_monolithic_scheme ||
84 process_id == _process_data.mechanics_process_id)
85 {
86 auto const& l = *_local_to_global_index_map;
87 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
88 &l.getGhostIndices(), &this->_sparsity_pattern};
89 }
90
91 // For staggered scheme and T process.
92 auto const& l = *_local_to_global_index_map_single_component;
93 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
94 &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
95}
96
97// TODO [WW]: remove if (_use_monolithic_scheme) during the refactoring of the
98// coupling part.
99template <int DisplacementDim>
101{
102 // Note: the heat conduction process and the mechanical process use the same
103 // order of shape functions.
104
105 if (_use_monolithic_scheme)
106 {
107 constructMonolithicProcessDofTable();
108 return;
109 }
110 constructDofTableOfSpecifiedProcessStaggeredScheme(
111 _process_data.mechanics_process_id);
112
113 // TODO move the two data members somewhere else.
114 // for extrapolation of secondary variables of stress or strain
115 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
116 *_mesh_subset_all_nodes};
117 _local_to_global_index_map_single_component.reset(
119 std::move(all_mesh_subsets_single_component),
120 // by location order is needed for output
122
123 if (!_use_monolithic_scheme)
124 {
125 _sparsity_pattern_with_single_component =
127 *_local_to_global_index_map_single_component, _mesh);
128 }
129}
130
131template <int DisplacementDim>
133 NumLib::LocalToGlobalIndexMap const& dof_table,
134 MeshLib::Mesh const& mesh,
135 unsigned const integration_order)
136{
138 DisplacementDim, ThermoMechanicsLocalAssembler>(
139 mesh.getElements(), dof_table, _local_assemblers,
140 mesh.isAxiallySymmetric(), integration_order, _process_data);
141
142 auto add_secondary_variable = [&](std::string const& name,
143 int const num_components,
144 auto get_ip_values_function)
145 {
146 _secondary_variables.addSecondaryVariable(
147 name,
148 makeExtrapolator(num_components, getExtrapolator(),
149 _local_assemblers,
150 std::move(get_ip_values_function)));
151 };
152
153 add_secondary_variable("sigma",
155 DisplacementDim>::RowsAtCompileTime,
156 &LocalAssemblerInterface::getIntPtSigma);
157
158 add_secondary_variable("epsilon",
160 DisplacementDim>::RowsAtCompileTime,
161 &LocalAssemblerInterface::getIntPtEpsilon);
162
163 //
164 // enable output of internal variables defined by material models
165 //
167 LocalAssemblerInterface>(_process_data.solid_materials,
168 add_secondary_variable);
169
170 // Set initial conditions for integration point data.
171 for (auto const& ip_writer : _integration_point_writer)
172 {
173 // Find the mesh property with integration point writer's name.
174 auto const& name = ip_writer->name();
175 if (!mesh.getProperties().existsPropertyVector<double>(name))
176 {
177 continue;
178 }
179 auto const& mesh_property =
180 *mesh.getProperties().template getPropertyVector<double>(name);
181
182 // The mesh property must be defined on integration points.
183 if (mesh_property.getMeshItemType() !=
185 {
186 continue;
187 }
188
189 auto const ip_meta_data =
191
192 // Check the number of components.
193 if (ip_meta_data.n_components !=
194 mesh_property.getNumberOfGlobalComponents())
195 {
196 OGS_FATAL(
197 "Different number of components in meta data ({:d}) than in "
198 "the integration point field data for '{:s}': {:d}.",
199 ip_meta_data.n_components, name,
200 mesh_property.getNumberOfGlobalComponents());
201 }
202
203 // Now we have a properly named vtk's field data array and the
204 // corresponding meta data.
205 std::size_t position = 0;
206 for (auto& local_asm : _local_assemblers)
207 {
208 std::size_t const integration_points_read =
209 local_asm->setIPDataInitialConditions(
210 name, &mesh_property[position],
211 ip_meta_data.integration_order);
212 if (integration_points_read == 0)
213 {
214 OGS_FATAL(
215 "No integration points read in the integration point "
216 "initial conditions set function.");
217 }
218 position += integration_points_read * ip_meta_data.n_components;
219 }
220 }
221
222 // Initialize local assemblers after all variables have been set.
224 &LocalAssemblerInterface::initialize, _local_assemblers,
225 *_local_to_global_index_map);
226}
227
228template <int DisplacementDim>
230{
231 if (_use_monolithic_scheme)
232 {
233 const int process_id_of_thermomechanics = 0;
234 initializeProcessBoundaryConditionsAndSourceTerms(
235 *_local_to_global_index_map, process_id_of_thermomechanics);
236 return;
237 }
238
239 // Staggered scheme:
240 // for the equations of heat conduction
241 initializeProcessBoundaryConditionsAndSourceTerms(
242 *_local_to_global_index_map_single_component,
243 _process_data.heat_conduction_process_id);
244
245 // for the equations of deformation.
246 initializeProcessBoundaryConditionsAndSourceTerms(
247 *_local_to_global_index_map, _process_data.mechanics_process_id);
248}
249
250template <int DisplacementDim>
252 const double t, double const dt, std::vector<GlobalVector*> const& x,
253 std::vector<GlobalVector*> const& xdot, int const process_id,
255{
256 DBUG("Assemble ThermoMechanicsProcess.");
257
258 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
259 dof_table = {std::ref(*_local_to_global_index_map)};
260 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
261
262 // Call global assembler for each local assembly item.
264 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
265 pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
266 b);
267}
268
269template <int DisplacementDim>
271 assembleWithJacobianConcreteProcess(const double t, double const dt,
272 std::vector<GlobalVector*> const& x,
273 std::vector<GlobalVector*> const& xdot,
274 int const process_id, GlobalMatrix& M,
276 GlobalMatrix& Jac)
277{
278 DBUG("AssembleJacobian ThermoMechanicsProcess.");
279
280 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
281 dof_tables;
282 // For the monolithic scheme
283 if (_use_monolithic_scheme)
284 {
285 DBUG(
286 "Assemble the Jacobian of ThermoMechanics for the monolithic"
287 " scheme.");
288 dof_tables.emplace_back(*_local_to_global_index_map);
289 }
290 else
291 {
292 // For the staggered scheme
293 if (process_id == _process_data.heat_conduction_process_id)
294 {
295 DBUG(
296 "Assemble the Jacobian equations of heat conduction process in "
297 "ThermoMechanics for the staggered scheme.");
298 }
299 else
300 {
301 DBUG(
302 "Assemble the Jacobian equations of mechanical process in "
303 "ThermoMechanics for the staggered scheme.");
304 }
305
306 // For the flexible appearance order of processes in the coupling.
307 if (_process_data.heat_conduction_process_id ==
308 0) // First: the heat conduction process
309 {
310 dof_tables.emplace_back(
311 *_local_to_global_index_map_single_component);
312 dof_tables.emplace_back(*_local_to_global_index_map);
313 }
314 else // vice versa
315 {
316 dof_tables.emplace_back(*_local_to_global_index_map);
317 dof_tables.emplace_back(
318 *_local_to_global_index_map_single_component);
319 }
320 }
321
322 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
323
326 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
327 process_id, M, K, b, Jac);
328
329 // TODO (naumov): Refactor the copy rhs part. This is copy from HM.
330 auto copyRhs = [&](int const variable_id, auto& output_vector)
331 {
332 if (_use_monolithic_scheme)
333 {
334 transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
335 output_vector,
336 std::negate<double>());
337 }
338 else
339 {
340 transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
341 output_vector,
342 std::negate<double>());
343 }
344 };
345 if (_use_monolithic_scheme ||
346 process_id == _process_data.heat_conduction_process_id)
347 {
348 copyRhs(0, *_heat_flux);
349 }
350 if (_use_monolithic_scheme ||
351 process_id == _process_data.mechanics_process_id)
352 {
353 copyRhs(1, *_nodal_forces);
354 }
355}
356
357template <int DisplacementDim>
359 std::vector<GlobalVector*> const& x, double const t, double const dt,
360 const int process_id)
361{
362 DBUG("PreTimestep ThermoMechanicsProcess.");
363
364 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
365
366 assert(process_id < 2);
367
368 if (process_id == _process_data.mechanics_process_id)
369 {
371 &LocalAssemblerInterface::preTimestep, _local_assemblers,
372 pv.getActiveElementIDs(), *_local_to_global_index_map,
373 *x[process_id], t, dt);
374 return;
375 }
376}
377
378template <int DisplacementDim>
380 std::vector<GlobalVector*> const& x, double const t, double const dt,
381 int const process_id)
382{
383 if (process_id != 0)
384 {
385 return;
386 }
387
388 DBUG("PostTimestep ThermoMechanicsProcess.");
389 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
390 auto const n_processes = x.size();
391 dof_tables.reserve(n_processes);
392 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
393 {
394 dof_tables.push_back(&getDOFTable(process_id));
395 }
396
397 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
399 &LocalAssemblerInterface::postTimestep, _local_assemblers,
400 pv.getActiveElementIDs(), dof_tables, x, t, dt);
401}
402
403template <int DisplacementDim>
406{
407 if (_process_data.mechanics_process_id == process_id)
408 {
409 return *_local_to_global_index_map;
410 }
411
412 // For the equation of pressure
413 return *_local_to_global_index_map_single_component;
414}
415
416template class ThermoMechanicsProcess<2>;
417template class ThermoMechanicsProcess<3>;
418
419} // namespace ThermoMechanics
420} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:131
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
Properties & getProperties()
Definition: Mesh.h:128
bool existsPropertyVector(std::string const &name) const
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
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)
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:352
Handles configuration of several secondary variables from the project file.
Local assembler of ThermoMechanics process.
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) 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
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
ThermoMechanicsProcess(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, ThermoMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const 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 initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
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
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 transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
Definition: DOFTableUtil.h:59
void solidMaterialInternalToSecondaryVariables(std::map< int, std::unique_ptr< MaterialLib::Solids::MechanicsBase< DisplacementDim > > > const &solid_materials, AddSecondaryVariableCallback const &add_secondary_variable)
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)
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Properties const &properties, std::string const &name)
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)
virtual std::vector< double > getEpsilonMechanical() const =0