Medusa  1.1
Coordinate Free Mehless Method implementation
DomainShape.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_DOMAINSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_DOMAINSHAPE_HPP_
3 
9 #include "DomainShape_fwd.hpp"
10 #include "ShapeDifference.hpp"
11 #include "ShapeUnion.hpp"
12 #include "TranslatedShape.hpp"
13 #include "RotatedShape.hpp"
16 #include <cmath>
17 
18 namespace mm {
19 
20 template <typename vec_t>
22  const vec_t& point, const vec_t& unit_normal) const {
23  assert_msg(unit_normal.norm() > 1e-9, "Normal %s given for point %s is zero.",
24  unit_normal, point);
25  if (dim == 1) return {false, point};
26  // Find point on the other side of the boundary with exponential search
27  vec_t start = point;
28  vec_t normal = unit_normal;
29  bool is_inside = contains(start);
30  scalar_t max_stretch = 100 * margin_;
31  while (true) {
32  if (contains(start + max_stretch * normal) != is_inside) break;
33  if (contains(start - max_stretch * normal) != is_inside) {
34  max_stretch *= -1;
35  break;
36  }
37  max_stretch *= 2;
38  if (std::isinf(max_stretch)) { // hint is bad
39  return {false, vec_t()};
40  }
41  }
42 
43  // Find the point on the boundary using bisection
44  scalar_t stretch = max_stretch;
45  while (std::abs(stretch) > margin_) {
46  stretch /= 2;
47  if (contains(start + stretch * normal) == is_inside) {
48  start += stretch * normal;
49  }
50  }
51 
52  // Make unit normal point outside
53  if (is_inside) {
54  normal *= signum(max_stretch);
55  } else {
56  normal *= -signum(max_stretch);
57  }
58 
59  // Make sure point is inside
60  while (!contains(start)) start -= margin_ * normal;
61 
62  return {true, start};
63 }
64 
65 template <typename vec_t>
67  return ShapeUnion<vec_t>(*this, other);
68 }
69 
70 template <typename vec_t>
72  return ShapeDifference<vec_t>(*this, other);
73 }
74 
75 template <typename vec_t>
77 DomainShape<vec_t>::discretizeWithDensity(const std::function<scalar_t(vec_t)>&, int, int) const {
78  assert_msg(false, "This domain does not support filling with density, "
79  "use a general fill engine instead.");
80  return DomainDiscretization<vec_t>(*this);
81 }
82 
83 template <typename vec_t>
85  return TranslatedShape<vec_t>(*this, a);
86 }
87 
88 template <typename vec_t>
90  assert_msg(dim == 2, "Angle rotation only available in 2D.");
91  scalar_t s = std::sin(angle);
92  scalar_t c = std::cos(angle);
93  Eigen::Matrix<double, dim, dim> Q; Q << c, -s, s, c;
94  return rotate(Q);
95 }
96 
97 template <typename vec_t>
98 RotatedShape<vec_t> DomainShape<vec_t>::rotate(const Eigen::Matrix<scalar_t, dim, dim>& Q) {
99  return RotatedShape<vec_t>(*this, Q);
100 }
101 
102 
103 } // namespace mm
104 
105 #endif // MEDUSA_BITS_DOMAINS_DOMAINSHAPE_HPP_
mm::DomainShape::add
ShapeUnion< vec_t > add(const DomainShape &other) const
Returns a shape representing a union of *this and other.
Definition: DomainShape.hpp:66
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::ShapeUnion
Class representing a union of two shapes.
Definition: ShapeUnion_fwd.hpp:28
mm::DomainShape::rotate
RotatedShape< vec_t > rotate(const Eigen::Matrix< scalar_t, dim, dim > &Q)
Transform the shape by given orthogonal matrix Q.
Definition: DomainShape.hpp:98
RotatedShape.hpp
mm::TranslatedShape
Class for working with translated domain shapes.
Definition: TranslatedShape_fwd.hpp:27
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::RotatedShape
Class for working with rotated (or mirrored) domain shapes.
Definition: RotatedShape_fwd.hpp:27
numutils.hpp
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::DomainShape
Base class for geometric shapes of domains.
Definition: DomainShape_fwd.hpp:52
mm::DomainShape::translate
TranslatedShape< vec_t > translate(const vec_t &a)
Translate the shape by given vector a.
Definition: DomainShape.hpp:84
ShapeUnion.hpp
TranslatedShape.hpp
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::DomainShape::subtract
ShapeDifference< vec_t > subtract(const DomainShape &other) const
Returns a shape representing a difference of *this and other.
Definition: DomainShape.hpp:71
mm::DomainShape::projectPointToBoundary
virtual std::pair< bool, vec_t > projectPointToBoundary(const vec_t &point, const vec_t &unit_normal) const
Project point to boundary using bisection along the line define by unit_normal.
Definition: DomainShape.hpp:21
mm::ShapeDifference
A class representing a set-difference of two shapes.
Definition: ShapeDifference_fwd.hpp:26
mm::signum
constexpr int signum(T x, std::false_type)
Signum overload for unsigned types.
Definition: numutils.hpp:79
mm::DomainShape::scalar_t
vec_t::Scalar scalar_t
Scalar data type used in computation.
Definition: DomainShape_fwd.hpp:55
DomainShape_fwd.hpp
ShapeDifference.hpp