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"
19 #include "DeactivatedSubdomain.h"
21 #include "MeshLib/Mesh.h"
23 #include "MeshLib/Node.h"
24 #include "ParameterLib/Parameter.h"
25 #include "ParameterLib/Utils.h"
26 
27 namespace ProcessLib
28 {
29 template <typename IsActive>
30 static std::pair<std::vector<std::size_t>, std::vector<std::size_t>>
32  MeshLib::Mesh const& sub_mesh,
33  IsActive is_active)
34 {
35  auto* const bulk_node_ids =
36  sub_mesh.getProperties().template getPropertyVector<std::size_t>(
37  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
38  if (bulk_node_ids == nullptr)
39  {
40  OGS_FATAL(
41  "Bulk node ids map is not available in the deactivate subdomain "
42  "mesh.");
43  }
44 
45  std::vector<std::size_t> inner_nodes;
46  // Reserve for more than needed, but not much, because almost all nodes will
47  // be the inner nodes.
48  inner_nodes.reserve(sub_mesh.getNumberOfNodes());
49 
50  std::vector<std::size_t> outer_nodes;
51 
52  ranges::partition_copy(
53  sub_mesh.getNodes() | MeshLib::views::ids,
54  back_inserter(inner_nodes),
55  back_inserter(outer_nodes),
56  [&](std::size_t const n)
57  {
58  auto const bulk_node = mesh.getNode((*bulk_node_ids)[n]);
59  const auto& connected_elements =
60  mesh.getElementsConnectedToNode(*bulk_node);
61  // Check whether this node is connected only to an active
62  // elements. Then it is an inner node, and outer otherwise.
63  return ranges::all_of(connected_elements | MeshLib::views::ids,
64  is_active);
65  });
66 
67  return {std::move(inner_nodes), std::move(outer_nodes)};
68 }
69 
70 static std::vector<std::vector<std::size_t>> extractElementsAlongOuterNodes(
71  MeshLib::Mesh const& mesh, MeshLib::Mesh const& sub_mesh,
72  std::vector<std::size_t> const& outer_nodes)
73 {
74  auto const& bulk_node_ids =
75  *sub_mesh.getProperties().template getPropertyVector<std::size_t>(
76  "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
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::vector<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 = [deactivated_subdomain_material_ids,
99  material_ids = *materialIDs(mesh)](std::size_t element_id)
100  {
101  int const element_material_id = material_ids[element_id];
102  return std::any_of(begin(deactivated_subdomain_material_ids),
103  end(deactivated_subdomain_material_ids),
104  [&](int const material_id)
105  { return material_id == element_material_id; });
106  };
107 
108  auto const& elements = mesh.getElements();
109  std::vector<MeshLib::Element*> deactivated_elements;
110  ranges::copy_if(elements, back_inserter(deactivated_elements), is_active,
111  [](auto const e) { return e->getID(); });
112 
113  auto bulk_element_ids =
114  deactivated_elements | MeshLib::views::ids | ranges::to<std::vector>;
115 
116  static int mesh_number = 0;
117  // Subdomain mesh consisting of deactivated elements.
119  "deactivate_subdomain_" + std::to_string(mesh_number++),
120  MeshLib::cloneElements(deactivated_elements));
121 
122  auto [inner_nodes, outer_nodes] =
123  extractInnerAndOuterNodes(mesh, *sub_mesh, is_active);
124 
125  auto outer_nodes_elements =
126  extractElementsAlongOuterNodes(mesh, *sub_mesh, outer_nodes);
127 
128  return {std::move(*sub_mesh), std::move(bulk_element_ids),
129  std::move(inner_nodes), std::move(outer_nodes),
130  std::move(outer_nodes_elements)};
131 }
132 
134  std::optional<BaseLib::ConfigTree> const& time_interval_config,
135  std::optional<std::string> const& curve_name,
136  std::map<std::string,
137  std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
138  curves)
139 {
140  // Check for the error case first: only one of the configs must be used.
141  if (time_interval_config && curve_name)
142  {
143  OGS_FATAL(
144  "In the deactivate subdomain either a time interval or a curve "
145  "must be given, not both.");
146  }
147 
148  // Parse time interval and return a curve
149  if (time_interval_config)
150  {
151  DBUG("Constructing time interval");
152  auto const start_time =
154  time_interval_config->getConfigParameter<double>("start");
155 
156  auto const end_time =
158  time_interval_config->getConfigParameter<double>("end");
159 
160  // Using very large value for the curve's value, s.t. for any time from
161  // start to end the whole subdomain is deactivated at once.
162  return {{start_time, end_time},
163  {std::numeric_limits<double>::max(),
164  std::numeric_limits<double>::max()},
165  false};
166  }
167 
168  // Try to find the curve.
169  if (curve_name)
170  {
171  DBUG("Using curve '{:s}'", *curve_name);
172  // Return a copy of the curve because the time interval in the other
173  // branch returns a temporary.
174  return *BaseLib::getOrError(curves, *curve_name,
175  "Could not find curve.");
176  }
177 
178  // If we got so far, there is an error: one of the configs must be
179  // available.
180  OGS_FATAL(
181  "In the deactivate subdomain neither a time interval nor a curve are "
182  "given. One of them must be specified.");
183 }
184 
186 static std::pair<Eigen::Vector3d, Eigen::Vector3d> parseLineSegment(
187  BaseLib::ConfigTree const& config)
188 {
189  DBUG("Constructing line segment");
190  auto const start =
192  config.getConfigParameter<std::vector<double>>("start");
193  if (start.size() != 3)
194  {
195  OGS_FATAL(
196  "For construction of a line segment the start point must be a 3D "
197  "point. Got a vector of size {}.",
198  start.size());
199  }
200 
201  auto const end =
203  config.getConfigParameter<std::vector<double>>("end");
204 
205  if (end.size() != 3)
206  {
207  OGS_FATAL(
208  "For construction of a line segment the end point must be a 3D "
209  "point. Got a vector of size {}.",
210  end.size());
211  }
212  return {Eigen::Vector3d{start[0], start[1], start[2]},
213  Eigen::Vector3d{end[0], end[1], end[2]}};
214 }
215 
217  BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh,
218  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
219  std::map<std::string,
220  std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
221  curves)
222 {
223  auto const& time_interval_config =
225  config.getConfigSubtreeOptional("time_interval");
226 
227  auto const& curve_name =
229  config.getConfigParameterOptional<std::string>("time_curve");
230  auto time_interval =
231  parseTimeIntervalOrCurve(time_interval_config, curve_name, curves);
232 
233  auto const line_segment_config =
235  config.getConfigSubtreeOptional("line_segment");
236 
237  if (time_interval_config && line_segment_config)
238  {
239  OGS_FATAL(
240  "When using time interval for subdomain deactivation a line "
241  "segment must not be specified.");
242  }
243 
244  if (curve_name && !line_segment_config)
245  {
246  OGS_FATAL(
247  "When using curve for subdomain deactivation a line segment must "
248  "be specified.");
249  }
250 
251  // If time interval was specified then an empty optional line segment is
252  // used *internally* because the whole selected material ids subdomain will
253  // be deactivated.
254  std::optional<std::pair<Eigen::Vector3d, Eigen::Vector3d>> line_segment;
255  if (curve_name)
256  {
257  line_segment = parseLineSegment(*line_segment_config);
258  }
259 
260  ParameterLib::Parameter<double>* boundary_value_parameter = nullptr;
261  auto boundary_value_parameter_name =
263  config.getConfigParameterOptional<std::string>("boundary_parameter");
264  if (boundary_value_parameter_name)
265  {
266  DBUG("Using parameter {:s}", *boundary_value_parameter_name);
267  boundary_value_parameter = &ParameterLib::findParameter<double>(
268  *boundary_value_parameter_name, parameters, 1, &mesh);
269  }
270 
271  auto deactivated_subdomain_material_ids =
273  config.getConfigParameter("material_ids", std::vector<int>{});
274 
275  if (deactivated_subdomain_material_ids.empty())
276  {
277  OGS_FATAL(
278  "The material IDs of the deactivated subdomains are not given. The "
279  "program terminates now.");
280  }
281 
282  std::sort(deactivated_subdomain_material_ids.begin(),
283  deactivated_subdomain_material_ids.end());
284 
285  if (materialIDs(mesh) == nullptr)
286  {
287  OGS_FATAL(
288  "The mesh doesn't contain materialIDs for subdomain deactivation. "
289  "The program terminates now.");
290  }
291 
292  auto deactivated_subdomain_mesh = createDeactivatedSubdomainMesh(
293  mesh, deactivated_subdomain_material_ids);
294 
295  return {std::move(time_interval), line_segment,
296  std::move(deactivated_subdomain_mesh), boundary_value_parameter};
297 }
298 
299 std::vector<DeactivatedSubdomain> createDeactivatedSubdomains(
300  BaseLib::ConfigTree const& config,
301  MeshLib::Mesh const& mesh,
302  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
303  std::map<std::string,
304  std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
305  curves)
306 {
307  std::vector<DeactivatedSubdomain> deactivated_subdomains;
308  // Deactivated subdomains
309  if (auto subdomains_config =
311  config.getConfigSubtreeOptional("deactivated_subdomains"))
312  {
313  INFO("There are subdomains being deactivated.");
314 
315  auto const deactivated_subdomain_configs =
317  subdomains_config->getConfigSubtreeList("deactivated_subdomain");
318  std::transform(std::begin(deactivated_subdomain_configs),
319  std::end(deactivated_subdomain_configs),
320  std::back_inserter(deactivated_subdomains),
321  [&](auto const& config) {
323  config, mesh, parameters, curves);
324  });
325  }
326  return deactivated_subdomains;
327 }
328 
329 } // 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:34
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
Definition of the Mesh class.
Definition of the Node class.
std::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
Definition: ConfigTree.cpp:159
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:100
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
Properties & getProperties()
Definition: Mesh.h:128
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:94
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:238
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Definition: Mesh.h:79
Map::mapped_type & getOrError(Map &map, Key const &key, std::string const &error_message)
Definition: Algorithm.h:147
constexpr bool all_of(List const &values)
Checks if all of the elements in the given list are true.
Definition: Algorithm.h:341
constexpr bool any_of(List const &values)
Checks if any of the elements in the given list is true.
Definition: Algorithm.h:325
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition: Mesh.h:313
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:264
std::unique_ptr< MeshLib::Mesh > createMeshFromElementSelection(std::string mesh_name, std::vector< MeshLib::Element * > const &elements)
Definition: Mesh.cpp:274
std::vector< Element * > cloneElements(std::vector< Element * > const &elements)
Clones a vector of elements using the Element::clone() function.
static DeactivatedSubdomainMesh createDeactivatedSubdomainMesh(MeshLib::Mesh const &mesh, std::vector< int > const &deactivated_subdomain_material_ids)
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< 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)
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 std::pair< Eigen::Vector3d, Eigen::Vector3d > parseLineSegment(BaseLib::ConfigTree const &config)
Returns a line segment represented by its begin and end points.
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)