OGS
ProcessVariable.cpp
Go to the documentation of this file.
1
11#include "ProcessVariable.h"
12
13#include <algorithm>
14#include <utility>
15
16#include "BaseLib/Algorithm.h"
17#include "BaseLib/Logging.h"
19#include "MeshLib/Mesh.h"
20#include "MeshLib/Node.h"
22#include "ParameterLib/Utils.h"
32
33namespace
34{
36 BaseLib::ConfigTree const& config,
37 std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
38{
39 //
40 // Get the mesh name from the config.
41 //
42 std::string mesh_name; // Either given directly in <mesh> or constructed
43 // from <geometrical_set>_<geometry>.
44
45#ifdef DOXYGEN_DOCU_ONLY
47 config.getConfigParameterOptional<std::string>("mesh");
48#endif // DOXYGEN_DOCU_ONLY
49
50 auto optional_mesh_name =
52 config.getConfigParameterOptional<std::string>("mesh");
53 if (optional_mesh_name)
54 {
55 mesh_name = *optional_mesh_name;
56 }
57 else
58 {
59#ifdef DOXYGEN_DOCU_ONLY
61 config.getConfigParameterOptional<std::string>("geometrical_set");
63 config.getConfigParameter<std::string>("geometry");
64#endif // DOXYGEN_DOCU_ONLY
65
66 // Looking for the mesh created before for the given geometry.
67 auto const geometrical_set_name =
69 config.getConfigParameter<std::string>("geometrical_set");
70 auto const geometry_name =
72 config.getConfigParameter<std::string>("geometry");
73
74 mesh_name = MeshGeoToolsLib::meshNameFromGeometry(geometrical_set_name,
75 geometry_name);
76 }
77
78 //
79 // Find and extract mesh from the list of meshes.
80 //
81 auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
82 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
83
84 return mesh;
85}
86} // namespace
87
88namespace ProcessLib
89{
91 BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
92 std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
93 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
94 std::map<std::string,
95 std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
96 curves)
97 :
98 _name(config.getConfigParameter<std::string>("name")),
99 _mesh(mesh),
101 _n_components(config.getConfigParameter<int>("components")),
103 _shapefunction_order(config.getConfigParameter<unsigned>("order")),
104 _deactivated_subdomains(
105 createDeactivatedSubdomains(config, mesh, parameters, curves)),
106 _initial_condition(ParameterLib::findParameter<double>(
108 config.getConfigParameter<std::string>("initial_condition"),
109 parameters, _n_components, &mesh)),
110 _compensate_non_equilibrium_initial_residuum(
112 config.getConfigParameter<bool>(
113 "compensate_non_equilibrium_initial_residuum", false))
114{
115 DBUG("Constructing process variable {:s}", _name);
116
118 {
119 OGS_FATAL("The given shape function order {:d} is not supported",
121 }
122
123 // Boundary conditions
124 if (auto bcs_config =
126 config.getConfigSubtreeOptional("boundary_conditions"))
127 {
128 for (
129 auto bc_config :
131 bcs_config->getConfigSubtreeList("boundary_condition"))
132 {
133 auto const& bc_mesh = findMeshInConfig(bc_config, meshes);
134 auto component_id =
136 bc_config.getConfigParameterOptional<int>("component");
137
138 if (!component_id && _n_components == 1)
139 {
140 // default value for single component vars.
141 component_id = 0;
142 }
143
144 _bc_configs.emplace_back(
145 std::move(bc_config), bc_mesh, component_id,
147 }
148 }
149 else
150 {
151 INFO("No boundary conditions for process variable '{:s}' found.",
152 _name);
153 }
154
155 // Source terms
157 if (auto sts_config = config.getConfigSubtreeOptional("source_terms"))
158 {
159 for (
160 auto st_config :
162 sts_config->getConfigSubtreeList("source_term"))
163 {
164 MeshLib::Mesh const& st_mesh = findMeshInConfig(st_config, meshes);
165 auto component_id =
167 st_config.getConfigParameterOptional<int>("component");
168 auto const type =
170 st_config.peekConfigParameter<std::string>("type");
171
172 if (!component_id)
173 {
174 if (_n_components == 1)
175 {
176 // default value for single component vars.
177 component_id = 0;
178 }
179 else if (type == "Anchor")
180 {
181 // dummy value
182 component_id = 0;
183 }
184 else if (type == "EmbeddedAnchor")
185 {
186 // dummy value
187 component_id = 0;
188 }
189 else
190 {
191 OGS_FATAL(
192 "Specifying the component id (<component>) for a "
193 "source term for a non-scalar process variable is "
194 "mandatory.");
195 }
196 }
197
198 _source_term_configs.emplace_back(std::move(st_config), st_mesh,
199 *component_id);
200 }
201 }
202 else
203 {
204 INFO("No source terms for process variable '{:s}' found.", _name);
205 }
206
207 if (!_deactivated_subdomains.empty())
208 {
209 _is_active = getOrCreateMeshProperty<unsigned char>(
210 _mesh, _name + "_active", MeshLib::MeshItemType::Cell, 1);
211 std::fill(std::begin(*_is_active), std::end(*_is_active), 1u);
212 }
213}
214
216
218{
219 return _mesh;
220}
221
222std::vector<std::unique_ptr<BoundaryCondition>>
224 const NumLib::LocalToGlobalIndexMap& dof_table,
225 const int variable_id,
226 unsigned const integration_order,
227 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
228 Process const& process,
229 std::vector<std::reference_wrapper<ProcessVariable>> const&
230 all_process_variables_for_this_process,
231 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
232{
233 std::vector<std::unique_ptr<BoundaryCondition>> bcs;
234 bcs.reserve(_bc_configs.size());
235
236 for (auto const& config : _bc_configs)
237 {
238 auto bc = createBoundaryCondition(
239 config, dof_table, _mesh, variable_id, integration_order,
240 _shapefunction_order, parameters, process,
241 all_process_variables_for_this_process, media);
242#ifdef USE_PETSC
243 if (bc == nullptr)
244 {
245 continue;
246 }
247#endif // USE_PETSC
248 bcs.push_back(std::move(bc));
249 }
250
252 parameters, bcs);
253
254 return bcs;
255}
256
258 const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
259 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
260 std::vector<std::unique_ptr<BoundaryCondition>>& bcs)
261{
262 for (auto const& deactivated_subdomain : _deactivated_subdomains)
263 {
264 auto const& deactivated_subdomain_mesh =
265 deactivated_subdomain.deactivated_subdomain_mesh;
266 auto const* parameter = &ParameterLib::findParameter<double>(
268 bool const set_outer_nodes_dirichlet_values =
269 deactivated_subdomain.boundary_value_parameter != nullptr;
270 if (set_outer_nodes_dirichlet_values)
271 {
272 parameter = deactivated_subdomain.boundary_value_parameter;
273 }
274
275 for (int component_id = 0;
276 component_id <
277 dof_table.getNumberOfVariableComponents(variable_id);
278 component_id++)
279 {
280 auto bc = std::make_unique<DeactivatedSubdomainDirichlet>(
281 *_is_active, deactivated_subdomain.time_interval, *parameter,
282 set_outer_nodes_dirichlet_values, deactivated_subdomain_mesh,
283 dof_table, variable_id, component_id);
284
285#ifdef USE_PETSC
286 // TODO: make it work under PETSc too.
287 if (bc == nullptr)
288 {
289 continue;
290 }
291#endif // USE_PETSC
292 bcs.push_back(std::move(bc));
293 }
294 }
295}
296
298{
299 if (_deactivated_subdomains.empty())
300 {
301 return;
302 }
303
305
306 // If none of the deactivated subdomains is active at current time, then the
307 // _ids_of_active_elements remain empty.
308 if (std::none_of(begin(_deactivated_subdomains),
309 end(_deactivated_subdomains), [&](auto const& ds)
310 { return ds.isInTimeSupportInterval(time); }))
311 {
312 // Also mark all of the elements as active.
313 assert(_is_active != nullptr); // guaranteed by constructor
314 std::fill(std::begin(*_is_active), std::end(*_is_active), 1u);
315
316 return;
317 }
318
319 auto is_active_in_subdomain = [&](std::size_t const element_id,
320 DeactivatedSubdomain const& ds) -> bool
321 {
322 return (!ds.isInTimeSupportInterval(time)) ||
323 !ds.isDeactivated(*_mesh.getElement(element_id), time);
324 };
325
326 auto is_active_in_all_subdomains = [&](std::size_t const element_id) -> bool
327 {
328 return std::all_of(begin(_deactivated_subdomains),
330 [&](auto const& ds)
331 { return is_active_in_subdomain(element_id, ds); });
332 };
333
334 auto const number_of_elements = _mesh.getNumberOfElements();
335 for (std::size_t element_id = 0; element_id < number_of_elements;
336 element_id++)
337 {
338 if (is_active_in_all_subdomains(element_id))
339 {
340 _ids_of_active_elements.push_back(element_id);
341 }
342 }
343
344 // all elements are deactivated
345 std::fill(std::begin(*_is_active), std::end(*_is_active), 0u);
346
347 for (auto const id : _ids_of_active_elements)
348 {
349 (*_is_active)[id] = 1u;
350 }
351}
352
353std::vector<std::unique_ptr<SourceTermBase>> ProcessVariable::createSourceTerms(
354 const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
355 unsigned const integration_order,
356 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
357 std::vector<std::reference_wrapper<ProcessVariable>> const&
358 all_process_variables_for_this_process,
359 const MeshLib::Mesh& bulk_mesh)
360{
361 std::vector<std::unique_ptr<SourceTermBase>> source_terms;
362
363 transform(cbegin(_source_term_configs), cend(_source_term_configs),
364 back_inserter(source_terms),
365 [&](auto const& config)
366 {
367 return createSourceTerm(
368 config, dof_table, config.mesh, variable_id,
369 integration_order, _shapefunction_order, parameters,
370 all_process_variables_for_this_process, bulk_mesh);
371 });
372
373 return source_terms;
374}
375
377
378} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the Mesh class.
Definition of the Node class.
std::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
std::optional< T > getConfigParameterOptional(std::string const &param) const
T getConfigParameter(std::string const &param) const
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:96
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
int getNumberOfVariableComponents(int variable_id) const
void updateDeactivatedSubdomains(double const time)
std::vector< std::unique_ptr< SourceTermBase > > createSourceTerms(const NumLib::LocalToGlobalIndexMap &dof_table, const int variable_id, unsigned const integration_order, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::vector< std::reference_wrapper< ProcessVariable > > const &all_process_variables_for_this_process, const MeshLib::Mesh &bulk_mesh)
std::vector< BoundaryConditionConfig > _bc_configs
ProcessVariable(BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh, std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::map< std::string, std::unique_ptr< MathLib::PiecewiseLinearInterpolation > > const &curves)
MeshLib::Mesh const & getMesh() const
Returns a mesh on which the process variable is defined.
std::vector< DeactivatedSubdomain > _deactivated_subdomains
const bool _compensate_non_equilibrium_initial_residuum
MeshLib::PropertyVector< unsigned char > * _is_active
std::vector< std::size_t > _ids_of_active_elements
std::vector< SourceTermConfig > _source_term_configs
std::vector< std::unique_ptr< BoundaryCondition > > createBoundaryConditions(const NumLib::LocalToGlobalIndexMap &dof_table, const int variable_id, unsigned const integration_order, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, Process const &process, std::vector< std::reference_wrapper< ProcessVariable > > const &all_process_variables_for_this_process, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
void createBoundaryConditionsForDeactivatedSubDomains(const NumLib::LocalToGlobalIndexMap &dof_table, const int variable_id, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::vector< std::unique_ptr< BoundaryCondition > > &bcs)
std::string meshNameFromGeometry(std::string const &geometrical_set_name, std::string const &geometry_name)
Mesh & findMeshByName(std::vector< std::unique_ptr< Mesh > > const &meshes, std::string_view const name)
Definition Mesh.cpp:365
OGS_NO_DANGLING Parameter< ParameterDataType > & findParameter(std::string const &parameter_name, std::vector< std::unique_ptr< ParameterBase > > const &parameters, int const num_components, MeshLib::Mesh const *const mesh=nullptr)
Definition Utils.h:102
std::vector< DeactivatedSubdomain > createDeactivatedSubdomains(BaseLib::ConfigTree const &config, MeshLib::Mesh const &mesh, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::map< std::string, std::unique_ptr< MathLib::PiecewiseLinearInterpolation > > const &curves)
std::unique_ptr< BoundaryCondition > createBoundaryCondition(const BoundaryConditionConfig &config, const NumLib::LocalToGlobalIndexMap &dof_table, const MeshLib::Mesh &bulk_mesh, const int variable_id, const unsigned integration_order, const unsigned shapefunction_order, const std::vector< std::unique_ptr< ParameterLib::ParameterBase > > &parameters, const Process &process, std::vector< std::reference_wrapper< ProcessVariable > > const &all_process_variables_for_this_process, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
std::unique_ptr< SourceTermBase > createSourceTerm(const SourceTermConfig &config, const NumLib::LocalToGlobalIndexMap &dof_table_bulk, const MeshLib::Mesh &source_term_mesh, const int variable_id, const unsigned integration_order, const unsigned shapefunction_order, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::vector< std::reference_wrapper< ProcessVariable > > const &all_process_variables_for_this_process, const MeshLib::Mesh &bulk_mesh)
MeshLib::Mesh const & findMeshInConfig(BaseLib::ConfigTree const &config, std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes)
static PROCESSLIB_EXPORT const std::string zero_parameter_name