Medusa  1.1
Coordinate Free Mehless Method implementation
BallShape.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_BALLSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_BALLSHAPE_HPP_
3 
4 #include "BallShape_fwd.hpp"
7 #include "GeneralFill.hpp"
10 
16 namespace mm {
17 
18 template <typename vec_t>
20  int type) const {
21  if (type == 0) type = -1;
22  DomainDiscretization<vec_t> domain(*this);
23  int num = iceil(2*Pi<scalar_t>::value*radius_/step);
25  radius_, num);
26  for (const auto& p : points) {
27  domain.addBoundaryNode(center_+p, type, p.normalized());
28  }
29  return domain;
30 }
31 
32 template <typename vec_t>
34  scalar_t step, int internal_type, int boundary_type) const {
35  auto d = discretizeBoundaryWithStep(step, boundary_type);
36  if (internal_type == 0) internal_type = 1;
37  if (dim == 1) {
38  int num_of_points = iceil(2*radius_/step) - 1;
39  for (const vec_t& v : linspace(bbox().first, bbox().second, {num_of_points}, false)) {
40  d.addInternalNode(v, internal_type);
41  }
42  } else if (dim == 2) {
43  // Use concentric circles with step dx.
44  int num_of_circles = iceil(radius_ / step); // not counting outer circle
45  scalar_t dr = radius_ / num_of_circles;
46  d.addInternalNode(center_, internal_type);
47  for (int i = 1; i < num_of_circles; ++i) {
48  scalar_t r = i*dr;
49  int num_of_points = iceil(Pi<scalar_t>::value / std::asin(step/2.0/r));
50  scalar_t dfi = 2*Pi<scalar_t>::value / num_of_points;
51  for (int j = 0; j < num_of_points; ++j) {
52  d.addInternalNode(center_ + r * vec_t({std::cos(j*dfi), std::sin(j*dfi)}),
53  internal_type);
54  }
55  }
56  } else if (dim >= 3) { // Call default fill engine (could do concentric n-dim spheres).
57  assert_msg(false, "This domain does not support filling with density, "
58  "use a GeneralSurfaceFill engine instead.");
59  }
60  return d;
61 }
62 
63 template <typename vec_t>
64 std::ostream& BallShape<vec_t>::print(std::ostream& os) const {
65  return os << "BallShape(" << center_.transpose() << ", " << radius_ << ")";
66 }
67 
68 template <typename vec_t>
71  int type) const {
72  if (type == 0) type = -1;
73  DomainDiscretization<vec_t> domain(*this);
74  if (dim == 1) { // always only two points
75  domain.addBoundaryNode({center_[0] - radius_}, type, vec_t(-1.0));
76  domain.addBoundaryNode({center_[0] + radius_}, type, vec_t(1.0));
77  } else if (dim == 2) {
78  scalar_t cur_alpha = 0.0;
79  scalar_t max_alpha = 2*Pi<scalar_t>::value;
80  Range<scalar_t> alphas = {cur_alpha};
81  while (cur_alpha < max_alpha) {
82  scalar_t r = dr(center_ + radius_ * vec_t({std::cos(cur_alpha), std::sin(cur_alpha)}));
83  cur_alpha += 2*std::asin(r/2.0/radius_);
84  alphas.push_back(cur_alpha);
85  }
86 
87  scalar_t shrink = 2*Pi<scalar_t>::value / cur_alpha;
88  for (int i = 0; i < alphas.size()-1; ++i) {
89  scalar_t a = alphas[i] * shrink;
90  vec_t normal = {std::cos(a), std::sin(a)};
91  domain.addBoundaryNode(center_ + radius_ * normal, type, normal);
92  }
93  } else if (dim == 3) {
94  // define modified density
95  auto dr2upper = [&](const Vec<scalar_t, 2>& v) {
96  double d = (1 + v.squaredNorm());
97  vec_t p = {2*v[0]/d, 2*v[1]/d, (1-v.squaredNorm())/d};
98  return d/2 * dr(center_ + radius_ * p) / radius_;
99  };
100 
101  BallShape<Vec<scalar_t, 2>> circ(0.0, 1.0);
103  fill.seed(0).proximityTolerance(0.99);
104  auto domain2d = circ.discretizeWithDensity(dr2upper, fill);
105 
106  for (int i : domain2d.interior()) {
107  auto v = domain2d.pos(i);
108  double d = (1 + v.squaredNorm());
109  vec_t p = {2*v[0]/d, 2*v[1]/d, (1-v.squaredNorm())/d};
110  domain.addBoundaryNode(center_ + radius_*p, type, p.normalized());
111  }
112 
113  for (int i : domain2d.boundary()) {
114  auto v = domain2d.pos(i);
115  vec_t p = {v[0], v[1], 0};
116  domain.addBoundaryNode(center_ + radius_*p, type, p.normalized());
117  }
118 
119  // lower hemisphere
120  auto dr2lower = [&](const Vec<scalar_t, 2>& v) {
121  double d = (1 + v.squaredNorm());
122  vec_t p = {2*v[0]/d, 2*v[1]/d, -(1-v.squaredNorm())/d};
123  return d/2 * dr(center_ + radius_ * p) / radius_;
124  };
125  domain2d = circ.discretizeWithDensity(dr2lower, fill);
126  for (int i : domain2d.interior()) {
127  auto v = domain2d.pos(i);
128  double d = (1 + v.squaredNorm());
129  vec_t p = {2*v[0]/d, 2*v[1]/d, -(1-v.squaredNorm())/d};
130  domain.addBoundaryNode(center_ + radius_*p, type, p.normalized());
131  }
132  } else {
133  assert_msg(false, "This domain does not support filling with density, "
134  "use a GeneralSurfaceFill engine instead.");
135  }
136  return domain;
137 }
138 
139 } // namespace mm
140 
141 #endif // MEDUSA_BITS_DOMAINS_BALLSHAPE_HPP_
mm::BallShape
Class for working with ball shaped domains.
Definition: BallShape_fwd.hpp:24
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::DomainDiscretization
Class representing domain discretization along with an associated shape.
Definition: DomainDiscretization_fwd.hpp:46
mm::BallShape::discretizeBoundaryWithStep
DomainDiscretization< vec_t > discretizeBoundaryWithStep(scalar_t step, int type) const override
Returns a discretization of the boundary of this shape with approximately uniform distance step betwe...
Definition: BallShape.hpp:19
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
DomainDiscretization.hpp
mm::BallShape::discretizeWithStep
DomainDiscretization< vec_t > discretizeWithStep(scalar_t step, int internal_type, int boundary_type) const override
Returns a discretization of this shape with approximately uniform distance step between nodes.
Definition: BallShape.hpp:33
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::GeneralFill::proximityTolerance
GeneralFill & proximityTolerance(scalar_t zeta)
Set proximity tolerance.
numutils.hpp
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::iceil
int iceil(T x)
Ceils a floating point to an integer.
Definition: numutils.hpp:27
mm::linspace
Range< Vec< scalar_t, dim > > linspace(const Vec< scalar_t, dim > &beg, const Vec< scalar_t, dim > &end, const Vec< int, dim > &counts, const Vec< bool, dim > include_boundary=true)
Multidimensional clone of Matlab's linspace function.
Definition: numutils.hpp:210
mm::BallShape::print
std::ostream & print(std::ostream &os) const override
Output information about this shape to given output stream os.
Definition: BallShape.hpp:64
mm::GeneralFill::seed
GeneralFill & seed(int seed)
Set custom seed for the random number generator.
Definition: GeneralFill_fwd.hpp:65
assert.hpp
mm::DomainShape::discretizeWithDensity
virtual DomainDiscretization< vec_t > discretizeWithDensity(const std::function< scalar_t(vec_t)> &dr, int internal_type, int boundary_type) const
Returns a discretization of the domain with spatially variable step.
Definition: DomainShape.hpp:77
mm::BallShape::discretizeBoundaryWithDensity
DomainDiscretization< vec_t > discretizeBoundaryWithDensity(const std::function< scalar_t(vec_t)> &dr, int type) const override
Discretizes boundary with given density and fill engine.
Definition: BallShape.hpp:70
discretization_helpers.hpp
BallShape_fwd.hpp
mm::Range::size
int size() const
Returns number of elements.
Definition: Range_fwd.hpp:185
mm::Pi
Value of Pi in type T. Usage:
Definition: Config.hpp:40
mm::discretization_helpers::SphereDiscretization::construct
static std::vector< Eigen::Matrix< scalar_t, dim, 1 >, Eigen::aligned_allocator< Eigen::Matrix< scalar_t, dim, 1 > > > construct(scalar_t radius, int num_samples, generator_t &generator)
Construct a randomized discretization.
Definition: discretization_helpers.hpp:85
mm::Range
An extension of std::vector<T> to support additional useful operations.
Definition: Range_fwd.hpp:30
GeneralFill.hpp
mm::DomainDiscretization::addBoundaryNode
int addBoundaryNode(const vec_t &point, int type, const vec_t &normal)
Adds a boundary node with given type and normal to the domain.
Definition: DomainDiscretization.hpp:56
mm::DomainShape::scalar_t
vec_t::Scalar scalar_t
Scalar data type used in computation.
Definition: DomainShape_fwd.hpp:55