OGS
AddProcessDataToMesh.cpp
Go to the documentation of this file.
1
12
13#include <span>
14
15#include "InfoLib/GitInfo.h"
21#ifdef USE_PETSC
23#endif
25
28static void addOgsVersion(MeshLib::Mesh& mesh)
29{
30 auto& ogs_version_field = *MeshLib::getOrCreateMeshProperty<char>(
33
34 ogs_version_field.assign(GitInfoLib::GitInfo::ogs_version);
35}
36
38 NumLib::Time const& t,
39 std::vector<GlobalVector*> const& x,
40 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
42 std::string const& output_name,
43 MeshLib::Mesh& mesh)
44{
45 DBUG(" secondary variable {:s}", output_name);
46
47 auto& nodal_values_mesh = *MeshLib::getOrCreateMeshProperty<double>(
48 mesh, output_name, MeshLib::MeshItemType::Node,
50 if (nodal_values_mesh.size() !=
52 {
54 "Nodal property `{:s}' does not have the right number of "
55 "components. Expected: {:d}, actual: {:d}",
56 output_name,
58 nodal_values_mesh.size());
59 }
60
61 std::unique_ptr<GlobalVector> result_cache;
62 auto const& nodal_values =
63 var.fcts.eval_field(t(), x, dof_tables, result_cache);
64
65#ifdef USE_PETSC
66 std::size_t const global_vector_size =
67 nodal_values.getLocalSize() + nodal_values.getGhostSize();
68#else
69 std::size_t const global_vector_size = nodal_values.size();
70#endif
71 if (nodal_values_mesh.size() != global_vector_size)
72 {
74 "Secondary variable `{:s}' did not evaluate to the right number of "
75 "components. Expected: {:d}, actual: {:d}.",
76 var.name, nodal_values_mesh.size(), global_vector_size);
77 }
78
79 // Copy result
80 nodal_values.copyValues(std::span{nodal_values_mesh});
81}
82
84 NumLib::Time const& t,
85 std::vector<GlobalVector*> const& x,
86 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
88 std::string const& output_name,
89 MeshLib::Mesh& mesh)
90{
91 if (!var.fcts.eval_residuals)
92 {
93 return;
94 }
95
96 DBUG(" secondary variable {:s} residual", output_name);
97 auto const& property_name_res = output_name + "_residual";
98
99 auto& residuals_mesh = *MeshLib::getOrCreateMeshProperty<double>(
100 mesh, property_name_res, MeshLib::MeshItemType::Cell,
101 var.fcts.num_components);
102 if (residuals_mesh.size() !=
104 {
105 OGS_FATAL(
106 "Cell property `{:s}' does not have the right number of "
107 "components. Expected: {:d}, actual: {:d}",
108 property_name_res,
110 residuals_mesh.size());
111 }
112
113 std::unique_ptr<GlobalVector> result_cache;
114 auto const& residuals =
115 var.fcts.eval_residuals(t(), x, dof_table, result_cache);
116#ifdef USE_PETSC
117 std::size_t const global_vector_size =
118 residuals.getLocalSize() + residuals.getGhostSize();
119#else
120 std::size_t const global_vector_size = residuals.size();
121#endif
122 if (residuals_mesh.size() != global_vector_size)
123 {
124 OGS_FATAL(
125 "The residual of secondary variable `{:s}' did not evaluate to the "
126 "right number of components. Expected: {:d}, actual: {:d}.",
127 var.name, residuals_mesh.size(), global_vector_size);
128 }
129
130 // Copy result
131 residuals.copyValues(std::span{residuals_mesh});
132}
133
134static std::vector<double> copySolutionVector(GlobalVector const& x)
135{
136 std::vector<double> x_copy;
137 x.copyValues(x_copy);
138 return x_copy;
139}
140
142 [[maybe_unused]] MeshLib::Mesh const& mesh,
143 [[maybe_unused]] NumLib::LocalToGlobalIndexMap const& dof_table,
144 [[maybe_unused]] NumLib::LocalToGlobalIndexMap const& bulk_mesh_dof_table)
145{
146#ifdef USE_PETSC
147
148 if (&bulk_mesh_dof_table != &dof_table)
149 {
150 auto const bulk_id_string =
152 if (!mesh.getProperties().existsPropertyVector<std::size_t>(
153 bulk_id_string))
154 {
155 OGS_FATAL(
156 "The required bulk node ids map does not exist in "
157 "the boundary mesh '{:s}' or has the wrong data "
158 "type (should be equivalent to C++ data type "
159 "std::size_t which is an unsigned integer of size "
160 "{:d} or UInt64 in vtk terminology).",
161 mesh.getName(), sizeof(std::size_t));
162 }
163 return mesh.getProperties().getPropertyVector<std::size_t>(
164 bulk_id_string);
165 }
166#endif
167
168 return nullptr;
169}
170
172 std::size_t const mesh_id, std::size_t const node_id,
173 [[maybe_unused]] bool const is_ghost_node, int const global_component_id,
174 GlobalVector const& x, NumLib::LocalToGlobalIndexMap const& dof_table,
175 [[maybe_unused]] NumLib::LocalToGlobalIndexMap const& bulk_mesh_dof_table,
176 [[maybe_unused]] MeshLib::PropertyVector<std::size_t> const* const
177 bulk_node_id_map)
178{
179#ifdef USE_PETSC
180 if (is_ghost_node && &bulk_mesh_dof_table != &dof_table)
181 {
182 auto const bulk_node_id = (*bulk_node_id_map)[node_id];
183 std::size_t const bulk_mesh_id = 0;
184
185 MeshLib::Location const loc(bulk_mesh_id, MeshLib::MeshItemType::Node,
186 bulk_node_id);
187
188 // early return!
189 return bulk_mesh_dof_table.getLocalIndex(
190 loc, global_component_id, x.getRangeBegin(), x.getRangeEnd());
191 }
192#endif
193
194 MeshLib::Location const loc(mesh_id, MeshLib::MeshItemType::Node, node_id);
195
196 return dof_table.getLocalIndex(loc, global_component_id, x.getRangeBegin(),
197 x.getRangeEnd());
198}
199
200static bool isGhostNode([[maybe_unused]] MeshLib::Mesh const& mesh,
201 [[maybe_unused]] std::size_t const node_id)
202{
203#ifndef USE_PETSC
204 return false;
205#else
206 return static_cast<MeshLib::NodePartitionedMesh const&>(mesh).isGhostNode(
207 node_id);
208#endif
209}
210
240static std::set<std::string> addPrimaryVariablesToMesh(
241 MeshLib::Mesh& mesh,
242 GlobalVector const& x,
243 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>> const&
244 process_variables,
245 std::set<std::string> const& output_variables,
246 NumLib::LocalToGlobalIndexMap const& dof_table,
247 NumLib::LocalToGlobalIndexMap const& bulk_mesh_dof_table)
248{
249 if (dof_table.getNumberOfVariables() !=
250 bulk_mesh_dof_table.getNumberOfVariables() ||
251 dof_table.getNumberOfGlobalComponents() !=
252 bulk_mesh_dof_table.getNumberOfGlobalComponents())
253 {
254 OGS_FATAL(
255 "The d.o.f. table for the passed mesh must have the same number of "
256 "variables and global components as the d.o.f. table for the full "
257 "simulation domain. But the values differ: {} != {} (variables) or "
258 "{} != {} (global components).",
259 dof_table.getNumberOfVariables(),
260 bulk_mesh_dof_table.getNumberOfVariables(),
261 dof_table.getNumberOfGlobalComponents(),
262 bulk_mesh_dof_table.getNumberOfGlobalComponents());
263 }
264
265 auto const number_of_dof_variables = dof_table.getNumberOfVariables();
266
267 if (number_of_dof_variables != static_cast<int>(process_variables.size()))
268 {
269 OGS_FATAL(
270 "The number of variables in the d.o.f. table differs from the "
271 "number of primary variables of the process {} != {}.",
272 number_of_dof_variables, process_variables.size());
273 }
274
275 auto const x_copy = copySolutionVector(x);
276 std::set<std::string> names_of_already_output_variables;
277
278 int global_component_offset_next = 0;
279
280 auto const* const bulk_node_id_map = getBulkNodeIdMapForPetscIfNecessary(
281 mesh, dof_table, bulk_mesh_dof_table);
282
283 for (int variable_id = 0; variable_id < number_of_dof_variables;
284 ++variable_id)
285 {
286 auto const& pv = process_variables[variable_id].get();
287 auto const n_components = pv.getNumberOfGlobalComponents();
288
289 // increase global component offset even if we do not add anything in
290 // this iteration
291 int global_component_offset = global_component_offset_next;
292 global_component_offset_next += n_components;
293
294 if (!output_variables.empty() &&
295 !output_variables.contains(pv.getName()))
296 {
297 continue;
298 }
299
300 names_of_already_output_variables.insert(pv.getName());
301
302 DBUG(" process variable {:s}", pv.getName());
303
304 auto& output_data = *MeshLib::getOrCreateMeshProperty<double>(
305 mesh, pv.getName(), MeshLib::MeshItemType::Node, n_components);
306
307 // mesh subsets are the same for all components
308 int const dummy_component_id = 0;
309 auto const& mesh_subset =
310 dof_table.getMeshSubset(variable_id, dummy_component_id);
311 auto const mesh_id = mesh_subset.getMeshID();
312
313 for (auto const* node : mesh_subset.getNodes())
314 {
315 auto const node_id = node->getID();
316 auto const is_ghost_node = isGhostNode(mesh, node_id);
317
318 for (int component_id = 0; component_id < n_components;
319 ++component_id)
320 {
321 auto const global_component_id =
322 global_component_offset + component_id;
323
324 auto const in_index = getIndexForComponentInSolutionVector(
325 mesh_id, node_id, is_ghost_node, global_component_id, x,
326 dof_table, bulk_mesh_dof_table, bulk_node_id_map);
327
328 // per node ordering of components
329 auto const out_index = node_id * n_components + component_id;
330
331 // request for index of linear quantities at higher order nodes
332 // results in returning NumLib::MeshComponentMap::nop
333 if (in_index == NumLib::MeshComponentMap::nop)
334 {
335 output_data[out_index] = 0;
336 continue;
337 }
338
339 output_data[out_index] = x_copy[in_index];
340 }
341 }
342 }
343
344 return names_of_already_output_variables;
345}
346
348 ProcessLib::SecondaryVariableCollection const& secondary_variables,
349 std::set<std::string>& names_of_already_output_variables,
350 NumLib::Time const& t, std::vector<GlobalVector*> const& xs,
351 MeshLib::Mesh& mesh,
352 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
353 bool const output_residuals)
354{
355 for (auto const& external_variable_name : secondary_variables)
356 {
357 auto const& name = external_variable_name.first;
358 if (!names_of_already_output_variables.insert(name).second)
359 {
360 // no insertion took place, output already done
361 continue;
362 }
363
364 addSecondaryVariableNodes(t, xs, dof_tables,
365 secondary_variables.get(name), name, mesh);
366
367 if (output_residuals)
368 {
370 t, xs, dof_tables, secondary_variables.get(name), name, mesh);
371 }
372 }
373}
374
375namespace ProcessLib
376{
378 std::vector<GlobalVector*> const& xs,
379 int const process_id,
380 ProcessOutputData const& process_output_data,
381 bool const output_secondary_variables,
382 OutputDataSpecification const& process_output)
383{
384 DBUG("Process output data.");
385
386 auto const& pod = process_output_data;
387 auto const& process_variables = pod.getProcessVariables(process_id);
388 auto const& secondary_variables = pod.getSecondaryVariables();
389 auto const* const integration_point_writers =
390 pod.getIntegrationPointWriters();
391 auto const& bulk_mesh_dof_table = pod.getBulkMeshDofTable(process_id);
392 auto const& output_mesh_dof_table = pod.getOutputMeshDofTable(process_id);
393 auto& output_mesh = pod.getOutputMesh();
394
395 auto const& output_variables = process_output.output_variables;
396
397 addOgsVersion(output_mesh);
398
399 auto names_of_already_output_variables = addPrimaryVariablesToMesh(
400 output_mesh, *xs[process_id], process_variables, output_variables,
401 output_mesh_dof_table, bulk_mesh_dof_table);
402
403 if (output_secondary_variables)
404 {
405 auto const& output_mesh_dof_tables =
406 pod.getOutputMeshDofTablesOfAllProcesses();
407
408 addSecondaryVariablesToMesh(secondary_variables,
409 names_of_already_output_variables, t, xs,
410 output_mesh, output_mesh_dof_tables,
411 process_output.output_residuals);
412 }
413
414 if (integration_point_writers)
415 {
416 addIntegrationPointDataToMesh(output_mesh, *integration_point_writers);
417 }
418}
419} // namespace ProcessLib
static void addSecondaryVariableNodes(NumLib::Time const &t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, ProcessLib::SecondaryVariable const &var, std::string const &output_name, MeshLib::Mesh &mesh)
static std::vector< double > copySolutionVector(GlobalVector const &x)
static GlobalIndexType getIndexForComponentInSolutionVector(std::size_t const mesh_id, std::size_t const node_id, bool const is_ghost_node, int const global_component_id, GlobalVector const &x, NumLib::LocalToGlobalIndexMap const &dof_table, NumLib::LocalToGlobalIndexMap const &bulk_mesh_dof_table, MeshLib::PropertyVector< std::size_t > const *const bulk_node_id_map)
MeshLib::PropertyVector< std::size_t > const * getBulkNodeIdMapForPetscIfNecessary(MeshLib::Mesh const &mesh, NumLib::LocalToGlobalIndexMap const &dof_table, NumLib::LocalToGlobalIndexMap const &bulk_mesh_dof_table)
static std::set< std::string > addPrimaryVariablesToMesh(MeshLib::Mesh &mesh, GlobalVector const &x, std::vector< std::reference_wrapper< ProcessLib::ProcessVariable > > const &process_variables, std::set< std::string > const &output_variables, NumLib::LocalToGlobalIndexMap const &dof_table, NumLib::LocalToGlobalIndexMap const &bulk_mesh_dof_table)
static void addSecondaryVariablesToMesh(ProcessLib::SecondaryVariableCollection const &secondary_variables, std::set< std::string > &names_of_already_output_variables, NumLib::Time const &t, std::vector< GlobalVector * > const &xs, MeshLib::Mesh &mesh, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, bool const output_residuals)
static void addSecondaryVariableResiduals(NumLib::Time const &t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, ProcessLib::SecondaryVariable const &var, std::string const &output_name, MeshLib::Mesh &mesh)
static bool isGhostNode(MeshLib::Mesh const &mesh, std::size_t const node_id)
static void addOgsVersion(MeshLib::Mesh &mesh)
#define OGS_FATAL(...)
Definition Error.h:26
Git information.
GlobalMatrix::IndexType GlobalIndexType
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of mesh class for partitioned mesh (by node) for parallel computing within the framework o...
Global vector based on Eigen vector.
Definition EigenVector.h:26
IndexType getRangeEnd() const
return an end index of the active data range
Definition EigenVector.h:49
void copyValues(std::vector< double > &u) const
static constexpr IndexType getRangeBegin()
return a start index of the active data range
Definition EigenVector.h:46
std::size_t getMeshID() const
return this mesh ID
Definition MeshSubset.h:76
const std::string getName() const
Get name of the mesh.
Definition Mesh.h:105
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:99
MeshLib::MeshSubset const & getMeshSubset(int const variable_id, int const component_id) const
GlobalIndexType getLocalIndex(MeshLib::Location const &l, std::size_t const comp_id, std::size_t const range_begin, std::size_t const range_end) const
static constexpr NUMLIB_EXPORT GlobalIndexType const nop
Holds all data of a process that are needed for output.
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(int const process_id) const
Handles configuration of several secondary variables from the project file.
SecondaryVariable const & get(std::string const &external_name) const
Returns the secondary variable with the given external name.
const std::string OGS_VERSION
Definition GitInfo.cpp:20
GITINFOLIB_EXPORT const std::string ogs_version
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition Properties.h:191
void addProcessDataToMesh(NumLib::Time const &t, std::vector< GlobalVector * > const &xs, int const process_id, ProcessOutputData const &process_output_data, bool const output_secondary_variables, OutputDataSpecification const &process_output)
Holds information about which variables to write to output files.
bool output_residuals
Tells if also to output extrapolation residuals.
std::set< std::string > output_variables
All variables that shall be output.
Function const eval_field
Computes the value of the field at every node of the underlying mesh.
const unsigned num_components
Number of components of the variable.
Stores information about a specific secondary variable.
SecondaryVariableFunctions fcts
Functions used for computing the secondary variable.
std::string const name
Name of the variable; used, e.g., for output.