GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <exploragram/optimal_transport/optimal_transport_3d.h>
41			#include <exploragram/optimal_transport/linear_least_squares.h>
42
43			#include <geogram/mesh/mesh_io.h>
44			#include <geogram/mesh/mesh_reorder.h>
45			#include <geogram/mesh/mesh_geometry.h>
46			#include <geogram/mesh/mesh_AABB.h>
47			#include <geogram/voronoi/CVT.h>
48			#include <geogram/voronoi/generic_RVD_vertex.h>
49			#include <geogram/voronoi/RVD_callback.h>
50			#include <geogram/voronoi/generic_RVD_cell.h>
51			#include <geogram/delaunay/delaunay_3d.h>
52			#include <geogram/points/nn_search.h>
53			#include <geogram/basic/stopwatch.h>
54			#include <geogram/basic/progress.h>
55
56			#include <stack>
57			#include <iterator>
58
59			namespace {
60			using namespace GEO;
61
62			// Implementation notes on OTMPolyhedronCallback:
63			// [Bruno Levy, Sept. 13th, 2017]
64			//
65			// OTMPolyhedronCallback computes the value, gradient
66			// and Hessian of the objective function of semi-discrete
67			// optimal transport.
68			//
69			// Its operator() is called for each pair of:
70			// (Laguerre cell of vertex v, background mesh tetrahedron t)
71			// that have a non-empty intersection. The intersection is
72			// represented as a GEOGen::ConvexCell object, that stores all
73			// the geometric and combinatorial information associated with
74			// the intersection.
75			//
76			// The objective function, gradient and Hessian depend on integrals
77			// computed on the Laguerre cells (objective function and gradient)
78			// and Laguerre cells boundaries (Hessian). The integrals are
79			// integrals of linear functions over the
80			// (Laguerre cell /\ background tetrahedron) intersections, evaluated in
81			// closed form (integral of linear function over a polytope). Volumetric
82			// integrals over the cells are evaluated by decomposing the cell into
83			// tetrahedra. Surfacic integrals over the cells boundary are evaluated
84			// by decomposing the cells facets into triangles.
85			//
86			// Each (Laguerre cell /\ background tetrahedron) is represented
87			// by a GEOGen::ConvexCell object, represented in dual form, that is its
88			// vertices correspond to facets, and its
89			// triangles to vertices. This is because computing the intersection
90			// between halfspaces can be done very efficiently with this form.
91			// The code looks more complicated than it should be, due to the dual
92			// representation of GEOGen::ConvexCell. Two nested loops iterate on the
93			// facets (i.e., the dual vertices) and around each facet (i.e. the loop
94			// of dual triangles incident to each dual vertex). It's a bit painful,
95			// but remember, computing the intersections in this form is much much
96			// faster.
97			//
98			// The value of the objective function is computed by
99			// eval_F() and eval_F_weighted() (if varying density is attached
100			// to the background mesh as a vertex attribute named "weight").
101			// The objective function is the most painful to compute, but it
102			// is not needed by the Newton solver (only the hierarchical BFGS
103			// solver needs it).
104			//
105			// The gradient of the objective function corrresponds to the
106			// difference between the (possibly weighted) mass (or volume) of
107			// each Laguerre cell and the prescribed masses. The mass (and
108			// optionally centers of mass) is computed by compute_m_and_mg().
109			// The callback operator() only computes the mass of each
110			// Laguerre cell. The prescribed mass is subtracted by the OptimalTransport
111			// class.
112			// The center of mass is used by some applications (for instance
113			// the Euler fluid simulator). The center of mass is weighted by the
114			// mass so that the contribution of different intersections can be
115			// summed (it is divided in the end by the total mass in the
116			// OptimalTransport class).
117			//
118			// The Hessian of the objective function has a non-zero coefficient
119			// hij each time the Laguerre cells of vertex i and vertex j have a
120			// non-empty intersection (a common facet). The value of the coefficient
121			// is the mass (area or weighted area) of the face divided by twice the
122			// distance between vertex i and vertex j. They
123			// are computed by update_Hessian(), that calls OpenNL functions to
124			// construct the (sparse) Hessian. The right-hand side of the linear system
125			// (minus the gradient) is computed by the OptimalTransport class (that
126			// copies it from the gradient).
127			//
128			// Important note/gotcha: Beware the signs of everything:
129			// F = \sum_i \nu_i \psi_i -
130			// \sum_i \int_{Lag_i} \|x - y_i\|^2 - \psi_i u(x) dx
131			// \frac{\partial F}{\partial \psi_i} = \nu_i - \int_{Lag_i} u(x)dx
132			// There is a minus sign in the mass of the Laguerre cells in the gradient
133			// component.
134			// ... but remember, we maximize F, and numerical code prefers to minimize
135			// things, so finally we minimize -F, and there is no minus sign.
136			// ... but again, when doing Newton, we solve \nabla^2 F p = -nabla F, so
137			// finally, there is a minus sign for the components of the gradient when
138			// assembling the RHS for the Newton solves.
139
140			/**
141			* \brief Computes the contribution of a polyhedron
142			* to the objective function minimized by a semi-discrete
143			* optimal transport map.
144			* \details Works in uniform and weighted mode. Works in
145			* Newton and BFGS mode.
146			*/
147			class OTMPolyhedronCallback :
148			public OptimalTransportMap::Callback,
149			public RVDPolyhedronCallback {
150
151			public:
152			/**
153			* \brief OTMPolyhedronCallback constructor.
154			* \param[in] OTM a pointer to the OptimalTransportMap3d
155			*/
156		✗	OTMPolyhedronCallback(OptimalTransportMap3d* OTM) :
157		✗	OptimalTransportMap::Callback(OTM) {
158		✗	}
159
160			/**
161			* \copydoc RVDPolyhedronCallback::operator()
162			*/
163		✗	void operator() (
164			index_t v,
165			index_t t,
166			const GEOGen::ConvexCell& C
167			) const override {
168			// v can be an air particle.
169		✗	if(v >= n_) {
170		✗	return;
171			}
172
173		✗	geo_argused(t);
174
175			double m, mgx, mgy, mgz;
176		✗	compute_m_and_mg(C, m, mgx, mgy, mgz);
177
178		✗	if(spinlocks_ != nullptr) {
179		✗	spinlocks_->acquire_spinlock(v);
180			}
181
182			// +m because we maximize F <=> minimize -F
183		✗	g_[v] += m;
184
185		✗	if(Newton_step_) {
186			// ... but here -m because Newton step =
187			// solve H p = -g (minus g in the RHS).
188		✗	OTM_->add_i_right_hand_side(v,-m);
189			}
190
191		✗	if(mg_ != nullptr) {
192		✗	mg_[3*v] += mgx;
193		✗	mg_[3*v+1] += mgy;
194		✗	mg_[3*v+2] += mgz;
195			}
196
197		✗	if(spinlocks_ != nullptr) {
198		✗	spinlocks_->release_spinlock(v);
199			}
200
201		✗	if(Newton_step_) {
202			// Spinlocks are managed internally by update_Hessian().
203		✗	update_Hessian(C, v);
204			}
205
206		✗	if(eval_F_) {
207		✗	Thread* thread = Thread::current();
208			index_t current_thread_id =
209		✗	(thread == nullptr) ? 0 : thread->id();
210		✗	double F = weighted_ ? eval_F_weighted(C, v) : eval_F(C, v);
211			const_cast<OTMPolyhedronCallback*>(this)->
212		✗	funcval_[current_thread_id] += F;
213			}
214			}
215
216			protected:
217
218			/**
219			* \brief Computes the mass and mass times centroid of the
220			* current ConvexCell.
221			* \details Weights are taken into account if present.
222			* \param[in] C a const reference to the current ConvexCell
223			* \param[out] m , mgx , mgy , mgz the mass and the mass times the
224			* centroid of the ConvexCell. mgx, mgy and mgz are not computed
225			* if mg_ is nullptr.
226			*/
227		✗	void compute_m_and_mg(
228			const GEOGen::ConvexCell& C,
229			double& m, double& mgx, double& mgy, double& mgz
230			) const {
231
232		✗	m = 0.0;
233		✗	mgx = 0.0;
234		✗	mgy = 0.0;
235		✗	mgz = 0.0;
236
237
238			// Find one vertex V0 of the Convex Cell. It will be then decomposed
239			// into tetrahedra radiating from V0.
240		✗	const GEOGen::Vertex* V0 = nullptr;
241		✗	for(index_t ct=0; ct < C.max_t(); ++ct) {
242		✗	if(C.triangle_is_used(ct)) {
243		✗	V0 = &C.triangle_dual(ct);
244		✗	break;
245			}
246			}
247		✗	if(V0 == nullptr) {
248		✗	return;
249			}
250
251
252			// Iterate on all the vertices of the convex cell.
253			// The vertices associated to no triangle are skipped.
254			// Note: the convex cell is in dual form, thus a
255			// vertex of the ConvexCell corresponds to a polygonal
256			// facet.
257		✗	for(index_t cv = 0; cv < C.max_v(); ++cv) {
258		✗	signed_index_t ct = C.vertex_triangle(cv);
259		✗	if(ct == -1) {
260		✗	continue;
261			}
262		✗	geo_debug_assert(C.triangle_is_used(index_t(ct)));
263
264
265			// Iterate around the facet vertices (that correspond
266			// to the current vertex in dual form). The polygonal
267			// facet is triangulated. The triangles that contain
268			// the vertex V0 are skipped (their volume is zero).
269			GEOGen::ConvexCell::Corner first(
270			index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
271		✗	);
272		✗	const GEOGen::Vertex* V1 = &C.triangle_dual(first.t);
273		✗	const GEOGen::Vertex* V2 = nullptr;
274		✗	const GEOGen::Vertex* V3 = nullptr;
275		✗	const double* p0 = V0->point();
276		✗	GEOGen::ConvexCell::Corner c = first;
277			do {
278		✗	V2 = V3;
279		✗	V3 = &C.triangle_dual(c.t);
280		✗	if(
281		✗	V2 != nullptr && V3 != V1 &&
282		✗	V1 != V0 && V2 != V0 && V3 != V0
283			) {
284		✗	const double* p1 = V1->point();
285		✗	const double* p2 = V2->point();
286		✗	const double* p3 = V3->point();
287		✗	double cur_m = GEO::Geom::tetra_volume<3>(p0,p1,p2,p3);
288		✗	if(weighted_) {
289		✗	double w0 = V0->weight();
290		✗	double w1 = V1->weight();
291		✗	double w2 = V2->weight();
292		✗	double w3 = V3->weight();
293		✗	double S = w0+w1+w2+w3;
294		✗	if(mg_ != nullptr) {
295		✗	mgx += 0.25*cur_m*(
296		✗	w0*p0[0]+w1*p1[0]+w2*p2[0]+w3*p3[0]
297			);
298		✗	mgy += 0.25*cur_m*(
299		✗	w0*p0[1]+w1*p1[1]+w2*p2[1]+w3*p3[1]
300			);
301		✗	mgz += 0.25*cur_m*(
302		✗	w0*p0[2]+w1*p1[2]+w2*p2[2]+w3*p3[2]
303			);
304			}
305		✗	cur_m *= (S/4.0);
306			} else {
307		✗	if(mg_ != nullptr) {
308		✗	mgx += 0.25*cur_m*(p0[0]+p1[0]+p2[0]+p3[0]);
309		✗	mgy += 0.25*cur_m*(p0[1]+p1[1]+p2[1]+p3[1]);
310		✗	mgz += 0.25*cur_m*(p0[2]+p1[2]+p2[2]+p3[2]);
311			}
312			}
313		✗	m += cur_m;
314			}
315		✗	C.move_to_next_around_vertex(c);
316		✗	} while(c != first);
317			}
318
319			}
320
321			/**
322			* \brief Updates the Hessian according to the current ConvexCell.
323			* \param[in] C a const reference to the current ConvexCell
324			* \param[in] v the current seed
325			*/
326		✗	void update_Hessian(
327			const GEOGen::ConvexCell& C, index_t v
328			) const {
329			// The coefficient of the Hessian associated to a pair of
330			// adjacent cells Lag(i),Lag(j) is :
331			// - mass(Lag(i) /\ Lag(j)) / (2*distance(pi,pj))
332
333		✗	const double* p0 = OTM_->point_ptr(v);
334
335			// Iterate on all the vertices of the convex cell.
336			// The vertices associated to no triangle are skipped.
337			// Note: the convex cell is in dual form, thus a
338			// vertex of the ConvexCell corresponds to a polygonal
339			// facet.
340
341		✗	for(index_t cv = 0; cv < C.max_v(); ++cv) {
342		✗	signed_index_t ct = C.vertex_triangle(cv);
343		✗	if(ct == -1) {
344		✗	continue;
345			}
346		✗	geo_debug_assert(C.triangle_is_used(index_t(ct)));
347
348			// Get index of adjacent seed if any.
349		✗	index_t v_adj = index_t(-1);
350		✗	signed_index_t adjacent = C.vertex_id(cv);
351		✗	if(adjacent > 0) {
352			// Positive adjacent indices correspond to
353			// Voronoi seed - Voronoi seed link
354		✗	v_adj = index_t(adjacent - 1);
355			}
356		✗	if(v_adj == index_t(-1)) {
357		✗	continue;
358			}
359
360		✗	double hij = 0;
361
362			// Iterate around the facet vertices (that correspond
363			// to the current vertex in dual form). The polygonal
364			// facet is triangulated.
365
366			GEOGen::ConvexCell::Corner first(
367			index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
368		✗	);
369		✗	const GEOGen::Vertex* V1 = &C.triangle_dual(first.t);
370		✗	const GEOGen::Vertex* V2 = nullptr;
371		✗	const GEOGen::Vertex* V3 = nullptr;
372		✗	GEOGen::ConvexCell::Corner c = first;
373			do {
374		✗	V2 = V3;
375		✗	V3 = &C.triangle_dual(c.t);
376		✗	if(V2 != nullptr && V3 != V1) {
377		✗	double cur_m = GEO::Geom::triangle_area_3d(
378			V1->point(),V2->point(),V3->point()
379			);
380		✗	if(weighted_) {
381		✗	cur_m *= (
382		✗	V1->weight()+V2->weight()+V3->weight()
383		✗	)/3.0;
384			}
385		✗	hij += cur_m;
386			}
387		✗	C.move_to_next_around_vertex(c);
388		✗	} while(c != first);
389
390		✗	const double* p1 = OTM_->point_ptr(v_adj);
391		✗	hij /= (2.0 * GEO::Geom::distance(p0,p1,3));
392
393		✗	if(spinlocks_ != nullptr) {
394		✗	spinlocks_->acquire_spinlock(v);
395			}
396
397			// Diagonal is positive, extra-diagonal
398			// coefficients are negative,
399			// this is a convex function.
400
401		✗	if(v_adj < n_) {
402		✗	OTM_->add_ij_coefficient(
403			v, v_adj, -hij
404			);
405			}
406		✗	OTM_->add_ij_coefficient(
407			v, v, hij
408			);
409
410		✗	if(spinlocks_ != nullptr) {
411		✗	spinlocks_->release_spinlock(v);
412			}
413			}
414		✗	}
415
416			/**
417			* \brief Computes the contribution of the current ConvexCell
418			* to the objective function, in the uniform (non-weighted)
419			* case.
420			* \param[in] C a const reference to the current ConvexCell
421			* \param[in] v the current seed
422			*/
423		✗	double eval_F(const GEOGen::ConvexCell& C, index_t v) const {
424
425			// Note: we compute F as a sum of integrals over (signed)
426			// tetrahedra; formed by vertex v and triplets of cell
427			// vertices. The contribution of F is simpler to compute
428			// if v is a vertex of the tetrahedron. Since we use
429			// signed tetrahedra, portions of tetrahedra that are
430			// outside the cell cancel-out ("a-la" Stokes theorem).
431
432		✗	geo_debug_assert(!weighted_);
433
434		✗	double F = 0.0;
435
436			// Iterate on all the vertices of the convex cell.
437			// The vertices associated to no triangle are skipped.
438			// Note: the convex cell is in dual form, thus a
439			// vertex of the ConvexCell corresponds to a polygonal
440			// facet.
441
442		✗	for(index_t cv = 0; cv < C.max_v(); ++cv) {
443		✗	signed_index_t ct = C.vertex_triangle(cv);
444		✗	if(ct == -1) {
445		✗	continue;
446			}
447		✗	geo_debug_assert(C.triangle_is_used(index_t(ct)));
448
449
450			// Iterate around the facet vertices (that correspond
451			// to the current vertex in dual form). The polygonal
452			// facet is triangulated.
453
454			GEOGen::ConvexCell::Corner first(
455			index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
456		✗	);
457		✗	const GEOGen::Vertex* V1 = &C.triangle_dual(first.t);
458		✗	const GEOGen::Vertex* V2 = nullptr;
459		✗	const GEOGen::Vertex* V3 = nullptr;
460		✗	GEOGen::ConvexCell::Corner c = first;
461			do {
462		✗	V2 = V3;
463		✗	V3 = &C.triangle_dual(c.t);
464		✗	if(V2 != nullptr && V3 != V1) {
465		✗	const double* p0 = OTM_->point_ptr(v);
466		✗	const double* p1 = V1->point();
467		✗	const double* p2 = V2->point();
468		✗	const double* p3 = V3->point();
469		✗	double m = Geom::tetra_signed_volume(p0, p1, p2, p3);
470		✗	double fT = 0.0;
471		✗	for(coord_index_t cc = 0; cc < 3; ++cc) {
472		✗	double Uc = p1[cc] - p0[cc];
473		✗	double Vc = p2[cc] - p0[cc];
474		✗	double Wc = p3[cc] - p0[cc];
475		✗	fT +=
476		✗	Uc * Uc +
477		✗	Vc * Vc +
478		✗	Wc * Wc +
479		✗	Uc * Vc +
480		✗	Vc * Wc +
481		✗	Wc * Uc;
482			}
483		✗	fT = m * (fT / 10.0 - w_[v]);
484		✗	F += fT;
485			}
486		✗	C.move_to_next_around_vertex(c);
487		✗	} while(c != first);
488			}
489			// -F because we maximize F <=> minimize -F
490		✗	return -F;
491			}
492
493
494			/**
495			* \brief Computes the contribution of the current ConvexCell
496			* to the objective function, in the weighted case.
497			* \param[in] C a const reference to the current ConvexCell
498			* \param[in] v the current seed
499			*/
500		✗	double eval_F_weighted(const GEOGen::ConvexCell& C, index_t v) const {
501		✗	double F = 0.0;
502
503
504			// Find one vertex V0 of the Convex Cell. It will be then decomposed
505			// into tetrahedra radiating from V0.
506		✗	const GEOGen::Vertex* V0 = nullptr;
507		✗	for(index_t ct=0; ct < C.max_t(); ++ct) {
508		✗	if(C.triangle_is_used(ct)) {
509		✗	V0 = &C.triangle_dual(ct);
510		✗	break;
511			}
512			}
513		✗	if(V0 == nullptr) {
514		✗	return F;
515			}
516
517			// Iterate on all the vertices of the convex cell.
518			// The vertices associated to no triangle are skipped.
519			// Note: the convex cell is in dual form, thus a
520			// vertex of the ConvexCell corresponds to a polygonal
521			// facet.
522		✗	for(index_t cv = 0; cv < C.max_v(); ++cv) {
523		✗	signed_index_t ct = C.vertex_triangle(cv);
524		✗	if(ct == -1) {
525		✗	continue;
526			}
527		✗	geo_debug_assert(C.triangle_is_used(index_t(ct)));
528
529
530			// Iterate around the facet vertices (that correspond
531			// to the current vertex in dual form). The polygonal
532			// facet is triangulated. The triangles that contain
533			// the vertex V0 are skipped (their volume is zero).
534			GEOGen::ConvexCell::Corner first(
535			index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
536		✗	);
537		✗	const GEOGen::Vertex* V1 = &C.triangle_dual(first.t);
538		✗	const GEOGen::Vertex* V2 = nullptr;
539		✗	const GEOGen::Vertex* V3 = nullptr;
540		✗	GEOGen::ConvexCell::Corner c = first;
541			do {
542		✗	V2 = V3;
543		✗	V3 = &C.triangle_dual(c.t);
544		✗	if(
545		✗	V2 != nullptr && V3 != V1 &&
546		✗	V1 != V0 && V2 != V0 && V3 != V0
547			) {
548
549		✗	const double* p0 = V0->point();
550		✗	const double* p1 = V1->point();
551		✗	const double* p2 = V2->point();
552		✗	const double* p3 = V3->point();
553
554		✗	double p0_mass = V0->weight();
555		✗	double p1_mass = V1->weight();
556		✗	double p2_mass = V2->weight();
557		✗	double p3_mass = V3->weight();
558
559		✗	const double* q = OTM_->point_ptr(v);
560
561		✗	double Tvol = GEO::Geom::tetra_volume<3>(p0,p1,p2,p3);
562		✗	double Sp = p0_mass + p1_mass + p2_mass + p3_mass;
563
564		✗	double m = (Tvol * Sp) / 4.0;
565			double rho[4], alpha[4];
566
567		✗	rho[0] = p0_mass;
568		✗	rho[1] = p1_mass;
569		✗	rho[2] = p2_mass;
570		✗	rho[3] = p3_mass;
571
572		✗	alpha[0] = Sp + rho[0];
573		✗	alpha[1] = Sp + rho[1];
574		✗	alpha[2] = Sp + rho[2];
575		✗	alpha[3] = Sp + rho[3];
576
577		✗	double dotprod_00 = 0.0;
578		✗	double dotprod_10 = 0.0;
579		✗	double dotprod_11 = 0.0;
580		✗	double dotprod_20 = 0.0;
581		✗	double dotprod_21 = 0.0;
582		✗	double dotprod_22 = 0.0;
583		✗	double dotprod_30 = 0.0;
584		✗	double dotprod_31 = 0.0;
585		✗	double dotprod_32 = 0.0;
586		✗	double dotprod_33 = 0.0;
587
588		✗	for(coord_index_t cc = 0; cc < 3; ++cc) {
589		✗	double sp0 = q[cc] - p0[cc];
590		✗	double sp1 = q[cc] - p1[cc];
591		✗	double sp2 = q[cc] - p2[cc];
592		✗	double sp3 = q[cc] - p3[cc];
593		✗	dotprod_00 += sp0 * sp0;
594		✗	dotprod_10 += sp1 * sp0;
595		✗	dotprod_11 += sp1 * sp1;
596		✗	dotprod_20 += sp2 * sp0;
597		✗	dotprod_21 += sp2 * sp1;
598		✗	dotprod_22 += sp2 * sp2;
599		✗	dotprod_30 += sp3 * sp0;
600		✗	dotprod_31 += sp3 * sp1;
601		✗	dotprod_32 += sp3 * sp2;
602		✗	dotprod_33 += sp3 * sp3;
603			}
604
605		✗	double fT = 0.0;
606		✗	fT += (alpha[0] + rho[0]) * dotprod_00;
607		✗	fT += (alpha[1] + rho[0]) * dotprod_10;
608		✗	fT += (alpha[1] + rho[1]) * dotprod_11;
609		✗	fT += (alpha[2] + rho[0]) * dotprod_20;
610		✗	fT += (alpha[2] + rho[1]) * dotprod_21;
611		✗	fT += (alpha[2] + rho[2]) * dotprod_22;
612		✗	fT += (alpha[3] + rho[0]) * dotprod_30;
613		✗	fT += (alpha[3] + rho[1]) * dotprod_31;
614		✗	fT += (alpha[3] + rho[2]) * dotprod_32;
615		✗	fT += (alpha[3] + rho[3]) * dotprod_33;
616
617		✗	fT = Tvol * fT / 60.0 - m * w_[v];
618
619		✗	F += fT;
620			}
621		✗	C.move_to_next_around_vertex(c);
622		✗	} while(c != first);
623			}
624			// -F because we maximize F <=> minimize -F
625		✗	return -F;
626			}
627
628			};
629
630			/**********************************************************************/
631
632		✗	void cell_shrink_to_animation(Mesh& mesh) {
633		✗	Attribute<index_t> cell_region;
634		✗	Attribute<double> W(mesh.vertices.attributes(),"w");
635		✗	cell_region.bind_if_is_defined(
636		✗	mesh.cells.attributes(), "region"
637			);
638		✗	if(!cell_region.is_bound()) {
639		✗	Logger::err("OTM") << "region: no such cell attribute"
640		✗	<< std::endl;
641		✗	return;
642			}
643		✗	index_t nb_RVD_cells = 0;
644		✗	for(index_t c: mesh.cells) {
645		✗	nb_RVD_cells = std::max(nb_RVD_cells, cell_region[c]);
646			}
647		✗	++nb_RVD_cells;
648		✗	vector<vec3> center(nb_RVD_cells, vec3(0.0, 0.0, 0.0));
649		✗	vector<index_t> nb(nb_RVD_cells, 0);
650		✗	for(index_t c: mesh.cells) {
651		✗	index_t rvc = cell_region[c];
652		✗	for(index_t lv=0; lv<4; ++lv) {
653		✗	index_t v = mesh.cells.vertex(c,lv);
654		✗	center[rvc] += mesh.vertices.point(v);
655		✗	++nb[rvc];
656			}
657			}
658		✗	for(index_t rvc=0; rvc<nb_RVD_cells; ++rvc) {
659		✗	if(nb[rvc] != 0) {
660		✗	center[rvc] /= double(nb[rvc]);
661			}
662			}
663		✗	for(index_t v: mesh.vertices) {
664		✗	W[v] = mesh.vertices.point_ptr(v)[3];
665			}
666		✗	mesh.vertices.set_dimension(6);
667		✗	for(index_t c: mesh.cells) {
668		✗	index_t rvc = cell_region[c];
669		✗	for(index_t lv=0; lv<4; ++lv) {
670		✗	index_t v = mesh.cells.vertex(c,lv);
671		✗	const vec3& g = center[rvc];
672		✗	double* p = mesh.vertices.point_ptr(v);
673		✗	p[3] = g[0];
674		✗	p[4] = g[1];
675		✗	p[5] = g[2];
676			}
677			}
678		✗	}
679
680			/**********************************************************************/
681
682			/**
683			* \brief Gets the Delaunay implementation that best fits
684			* user desire.
685			* \param[in] user_delaunay the desired implementation.
686			* \return the available implementation nearest to user desire.
687			*/
688		✗	std::string default_delaunay(const std::string& user_delaunay) {
689		✗	std::string result = "PDEL";
690		✗	if(user_delaunay != "default") {
691		✗	result = user_delaunay;
692			}
693			#ifndef GEOGRAM_WITH_PDEL
694		✗	if(result == "PDEL") {
695		✗	result = "BPOW";
696			}
697			#endif
698		✗	return result;
699		✗	}
700
701
702			}
703
704
705
706
707			namespace GEO {
708
709		✗	OptimalTransportMap3d::OptimalTransportMap3d(
710			Mesh* mesh, const std::string& delaunay, bool BRIO
711		✗	) : OptimalTransportMap(
712			3,
713			mesh,
714		✗	default_delaunay(delaunay),
715			BRIO
716		✗	) {
717		✗	callback_ = new OTMPolyhedronCallback(this);
718		✗	total_mass_ = total_mesh_mass();
719		✗	}
720
721		✗	double OptimalTransportMap3d::total_mesh_mass() const {
722		✗	double result = 0.0;
723			// This is terribly confusing, the parameters for
724			// a power diagram are called "weights", and the
725			// standard attribute name for vertices density is
726			// also called "weight" (and is unrelated).
727			// In this program, what is called weight corresponds
728			// to the parameters of the power diagram (except the
729			// name of the attribute), and everything that corresponds
730			// to mass/density is called mass.
731		✗	Attribute<double> vertex_mass;
732		✗	vertex_mass.bind_if_is_defined(
733		✗	mesh_->vertices.attributes(), "weight"
734			);
735
736		✗	for(index_t t: mesh_->cells) {
737		✗	double tet_mass = GEO::Geom::tetra_volume<3>(
738		✗	mesh_->vertices.point_ptr(mesh_->cells.tet_vertex(t, 0)),
739		✗	mesh_->vertices.point_ptr(mesh_->cells.tet_vertex(t, 1)),
740		✗	mesh_->vertices.point_ptr(mesh_->cells.tet_vertex(t, 2)),
741		✗	mesh_->vertices.point_ptr(mesh_->cells.tet_vertex(t, 3))
742			);
743		✗	if(vertex_mass.is_bound()) {
744		✗	tet_mass *= (
745		✗	vertex_mass[mesh_->cells.tet_vertex(t, 0)] +
746		✗	vertex_mass[mesh_->cells.tet_vertex(t, 1)] +
747		✗	vertex_mass[mesh_->cells.tet_vertex(t, 2)] +
748		✗	vertex_mass[mesh_->cells.tet_vertex(t, 3)]
749		✗	) / 4.0;
750			}
751		✗	result += tet_mass;
752			}
753		✗	return result;
754		✗	}
755
756		✗	OptimalTransportMap3d::~OptimalTransportMap3d() {
757		✗	}
758
759		✗	void OptimalTransportMap3d::get_RVD(Mesh& RVD_mesh) {
760		✗	RVD_mesh.clear();
761		✗	Attribute<index_t> tet_region(RVD_mesh.cells.attributes(),"region");
762		✗	RVD()->compute_RVD(
763			RVD_mesh,
764			0, // dim (0 means use default)
765			false, // borders_only
766		✗	show_RVD_seed_ // integration_simplices
767			);
768		✗	if(!show_RVD_seed_) {
769		✗	cell_shrink_to_animation(RVD_mesh);
770			}
771		✗	}
772
773		✗	void OptimalTransportMap3d::compute_Laguerre_centroids(double* centroids) {
774		✗	vector<double> g(nb_points(), 0.0);
775		✗	Memory::clear(centroids, nb_points()*sizeof(double)*3);
776
777		✗	callback_->set_Laguerre_centroids(centroids);
778		✗	callback_->set_g(g.data());
779			{
780		✗	Stopwatch* W = nullptr;
781		✗	if(newton_) {
782		✗	W = new Stopwatch("RVD");
783		✗	Logger::out("OTM") << "In RVD (centroids)..." << std::endl;
784			}
785		✗	RVD_->for_each_polyhedron(
786		✗	*dynamic_cast<RVDPolyhedronCallback*>(callback_),
787			false,false,true
788			);
789		✗	if(newton_) {
790		✗	delete W;
791			}
792			}
793		✗	callback_->set_Laguerre_centroids(nullptr);
794
795		✗	for(index_t v=0; v<nb_points(); ++v) {
796		✗	centroids[3*v ] /= g[v];
797		✗	centroids[3*v+1] /= g[v];
798		✗	centroids[3*v+2] /= g[v];
799			}
800		✗	}
801
802		✗	void OptimalTransportMap3d::call_callback_on_RVD() {
803		✗	RVD_->for_each_polyhedron(
804		✗	*dynamic_cast<RVDPolyhedronCallback*>(callback_),
805			false, // symbolic
806			false, // connected components priority
807		✗	!clip_by_balls_ // parallel
808			);
809			// clip_by_balls deactivates parallel mode, because it needs to access
810			// the PointAllocator of the current thread, and we do not have any
811			// access (for now).
812		✗	}
813
814			/**********************************************************************/
815
816		✗	void compute_Laguerre_centroids_3d(
817			Mesh* omega,
818			index_t nb_points,
819			const double* points,
820			double* centroids,
821			RVDPolyhedronCallback* cb,
822			bool verbose,
823			index_t nb_iter
824			) {
825		✗	omega->vertices.set_dimension(4);
826
827			// false = no BRIO
828			// (OTM does not use multilevel and lets Delaunay
829			// reorder the vertices)
830			OptimalTransportMap3d OTM(
831			omega,
832		✗	"PDEL",
833			false
834		✗	);
835
836		✗	OTM.set_regularization(1e-3);
837		✗	OTM.set_Newton(true);
838		✗	OTM.set_points(nb_points, points);
839		✗	OTM.set_epsilon(0.01);
840		✗	OTM.set_Laguerre_centroids(centroids);
841		✗	OTM.set_verbose(verbose);
842		✗	OTM.optimize(nb_iter);
843
844		✗	if(cb != nullptr) {
845		✗	OTM.RVD()->for_each_polyhedron(*cb,false,false,false);
846			}
847
848		✗	omega->vertices.set_dimension(3);
849		✗	}
850			}
851
852			/****************************************************************************/
853
854			namespace {
855			using namespace GEO;
856
857
858			/**
859			* \brief Gets the number of connected components
860			* of each tetrahedra regions in a mesh.
861			* \param[in] RVD a const reference to the mesh
862			* \param[out] nb_cnx_comps on exit, nb_cnx_comps[r]
863			* contains the number of connected components of
864			* region r.
865			*/
866		✗	void get_nb_connected_components(
867			const Mesh& RVD, vector<index_t>& nb_cnx_comps
868			) {
869		✗	Attribute<index_t> tet_region(RVD.cells.attributes(),"region");
870		✗	vector<bool> marked(RVD.cells.nb(), false);
871		✗	std::stack<index_t> S;
872		✗	for(index_t t = 0; t < RVD.cells.nb(); ++t) {
873		✗	if(!marked[t]) {
874		✗	index_t cur_v = index_t(tet_region[t]);
875		✗	marked[t] = true;
876		✗	S.push(t);
877		✗	while(!S.empty()) {
878		✗	index_t cur_t = S.top();
879		✗	S.pop();
880		✗	for(index_t lf = 0; lf < 4; ++lf) {
881		✗	index_t neigh = RVD.cells.tet_adjacent(cur_t, lf);
882		✗	if(neigh != NO_CELL) {
883		✗	if(
884		✗	tet_region[neigh] == cur_v &&
885		✗	!marked[neigh]
886			) {
887		✗	marked[neigh] = true;
888		✗	S.push(neigh);
889			}
890			}
891			}
892			}
893		✗	if(cur_v >= nb_cnx_comps.size()) {
894		✗	nb_cnx_comps.resize(cur_v,0);
895			}
896		✗	++nb_cnx_comps[cur_v];
897			}
898			}
899		✗	}
900
901
902			/**
903			* \brief Computes the original and final vertices of
904			* the optimal transport.
905			* \param[in] CVT the Centroidal Voronoi Tesselation that
906			* samples the target mesh M2
907			* \param[in] OTM the Optimal Transport Map that back-projects
908			* the samples of the target mesh M2 onto the source mesh M1
909			* \param[out] M1_vertices the source vertices, computed from
910			* the centroids of the power cells restricted to M1
911			* \param[out] M2_vertices the destination vertices, copied
912			* from the sampling of M2
913			* \param[out] nb_cnx_comps gives for each vertex the number of
914			* connected components of the power cell restricted to M1
915			*/
916		✗	void get_vertices(
917			CentroidalVoronoiTesselation& CVT,
918			OptimalTransportMap3d& OTM,
919			vector<vec3>& M1_vertices,
920			vector<vec3>& M2_vertices,
921			vector<index_t>& nb_cnx_comps
922			) {
923		✗	index_t nb_vertices = OTM.RVD()->delaunay()->nb_vertices();
924		✗	M1_vertices.resize(nb_vertices);
925		✗	M2_vertices.resize(nb_vertices);
926		✗	nb_cnx_comps.assign(nb_vertices, 0);
927			{
928		✗	for(index_t v = 0; v < nb_vertices; ++v) {
929		✗	const double* p = CVT.RVD()->delaunay()->vertex_ptr(v);
930		✗	M2_vertices[v] = vec3(p[0], p[1], p[2]);
931			}
932			}
933
934			{
935		✗	vector<vec3> mg(nb_vertices, vec3(0.0, 0.0, 0.0));
936		✗	vector<double> m(nb_vertices, 0.0);
937
938		✗	Mesh RVD;
939		✗	Attribute<index_t> tet_region(RVD.cells.attributes(),"region");
940		✗	OTM.RVD()->compute_RVD(
941			RVD,
942			0, // dim (0 means use default)
943			false, // cells_borders_only
944			true // integration_simplices
945			);
946		✗	RVD.vertices.set_dimension(3);
947		✗	RVD.cells.connect();
948
949			/*
950			if(CmdLine::get_arg_bool("RVD")) {
951			MeshIOFlags flags;
952			flags.set_element(MESH_CELLS);
953			flags.set_attribute(MESH_CELL_REGION);
954			mesh_save(RVD, "RVD.meshb", flags);
955			}
956			*/
957
958		✗	for(index_t t = 0; t < RVD.cells.nb(); ++t) {
959		✗	index_t v = tet_region[t];
960		✗	index_t v0 = RVD.cells.tet_vertex(t, 0);
961		✗	index_t v1 = RVD.cells.tet_vertex(t, 1);
962		✗	index_t v2 = RVD.cells.tet_vertex(t, 2);
963		✗	index_t v3 = RVD.cells.tet_vertex(t, 3);
964		✗	vec3 p0(RVD.vertices.point_ptr(v0));
965		✗	vec3 p1(RVD.vertices.point_ptr(v1));
966		✗	vec3 p2(RVD.vertices.point_ptr(v2));
967		✗	vec3 p3(RVD.vertices.point_ptr(v3));
968		✗	double mt = GEO::Geom::tetra_signed_volume(p0, p1, p2, p3);
969		✗	mg[v] += (mt / 4.0) * (p0 + p1 + p2 + p3);
970		✗	m[v] += mt;
971			}
972		✗	for(index_t v = 0; v < nb_vertices; ++v) {
973		✗	double s = ::fabs(m[v]);
974		✗	if(s != 0.0) {
975		✗	s = 1.0 / s;
976			}
977		✗	M1_vertices[v] = s * mg[v];
978			}
979		✗	get_nb_connected_components(RVD, nb_cnx_comps);
980		✗	}
981		✗	}
982
983			/**
984			* \brief Tests whether a tetrahedron is included inside
985			* a given tetrahedral mesh.
986			* \details Implemented by sampling the tetrahedron and
987			* testing whether each sample is inside the mesh.
988			* \param[in] AABB a const reference to a MeshCellsAABB
989			* \param[in] p1 a const reference to the first vertex
990			* of the tetrahedron
991			* \param[in] p2 a const reference to the second vertex
992			* of the tetrahedron
993			* \param[in] p3 a const reference to the third vertex
994			* of the tetrahedron
995			* \param[in] p4 a const reference to the fourth vertex
996			* of the tetrahedron
997			*/
998		✗	bool mesh_contains_tet(
999			const MeshCellsAABB& AABB,
1000			const vec3& p1,
1001			const vec3& p2,
1002			const vec3& p3,
1003			const vec3& p4
1004			) {
1005		✗	const index_t NB = 5;
1006		✗	for(index_t u1=0; u1<=NB; ++u1) {
1007		✗	double s1 = double(u1)/double(NB);
1008		✗	for(index_t u2=0; u1+u2<=NB; ++u2) {
1009		✗	double s2 = double(u2)/double(NB);
1010		✗	for(index_t u3=0; u1+u2+u3<=NB; ++u3) {
1011		✗	double s3 = double(u3)/double(NB);
1012		✗	index_t u4=NB-u1-u2-u3;
1013		✗	double s4 = double(u4)/double(NB);
1014
1015			// Skip the four vertices of the tetrahedron
1016		✗	if(
1017		✗	(u1 == NB) || (u2 == NB) ||
1018		✗	(u3 == NB) || (u4 == NB)
1019			) {
1020		✗	continue;
1021			}
1022
1023		✗	vec3 g = s1*p1+s2*p2+s3*p3+s4*p4;
1024		✗	if(AABB.containing_tet(g) == MeshCellsAABB::NO_TET) {
1025		✗	return false;
1026			}
1027			}
1028			}
1029			}
1030		✗	return true;
1031			}
1032			}
1033
1034			/****************************************************************************/
1035
1036			namespace GEO {
1037
1038		✗	void compute_morph(
1039			CentroidalVoronoiTesselation& CVT,
1040			OptimalTransportMap3d& OTM,
1041			Mesh& morph,
1042			bool filter_tets
1043			) {
1044		✗	geo_assert(CVT.volumetric());
1045
1046		✗	Logger::out("OTM")
1047		✗	<< "Computing coherent tet mesh"
1048		✗	<< std::endl;
1049
1050		✗	morph.clear();
1051		✗	morph.vertices.set_dimension(6);
1052
1053			// Step 1: Compute the candidate tets from the Delaunay triangulation
1054			// of the samples restricted to M2.
1055		✗	vector<index_t> morph_tets;
1056		✗	morph_tets.clear();
1057			{
1058		✗	vector<double> embedding;
1059		✗	CVT.RVD()->compute_RDT(
1060			morph_tets,
1061			embedding,
1062			RestrictedVoronoiDiagram::RDT_SEEDS_ALWAYS
1063			);
1064		✗	}
1065
1066			// Step 2: Compute the original vertices location (centroids
1067			// of power cells restricted to M1) and get the final vertices
1068			// locations (sampling of M2).
1069		✗	vector<vec3> M1_vertices;
1070		✗	vector<vec3> M2_vertices;
1071		✗	vector<index_t> nb_cnx_comps;
1072		✗	get_vertices(CVT, OTM, M1_vertices, M2_vertices, nb_cnx_comps);
1073		✗	index_t nb_vertices = M1_vertices.size();
1074
1075			// Gather all the points coordinates in a vector of
1076			// 6d points.
1077		✗	vector<double> morph_vertices(nb_vertices*6);
1078		✗	for(index_t v=0; v<nb_vertices; ++v) {
1079		✗	morph_vertices[v*6 ] = M2_vertices[v].x;
1080		✗	morph_vertices[v*6+1] = M2_vertices[v].y;
1081		✗	morph_vertices[v*6+2] = M2_vertices[v].z;
1082		✗	morph_vertices[v*6+3] = M1_vertices[v].x;
1083		✗	morph_vertices[v*6+4] = M1_vertices[v].y;
1084		✗	morph_vertices[v*6+5] = M1_vertices[v].z;
1085			}
1086
1087			// Step 3: Filter-out the tets incident to a vertex
1088			// that splits during transport.
1089		✗	index_t nb_tets = morph_tets.size()/4;
1090		✗	vector<bool> tet_to_remove(nb_tets, false);
1091		✗	for(index_t t=0; t<nb_tets; ++t) {
1092		✗	for(index_t lv=0; lv<4; ++lv) {
1093		✗	index_t v = morph_tets[4*t+lv];
1094		✗	if(nb_cnx_comps[v] > 1) {
1095		✗	tet_to_remove[t] = true;
1096		✗	break;
1097			}
1098			}
1099			}
1100
1101			// Step 4: Filter-out the tets that are not contained by
1102			// the initial mesh M1.
1103		✗	if(filter_tets) {
1104		✗	MeshCellsAABB AABB(*OTM.RVD()->mesh());
1105			try {
1106		✗	ProgressTask progress("Classifying", 100);
1107		✗	for(index_t t=0; t<nb_tets; ++t) {
1108		✗	progress.progress(t * 100 / nb_tets);
1109		✗	if(!tet_to_remove[t]) {
1110		✗	vec3 p[4];
1111		✗	for(index_t lv=0; lv<4; ++lv) {
1112		✗	for(coord_index_t c=0; c<3; ++c) {
1113		✗	index_t v = morph_tets[4*t+lv];
1114		✗	p[lv][c] = morph_vertices[v*6+3+c];
1115			}
1116			}
1117		✗	if(!mesh_contains_tet(AABB, p[0], p[1], p[2], p[3])) {
1118		✗	tet_to_remove[t] = true;
1119			}
1120			}
1121			}
1122		✗	} catch(...) {
1123		✗	}
1124		✗	}
1125
1126			// Step 5: create the output mesh.
1127		✗	vector<index_t> filtered_morph_tets;
1128		✗	for(index_t t=0; t<nb_tets; ++t) {
1129		✗	if(!tet_to_remove[t]) {
1130		✗	filtered_morph_tets.push_back(morph_tets[4*t]);
1131		✗	filtered_morph_tets.push_back(morph_tets[4*t+1]);
1132		✗	filtered_morph_tets.push_back(morph_tets[4*t+2]);
1133		✗	filtered_morph_tets.push_back(morph_tets[4*t+3]);
1134			}
1135			}
1136
1137		✗	morph.cells.assign_tet_mesh(
1138			6, morph_vertices, filtered_morph_tets, true
1139			);
1140		✗	morph.cells.connect();
1141		✗	morph.cells.compute_borders();
1142		✗	}
1143
1144		✗	void compute_singular_surface(
1145			CentroidalVoronoiTesselation& CVT,
1146			OptimalTransportMap3d& OTM,
1147			Mesh& singular
1148			) {
1149
1150		✗	std::set<bindex> edges;
1151			{
1152		✗	vector<index_t> simplices;
1153		✗	vector<double> embedding;
1154		✗	CVT.RVD()->compute_RDT(
1155			simplices,
1156			embedding,
1157			RestrictedVoronoiDiagram::RDT_SEEDS_ALWAYS
1158			);
1159		✗	for(index_t t=0; t*4<simplices.size(); ++t) {
1160		✗	index_t v1 = simplices[t*4];
1161		✗	index_t v2 = simplices[t*4+1];
1162		✗	index_t v3 = simplices[t*4+2];
1163		✗	index_t v4 = simplices[t*4+3];
1164		✗	edges.insert(bindex(v1,v2));
1165		✗	edges.insert(bindex(v1,v3));
1166		✗	edges.insert(bindex(v1,v4));
1167		✗	edges.insert(bindex(v2,v3));
1168		✗	edges.insert(bindex(v2,v4));
1169		✗	edges.insert(bindex(v3,v4));
1170			}
1171		✗	}
1172
1173		✗	Mesh RVD;
1174		✗	MeshIOFlags flags;
1175		✗	flags.set_element(MESH_CELLS);
1176		✗	flags.set_attribute(MESH_CELL_REGION);
1177		✗	OTM.RVD()->compute_RVD(
1178			RVD,
1179			0,
1180			false, // cells_borders_only
1181			true // integration_simplices
1182			);
1183		✗	RVD.vertices.set_dimension(3);
1184		✗	RVD.cells.connect();
1185
1186		✗	singular.clear();
1187		✗	singular.vertices.set_dimension(3);
1188
1189		✗	vector<index_t> triangles;
1190
1191		✗	Attribute<index_t> tet_region(RVD.cells.attributes(),"region");
1192
1193		✗	for(index_t t=0; t<RVD.cells.nb(); ++t) {
1194		✗	index_t v1 = tet_region[t];
1195		✗	for(index_t f=0; f<4; ++f) {
1196		✗	index_t nt = RVD.cells.tet_adjacent(t,f);
1197		✗	if(nt != NO_CELL) {
1198		✗	index_t v2 = tet_region[nt];
1199		✗	if(v1 != v2 && edges.find(bindex(v1,v2)) == edges.end()) {
1200		✗	for(index_t i=0; i<3; ++i) {
1201			index_t lv =
1202		✗	RVD.cells.local_tet_facet_vertex_index(f,i);
1203		✗	index_t v = RVD.cells.tet_vertex(t,lv);
1204		✗	triangles.push_back(v);
1205			}
1206			}
1207			}
1208			}
1209			}
1210
1211		✗	singular.vertices.assign_points(
1212		✗	RVD.vertices.point_ptr(0),
1213			RVD.vertices.dimension(), RVD.vertices.nb()
1214			);
1215		✗	singular.facets.assign_triangle_mesh(
1216			triangles, true
1217			);
1218		✗	}
1219			}
1220