OGS
CreateDeactivatedSubdomain.cpp
Go to the documentation of this file.
1
10
11#include <Eigen/Core>
12#include <range/v3/algorithm/all_of.hpp>
13#include <range/v3/algorithm/copy_if.hpp>
14#include <range/v3/algorithm/partition_copy.hpp>
15#include <range/v3/range/conversion.hpp>
16
17#include "BaseLib/ConfigTree.h"
18#include "BaseLib/Error.h"
21#include "MeshLib/Mesh.h"
22#include "MeshLib/Node.h"
26#include "ParameterLib/Utils.h"
27
28namespace ProcessLib
29{
30template <typename IsActive>
31static std::pair<std::vector<std::size_t>, std::vector<std::size_t>>
33 MeshLib::Mesh const& sub_mesh,
34 IsActive is_active)
35{
36 auto* const bulk_node_ids = MeshLib::bulkNodeIDs(sub_mesh);
37 if (bulk_node_ids == nullptr)
38 {
40 "Bulk node ids map is not available in the deactivate subdomain "
41 "mesh.");
42 }
43
44 std::vector<std::size_t> inner_nodes;
45 // Reserve for more than needed, but not much, because almost all nodes will
46 // be the inner nodes.
47 inner_nodes.reserve(sub_mesh.getNumberOfNodes());
48
49 std::vector<std::size_t> outer_nodes;
50
51 ranges::partition_copy(
52 sub_mesh.getNodes() | MeshLib::views::ids,
53 back_inserter(inner_nodes),
54 back_inserter(outer_nodes),
55 [&](std::size_t const n)
56 {
57 auto const bulk_node = mesh.getNode((*bulk_node_ids)[n]);
58 const auto& connected_elements =
59 mesh.getElementsConnectedToNode(*bulk_node);
60 // Check whether this node is connected only to an active
61 // elements. Then it is an inner node, and outer otherwise.
62 return ranges::all_of(connected_elements | MeshLib::views::ids,
63 is_active);
64 });
65
66 return {std::move(inner_nodes), std::move(outer_nodes)};
67}
68
69static std::vector<std::vector<std::size_t>> extractElementsAlongOuterNodes(
70 MeshLib::Mesh const& mesh, MeshLib::Mesh const& sub_mesh,
71 std::vector<std::size_t> const& outer_nodes)
72{
73 auto const& bulk_node_ids =
74 *sub_mesh.getProperties().template getPropertyVector<std::size_t>(
77
78 auto to_bulk_node_id =
79 ranges::views::transform([&bulk_node_ids](std::size_t const node_id)
80 { return bulk_node_ids[node_id]; });
81 auto connected_elements = ranges::views::transform(
82 [&mesh](std::size_t const node_id)
83 {
84 return mesh.getElementsConnectedToNode(node_id) |
85 MeshLib::views::ids | ranges::to<std::vector>;
86 });
87
88 return outer_nodes | to_bulk_node_id | connected_elements |
89 ranges::to<std::vector>;
90}
91
93 MeshLib::Mesh const& mesh,
94 std::unordered_set<int> const& deactivated_subdomain_material_ids)
95{
96 // An element is active if its material id matches one of the deactivated
97 // subdomain material ids.
98 auto is_active =
99 [&deactivated_subdomain_material_ids,
100 &material_ids = *materialIDs(mesh)](std::size_t element_id)
101 {
102 return deactivated_subdomain_material_ids.contains(
103 material_ids[element_id]);
104 };
105
106 auto const& elements = mesh.getElements();
107 std::vector<MeshLib::Element*> deactivated_elements;
108 ranges::copy_if(elements, back_inserter(deactivated_elements), is_active,
109 [](auto const e) { return e->getID(); });
110
111 auto bulk_element_ids = deactivated_elements | MeshLib::views::ids |
112 ranges::to<std::unordered_set>;
113
114 static int mesh_number = 0;
115 // Subdomain mesh consisting of deactivated elements.
117 "deactivate_subdomain_" + std::to_string(mesh_number++),
118 MeshLib::cloneElements(deactivated_elements));
119
120 auto [inner_nodes, outer_nodes] =
121 extractInnerAndOuterNodes(mesh, *sub_mesh, is_active);
122
123 auto outer_nodes_elements =
124 extractElementsAlongOuterNodes(mesh, *sub_mesh, outer_nodes);
125
126 return {std::move(*sub_mesh), std::move(bulk_element_ids),
127 std::move(inner_nodes), std::move(outer_nodes),
128 std::move(outer_nodes_elements)};
129}
130
132 std::optional<BaseLib::ConfigTree> const& time_interval_config,
133 std::optional<std::string> const& curve_name,
134 std::map<std::string,
135 std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
136 curves)
137{
138 // Check for the error case first: only one of the configs must be used.
139 if (time_interval_config && curve_name)
140 {
141 OGS_FATAL(
142 "In the deactivate subdomain either a time interval or a curve "
143 "must be given, not both.");
144 }
145
146 // Parse time interval and return a curve
147 if (time_interval_config)
148 {
149 DBUG("Constructing time interval");
150 auto const start_time =
152 time_interval_config->getConfigParameter<double>("start");
153
154 auto const end_time =
156 time_interval_config->getConfigParameter<double>("end");
157
158 // Using very large value for the curve's value, s.t. for any time from
159 // start to end the whole subdomain is deactivated at once.
160 return {{start_time, end_time},
161 {std::numeric_limits<double>::max(),
162 std::numeric_limits<double>::max()},
163 false};
164 }
165
166 // Try to find the curve.
167 if (curve_name)
168 {
169 DBUG("Using curve '{:s}'", *curve_name);
170 // Return a copy of the curve because the time interval in the other
171 // branch returns a temporary.
172 return *BaseLib::getOrError(curves, *curve_name,
173 "Could not find curve.");
174 }
175
176 // If we got so far, there is an error: one of the configs must be
177 // available.
178 OGS_FATAL(
179 "In the deactivate subdomain neither a time interval nor a curve are "
180 "given. One of them must be specified.");
181}
182
184static std::pair<Eigen::Vector3d, Eigen::Vector3d> parseLineSegment(
185 BaseLib::ConfigTree const& config)
186{
187 DBUG("Constructing line segment");
188 auto const start =
190 config.getConfigParameter<std::vector<double>>("start");
191 if (start.size() != 3)
192 {
193 OGS_FATAL(
194 "For construction of a line segment the start point must be a 3D "
195 "point. Got a vector of size {}.",
196 start.size());
197 }
198
199 auto const end =
201 config.getConfigParameter<std::vector<double>>("end");
202
203 if (end.size() != 3)
204 {
205 OGS_FATAL(
206 "For construction of a line segment the end point must be a 3D "
207 "point. Got a vector of size {}.",
208 end.size());
209 }
210 return {Eigen::Vector3d{start[0], start[1], start[2]},
211 Eigen::Vector3d{end[0], end[1], end[2]}};
212}
213
216{
217 DBUG("Constructing a ball");
218 auto const center =
220 config.getConfigParameter<std::vector<double>>("center");
221 if (center.size() != 3)
222 {
223 OGS_FATAL(
224 "For construction of the center of a ball must be a 3D point. Got "
225 "a vector of size {}.",
226 center.size());
227 }
228
229 auto const radius =
231 config.getConfigParameter<double>("radius");
232
233 return {Eigen::Vector3d{center[0], center[1], center[2]}, radius};
234}
235
237 BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh,
238 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
239 std::map<std::string,
240 std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
241 curves)
242{
243 auto const& time_interval_config =
245 config.getConfigSubtreeOptional("time_interval");
246
247 auto const& curve_name =
249 config.getConfigParameterOptional<std::string>("time_curve");
250 auto time_interval =
251 parseTimeIntervalOrCurve(time_interval_config, curve_name, curves);
252
253 auto const line_segment_config =
255 config.getConfigSubtreeOptional("line_segment");
256
257 auto const ball_config =
259 config.getConfigSubtreeOptional("ball");
260
261 if (line_segment_config && ball_config)
262 {
263 OGS_FATAL(
264 "line_segment and ball can not be defined "
265 "simultaneously.");
266 }
267
268 if (time_interval_config && (line_segment_config || ball_config))
269 {
270 OGS_FATAL(
271 "When using time interval for subdomain deactivation a line "
272 "segment or a ball must not be specified.");
273 }
274
275 if (curve_name && !(line_segment_config || ball_config))
276 {
277 OGS_FATAL(
278 "When using curve for subdomain deactivation a line segment or a "
279 "ball must be specified.");
280 }
281
282 // If time interval was specified then an empty optional line segment is
283 // used *internally* because the whole selected material ids subdomain will
284 // be deactivated.
285 std::optional<std::pair<Eigen::Vector3d, Eigen::Vector3d>> line_segment;
286 if (line_segment_config)
287 {
288 line_segment = parseLineSegment(*line_segment_config);
289 }
290
291 std::optional<detail::Ball> ball;
292 if (ball_config)
293 {
294 ball = parseBall(*ball_config);
295 }
296
297 ParameterLib::Parameter<double>* boundary_value_parameter = nullptr;
298 auto boundary_value_parameter_name =
300 config.getConfigParameterOptional<std::string>("boundary_parameter");
301 if (boundary_value_parameter_name)
302 {
303 DBUG("Using parameter {:s}", *boundary_value_parameter_name);
304 boundary_value_parameter = &ParameterLib::findParameter<double>(
305 *boundary_value_parameter_name, parameters, 1, &mesh);
306 }
307
308 auto deactivated_subdomain_material_ids =
310 config.getConfigParameter("material_ids", std::vector<int>{}) |
311 ranges::to<std::unordered_set>();
312
313 if (deactivated_subdomain_material_ids.empty())
314 {
315 OGS_FATAL(
316 "The material IDs of the deactivated subdomains are not given. The "
317 "program terminates now.");
318 }
319
320 if (materialIDs(mesh) == nullptr)
321 {
322 OGS_FATAL(
323 "The mesh doesn't contain materialIDs for subdomain deactivation. "
324 "The program terminates now.");
325 }
326
327 auto deactivated_subdomain_mesh = createDeactivatedSubdomainMesh(
328 mesh, deactivated_subdomain_material_ids);
329
330 return {std::move(time_interval), line_segment, ball,
331 std::move(deactivated_subdomain_mesh), boundary_value_parameter};
332}
333
334std::vector<DeactivatedSubdomain> createDeactivatedSubdomains(
335 BaseLib::ConfigTree const& config,
336 MeshLib::Mesh const& mesh,
337 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
338 std::map<std::string,
339 std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
340 curves)
341{
342 std::vector<DeactivatedSubdomain> deactivated_subdomains;
343 // Deactivated subdomains
344 if (auto subdomains_config =
346 config.getConfigSubtreeOptional("deactivated_subdomains"))
347 {
348 INFO("There are subdomains being deactivated.");
349
350 auto const deactivated_subdomain_configs =
352 subdomains_config->getConfigSubtreeList("deactivated_subdomain");
353 std::transform(std::begin(deactivated_subdomain_configs),
354 std::end(deactivated_subdomain_configs),
355 std::back_inserter(deactivated_subdomains),
356 [&](auto const& config) {
358 config, mesh, parameters, curves);
359 });
360 }
361 return deactivated_subdomains;
362}
363
364} // namespace ProcessLib
Definition of Duplicate functions.
Definition of the Element class.
#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
std::optional< T > getConfigParameterOptional(std::string const &param) const
T getConfigParameter(std::string const &param) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition Mesh.h:93
Properties & getProperties()
Definition Mesh.h:136
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:102
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:256
OGS_NO_DANGLING Map::mapped_type & getOrError(Map &map, Key const &key, std::string const &error_message)
Definition Algorithm.h:118
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
std::unique_ptr< MeshLib::Mesh > createMeshFromElementSelection(std::string mesh_name, std::vector< MeshLib::Element * > const &elements)
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition Properties.h:185
std::vector< Element * > cloneElements(std::vector< Element * > const &elements)
Clones a vector of elements using the Element::clone() function.
PropertyVector< std::size_t > const * bulkNodeIDs(Mesh const &mesh)
Definition Mesh.cpp:292
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
static DeactivatedSubdomainMesh createDeactivatedSubdomainMesh(MeshLib::Mesh const &mesh, std::unordered_set< int > const &deactivated_subdomain_material_ids)
static detail::Ball parseBall(BaseLib::ConfigTree const &config)
Returns a ball represented by its center and its radius.
static std::pair< std::vector< std::size_t >, std::vector< std::size_t > > extractInnerAndOuterNodes(MeshLib::Mesh const &mesh, MeshLib::Mesh const &sub_mesh, IsActive is_active)
DeactivatedSubdomain createDeactivatedSubdomain(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::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)
static std::pair< Eigen::Vector3d, Eigen::Vector3d > parseLineSegment(BaseLib::ConfigTree const &config)
Returns a line segment represented by its begin and end points.
static std::vector< std::vector< std::size_t > > extractElementsAlongOuterNodes(MeshLib::Mesh const &mesh, MeshLib::Mesh const &sub_mesh, std::vector< std::size_t > const &outer_nodes)
static MathLib::PiecewiseLinearInterpolation parseTimeIntervalOrCurve(std::optional< BaseLib::ConfigTree > const &time_interval_config, std::optional< std::string > const &curve_name, std::map< std::string, std::unique_ptr< MathLib::PiecewiseLinearInterpolation > > const &curves)