Caribou
Mesh.h
1 #pragma once
2 
3 #include <Caribou/config.h>
4 #include <Caribou/Topology/BaseMesh.h>
5 #include <Caribou/Topology/Domain.h>
6 
7 #include <Eigen/Dense>
8 #include <memory>
9 #include <vector>
10 
11 namespace caribou::topology {
12 
13  template <
14  UNSIGNED_INTEGER_TYPE WorldDimension,
15  typename MatrixType = Eigen::Matrix<FLOATING_POINT_TYPE, Eigen::Dynamic, WorldDimension, (WorldDimension>1?Eigen::RowMajor:Eigen::ColMajor)>>
16  class Mesh : public BaseMesh {
17  public:
18  using PositionsContainer = MatrixType;
19 
20  static_assert(std::is_base_of_v<Eigen::MatrixBase<PositionsContainer>, PositionsContainer>, "The matrix type must inherit Eigen::MatrixBase");
21  static_assert(Eigen::MatrixBase<PositionsContainer>::ColsAtCompileTime == WorldDimension, "The matrix type must have N columns at compile time for a mesh of dimension N.");
22  static_assert(WorldDimension == 1 or WorldDimension == 2 or WorldDimension == 3, "The world dimension must be 1, 2 or 3.");
23 
24  static constexpr INTEGER_TYPE Dimension = WorldDimension;
25  using Real = typename PositionsContainer::Scalar;
26  using WorldCoordinates = Eigen::Matrix<Real, 1, Dimension>;
27 
28  template <typename Element, typename NodeIndex = UNSIGNED_INTEGER_TYPE>
29  using Domain = Domain<Mesh<Dimension, MatrixType>, Element, NodeIndex>;
30 
34  Mesh() = default;
35 
39  virtual ~Mesh() = default;
40 
46  explicit Mesh(const std::vector<WorldCoordinates> & positions) {
47  p_positions.resize(positions.size(), Dimension);
48  for (std::size_t i = 0; i < positions.size(); ++i) {
49  p_positions.row(static_cast<Eigen::Index>(i)) = positions[i];
50  }
51  }
52 
58  explicit Mesh(const MatrixType & positions) : p_positions(positions) {}
59 
65  template <typename EigenDerived>
66  explicit Mesh(const EigenDerived * positions) : p_positions(*positions) {}
67 
69  Mesh(const Mesh & other)
70  : p_positions(other.p_positions) {}
71 
73  Mesh(Mesh && other) noexcept {
74  swap(*this, other);
75  }
76 
78  auto operator=(Mesh other) noexcept -> Mesh &
79  {
80  swap(*this, other);
81  return *this;
82  }
83 
87  [[nodiscard]]
88  inline auto dimension() const -> UNSIGNED_INTEGER_TYPE final {
89  return Dimension;
90  };
91 
95  [[nodiscard]]
96  inline auto number_of_domains() const -> UNSIGNED_INTEGER_TYPE final {
97  return p_domains.size();
98  }
99 
103  [[nodiscard]]
104  inline auto number_of_nodes() const -> UNSIGNED_INTEGER_TYPE final {return p_positions.rows();};
105 
109  [[nodiscard]]
110  inline auto domains() const -> std::vector<std::pair<std::string, const BaseDomain *>> {
111  std::vector<std::pair<std::string, const BaseDomain *>> domains;
112  domains.reserve(p_domains.size());
113  for (const auto & d : p_domains) {
114  domains.emplace_back(d.first, d.second.get());
115  }
116  return domains;
117  }
118 
123  [[nodiscard]]
124  inline auto domain(const UNSIGNED_INTEGER_TYPE & i) const -> const BaseDomain * {
125  caribou_assert(
126  (i < p_domains.size()) &&
127  "Trying to get a domain outside from a domain id larger than the number of domains."
128  );
129 
130  return p_domains[i].second.get();
131  }
132 
134  [[nodiscard]]
135  inline auto domain(const std::string & name) const -> const BaseDomain * {
136  for (const auto & d : p_domains) {
137  if (d.first == name) {
138  return d.second.get();
139  }
140  }
141 
142  return nullptr;
143  }
144 
153  template<typename Element, typename NodeIndex, typename... Args>
154  inline
155  typename std::enable_if_t<std::is_integral_v<NodeIndex>, Domain<Element, NodeIndex> *>
156  add_domain(const std::string & name, Args ...args) {
157  static_assert(geometry::traits<Element>::Dimension == Dimension, "The dimension of the mesh doesn't match the dimension of the element type.");
158  for (const auto & d : p_domains) {
159  if (d.first == name) {
160  throw std::domain_error(std::string("A domain with the name ") + name + " already exists.");
161  }
162  }
163 
164  auto domain_ptr = new Domain<Element, NodeIndex>(this, std::forward<Args>(args)...);
165  auto res = p_domains.emplace(p_domains.end(), name, std::unique_ptr<BaseDomain>(static_cast<BaseDomain *>(domain_ptr)));
166  if (res->second) {
167  return domain_ptr;
168  } else {
169  delete domain_ptr;
170  return static_cast<Domain<Element, NodeIndex>*> (nullptr);
171  }
172  }
173 
182  template<typename Element, typename NodeIndex, typename... Args>
183  inline
184  typename std::enable_if_t<std::is_integral_v<NodeIndex>, Domain<Element, NodeIndex> *>
185  add_domain(Args ...args) {
186  static_assert(geometry::traits<Element>::Dimension == Dimension, "The dimension of the mesh doesn't match the dimension of the element type.");
187 
188  auto domain_ptr = new Domain<Element, NodeIndex>(this, std::forward<Args>(args)...);
189 
190  std::ostringstream ss;
191  ss << "domain_" << domain_ptr;
192 
193  auto res = p_domains.emplace(p_domains.end(), ss.str(), std::unique_ptr<BaseDomain>(static_cast<BaseDomain *>(domain_ptr)));
194  if (res->second) {
195  return domain_ptr;
196  } else {
197  delete domain_ptr;
198  return static_cast<Domain<Element, NodeIndex>*> (nullptr);
199  }
200  }
201 
210  template<typename Element, typename... Args>
211  inline auto add_domain(const std::string & name, Args ...args) -> Domain<Element> * {
212  static_assert(geometry::traits<Element>::Dimension == Dimension, "The dimension of the mesh doesn't match the dimension of the element type.");
213  return add_domain<Element, UNSIGNED_INTEGER_TYPE, Args...>(name, std::forward<Args>(args)...);
214  }
215 
224  template<typename Element, typename... Args>
225  inline auto add_domain(const char * name, Args ...args) -> Domain<Element> * {
226  return add_domain<Element, Args...>(std::string(name), std::forward<Args>(args)...);
227  }
228 
236  template<typename Element, typename... Args>
237  inline auto add_domain(Args ...args) -> Domain<Element> * {
238  static_assert(geometry::traits<Element>::Dimension == Dimension, "The dimension of the mesh doesn't match the dimension of the element type.");
239  auto domain_ptr = new Domain<Element>(this, std::forward<Args>(args)...);
240 
241  std::ostringstream ss;
242  ss << "domain_" << domain_ptr;
243 
244  auto res = p_domains.emplace(p_domains.end(), ss.str(), std::unique_ptr<BaseDomain>(static_cast<BaseDomain *>(domain_ptr)));
245  if (res->second) {
246  return domain_ptr;
247  } else {
248  delete domain_ptr;
249  return static_cast<Domain<Element>*> (nullptr);
250  }
251  }
252 
258  inline auto remove_domain(const std::string & name) -> bool {
259  auto it = p_domains.begin();
260  while (it != p_domains.end() and it->first != name);
261  if (it != p_domains.end()) {
262  p_domains.erase(it);
263  return true;
264  }
265  return false;
266  }
267 
273  inline auto remove(const BaseDomain * domain) -> bool {
274  auto it = p_domains.begin();
275  while (it != p_domains.end() and it->second.get() != domain);
276  if (it != p_domains.end()) {
277  p_domains.erase(it);
278  return true;
279  }
280  return false;
281  }
282 
286  [[nodiscard]]
287  inline auto position(UNSIGNED_INTEGER_TYPE index) const {
288  caribou_assert(static_cast<Eigen::Index>(index) < p_positions.rows()
289  && "Trying to access the position vector at a node index greater "
290  "than the number of nodes in the mesh.");
291  return p_positions.row(index);
292  }
293 
308  template <typename IntegerType, std::size_t N>
309  inline auto positions(const IntegerType(&indices)[N]) const {
310  Eigen::Matrix<Real, N, Dimension, (Dimension>1?Eigen::RowMajor:Eigen::ColMajor)> positions;
311  for (std::size_t i = 0; i < N; ++i) {
312  positions.row(i) = p_positions.row(indices[i]);
313  }
314  return positions;
315  }
316 
332  template <typename IntegerType, std::size_t N>
333  inline auto positions(const std::array<IntegerType, N> & indices) const {
334  return positions<IntegerType, N>(indices.data());
335  }
336 
352  template <typename IntegerType>
353  inline auto positions(const std::vector<IntegerType> & indices) const {
354  const auto number_of_indices = indices.size();
355  Eigen::Matrix<Real, Eigen::Dynamic, Dimension, (Dimension == 1) ? Eigen::ColMajor : Eigen::RowMajor> positions;
356  positions.resize(number_of_indices, Dimension);
357  for (std::size_t i = 0; i < number_of_indices; ++i) {
358  positions.row(i) = p_positions.row(indices[i]);
359  }
360  return positions;
361  }
362 
379  template <typename EigenDerived>
380  inline auto positions(const Eigen::EigenBase<EigenDerived> & indices) const {
381  static_assert(EigenDerived::RowsAtCompileTime == 1 or EigenDerived::ColsAtCompileTime == 1,
382  "Only vector type matrix can be used as the indices container.");
383  const auto number_of_indices = indices.size();
384  Eigen::Matrix<Real,
385  EigenDerived::RowsAtCompileTime == 1 ? EigenDerived::ColsAtCompileTime : EigenDerived::RowsAtCompileTime,
386  Dimension,
387  (Dimension>1 ? Eigen::RowMajor : Eigen::ColMajor)> positions;
388  positions.resize(number_of_indices, Dimension);
389  for (Eigen::Index i = 0; i < number_of_indices; ++i) {
390  positions.row(i) = p_positions.row(indices.derived()[i]);
391  }
392  return positions;
393  }
394 
396  friend void swap(Mesh & first, Mesh& second) noexcept
397  {
398  // enable ADL
399  using std::swap;
400  swap(first.p_positions, second.p_positions);
401  }
402 
403  private:
404  PositionsContainer p_positions;
405  std::vector<std::pair<std::string, std::unique_ptr<BaseDomain>>> p_domains;
406  };
407 
408  // Template deduction guides
409  template <typename MatrixType>
410  Mesh(const MatrixType &) -> Mesh<MatrixType::ColsAtCompileTime, MatrixType>;
411 
412  template <typename MatrixType>
413  Mesh(const MatrixType *) ->
414  Mesh<
415  MatrixType::ColsAtCompileTime,
416  Eigen::Ref<
417  const std::decay_t<MatrixType>,
418  Eigen::internal::traits<MatrixType>::Alignment,
419  Eigen::Stride<MatrixType::OuterStrideAtCompileTime,MatrixType::InnerStrideAtCompileTime>
420  >
421  >;
422 }
caribou::topology::Mesh::number_of_nodes
auto number_of_nodes() const -> UNSIGNED_INTEGER_TYPE final
Definition: Mesh.h:104
caribou::topology::Mesh::add_domain
auto add_domain(const char *name, Args ...args) -> Domain< Element > *
Definition: Mesh.h:225
caribou::topology::Mesh::add_domain
auto add_domain(const std::string &name, Args ...args) -> Domain< Element > *
Definition: Mesh.h:211
caribou::topology::Mesh::~Mesh
virtual ~Mesh()=default
caribou::topology::Mesh::position
auto position(UNSIGNED_INTEGER_TYPE index) const
Definition: Mesh.h:287
caribou::topology::Mesh::swap
friend void swap(Mesh &first, Mesh &second) noexcept
Definition: Mesh.h:396
caribou::topology::Domain
Definition: Domain.h:75
caribou::topology::Mesh::positions
auto positions(const IntegerType(&indices)[N]) const
Definition: Mesh.h:309
caribou::topology::Mesh::positions
auto positions(const Eigen::EigenBase< EigenDerived > &indices) const
Definition: Mesh.h:380
caribou::topology::Mesh::Mesh
Mesh(const std::vector< WorldCoordinates > &positions)
Construct the unstructured mesh with a set of point positions from a position vector.
Definition: Mesh.h:46
caribou::topology::Mesh::Mesh
Mesh(const Mesh &other)
Definition: Mesh.h:69
caribou::topology::Mesh::add_domain
std::enable_if_t< std::is_integral_v< NodeIndex >, Domain< Element, NodeIndex > * > add_domain(Args ...args)
Definition: Mesh.h:185
caribou::topology::Mesh::Mesh
Mesh()=default
caribou::topology::Mesh::dimension
auto dimension() const -> UNSIGNED_INTEGER_TYPE final
Definition: Mesh.h:88
caribou::topology::Mesh::add_domain
auto add_domain(Args ...args) -> Domain< Element > *
Definition: Mesh.h:237
caribou::geometry::traits
Definition: Element.h:14
caribou::topology::Mesh::operator=
auto operator=(Mesh other) noexcept -> Mesh &
Definition: Mesh.h:78
caribou::topology::Mesh::remove
auto remove(const BaseDomain *domain) -> bool
Definition: Mesh.h:273
caribou::topology::Mesh::positions
auto positions(const std::array< IntegerType, N > &indices) const
Definition: Mesh.h:333
caribou::topology::Mesh::domains
auto domains() const -> std::vector< std::pair< std::string, const BaseDomain * >>
Definition: Mesh.h:110
caribou::topology::Mesh::Mesh
Mesh(Mesh &&other) noexcept
Definition: Mesh.h:73
caribou::topology::Mesh
Definition: Mesh.h:16
caribou::topology::Mesh::Mesh
Mesh(const MatrixType &positions)
Construct the unstructured mesh with an Eigen matrix containing the position vector (NxD with N nodes...
Definition: Mesh.h:58
caribou::topology::BaseDomain
A domain is subspace of a mesh containing a set of points and the topological relation between them.
Definition: BaseDomain.h:13
caribou::topology::Mesh::add_domain
std::enable_if_t< std::is_integral_v< NodeIndex >, Domain< Element, NodeIndex > * > add_domain(const std::string &name, Args ...args)
Definition: Mesh.h:156
caribou::topology::Mesh::number_of_domains
auto number_of_domains() const -> UNSIGNED_INTEGER_TYPE final
Definition: Mesh.h:96
caribou::topology::Mesh::domain
auto domain(const std::string &name) const -> const BaseDomain *
Definition: Mesh.h:135
caribou::topology::BaseMesh
A mesh is a discrete view of a space by the mean of nodes.
Definition: BaseMesh.h:16
caribou::topology::Mesh::remove_domain
auto remove_domain(const std::string &name) -> bool
Definition: Mesh.h:258
caribou::topology::Mesh::domain
auto domain(const UNSIGNED_INTEGER_TYPE &i) const -> const BaseDomain *
Definition: Mesh.h:124
caribou::topology::Mesh::Mesh
Mesh(const EigenDerived *positions)
Construct the unstructured mesh with an Eigen matrix containing the position vector (NxD with N nodes...
Definition: Mesh.h:66
caribou::topology::Mesh::positions
auto positions(const std::vector< IntegerType > &indices) const
Definition: Mesh.h:353