Medusa  1.1
Coordinate Free Mehless Method implementation
BoxShape.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_BOXSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_BOXSHAPE_HPP_
3 
9 #include "BoxShape_fwd.hpp"
10 #include <medusa/bits/utils/assert.hpp>
11 #include "DomainDiscretization.hpp"
12 #include "discretization_helpers.hpp"
13 #include "DomainShape.hpp"
14 #include "GeneralFill.hpp"
15 #include <medusa/bits/utils/numutils.hpp>
16 #include <medusa/bits/types/Vec.hpp>
17 
18 namespace mm {
19 
20 template <typename vec_t>
21 BoxShape<vec_t>::BoxShape(vec_t beg, vec_t end) : beg_(beg), end_(end) {
22  for (int i = 0; i < dim; ++i) {
23  assert_msg(std::abs(beg_[i] - end_[i]) >= margin_,
24  "Box has no volume, sides in dimension %d have same coordinate %g.",
25  i, beg_[i]);
26  if (beg_[i] > end_[i]) std::swap(beg_[i], end_[i]);
27  }
28 }
29 
30 template <typename vec_t>
31 bool BoxShape<vec_t>::contains(const vec_t& point) const {
32  for (int i = 0; i < dim; ++i) {
33  if (!(beg_[i] - margin_ <= point[i] && point[i] <= end_[i] + margin_)) {
34  return false;
35  }
36  }
37  return true;
38 }
39 
40 template <typename vec_t>
41 DomainDiscretization<vec_t> BoxShape<vec_t>::discretizeBoundaryWithStep(scalar_t step,
42  int type) const {
43  Vec<int, dim> counts;
44  for (int i = 0; i < dim; ++i) counts[i] = iceil((end_[i] - beg_[i]) / step) + 1;
45  return discretizeBoundary(counts, type);
46 }
47 
48 template <typename vec_t>
49 DomainDiscretization<vec_t> BoxShape<vec_t>::discretizeBoundary(const Vec<int, dim>& counts,
50  int type) const {
51  for (int x : counts) {
52  assert_msg(x >= 2, "All counts must be greater than 2, got counts = %s", counts);
53  }
54  DomainDiscretization<vec_t> domain(*this);
55  vec_t normal, point;
56  if (dim == 1) {
57  domain.addBoundaryNode(beg_, type == 0 ? -1 : type, vec_t(-1.0));
58  domain.addBoundaryNode(end_, type == 0 ? -2 : type, vec_t(1.0));
59  } else if (dim == 2) {
60  const int TOP = (type == 0) ? -4 : type;
61  const int RIGHT = (type == 0) ? -2 : type;
62  const int BOTTOM = (type == 0) ? -3 : type;
63  const int LEFT = (type == 0) ? -1 : type;
64 
65  scalar_t step = (end_[0] - beg_[0]) / (counts[0] - 1);
66  for (int i = 1; i < counts[0]-1; ++i) {
67  point[0] = beg_[0] + i*step;
68  point[1] = beg_[1];
69  normal[0] = 0; normal[1] = -1;
70  domain.addBoundaryNode(point, BOTTOM, normal);
71  point[1] = end_[1];
72  normal[0] = 0; normal[1] = 1;
73  domain.addBoundaryNode(point, TOP, normal);
74  }
75  step = (end_[1] - beg_[1]) / (counts[1] - 1);
76  for (int i = 1; i < counts[1]-1; ++i) {
77  point[1] = beg_[1] + i*step;
78  point[0] = beg_[0];
79  normal[0] = -1; normal[1] = 0;
80  domain.addBoundaryNode(point, LEFT, normal);
81  point[0] = end_[0];
82  normal[0] = 1; normal[1] = 0;
83  domain.addBoundaryNode(point, RIGHT, normal);
84  }
85  domain.addBoundaryNode(beg_, LEFT, vec_t({-1, -1}).normalized());
86  domain.addBoundaryNode(end_, RIGHT, vec_t({1, 1}).normalized());
87  domain.addBoundaryNode(vec_t({beg_[0], end_[1]}), TOP, vec_t({-1, 1}).normalized());
88  domain.addBoundaryNode(vec_t({end_[0], beg_[1]}), BOTTOM, vec_t({1, -1}).normalized());
89  } else if (dim == 3) {
90  int LEFT = -1, RIGHT = -2, FRONT = -3, BACK = -4, BOTTOM = -5, TOP = -6;
91  if (type != 0) LEFT = RIGHT = FRONT = BACK = BOTTOM = TOP = type;
92  scalar_t step = (end_[0] - beg_[0]) / (counts[0] - 1); // x
93  for (int i = 1; i < counts[0]-1; ++i) {
94  point[0] = beg_[0] + i*step;
95  point[1] = beg_[1];
96  point[2] = beg_[2];
97  normal = vec_t({0, -1, -1}).normalized();
98  domain.addBoundaryNode(point, BOTTOM, normal);
99  point[1] = beg_[1];
100  point[2] = end_[2];
101  normal = vec_t({0, -1, 1}).normalized();
102  domain.addBoundaryNode(point, FRONT, normal);
103  point[1] = end_[1];
104  point[2] = beg_[2];
105  normal = vec_t({0, 1, -1}).normalized();
106  domain.addBoundaryNode(point, BACK, normal);
107  point[1] = end_[1];
108  point[2] = end_[2];
109  normal = vec_t({0, 1, 1}).normalized();
110  domain.addBoundaryNode(point, TOP, normal);
111 
112  scalar_t inner_step = (end_[1] - beg_[1]) / (counts[1] - 1); // y
113  for (int j = 1; j < counts[1]-1; ++j) {
114  point[1] = beg_[1] + j*inner_step;
115  point[2] = beg_[2];
116  normal = vec_t({0, 0, -1});
117  domain.addBoundaryNode(point, BOTTOM, normal);
118  point[2] = end_[2];
119  normal = vec_t({0, 0, 1});
120  domain.addBoundaryNode(point, TOP, normal);
121  }
122  inner_step = (end_[2] - beg_[2]) / (counts[2] - 1); // z
123  for (int j = 1; j < counts[2]-1; ++j) {
124  point[2] = beg_[2] + j*inner_step;
125  point[1] = beg_[1];
126  normal = vec_t({0, -1, 0});
127  domain.addBoundaryNode(point, FRONT, normal);
128  point[1] = end_[1];
129  normal = vec_t({0, 1, 0});
130  domain.addBoundaryNode(point, BACK, normal);
131  }
132  }
133  step = (end_[1] - beg_[1]) / (counts[1] - 1); // y
134  for (int i = 1; i < counts[1]-1; ++i) {
135  point[1] = beg_[1] + i*step;
136  point[0] = beg_[0];
137  point[2] = beg_[2];
138  normal = vec_t({-1, 0, -1}).normalized();
139  domain.addBoundaryNode(point, LEFT, normal);
140  point[0] = beg_[0];
141  point[2] = end_[2];
142  normal = vec_t({-1, 0, 1}).normalized();
143  domain.addBoundaryNode(point, TOP, normal);
144  point[0] = end_[0];
145  point[2] = beg_[2];
146  normal = vec_t({1, 0, -1}).normalized();
147  domain.addBoundaryNode(point, BOTTOM, normal);
148  point[0] = end_[0];
149  point[2] = end_[2];
150  normal = vec_t({1, 0, 1}).normalized();
151  domain.addBoundaryNode(point, RIGHT, normal);
152 
153  scalar_t inner_step = (end_[2] - beg_[2]) / (counts[2] - 1); // z
154  for (int j = 1; j < counts[2]-1; ++j) {
155  point[2] = beg_[2] + j*inner_step;
156  point[0] = beg_[0];
157  normal = vec_t({-1, 0, 0});
158  domain.addBoundaryNode(point, LEFT, normal);
159  point[0] = end_[0];
160  normal = vec_t({1, 0, 0});
161  domain.addBoundaryNode(point, RIGHT, normal);
162  }
163  }
164  step = (end_[2] - beg_[2]) / (counts[2] - 1); // z
165  for (int i = 1; i < counts[2]-1; ++i) {
166  point[2] = beg_[2] + i*step;
167  point[0] = beg_[0];
168  point[1] = beg_[1];
169  normal = vec_t({-1, -1, 0}).normalized();
170  domain.addBoundaryNode(point, FRONT, normal);
171  point[0] = beg_[0];
172  point[1] = end_[1];
173  normal = vec_t({-1, 1, 0}).normalized();
174  domain.addBoundaryNode(point, LEFT, normal);
175  point[0] = end_[0];
176  point[1] = beg_[1];
177  normal = vec_t({1, -1, 0}).normalized();
178  domain.addBoundaryNode(point, RIGHT, normal);
179  point[0] = end_[0];
180  point[1] = end_[1];
181  normal = vec_t({1, 1, 0}).normalized();
182  domain.addBoundaryNode(point, BACK, normal);
183  }
184  domain.addBoundaryNode(beg_, BOTTOM, vec_t({-1, -1, -1}).normalized());
185  domain.addBoundaryNode(end_, TOP, vec_t({1, 1, 1}).normalized());
186  domain.addBoundaryNode(vec_t({end_[0], beg_[1], beg_[2]}), BOTTOM,
187  vec_t({1, -1, -1}).normalized());
188  domain.addBoundaryNode(vec_t({beg_[0], end_[1], end_[2]}), TOP,
189  vec_t({-1, 1, 1}).normalized());
190  domain.addBoundaryNode(vec_t({end_[0], end_[1], beg_[2]}), BACK,
191  vec_t({1, 1, -1}).normalized());
192  domain.addBoundaryNode(vec_t({beg_[0], end_[1], beg_[2]}), LEFT,
193  vec_t({-1, 1, -1}).normalized());
194  domain.addBoundaryNode(vec_t({beg_[0], beg_[1], end_[2]}), FRONT,
195  vec_t({-1, -1, 1}).normalized());
196  domain.addBoundaryNode(vec_t({end_[0], beg_[1], end_[2]}), RIGHT,
197  vec_t({1, -1, 1}).normalized());
198  }
199  return domain;
200 }
201 
202 template <typename vec_t>
203 std::ostream& BoxShape<vec_t>::print(std::ostream& os) const {
204  return os << "BoxShape(" << beg_.transpose() << ", " << end_.transpose() << ")";
205 }
206 
207 template <typename vec_t>
208 DomainDiscretization<vec_t>
209 BoxShape<vec_t>::discretizeBoundaryWithDensity(const std::function<scalar_t(vec_t)>& dr,
210  int type) const {
211  static_assert(0 <= dim && dim <= 3, "Only dimensions up to 3 are supported.");
212  DomainDiscretization<vec_t> domain(*this);
213  if (dim == 1) {
214  domain.addBoundaryNode(beg_, type == 0 ? -1 : type, vec_t(-1.0));
215  domain.addBoundaryNode(end_, type == 0 ? -2 : type, vec_t(1.0));
216  } else if (dim == 2) {
217  const int TOP = (type == 0) ? -4 : type;
218  const int RIGHT = (type == 0) ? -2 : type;
219  const int BOTTOM = (type == 0) ? -3 : type;
220  const int LEFT = (type == 0) ? -1 : type;
221 
222  vec_t p1 = beg_;
223  vec_t p2 = {beg_[0], end_[1]};
224  vec_t p3 = end_;
225  vec_t p4 = {end_[0], beg_[1]};
226 
227  for (const auto& p : discretization_helpers::discretizeLineWithDensity(p1, p2, dr)) {
228  domain.addBoundaryNode(p, LEFT, {-1, 0});
229  }
230  for (const auto& p : discretization_helpers::discretizeLineWithDensity(p2, p3, dr)) {
231  domain.addBoundaryNode(p, TOP, {0, 1});
232  }
233  for (const auto& p : discretization_helpers::discretizeLineWithDensity(p3, p4, dr)) {
234  domain.addBoundaryNode(p, RIGHT, {1, 0});
235  }
236  for (const auto& p : discretization_helpers::discretizeLineWithDensity(p4, p1, dr)) {
237  domain.addBoundaryNode(p, BOTTOM, {0, -1});
238  }
239 
240  domain.addBoundaryNode(p1, LEFT, vec_t({-1, -1}).normalized());
241  domain.addBoundaryNode(p2, TOP, vec_t({-1, 1}).normalized());
242  domain.addBoundaryNode(p3, RIGHT, vec_t({1, 1}).normalized());
243  domain.addBoundaryNode(p4, BOTTOM, vec_t({1, -1}).normalized());
244  } else if (dim == 3) {
245  int LEFT = -1, RIGHT = -2, FRONT = -3, BACK = -4, BOTTOM = -5, TOP = -6;
246  if (type != 0) LEFT = RIGHT = FRONT = BACK = BOTTOM = TOP = type;
247  Vec<scalar_t, 2> side_beg, side_end;
248  GeneralFill<Vec<scalar_t, 2>> fill;
249  fill.seed(0).proximityTolerance(0.99);
250  {
251  side_beg[0] = beg_[1]; side_beg[1] = beg_[2];
252  side_end[0] = end_[1]; side_end[1] = end_[2];
253  BoxShape<Vec<scalar_t, 2>> side(side_beg, side_end);
254  auto dr2 = [&](const Vec<scalar_t, 2>& p) { return dr({beg_[0], p[0], p[1]}); };
255  auto domain2d = side.discretizeWithDensity(dr2, fill);
256  vec_t pos; pos[0] = beg_[0];
257  for (int i : domain2d.interior()) {
258  pos[1] = domain2d.pos(i, 0); pos[2] = domain2d.pos(i, 1);
259  domain.addBoundaryNode(pos, LEFT, {-1, 0, 0});
260  }
261  auto dr22 = [&](const Vec<scalar_t, 2>& p) { return dr({end_[0], p[0], p[1]}); };
262  domain2d = side.discretizeWithDensity(dr22, fill);
263  pos[0] = end_[0];
264  for (int i : domain2d.interior()) {
265  pos[1] = domain2d.pos(i, 0); pos[2] = domain2d.pos(i, 1);
266  domain.addBoundaryNode(pos, RIGHT, {1, 0, 0});
267  }
268  }
269  {
270  side_beg[0] = beg_[0]; side_beg[1] = beg_[2];
271  side_end[0] = end_[0]; side_end[1] = end_[2];
272  BoxShape<Vec<scalar_t, 2>> side(side_beg, side_end);
273  auto dr2 = [&](const Vec<scalar_t, 2>& p) { return dr({p[0], beg_[1], p[1]}); };
274  auto domain2d = side.discretizeWithDensity(dr2, fill);
275  vec_t pos; pos[1] = beg_[1];
276  for (int i : domain2d.interior()) {
277  pos[0] = domain2d.pos(i, 0); pos[2] = domain2d.pos(i, 1);
278  domain.addBoundaryNode(pos, FRONT, {0, -1, 0});
279  }
280  auto dr22 = [&](const Vec<scalar_t, 2>& p) { return dr({p[0], end_[1], p[1]}); };
281  domain2d = side.discretizeWithDensity(dr22, fill);
282  pos[1] = end_[1];
283  for (int i : domain2d.interior()) {
284  pos[0] = domain2d.pos(i, 0); pos[2] = domain2d.pos(i, 1);
285  domain.addBoundaryNode(pos, BACK, {0, 1, 0});
286  }
287  }
288  {
289  side_beg[0] = beg_[0]; side_beg[1] = beg_[1];
290  side_end[0] = end_[0]; side_end[1] = end_[1];
291  BoxShape<Vec<scalar_t, 2>> side(side_beg, side_end);
292  auto dr2 = [&](const Vec<scalar_t, 2>& p) { return dr({p[0], p[1], beg_[2]}); };
293  auto domain2d = side.discretizeWithDensity(dr2, fill);
294  vec_t pos; pos[2] = beg_[2];
295  for (int i : domain2d.interior()) {
296  pos[0] = domain2d.pos(i, 0); pos[1] = domain2d.pos(i, 1);
297  domain.addBoundaryNode(pos, BOTTOM, {0, 0, -1});
298  }
299  auto dr22 = [&](const Vec<scalar_t, 2>& p) { return dr({p[0], p[1], end_[2]}); };
300  domain2d = side.discretizeWithDensity(dr22, fill);
301  pos[2] = end_[2];
302  for (int i : domain2d.interior()) {
303  pos[0] = domain2d.pos(i, 0); pos[1] = domain2d.pos(i, 1);
304  domain.addBoundaryNode(pos, TOP, {0, 0, 1});
305  }
306  }
307  }
308  return domain;
309 }
310 
311 template <typename vec_t>
312 DomainDiscretization<vec_t> BoxShape<vec_t>::discretizeWithStep(
313  scalar_t step, int internal_type, int boundary_type) const {
314  Vec<int, dim> counts;
315  for (int i = 0; i < dim; ++i) counts[i] = iceil((end_[i] - beg_[i]) / step) + 1;
316  return discretize(counts, internal_type, boundary_type);
317 }
318 
319 template <typename vec_t>
320 DomainDiscretization<vec_t> BoxShape<vec_t>::discretize(
321  const Vec<int, dim>& counts, int internal_type, int boundary_type) const {
322  auto domain = discretizeBoundary(counts, boundary_type);
323  Vec<int, dim> counts_internal = counts - Vec<int, dim>::Constant(2);
324  if (internal_type == 0) internal_type = 1;
325  for (const vec_t& p : linspace(beg_, end_, counts_internal, false)) {
326  domain.addInternalNode(p, internal_type);
327  }
328  return domain;
329 }
330 
331 } // namespace mm
332 
333 #endif // MEDUSA_BITS_DOMAINS_BOXSHAPE_HPP_
mm::DomainShape::margin_
scalar_t margin_
Tolerance for the geometric operation of the domain.
Definition: DomainShape_fwd.hpp:64
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::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
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.
mm::DomainShape::dim
@ dim
Dimensionality of the domain.
Definition: DomainShape_fwd.hpp:57
numutils.hpp
BoxShape_fwd.hpp
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::BoxShape::discretizeBoundary
DomainDiscretization< vec_t > discretizeBoundary(const Vec< int, dim > &counts, int type=0) const
Uniformly discretizes underlying domain boundary.
Definition: BoxShape.hpp:49
DomainShape.hpp
mm::BoxShape::BoxShape
BoxShape(vec_t beg, vec_t end)
Constructs a n dimensional box shaped domain defined by two n dimensional points - vertices at ends o...
Definition: BoxShape.hpp:21
mm::BoxShape::end_
vec_t end_
Larger of the edge points.
Definition: BoxShape_fwd.hpp:30
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::BoxShape::contains
bool contains(const vec_t &point) const override
Return true if point is not more than margin() outside the domain.
Definition: BoxShape.hpp:31
mm::BoxShape::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: BoxShape.hpp:209
mm::BoxShape
Class for working with box shaped domains.
Definition: BoxShape_fwd.hpp:28
mm::BoxShape::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: BoxShape.hpp:312
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::discretization_helpers::discretizeLineWithDensity
Range< vec_t > discretizeLineWithDensity(const vec_t &p, const vec_t &q, const func_t &delta_r)
Returns nodes lying on the line segment pq with approximate distances delta_r.
Definition: discretization_helpers.hpp:35
mm::BoxShape::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: BoxShape.hpp:41
discretization_helpers.hpp
mm::BoxShape::discretize
DomainDiscretization< vec_t > discretize(const Vec< int, dim > &counts, int internal_type=0, int boundary_type=0) const
Uniformly discretizes underlying domain.
Definition: BoxShape.hpp:320
Vec.hpp
mm::BoxShape::print
std::ostream & print(std::ostream &os) const override
Output information about this shape to given output stream os.
Definition: BoxShape.hpp:203
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
mm::BoxShape::beg_
vec_t beg_
Smaller of the edge points.
Definition: BoxShape_fwd.hpp:29