Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_locate.cc
Go to the documentation of this file.
1
21//
22// given x, search K in mesh such that x in K
23//
24// author: Pierre.Saramito@imag.fr
25//
26// date: 12 march 2012
27//
28// implementation note:
29// use CGAL::Segment_tree for element location
30//
31#include "rheolef/geo_locate.h"
32#include "rheolef/geo.h"
33#include "rheolef/geo_element_contains.h"
34#include "rheolef/point_util.h"
35
36// internal includes:
37#ifdef _RHEOLEF_HAVE_CGAL
38#include "rheolef/cgal_traits.h"
39#include <CGAL/Segment_tree_k.h>
40#include <CGAL/Range_segment_tree_traits.h>
41#endif // _RHEOLEF_HAVE_CGAL
42
43namespace rheolef {
44
45// ---------------------------------------------------------------------
46// 1) utility: bbox of a geo_element
47// ---------------------------------------------------------------------
48template <class T, class M>
49void
51{
52 // TODO: omega.order == 1 only: K.size == K.n_node
54 size_type d = omega.dimension();
55 for (size_type j = 0; j < d; j++) {
56 xmin[j] = std::numeric_limits<T>::max();
57 xmax[j] = -std::numeric_limits<T>::max();
58 }
59 for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
60 size_type dis_inod = K[iloc];
61 const point_basic<T>& x = omega.dis_node(dis_inod);
62 for (size_type j = 0; j < d; j++) {
63 xmin[j] = std::min(x[j], xmin[j]);
64 xmax[j] = std::max(x[j], xmax[j]);
65 }
66 }
67}
68#ifdef _RHEOLEF_HAVE_CGAL
69// -----------------------------------------------------------------------------------
70// 2) missing in cgal: inspirated from "examples/RangeSegmentTrees/include/Tree_Traits.h"
71// -----------------------------------------------------------------------------------
72template <class Kernel, class Val>
73class Segment_tree_map_traits_1 {
74 public:
75 typedef Val Value;
76 typedef typename Kernel::FT Point_1 ;
77 typedef Point_1 Key;
78 typedef Point_1 Key_1;
79 typedef std::pair<Key,Key> Pure_interval;
80 typedef std::pair<Pure_interval, Val> Interval;
81
82 class C_Low_1{
83 public:
84 Key_1 operator()(const Interval& i)
85 { return i.first.first;}
86 };
87
88 class C_High_1{
89 public:
90 Key_1 operator()(const Interval& i)
91 { return i.first.second;}
92 };
93
94 class C_Key_1{
95 public:
96 Key_1 operator()(const Key& k)
97 { return k;}
98 };
99
100 class C_Compare_1{
101 public:
102 bool operator()(Key_1 k1, Key_1 k2)
103 {
104 return std::less<Float>()(k1,k2);
105 }
106 };
107
108 typedef C_Compare_1 compare_1;
109 typedef C_Low_1 low_1;
110 typedef C_High_1 high_1;
111 typedef C_Key_1 key_1;
112};
113// -----------------------------------------------------------------------------------
114// 3) helper a generic dimension=D implementation with cgal
115// -----------------------------------------------------------------------------------
116
117template <class T, size_t D> struct cgal_locate_traits {};
118
119template <class T>
120struct cgal_locate_traits<T,1> {
121
122// typedef:
123
124 typedef typename geo_base_rep<T,sequential>::size_type size_type;
125
126 typedef typename geo_cgal_traits<T,1>::Kernel Kernel;
127 typedef Segment_tree_map_traits_1<Kernel,size_type> Traits;
128 typedef ::CGAL::Segment_tree_1<Traits> Segment_tree_type;
129 typedef typename Traits::Interval Interval;
130 typedef typename Traits::Pure_interval Pure_interval;
131 typedef typename Traits::Key Key;
132
133 static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
134 {
135 return Pure_interval(Key(xmin[0]),
136 Key(xmax[0]));
137 }
138 static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
139 {
140 return Pure_interval (Key(x[0]-eps),
141 Key(x[0]+eps));
142 }
143 static void put (std::ostream& out, const Pure_interval& b) {
144 out << "[" << b.first << "," << b.second << "[";
145 }
146};
147template <class T>
148struct cgal_locate_traits<T,2> {
149
150// typedef:
151
152 typedef typename geo_base_rep<T,sequential>::size_type size_type;
153
154 typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
155 typedef ::CGAL::Segment_tree_map_traits_2<Kernel,size_type> Traits;
156 typedef ::CGAL::Segment_tree_2<Traits> Segment_tree_type;
157 typedef typename Traits::Interval Interval;
158 typedef typename Traits::Pure_interval Pure_interval;
159 typedef typename Traits::Key Key;
160
161 static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
162 {
163 return Pure_interval(Key(xmin[0], xmin[1]),
164 Key(xmax[0], xmax[1]));
165 }
166 static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
167 {
168 return Pure_interval (Key(x[0]-eps, x[1]-eps),
169 Key(x[0]+eps, x[1]+eps));
170 }
171 static void put (std::ostream& out, const Pure_interval& b) {
172 out << "[" << b.first.x() << "," << b.second.x() << "[x["
173 << b.first.y() << "," << b.second.y() << "[";
174 }
175};
176template <class T>
177struct cgal_locate_traits<T,3> {
178
179// typedef:
180
181 typedef typename geo_base_rep<T,sequential>::size_type size_type;
182
183 typedef typename geo_cgal_traits<T,3>::Kernel Kernel;
184 typedef ::CGAL::Segment_tree_map_traits_3<Kernel,size_type> Traits;
185 typedef ::CGAL::Segment_tree_3<Traits> Segment_tree_type;
186 typedef typename Traits::Interval Interval;
187 typedef typename Traits::Pure_interval Pure_interval;
188 typedef typename Traits::Key Key;
189
190 static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
191 {
192 return Pure_interval(Key(xmin[0], xmin[1], xmin[2]),
193 Key(xmax[0], xmax[1], xmax[2]));
194 }
195 static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
196 {
197 return Pure_interval (Key(x[0]-eps, x[1]-eps, x[2]-eps),
198 Key(x[0]+eps, x[1]+eps, x[2]+eps));
199 }
200 static void put (std::ostream& out, const Pure_interval& b) {
201 out << "[" << b.first.x() << "," << b.second.x() << "[x["
202 << b.first.y() << "," << b.second.y() << "[x["
203 << b.first.z() << "," << b.second.z() << "[";
204 }
205};
206// ---------------------------------------------------------------------
207// 4) the geo_locate interface
208// ---------------------------------------------------------------------
209template <class T, class M>
211{
212 if (_ptr != 0) {
213 delete_macro(_ptr);
214 }
215}
216template <class T, class M>
219 const geo_base_rep<T,M>& omega,
220 const point_basic<T>& x,
221 size_type dis_ie_guest) const
222{
223 if (_ptr == 0) { _ptr = make_ptr(omega); }
224 return _ptr->seq_locate (omega, x, dis_ie_guest);
225}
226template <class T, class M>
229 const geo_base_rep<T,M>& omega,
230 const point_basic<T>& x,
231 size_type dis_ie_guest) const
232{
233 if (_ptr == 0) { _ptr = make_ptr(omega); }
234 return _ptr->dis_locate (omega, x, dis_ie_guest);
235}
236// ---------------------------------------------------------------------
237// 5) the tree box abstract data structure
238// ---------------------------------------------------------------------
239template <class T, class M>
240class geo_locate_abstract_rep {
241public:
242 typedef typename disarray<T,M>::size_type size_type;
243 virtual ~geo_locate_abstract_rep() {}
244 virtual void initialize (const geo_base_rep<T,M>& omega) const = 0;
245 virtual size_type seq_locate (
246 const geo_base_rep<T,M>& omega,
247 const point_basic<T>& x,
248 size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
249 virtual size_type dis_locate (
250 const geo_base_rep<T,M>& omega,
251 const point_basic<T>& x,
252 size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
253};
254// ---------------------------------------------------------------------
255// 5) the tree box concrete, and dimension dependent, data structure
256// = cgal "segment" tree
257// ---------------------------------------------------------------------
258template <class T, class M, size_t D>
259class geo_locate_rep : public geo_locate_abstract_rep<T,M> {
260public:
261
262// typedef:
263
265
266 typedef typename cgal_locate_traits<T,D>::Segment_tree_type Segment_tree_type;
267 typedef typename cgal_locate_traits<T,D>::Interval Interval;
268
269// allocators:
270
271 geo_locate_rep() : _tree() {}
272 geo_locate_rep(const geo_base_rep<T,M>& omega) : _tree() { initialize(omega); }
273 ~geo_locate_rep() {}
274 void initialize (const geo_base_rep<T,M>& omega) const;
275
276// accessors:
277
278 size_type seq_locate (
279 const geo_base_rep<T,M>& omega,
280 const point_basic<T>& x,
281 size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const;
282 size_type dis_locate (
283 const geo_base_rep<T,M>& omega,
284 const point_basic<T>& x,
285 size_type dis_ie_guest = std::numeric_limits<size_type>::max()) const;
286
287// data:
288protected:
289 mutable Segment_tree_type _tree; // cgal window query requires a mutable _tree (why?)
290};
291// ---------------------------------------------------------------------
292// 5a) allocator
293// ---------------------------------------------------------------------
294template <class T, class M>
295geo_locate_abstract_rep<T,M>*
297{
298 check_macro (omega.dimension() == omega.map_dimension(), "geo_locate: map_dim < dim: not supported");
299 switch (omega.dimension()) {
300 case 1: return new_macro((geo_locate_rep<T,M,1>)(omega));
301 case 2: return new_macro((geo_locate_rep<T,M,2>)(omega));
302 case 3: return new_macro((geo_locate_rep<T,M,3>)(omega));
303 default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
304 }
305}
306// ---------------------------------------------------------------------
307// 5b) initialize
308// ---------------------------------------------------------------------
309template <class T, class M, size_t D>
310void
311geo_locate_rep<T,M,D>::initialize (const geo_base_rep<T,M>& omega) const
312{
313 // create the corresponding vector of bounding boxes
314 trace_macro ("geo::locate initialize...");
315 std::list<Interval> boxes;
316 point_basic<T> xmin, xmax;
317 for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
318 const geo_element& K = omega[ie];
319 compute_bbox (omega, K, xmin, xmax);
320 boxes.push_back (Interval(cgal_locate_traits<T,D>::make_cgal_bbox (xmin,xmax),ie));
321 }
322 // create tree:
323 // _tree.clear(); // TODO : when e.g. nodes update, rebuild boxes & tree
324 _tree.make_tree (boxes.begin(), boxes.end());
325 trace_macro ("geo::locate initialize done");
326}
327// ---------------------------------------------------------------------
328// 5c) locate
329// ---------------------------------------------------------------------
330template <class T, class M, size_t D>
331typename geo_locate_rep<T,M,D>::size_type
332geo_locate_rep<T,M,D>::seq_locate (
333 const geo_base_rep<T,M>& omega,
334 const point_basic<T>& x,
335 size_type dis_ie_guest) const
336{
337 if (dis_ie_guest != std::numeric_limits<size_type>::max()) {
338 // have a guest:
339 const distributor& ownership = omega.ownership();
340 if (ownership.is_owned (dis_ie_guest)) {
341 size_type first_dis_ie = ownership.first_index();
342 size_type ie_guest = dis_ie_guest - first_dis_ie;
343 const geo_element& K_guest = omega[ie_guest];
344 bool intersect = details::contains (K_guest, omega.get_nodes(), x);
345 if (intersect) {
346 return K_guest.dis_ie();
347 }
348 }
349 }
351 size_type dis_ie = std::numeric_limits<size_type>::max();
352 // epsilon is only for bbox: then, predicate on element K is exact
353 // TODO: compute epsilon with omega.hmin scale ?
354 // static const T eps = 1e5*std::numeric_limits<T>::epsilon();
355 // static const T eps = 1000*std::numeric_limits<T>::epsilon();
356 static const T eps = sqrt(std::numeric_limits<T>::epsilon());
357 Interval xe = Interval (cgal_locate_traits<T,D>::make_cgal_point_window (x, eps), 0);
358 std::list<Interval> intersected_boxes;
359 // point query = inverse range query ; from ::CGAL documentation:
360 // "In order to perform an inverse range query, a range query of epsilon width has to be performed.
361 // We preferred not to offer an extra function for this sort of query, since the inverse range
362 // query is a special case of the range query (window_query)"
363 _tree.window_query (xe, std::back_inserter(intersected_boxes));
364 for (typename std::list<Interval>::iterator j = intersected_boxes.begin(); j != intersected_boxes.end(); j++) {
365 size_type ie = (*j).second;
366 const geo_element& K = omega[ie];
367 bool intersect = details::contains (K, omega.get_nodes(), x);
368 if (intersect) {
369 dis_ie = K.dis_ie();
370 break;
371 }
372 }
373 return dis_ie;
374}
375// ---------------------------------------------------------------------
376// 5c) dis_locate = one distributed computation
377// ---------------------------------------------------------------------
378template <class T, class M, size_t D>
379typename geo_locate_rep<T,M,D>::size_type
380geo_locate_rep<T,M,D>::dis_locate (
381 const geo_base_rep<T,M>& omega,
382 const point_basic<T>& x,
383 size_type dis_ie_guest) const
384{
385 // mpi::reduce with std::minimum handles unsigned=size_t as signed=long
386 // so use here long to reduce, otherwise leads to a bug
387 using signed_size_type = long;
388 const signed_size_type signed_large = -1;
389 size_type dis_ie = seq_locate (omega, x, dis_ie_guest);
390#ifdef _RHEOLEF_HAVE_MPI
391 if (omega.comm().size() > 1 && is_distributed<M>::value) {
392 // fix for numeric_limits<size_t> == -1 with long int
393 dis_ie = mpi::all_reduce (omega.comm(), signed_size_type(dis_ie), mpi::maximum<signed_size_type>());
394 }
395#endif // _RHEOLEF_HAVE_MPI
396 return dis_ie;
397}
398// ---------------------------------------------------------------------
399// 6) geo::locate()
400// ---------------------------------------------------------------------
401// 6a) scalar case:
402// ----------------
403template <class T, class M>
406 const point_basic<T>& x,
407 size_type dis_ie_guest) const
408{
409 return _locator.seq_locate (*this, x, dis_ie_guest);
410}
411template <class T, class M>
414 const point_basic<T>& x,
415 size_type dis_ie_guest) const
416{
417 return _locator.dis_locate (*this, x, dis_ie_guest);
418}
419// ----------------
420// 6b) disarray case:
421// ----------------
422template <class T>
423void
427 bool do_check) const
428{
429 dis_ie.resize (x.ownership());
430 for (size_type i = 0, n = x.size(); i < n; i++) {
431 dis_ie[i] = base::_locator.seq_locate (*this, x[i], dis_ie[i]);
432 if (do_check) {
433 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
434 "locate: failed at x="<<ptos(x[i],base::_dimension));
435 }
436 }
437}
438#ifdef _RHEOLEF_HAVE_MPI
439template <class T>
440void
444 // dis_ie(ix) can contains a guest for x in K=omega[dis_ie(ix)]
445 bool do_check) const
446{
447 // mpi::reduce with std::minimum handles unsigned=size_t as signed=long
448 // so use here long to reduce, otherwise leads to a bug
449 using signed_size_type = long;
450 const signed_size_type signed_large = -1;
451 const size_type large = std::numeric_limits<size_type>::max();
452 const T infty = std::numeric_limits<T>::max();
453 // 1) scan x and locate into the local mesh partition ; list when failed
454 distributor ownership = x.ownership();
455 if (dis_ie.ownership() != ownership) dis_ie.resize (ownership);
456 size_type first_dis_i = ownership.first_index();
457 std::list<id_pt_t<T> > failed;
458 for (size_type i = 0, n = x.size(); i < n; i++) {
459 dis_ie[i] = base::_locator.seq_locate (*this, x[i], dis_ie[i]);
460 if (dis_ie[i] == large) {
461 size_type dis_i = first_dis_i + i;
462 failed.push_back (id_pt_t<T>(dis_i, x[i]));
463 }
464 }
465 // 2) merge the failed list into a massive disarray, then distributed to all
466 communicator comm = ownership.comm();
467 distributor fld_ownership (distributor::decide, comm, failed.size());
468 size_type fld_dis_size = fld_ownership.dis_size();
469 if (comm.size() == 1 || fld_dis_size == 0) {
470 // no unsolved, on any procs: nothing to do !
471trace_macro ("locate done(1)");
472 return;
473 }
474 size_type first_fld_dis_i = fld_ownership.first_index();
475 size_type last_fld_dis_i = fld_ownership. last_index();
476 id_pt_t<T> unset (large, point_basic<T>(infty,infty,infty));
477 std::vector<id_pt_t<T> > massive_failed (fld_dis_size, unset);
478 typename std::list<id_pt_t<T> >::iterator iter = failed.begin();
479 for (size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
480 massive_failed [fld_dis_i] = *iter;
481 }
482 std::vector<id_pt_t<T> > massive_query (fld_dis_size, unset);
483 mpi::all_reduce (
484 comm,
485 massive_failed.begin().operator->(),
486 massive_failed.size(),
487 massive_query.begin().operator->(),
489
490 // 3) run the locator on all failed points ON ALL PROCS, skipping local queries (already failed)
491 std::vector<signed_size_type> massive_result (fld_dis_size, signed_large);
492 // 3a) range [0:first[
493 for (size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
494 massive_result [fld_dis_i] = base::_locator.seq_locate (*this, massive_query[fld_dis_i].second);
495 }
496 // 3b) range [last,dis_size[
497 for (size_type fld_dis_i = last_fld_dis_i; fld_dis_i < fld_dis_size; ++fld_dis_i) {
498 massive_result [fld_dis_i] = base::_locator.seq_locate (*this, massive_query[fld_dis_i].second);
499 }
500 // 4) send & merge the results to all
501 std::vector<signed_size_type> massive_merged (fld_dis_size, signed_large);
502 mpi::all_reduce (
503 comm,
504 massive_result.begin().operator->(),
505 massive_result.size(),
506 massive_merged.begin().operator->(),
507 mpi::maximum<signed_size_type>());
508
509 // 5) store the local range into the distributed disarray dis_ie:
510 for (size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
511 size_type dis_i = massive_query [fld_dis_i].first;
512 check_macro (dis_i >= first_dis_i, "invalid index");
513 size_type i = dis_i - first_dis_i;
514 dis_ie[i] = massive_merged [fld_dis_i];
515 if (do_check) {
516 check_macro (dis_ie[i] != large,
517 "dis_locate: failed at x="<<ptos(x[i],base::_dimension));
518 }
519 }
520trace_macro ("locate done(2)");
521}
522#endif // _RHEOLEF_HAVE_MPI
523// ----------------------------------------------------------------------------
524// no CGAL: no locators
525// ----------------------------------------------------------------------------
526#else // _RHEOLEF_HAVE_CGAL
527template <class T, class M>
529{
530}
531template <class T, class M>
534 const geo_base_rep<T,M>& omega,
535 const point_basic<T>& x,
536 size_type dis_ie_guest) const
537{
538 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
539 return 0;
540}
541template <class T, class M>
544 const geo_base_rep<T,M>& omega,
545 const point_basic<T>& x,
546 size_type dis_ie_guest) const
547{
548 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
549 return 0;
550}
551template <class T, class M>
554 const point_basic<T>& x,
555 size_type dis_ie_guest) const
556{
557 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
558 return 0;
559}
560template <class T, class M>
563 const point_basic<T>& x,
564 size_type dis_ie_guest) const
565{
566 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
567 return 0;
568}
569template <class T>
570void
571geo_rep<T,sequential>::locate (
572 const disarray<point_basic<T>, sequential>& x,
573 disarray<size_type, sequential>& dis_ie,
574 bool do_check) const
575{
576 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
577}
578#ifdef _RHEOLEF_HAVE_MPI
579template <class T>
580void
581geo_rep<T,distributed>::locate (
582 const disarray<point_basic<T>, distributed>& x,
583 disarray<size_type, distributed>& dis_ie,
584 // dis_ie(ix) can contains a guest for x in K=omega[dis_ie(ix)]
585 bool do_check) const
586{
587 fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
588}
589#endif // _RHEOLEF_HAVE_MPI
590#endif // _RHEOLEF_HAVE_CGAL
591// ----------------------------------------------------------------------------
592// instanciation in library
593// ----------------------------------------------------------------------------
594#define _RHEOLEF_instanciation(T,M) \
595template class geo_locate<Float,M>; \
596template class geo_base_rep<Float,M>; \
597template class geo_rep<Float,M>;
598
600#ifdef _RHEOLEF_HAVE_MPI
602#endif // _RHEOLEF_HAVE_MPI
603
604} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
rep::base::size_type size_type
Definition disarray.h:501
see the distributor page for the full documentation
Definition distributor.h:69
size_type dis_size() const
global and local sizes
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
static const size_type decide
Definition distributor.h:83
const communicator_type & comm() const
base class for M=sequential or distributed meshes representations
Definition geo.h:528
size_type dis_locate(const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
const node_type & dis_node(size_type dis_inod) const
Definition geo.h:589
size_type map_dimension() const
Definition geo.h:564
size_type seq_locate(const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
base::size_type size_type
Definition geo.h:533
size_type dimension() const
Definition geo.h:563
size_type size(size_type dim) const
Definition geo.cc:456
see the geo_element page for the full documentation
size_type size() const
size_type seq_locate(const geo_base_rep< T, M > &omega, const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
size_type dis_locate(const geo_base_rep< T, M > &omega, const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
static geo_locate_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
disarray< T, distributed >::size_type size_type
Definition geo_locate.h:42
base::size_type size_type
Definition geo.h:934
base::size_type size_type
Definition geo.h:791
sequential mesh representation
Definition geo.h:778
#define trace_macro(message)
Definition dis_macros.h:111
#define error_macro(message)
Definition dis_macros.h:49
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
bool contains(const geo_element &K, const disarray< point_basic< T >, M > &node, const point_basic< T > &x)
This file is part of Rheolef.
void compute_bbox(const geo_base_rep< T, M > &omega, const geo_element &K, point_basic< T > &xmin, point_basic< T > &xmax)
Definition geo_locate.cc:50
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
t operator()(const t &a, const t &b)
Definition space.cc:386
static const bool value