Medusa  1.1
Coordinate Free Mehless Method implementation
discretization_helpers_advanced.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_ADVANCED_HPP_
2 #define MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_ADVANCED_HPP_
3 
12 #include <medusa/Config.hpp>
13 #include <medusa/bits/domains/PolygonShape.hpp>
14 #include <medusa/bits/domains/GeneralFill.hpp>
15 #include <medusa/bits/types/Range.hpp>
16 
17 #include <Eigen/Geometry>
18 
19 #include <vector>
20 
21 namespace mm {
22 namespace discretization_helpers {
23 
44 template <typename vec_t, typename func_t>
45 Range<vec_t> discretizeTriangleWithDensity(const vec_t& p1, const vec_t& p2, const vec_t& p3,
46  const vec_t& normal, const func_t& h,
47  bool only_interior = true) {
48  static_assert(vec_t::dim == 3, "This method is intended for 3D vectors only.");
49  using scalar_t = typename vec_t::scalar_t;
50  using V2 = mm::Vec<scalar_t, 2>;
51  static constexpr int dim = vec_t::dim;
52 
53  // Rotate to flat
54  vec_t zunit = {0, 0, 1};
55  scalar_t c = normal.dot(zunit);
56  if (c < 0) {
57  zunit = -zunit;
58  c = -c;
59  }
60  vec_t vp = normal.cross(zunit);
61 
62  Eigen::Matrix<scalar_t, dim, dim> R;
63  R.setIdentity();
64  Eigen::Matrix<scalar_t, dim, dim> VP;
65  VP << 0, -vp[2], vp[1], vp[2], 0, -vp[0], -vp[1], vp[0], 0;
66  R += VP + 1 / (1 + c) * VP * VP;
67 
68  vec_t np1 = R * p1;
69  vec_t np2 = R * p2;
70  vec_t np3 = R * p3;
71  // All z coordinates are the same, but still take the average.
72  scalar_t z = 1.0 / 3.0 * (np1[2] + np2[2] + np3[2]);
73 
74 
75  // Fill the flat triangle.
76  std::vector<V2> pts = {np1.template head<2>(), np2.template head<2>(),
77  np3.template head<2>()};
78  mm::PolygonShape<V2> triangle(pts);
79  auto new_h = [&](const V2& p) { return h(R.transpose() * vec_t(p[0], p[1], z)); };
80  auto domain = triangle.discretizeBoundaryWithDensity(new_h);
81  mm::GeneralFill<V2> fill;
82  fill.seed(0);
83  fill(domain, new_h);
84 
85  auto idxs = only_interior ? domain.interior() : domain.all();
86  std::vector<vec_t> result;
87  result.reserve(idxs.size());
88  for (int i : idxs) {
89  auto ip = domain.pos(i);
90  result.emplace_back(R.transpose() * vec_t(ip[0], ip[1], z));
91  }
92  return result;
93 }
94 
95 } // namespace discretization_helpers
96 } // namespace mm
97 
98 #endif // MEDUSA_BITS_DOMAINS_DISCRETIZATION_HELPERS_ADVANCED_HPP_
mm::PolygonShape::discretizeBoundaryWithDensity
DomainDiscretization< vec_t > discretizeBoundaryWithDensity(const std::function< scalar_t(vec_t)> &dr, int type) const override
Returns a discretization of the boundary of this shape with approximately uniform distance step betwe...
Definition: PolygonShape.hpp:118
PolygonShape.hpp
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::Vec
Eigen::Matrix< scalar_t, dim, 1, Eigen::ColMajor|Eigen::AutoAlign, dim, 1 > Vec
Fixed size vector type, representing a mathematical 1d/2d/3d vector.
Definition: Vec_fwd.hpp:31
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::GeneralFill
Implements general n-d node placing algorithm, as described in https://arxiv.org/abs/1812....
Definition: GeneralFill_fwd.hpp:43
mm::PolygonShape
Shape representing a simple (i.e. non self intersecting) nonempty polygon in 2D, which is given as a ...
Definition: PolygonShape_fwd.hpp:27
mm::discretization_helpers::discretizeTriangleWithDensity
Range< vec_t > discretizeTriangleWithDensity(const vec_t &p1, const vec_t &p2, const vec_t &p3, const vec_t &normal, const func_t &h, bool only_interior=true)
Discretize a triangle in 3D space with given density.
Definition: discretization_helpers_advanced.hpp:45
Config.hpp
mm::GeneralFill::seed
GeneralFill & seed(int seed)
Set custom seed for the random number generator.
Definition: GeneralFill_fwd.hpp:65
mm::Range< vec_t >
Range.hpp
GeneralFill.hpp