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"
21#include "ParameterLib/Utils.h"
31
32namespace
33{
35 BaseLib::ConfigTree const& config,
36 std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
37{
38 //
39 // Get the mesh name from the config.
40 //
41 std::string mesh_name; // Either given directly in <mesh> or constructed
42 // from <geometrical_set>_<geometry>.
43
44#ifdef DOXYGEN_DOCU_ONLY
46 config.getConfigParameterOptional<std::string>("mesh");
47#endif // DOXYGEN_DOCU_ONLY
48
49 auto optional_mesh_name =
51 config.getConfigParameterOptional<std::string>("mesh");
52 if (optional_mesh_name)
53 {
54 mesh_name = *optional_mesh_name;
55 }
56 else
57 {
58#ifdef DOXYGEN_DOCU_ONLY
60 config.getConfigParameterOptional<std::string>("geometrical_set");
62 config.getConfigParameter<std::string>("geometry");
63#endif // DOXYGEN_DOCU_ONLY
64
65 // Looking for the mesh created before for the given geometry.
66 auto const geometrical_set_name =
68 config.getConfigParameter<std::string>("geometrical_set");
69 auto const geometry_name =
71 config.getConfigParameter<std::string>("geometry");
72
73 mesh_name = MeshGeoToolsLib::meshNameFromGeometry(geometrical_set_name,
74 geometry_name);
75 }
76
77 //
78 // Find and extract mesh from the list of meshes.
79 //
80 auto const& mesh = *BaseLib::findElementOrError(
81 begin(meshes), end(meshes),
82 [&mesh_name](auto const& mesh)
83 {
84 assert(mesh != nullptr);
85 return mesh->getName() == mesh_name;
86 },
87 "Required mesh with name '" + mesh_name + "' not found.");
88 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
89
90 return mesh;
91}
92} // namespace
93
94namespace ProcessLib
95{
97 BaseLib::ConfigTree const& config)
98{
99 auto const compensate_non_equilibrium_initial_residuum_ptr =
101 config.getConfigParameterOptional<bool>(
102 "compensate_non_equilibrium_initial_residuum");
103
104 return (compensate_non_equilibrium_initial_residuum_ptr &&
105 *compensate_non_equilibrium_initial_residuum_ptr);
106}
107
109 BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
110 std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
111 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
112 std::map<std::string,
113 std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
114 curves)
115 :
116 _name(config.getConfigParameter<std::string>("name")),
117 _mesh(mesh),
119 _n_components(config.getConfigParameter<int>("components")),
121 _shapefunction_order(config.getConfigParameter<unsigned>("order")),
122 _deactivated_subdomains(
123 createDeactivatedSubdomains(config, mesh, parameters, curves)),
124 _initial_condition(ParameterLib::findParameter<double>(
126 config.getConfigParameter<std::string>("initial_condition"),
127 parameters, _n_components, &mesh)),
128 _compensate_non_equilibrium_initial_residuum(
130{
131 DBUG("Constructing process variable {:s}", _name);
132
134 {
135 OGS_FATAL("The given shape function order {:d} is not supported",
137 }
138
139 // Boundary conditions
140 if (auto bcs_config =
142 config.getConfigSubtreeOptional("boundary_conditions"))
143 {
144 for (
145 auto bc_config :
147 bcs_config->getConfigSubtreeList("boundary_condition"))
148 {
149 auto const& bc_mesh = findMeshInConfig(bc_config, meshes);
150 auto component_id =
152 bc_config.getConfigParameterOptional<int>("component");
153
154 if (!component_id && _n_components == 1)
155 {
156 // default value for single component vars.
157 component_id = 0;
158 }
159
160 _bc_configs.emplace_back(std::move(bc_config), bc_mesh,
161 component_id);
162 }
163 }
164 else
165 {
166 INFO("No boundary conditions for process variable '{:s}' found.",
167 _name);
168 }
169
170 // Source terms
172 if (auto sts_config = config.getConfigSubtreeOptional("source_terms"))
173 {
174 for (
175 auto st_config :
177 sts_config->getConfigSubtreeList("source_term"))
178 {
179 MeshLib::Mesh const& st_mesh = findMeshInConfig(st_config, meshes);
180 auto component_id =
182 st_config.getConfigParameterOptional<int>("component");
183
184 if (!component_id && _n_components == 1)
185 {
186 // default value for single component vars.
187 component_id = 0;
188 }
189
190 _source_term_configs.emplace_back(std::move(st_config), st_mesh,
191 component_id);
192 }
193 }
194 else
195 {
196 INFO("No source terms for process variable '{:s}' found.", _name);
197 }
198
199 if (!_deactivated_subdomains.empty())
200 {
201 _is_active = getOrCreateMeshProperty<unsigned char>(
202 _mesh, _name + "_active", MeshLib::MeshItemType::Cell, 1);
203 std::fill(std::begin(*_is_active), std::end(*_is_active), 1u);
204 }
205}
206
208
209std::string const& ProcessVariable::getName() const
210{
211 return _name;
212}
213
215{
216 return _mesh;
217}
218
219std::vector<std::unique_ptr<BoundaryCondition>>
221 const NumLib::LocalToGlobalIndexMap& dof_table,
222 const int variable_id,
223 unsigned const integration_order,
224 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
225 Process const& process,
226 std::vector<std::reference_wrapper<ProcessVariable>> const&
227 all_process_variables_for_this_process)
228{
229 std::vector<std::unique_ptr<BoundaryCondition>> bcs;
230 bcs.reserve(_bc_configs.size());
231
232 for (auto& config : _bc_configs)
233 {
234 auto bc = createBoundaryCondition(
235 config, dof_table, _mesh, variable_id, integration_order,
236 _shapefunction_order, parameters, process,
237 all_process_variables_for_this_process);
238#ifdef USE_PETSC
239 if (bc == nullptr)
240 {
241 continue;
242 }
243#endif // USE_PETSC
244 bcs.push_back(std::move(bc));
245 }
246
248 parameters, bcs);
249
250 return bcs;
251}
252
254 const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
255 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
256 std::vector<std::unique_ptr<BoundaryCondition>>& bcs)
257{
258 for (auto const& deactivated_subdomain : _deactivated_subdomains)
259 {
260 auto const& deactivated_subdomain_mesh =
261 deactivated_subdomain.deactivated_subdomain_mesh;
262 auto const* parameter = &ParameterLib::findParameter<double>(
264 bool const set_outer_nodes_dirichlet_values =
265 deactivated_subdomain.boundary_value_parameter != nullptr;
266 if (set_outer_nodes_dirichlet_values)
267 {
268 parameter = deactivated_subdomain.boundary_value_parameter;
269 }
270
271 for (int component_id = 0;
272 component_id <
273 dof_table.getNumberOfVariableComponents(variable_id);
274 component_id++)
275 {
276 auto bc = std::make_unique<DeactivatedSubdomainDirichlet>(
277 *_is_active, deactivated_subdomain.time_interval, *parameter,
278 set_outer_nodes_dirichlet_values, deactivated_subdomain_mesh,
279 dof_table, variable_id, component_id);
280
281#ifdef USE_PETSC
282 // TODO: make it work under PETSc too.
283 if (bc == nullptr)
284 {
285 continue;
286 }
287#endif // USE_PETSC
288 bcs.push_back(std::move(bc));
289 }
290 }
291}
292
294{
295 if (_deactivated_subdomains.empty())
296 {
297 return;
298 }
299
301
302 // If none of the deactivated subdomains is active at current time, then the
303 // _ids_of_active_elements remain empty.
304 if (std::none_of(
306 [&](auto const& ds) { return ds.isInTimeSupportInterval(time); }))
307 {
308 // Also mark all of the elements as active.
309 assert(_is_active != nullptr); // guaranteed by constructor
310 std::fill(std::begin(*_is_active), std::end(*_is_active), 1u);
311
312 return;
313 }
314
315 auto is_active_in_subdomain = [&](std::size_t const element_id,
316 DeactivatedSubdomain const& ds) -> bool
317 {
318 return (!ds.isInTimeSupportInterval(time)) ||
319 !ds.isDeactivated(*_mesh.getElement(element_id), time);
320 };
321
322 auto is_active_in_all_subdomains = [&](std::size_t const element_id) -> bool
323 {
324 return std::all_of(begin(_deactivated_subdomains),
326 [&](auto const& ds)
327 { return is_active_in_subdomain(element_id, ds); });
328 };
329
330 auto const number_of_elements = _mesh.getNumberOfElements();
331 for (std::size_t element_id = 0; element_id < number_of_elements;
332 element_id++)
333 {
334 if (is_active_in_all_subdomains(element_id))
335 {
336 _ids_of_active_elements.push_back(element_id);
337 }
338 }
339
340 // all elements are deactivated
341 std::fill(std::begin(*_is_active), std::end(*_is_active), 0u);
342
343 for (auto const id : _ids_of_active_elements)
344 {
345 (*_is_active)[id] = 1u;
346 }
347}
348
349std::vector<std::unique_ptr<SourceTerm>> ProcessVariable::createSourceTerms(
350 const NumLib::LocalToGlobalIndexMap& dof_table,
351 const int variable_id,
352 unsigned const integration_order,
353 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
354 std::vector<std::reference_wrapper<ProcessVariable>> const&
355 all_process_variables_for_this_process)
356{
357 std::vector<std::unique_ptr<SourceTerm>> source_terms;
358
359 transform(cbegin(_source_term_configs), cend(_source_term_configs),
360 back_inserter(source_terms),
361 [&](auto const& config)
362 {
363 return createSourceTerm(
364 config, dof_table, config.mesh, variable_id,
365 integration_order, _shapefunction_order, parameters,
366 all_process_variables_for_this_process);
367 });
368
369 return source_terms;
370}
371
373
374} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:35
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
Definition: ConfigTree.cpp:160
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:89
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:92
int getNumberOfVariableComponents(int variable_id) const
void updateDeactivatedSubdomains(double const time)
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< std::unique_ptr< SourceTerm > > 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)
std::vector< DeactivatedSubdomain > _deactivated_subdomains
std::string const & getName() const
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)
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::iterator_traits< InputIt >::reference findElementOrError(InputIt begin, InputIt end, Predicate predicate, std::string const &error="")
Definition: Algorithm.h:69
std::string meshNameFromGeometry(std::string const &geometrical_set_name, std::string const &geometry_name)
bool parseCompensateNonEquilibriumInitialResiduum(BaseLib::ConfigTree const &config)
std::unique_ptr< SourceTerm > 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)
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)
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