Medusa  1.1
Coordinate Free Mehless Method implementation
NURBSShape.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
2 #define MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
3 
9 #include "NURBSShape_fwd.hpp"
10 #include "GeneralSurfaceFill.hpp"
11 
12 #include <medusa/bits/domains/DomainDiscretization.hpp>
13 #include <medusa/bits/utils/assert.hpp>
14 #include <medusa/bits/spatial_search/KDTreeMutable.hpp>
15 #include <medusa/bits/domains/NURBSPatch.hpp>
16 #include <medusa/bits/domains/BoxShape.hpp>
17 #include <medusa/bits/domains/compute_normal.hpp>
18 
19 namespace mm {
20 
21 template <typename vec_t, typename param_vec_t>
22 NURBSShape<vec_t, param_vec_t>::NURBSShape(const Range<NURBSPatch<vec_t, param_vec_t>>& patches_in)
23  : patches { patches_in }, seed_(get_seed()) {
24  if (param_dim == 1) {
25  n_samples = 2;
26  }
27  for (NURBSPatch<vec_t, param_vec_t>& patch : patches) {
28  patch.computeDerivativeStructure();
29  }
30 }
31 
32 template <typename vec_t, typename param_vec_t>
33 NURBSShape<vec_t, param_vec_t>::NURBSShape(Range<NURBSPatch<vec_t, param_vec_t>>&& patches_in)
34  : patches { patches_in }, seed_(get_seed()) {
35  if (param_dim == 1) {
36  n_samples = 2;
37  }
38  for (NURBSPatch<vec_t, param_vec_t> &patch : patches) {
39  patch.computeDerivativeStructure();
40  }
41 }
42 
43 template <typename vec_t, typename param_vec_t>
44 DomainDiscretization<vec_t> NURBSShape<vec_t, param_vec_t>::discretizeBoundaryWithDensity(
45  const std::function<scalar_t(vec_t)>& h, int type) const {
46  return nurbs_shape_internal::NURBSShapeHelper<vec_t, param_vec_t>::
47  discretizeBoundaryWithDensity(*this, h, type);
48 }
49 
50 template <typename vec_t, typename param_vec_t>
51 NURBSShape<vec_t, param_vec_t>& NURBSShape<vec_t, param_vec_t>::proximityTolerance(
52  scalar_t zeta) {
53  assert_msg((0 < zeta && zeta < 1), "Zeta must be between 0 and 1, got %f.", zeta);
54  this->zeta = zeta;
55  return *this;
56 }
57 
58 template <typename vec_t, typename param_vec_t>
59 NURBSShape<vec_t, param_vec_t>& NURBSShape<vec_t, param_vec_t>::boundaryProximity(
60  scalar_t epsilon) {
61  assert_msg((0 < epsilon && epsilon < 1), "Epsilon must be between 0 and 1, got %f.", epsilon);
62  this->epsilon = epsilon;
63  return *this;
64 }
65 
67 template <typename vec_t>
68 struct nurbs_shape_internal::NURBSShapeHelper<vec_t, Vec<typename vec_t::scalar_t, 2>>{
69  typedef Vec<typename vec_t::scalar_t, 2> param_vec_t;
70  typedef typename vec_t::scalar_t scalar_t;
71  enum { dim = vec_t::dim,
72  proj_dim = dim + 1,
73  param_dim = param_vec_t::dim};
74  typedef Vec<scalar_t, proj_dim> proj_vec_t;
75 
76  static DomainDiscretization<vec_t> discretizeBoundaryWithDensity(
77  const NURBSShape<vec_t, param_vec_t>& shape,
78  const std::function<scalar_t(vec_t)>& h, int type) {
79  if (type == 0) type = -1;
80  std::mt19937 gen(shape.seed_);
81 
82  GeneralSurfaceFill<vec_t, param_vec_t> gsf;
83  gsf.seed(shape.seed_).maxPoints(shape.max_points).proximityTolerance(shape.zeta)
84  .numSamples(shape.n_samples);
85 
86  // Global domain and KDTree.
87  DomainDiscretization<vec_t> domain(shape);
88  KDTreeMutable<vec_t> tree_global;
89 
90  // Discretize every patch.
91  for (const NURBSPatch<vec_t, param_vec_t>& patch : shape.patches) {
92  BoxShape<param_vec_t> bs = patch.getDomain();
93  DomainDiscretization<param_vec_t> local_param_d(bs);
94 
95  // Discretize boundaries first.
96  Range<NURBSPatch<vec_t, Vec<scalar_t, 1>>> boundaries = patch.getBoundaries();
97  int i = -1;
98  for (NURBSPatch<vec_t, Vec<scalar_t, 1>>& boundary : boundaries) {
99  i++;
100  boundary.computeDerivativeStructure();
101 
102  KDTreeMutable<vec_t> tree_local; // Local KD-Tree should be per-boundary.
103  // Generate seed node.
104  BoxShape<Vec<scalar_t, 1>> boundary_bs = boundary.getDomain();
105  Vec<scalar_t, 1> t_min = boundary_bs.beg();
106  Vec<scalar_t, 1> t_max = boundary_bs.end();
107 
108  std::uniform_real_distribution<> dis(t_min(0), t_max(0));
109  Vec<scalar_t, 1> t_seed = Vec<scalar_t, 1>(dis(gen));
110 
111  // Insert into local structures.
112  param_vec_t param_seed = patch.getPatchParameterFromBoundaryParameter(t_seed,
113  i, shape.epsilon);
114  scalar_t w_seed;
115  vec_t pt_seed = boundary.evaluate(t_seed, &w_seed);
116  local_param_d.addInternalNode(param_seed, 1);
117  tree_local.insert(pt_seed);
118 
119  // Insert into global structures.
120  scalar_t check_radius_seed = h(pt_seed);
121  scalar_t d_sq = tree_global.size() == 0 ? 10 * check_radius_seed * check_radius_seed
122  : tree_global.query(pt_seed).second[0];
123  if (d_sq >= (shape.zeta * check_radius_seed) * (shape.zeta * check_radius_seed)) {
124  vec_t normal_seed = surface_fill_internal::
125  compute_normal(patch.jacobian(param_seed, pt_seed, w_seed));
126  domain.addBoundaryNode(pt_seed, type, normal_seed);
127  tree_global.insert(pt_seed);
128  }
129 
130  // General surface fill the boundary.
131  int cur_node = local_param_d.size() - 1;
132  int end_node = local_param_d.size();
133  while (cur_node < end_node && end_node < shape.max_points) {
134  param_vec_t param = local_param_d.pos(cur_node);
135  Vec<scalar_t, 1> t = patch.getBoundaryParameterFromPatchParameter(param, i);
136  vec_t pt, der;
137  std::tie(pt, der) = boundary.evaluatePointAndJacobian(t);
138 
139  std::array<Vec<scalar_t, 1>, 2> candidates{Vec<scalar_t, 1>(1),
140  Vec<scalar_t, 1>(-1)};
141 
142  // Filter candidates.
143  scalar_t alpha = h(pt) / der.norm();
144  for (const auto& u_cand : candidates) {
145  Vec<scalar_t, 1> t_new = t + alpha * u_cand; // Generate candidate.
146  param_vec_t param_new = patch.getPatchParameterFromBoundaryParameter(
147  t_new, i, shape.epsilon);
148  if (!local_param_d.contains(param_new)) continue;
149 
150  scalar_t w_new;
151  vec_t pt_new = boundary.evaluate(t_new, &w_new);
152 
153  // Check radius must be new radius.
154  scalar_t check_radius = (pt - pt_new).norm();
155 
156  // Check for local insertion.
157  scalar_t d_sq_loc = tree_local.query(pt_new).second[0];
158  if (d_sq_loc >= (shape.zeta * check_radius) *
159  (shape.zeta * check_radius)) {
160  local_param_d.addInternalNode(param_new, 1);
161  tree_local.insert(pt_new);
162  end_node++;
163 
164  // Check for global insertion.
165  scalar_t d_sq_global = tree_global.query(pt_new).second[0];
166  if (d_sq_global >= (shape.zeta * check_radius) *
167  (shape.zeta * check_radius)) {
168  vec_t normal_new = surface_fill_internal::
169  compute_normal(patch.jacobian(param_new, pt_new, w_new));
170  domain.addBoundaryNode(pt_new, type, normal_new);
171  tree_global.insert(pt_new);
172  }
173  }
174  }
175  cur_node++;
176  }
177  }
178 
179  // Random another parameter seed around the middle of the domain, so that fill is
180  // successful even if patches overlap.
181  std::uniform_real_distribution<> dis(0.0, 1.0);
182  param_vec_t param_another = bs.beg() + 0.25 * (bs.end() - bs.beg()) +
183  dis(gen) * 0.25 * (bs.end() - bs.beg());
184  local_param_d.addInternalNode(param_another, 1);
185 
186  scalar_t w_another;
187  vec_t pt_another = patch.evaluate(param_another, &w_another);
188  scalar_t check_radius_another = h(pt_another);
189  scalar_t d_sq_another = tree_global.query(pt_another).second[0];
190  if (d_sq_another >= (shape.zeta * check_radius_another) *
191  (shape.zeta * check_radius_another)) {
192  vec_t normal_another = surface_fill_internal::
193  compute_normal(patch.jacobian(param_another, pt_another, w_another));
194  domain.addBoundaryNode(pt_another, type, normal_another);
195  tree_global.insert(pt_another);
196  }
197 
198  // Surface fill.
199  auto param = [&patch](const param_vec_t& t) {
200  return patch.evaluatePointAndJacobian(t);
201  };
202  gsf.fillParametrization(domain, local_param_d, param, h, tree_global, type);
203  }
204 
205  return domain;
206  }
207 };
208 
209 template <typename vec_t>
210 struct nurbs_shape_internal::NURBSShapeHelper<vec_t, Vec<typename vec_t::scalar_t, 1>>{
211  typedef Vec<typename vec_t::scalar_t, 1> param_vec_t;
212  typedef typename vec_t::scalar_t scalar_t;
213  enum { dim = vec_t::dim,
214  proj_dim = dim + 1,
215  param_dim = param_vec_t::dim};
216  typedef Vec<scalar_t, proj_dim> proj_vec_t;
217 
218  static DomainDiscretization<vec_t> discretizeBoundaryWithDensity(
219  const NURBSShape<vec_t, param_vec_t>& shape,
220  const std::function<scalar_t(vec_t)>& h, int type) {
221  if (type == 0) type = -1;
222  std::mt19937 gen(shape.seed_);
223 
224  GeneralSurfaceFill<vec_t, param_vec_t> gsf;
225  gsf.seed(shape.seed_).maxPoints(shape.max_points).proximityTolerance(shape.zeta)
226  .numSamples(shape.n_samples);
227 
228  // Global domain and KDTree.
229  DomainDiscretization<vec_t> domain(shape);
230  KDTreeMutable<vec_t> tree_global;
231 
232  // Discretize every patch.
233  for (const NURBSPatch<vec_t, param_vec_t>& patch : shape.patches) {
234  BoxShape<param_vec_t> bs = patch.getDomain();
235  DomainDiscretization<param_vec_t> local_param_d(bs);
236 
237  // Insert boundary edges and random parameter as seeds.
238  std::array<param_vec_t, 2> param_seeds{
239  patch.getPatchParameterFromBoundaryParameter({}, 0, shape.epsilon),
240  patch.getPatchParameterFromBoundaryParameter({}, 1, shape.epsilon)};
241  // Another random seed could be generated here to make this algorithm also work on
242  // overlapping curves. This would, however, create two "gaps" of size between h
243  // and 2h where the fronts meet, instead of just one.
244 
245  for (const Vec1d& param_seed : param_seeds) {
246  scalar_t w_seed;
247  vec_t pt_seed = patch.evaluate(param_seed, &w_seed);
248  local_param_d.addInternalNode(param_seed, 1);
249 
250  // Insert into global structure.
251  scalar_t check_radius_seed = h(pt_seed);
252  scalar_t d_sq = tree_global.size() == 0 ? 10 * check_radius_seed * check_radius_seed
253  : tree_global.query(pt_seed).second[0];
254  if (d_sq >= (shape.zeta * check_radius_seed) * (shape.zeta * check_radius_seed)) {
255  vec_t normal_seed = surface_fill_internal::
256  compute_normal(patch.jacobian(param_seed, pt_seed, w_seed));
257  domain.addBoundaryNode(pt_seed, type, normal_seed);
258  tree_global.insert(pt_seed);
259  }
260  }
261 
262  // Surface fill.
263  auto param = [&patch](const param_vec_t& t) {
264  return patch.evaluatePointAndJacobian(t);
265  };
266  gsf.fillParametrization(domain, local_param_d, param, h, tree_global, type);
267  }
268 
269  return domain;
270  }
271 };
273 
274 } // namespace mm
275 
276 #endif // MEDUSA_BITS_DOMAINS_NURBSSHAPE_HPP_
compute_normal.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::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
mm::NURBSShape::proximityTolerance
NURBSShape & proximityTolerance(scalar_t zeta)
Set proximity tolerance in surface fill.
Definition: NURBSShape.hpp:51
mm::nurbs_shape_internal::NURBSShapeHelper
Internal structure of NURBSShape that helps with partial class specialization.
Definition: NURBSShape_fwd.hpp:35
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
GeneralSurfaceFill.hpp
mm::NURBSShape::discretizeBoundaryWithDensity
DomainDiscretization< vec_t > discretizeBoundaryWithDensity(const std::function< scalar_t(vec_t)> &, int) const override
Discretizes boundary with given density and fill engine.
Definition: NURBSShape.hpp:44
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
NURBSShape_fwd.hpp
mm::NURBSPatch
Class representing a single NURBS patch in an arbitrary dimensional space, defined on an arbitrary di...
Definition: NURBSPatch_fwd.hpp:72
KDTreeMutable.hpp
NURBSPatch.hpp
mm::Vec1d
Vec< double, 1 > Vec1d
Convenience typedef for 1d vector of doubles.
Definition: Vec_fwd.hpp:33
mm::get_seed
unsigned int get_seed()
Return a random seed.
Definition: randutils.hpp:25
assert.hpp
mm::NURBSShape
Class representing a shape made out of NURBS patches in an arbitrary dimensional space.
Definition: NURBSShape_fwd.hpp:52
mm::surface_fill_internal::compute_normal
Vec< scalar_t, dim_to > compute_normal(Eigen::Matrix< scalar_t, dim_to, dim_from > J)
Compute normal from given Jacobian.
mm::NURBSShape::NURBSShape
NURBSShape(const Range< NURBSPatch< vec_t, param_vec_t >> &patches_in)
Construct NURBSShape from a Range of patches.
Definition: NURBSShape.hpp:22
mm::NURBSShape::boundaryProximity
NURBSShape & boundaryProximity(scalar_t epsilon)
Calculate boundary normals epsilon times the size of domain away from the boundary.
Definition: NURBSShape.hpp:59
mm::Range
An extension of std::vector<T> to support additional useful operations.
Definition: Range_fwd.hpp:30
BoxShape.hpp
mm::DomainShape::scalar_t
vec_t::Scalar scalar_t
Scalar data type used in computation.
Definition: DomainShape_fwd.hpp:55