subdivision3.hpp 18 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
* version 0.1                                                                  *
* Copyright (C) 2009, IGG Team, LSIIT, University of Strasbourg                *
*                                                                              *
* This library is free software; you can redistribute it and/or modify it      *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 of the License, or (at your     *
* option) any later version.                                                   *
*                                                                              *
* This library is distributed in the hope that it will be useful, but WITHOUT  *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
* for more details.                                                            *
*                                                                              *
* You should have received a copy of the GNU Lesser General Public License     *
* along with this library; if not, write to the Free Software Foundation,      *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
*                                                                              *
* Web site: https://iggservis.u-strasbg.fr/CGoGN/                              *
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#include "Algo/Geometry/centroid.h"
#include "Algo/Modelisation/subdivision.h"

namespace CGoGN
{

namespace Algo
{

untereiner's avatar
untereiner committed
34
namespace IHM
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
{

template <typename PFP>
void subdivideEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
	assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
	assert(!map.edgeIsSubdivided(d) || !"Trying to subdivide an already subdivided edge") ;

	unsigned int eLevel = map.edgeLevel(d) ;

	unsigned int cur = map.getCurrentLevel() ;
	map.setCurrentLevel(eLevel) ;

	Dart dd = map.phi2(d) ;
	typename PFP::VEC3 p1 = position[d] ;
	typename PFP::VEC3 p2 = position[map.phi1(d)] ;

	map.setCurrentLevel(eLevel + 1) ;

	map.cutEdge(d) ;
	unsigned int eId = map.getEdgeId(d) ;
56 57 58
	map.setEdgeId(map.phi1(d), eId, EDGE_ORBIT) ; //mise a jour de l'id d'arrete sur chaque moitie d'arete
	map.setEdgeId(map.phi1(dd), eId, EDGE_ORBIT) ;

59

60 61
	map.setFaceId(EDGE_ORBIT, d) ; //mise a jour de l'id de face sur chaque brin de chaque moitie d'arete
	map.setFaceId(EDGE_ORBIT, dd) ;
62 63 64 65 66 67 68

	position[map.phi1(d)] = (p1 + p2) * typename PFP::REAL(0.5) ;

	map.setCurrentLevel(cur) ;
}

template <typename PFP>
69
void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType)
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
{
	assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
	assert(!map.faceIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ;

	unsigned int fLevel = map.faceLevel(d) ;
	Dart old = map.faceOldestDart(d) ;

	unsigned int cur = map.getCurrentLevel() ;
	map.setCurrentLevel(fLevel) ;		// go to the level of the face to subdivide its edges

	unsigned int degree = 0 ;
	typename PFP::VEC3 p ;
	Dart it = old ;
	do
	{
		++degree;
		p += position[it] ;
		if(!map.edgeIsSubdivided(it))							// first cut the edges (if they are not already)
untereiner's avatar
untereiner committed
88
			Algo::IHM::subdivideEdge<PFP>(map, it, position) ;	// and compute the degree of the face
89 90 91 92 93 94
		it = map.phi1(it) ;
	} while(it != old) ;
	p /= typename PFP::REAL(degree) ;

	map.setCurrentLevel(fLevel + 1) ;			// go to the next level to perform face subdivision

95
	Dart res;
96

97 98 99 100 101 102 103
	if(degree == 3 && sType ==  S_TRI)		//subdiviser une face triangulaire
	{
		Dart dd = map.phi1(old) ;
		Dart e = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, e) ;					// insert a new edge
		unsigned int id = map.getNewEdgeId() ;
		map.setEdgeId(map.phi_1(dd), id, EDGE_ORBIT) ;		// set the edge id of the inserted edge to the next available id
104

105
		//res= map.phi_1(e); //un brin de la nouvelle face du milieu
106

107 108 109
		unsigned int idface = map.getFaceId(old);
		map.setFaceId(dd, idface, FACE_ORBIT) ;
		map.setFaceId(e, idface, FACE_ORBIT) ;
110

111 112 113 114 115 116 117
		dd = e ;
		e = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, e) ;
		id = map.getNewEdgeId() ;
		map.setEdgeId(map.phi_1(dd), id, EDGE_ORBIT) ;

		map.setFaceId(e, idface, FACE_ORBIT) ;
118

119 120 121
		dd = e ;
		e = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, e) ;
122
		id = map.getNewEdgeId() ;
123
		map.setEdgeId(map.phi_1(dd), id, EDGE_ORBIT) ;
124

125
		map.setFaceId(e, idface, FACE_ORBIT) ;
126
	}
127 128 129 130 131 132 133
	else
	{
		Dart dd = map.phi1(old) ;
		map.splitFace(dd, map.phi1(map.phi1(dd))) ;

		Dart ne = map.phi2(map.phi_1(dd));
		Dart ne2 = map.phi2(ne);
134

135 136 137 138 139
		map.cutEdge(ne) ;
		unsigned int id = map.getNewEdgeId() ;
		map.setEdgeId(ne, id, EDGE_ORBIT) ;
		id = map.getNewEdgeId() ;
		map.setEdgeId(ne2, id, EDGE_ORBIT) ;
140

141
		position[map.phi2(ne)] = p ;
142

143 144 145 146 147 148
		dd = map.phi1(map.phi1(map.phi1(map.phi1(ne)))) ;
		while(dd != ne)
		{
			Dart next = map.phi1(map.phi1(dd)) ;
			map.splitFace(map.phi1(ne), dd) ;
			Dart nne = map.phi2(map.phi_1(dd)) ;
149

150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
			id = map.getNewEdgeId() ;
			map.setEdgeId(nne, id, EDGE_ORBIT) ;

			dd = next ;
		}

		unsigned int idface = map.getFaceId(old);
		Dart e = dd;
		do
		{
			map.setFaceId(dd, idface, DART_ORBIT) ;
			map.setFaceId(map.phi2(dd), idface, DART_ORBIT) ;
			dd = map.phi2(map.phi1(dd));
		}
		while(dd != ne);
	}


	map.setCurrentLevel(cur) ;
169 170 171
}

template <typename PFP>
172
void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType)
173 174 175 176 177 178 179 180 181 182 183 184
{
	assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
	assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ;

	unsigned int vLevel = map.volumeLevel(d) ;
	Dart old = map.volumeOldestDart(d) ;

	unsigned int cur = map.getCurrentLevel() ;
	map.setCurrentLevel(vLevel) ;		// go to the level of the face to subdivide its edges


	/*
185
	 * au niveau du volume courant i
186
	 * stockage d'un brin de chaque face de celui-ci
187
	 * avec calcul du centroid
188 189
	 */

190 191 192 193 194
	DartMarkerStore mf(map);		// Lock a face marker to save one dart per face
	DartMarkerStore mv(map);		// Lock a vertex marker to compute volume center

	typename PFP::VEC3 volCenter;
	unsigned count = 0 ;
195 196 197 198 199 200 201 202 203 204 205

	//Store faces that are traversed and start with the face of d
	std::vector<Dart> visitedFaces;
	visitedFaces.reserve(20);
	visitedFaces.push_back(old);

	//Store the edges before the cutEdge
	std::vector<Dart> oldEdges;
	oldEdges.reserve(20);

	mf.markOrbit(FACE_ORBIT, old) ;
206

207 208 209 210
	for(std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
	{
		Dart e = *face ;

211
		do
212
		{
213 214 215 216 217 218 219 220 221 222 223 224
			// add all face neighbours to the table
			//oldEdges.push_back(e);

			//compute volume centroid
			if(!mv.isMarked(e))
			{
				volCenter += position[e];
				++count;
				mv.markOrbit(VERTEX_ORBIT, e);

				oldEdges.push_back(e);
			}
untereiner's avatar
untereiner committed
225 226


227 228 229 230 231 232
			Dart ee = map.phi2(e) ;
			if(!mf.isMarked(ee)) // not already marked
			{
				visitedFaces.push_back(ee) ;
				mf.markOrbit(FACE_ORBIT, ee) ;
			}
233

234 235 236 237
			e = map.phi1(e) ;
		} while(e != *face) ;
	}

238 239
	volCenter /= typename PFP::REAL(count) ;

untereiner's avatar
untereiner committed
240
	/*
241
	 * Subdivision
untereiner's avatar
untereiner committed
242
	 */
243

untereiner's avatar
untereiner committed
244 245 246
	//Store the darts from quadrangulated faces
	std::vector<Dart> quadfaces;
	quadfaces.reserve(100);
247 248 249 250 251 252 253

	//First step : subdivide edges and faces
	//creates a i+1 edge level and i+1 face level
	for (std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
	{
		Dart d = *face;

254
		//if needed subdivide face
untereiner's avatar
untereiner committed
255
		if(!map.faceIsSubdivided(d))
256
			Algo::IHM::subdivideFace<PFP>(map, d, position, sType);
257

258 259
		//save a dart from the subdivided face
		unsigned int cur = map.getCurrentLevel() ;
untereiner's avatar
untereiner committed
260

261 262
		unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivise
		map.setCurrentLevel(fLevel) ;
untereiner's avatar
untereiner committed
263

264 265 266 267 268 269 270 271 272 273
		Dart cf = map.phi1(d);
		Dart e = cf;
		do
		{
			quadfaces.push_back(e);
			quadfaces.push_back(map.phi2(e));
			e = map.phi2(map.phi1(e));
		}while (e != cf);
		map.setCurrentLevel(cur);
	}
untereiner's avatar
untereiner committed
274

275
	map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision
untereiner's avatar
untereiner committed
276

277 278
	std::vector<Dart> newEdges;	//save darts from inner edges
	newEdges.reserve(100);
279

280 281 282 283
	//Second step : deconnect each corner, close each hole, subdivide each new face
	for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
	{
		Dart e = *edge;
284

285 286 287 288 289 290
		do
		{
			map.unsewFaces(map.phi1(map.phi1(e)));
			map.unsewFaces(map.phi1(e));
			e = map.phi2(map.phi_1(e));
		}while(e != *edge);
untereiner's avatar
untereiner committed
291

292
		map.closeHole(map.phi1(e));
untereiner's avatar
untereiner committed
293

294 295 296
		Dart old = map.phi1(map.phi2(map.phi1(e)));
		Dart dd = map.phi1(old) ;
		map.splitFace(dd, map.phi1(map.phi1(dd))) ;
untereiner's avatar
untereiner committed
297

298 299
		Dart ne = map.phi2(map.phi_1(dd));
		Dart ne2 = map.phi2(ne);
300

301 302
		map.cutEdge(ne) ;
		position[map.phi2(ne)] = volCenter;
303

304 305
		newEdges.push_back(ne);
		newEdges.push_back(map.phi1(ne));
untereiner's avatar
untereiner committed
306

307 308 309 310 311
		dd = map.phi1(map.phi1(map.phi1(map.phi1(ne)))) ;
		while(dd != ne)
		{
			Dart next = map.phi1(map.phi1(dd)) ;
			map.splitFace(map.phi1(ne), dd) ;
untereiner's avatar
untereiner committed
312

313 314
			newEdges.push_back(ne);
			newEdges.push_back(map.phi1(ne));
untereiner's avatar
untereiner committed
315

316 317 318
			dd = next ;
		}
	}
untereiner's avatar
untereiner committed
319

320 321 322 323 324 325 326 327 328 329
	//Third step : 3-sew internal faces
	for (std::vector<Dart>::iterator nvol = quadfaces.begin(); nvol != quadfaces.end(); nvol = nvol + 2)
	{
		if(map.phi3(map.phi2(*nvol)) == map.phi2(*nvol) && map.phi3(map.phi2(*(nvol+1))) == map.phi2(*(nvol+1)))
		{
			//id pour toutes les faces interieures
			map.sewVolumes(map.phi2(*nvol), map.phi2(*(nvol+1)));
			unsigned int idface = map.getNewFaceId();
			map.setFaceId(map.phi2(*nvol),idface, FACE_ORBIT);
		}
untereiner's avatar
untereiner committed
330

331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
		//id pour toutes les aretes exterieurs des faces quadrangulees
		unsigned int idedge = map.getEdgeId(*nvol);
		map.setEdgeId(map.phi2(*nvol), idedge, DART_ORBIT);
		map.setEdgeId( map.phi2(*(nvol+1)), idedge, DART_ORBIT);
	}

	//manque id pour les aretes interieurs : (i.e. 6 pour un hexa)
	DartMarker mne(map);
	for(std::vector<Dart>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it)
	{
		if(!mne.isMarked(*it))
		{
			unsigned int idedge = map.getNewEdgeId();
			map.setEdgeId(*it, idedge, EDGE_ORBIT);
			mne.markOrbit(EDGE_ORBIT,*it);
		}
	}
untereiner's avatar
untereiner committed
348

349 350
	map.setCurrentLevel(cur) ;
}
untereiner's avatar
untereiner committed
351

352 353 354
template <typename PFP>
void catmullClarck(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
untereiner's avatar
untereiner committed
355 356


357
}
untereiner's avatar
untereiner committed
358 359 360



361 362 363 364 365
template <typename PFP>
Dart subdivideFaceTri(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
	assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
	assert(!map.faceIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ;
untereiner's avatar
untereiner committed
366

367 368
	unsigned int fLevel = map.faceLevel(d) ;
	Dart old = map.faceOldestDart(d) ;
untereiner's avatar
untereiner committed
369

370 371
	unsigned int cur = map.getCurrentLevel() ;
	map.setCurrentLevel(fLevel) ;		// go to the level of the face to subdivide its edges
untereiner's avatar
untereiner committed
372

373
	Dart d1 = map.phi1(d);
untereiner's avatar
untereiner committed
374

375
	typename PFP::VEC3 center = Algo::Geometry::faceCentroidGen<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, d, position);	// compute center
untereiner's avatar
untereiner committed
376 377


378 379 380 381 382 383 384 385 386
	map.splitFace(d, d1) ;
	map.cutEdge(map.phi_1(d)) ;
	Dart x = map.phi2(map.phi_1(d)) ;
	Dart dd = map.template phi<111>(x) ;
	while(dd != x)
	{
		Dart next = map.phi1(dd) ;
		map.splitFace(dd, map.phi1(x)) ;
		dd = next ;
387 388
	}

389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
	position[map.phi2(x)] = center;

	map.setCurrentLevel(cur) ;

	return map.phi2(x);	// Return a dart of the central vertex
}

template <typename PFP>
void subdivideVolumeTri(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
	assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
	assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ;

	unsigned int vLevel = map.volumeLevel(d) ;
	Dart old = map.volumeOldestDart(d) ;


	unsigned int cur = map.getCurrentLevel() ;
	map.setCurrentLevel(vLevel) ;		// go to the level of the face to subdivide its edges

untereiner's avatar
untereiner committed
409
	/*
410 411
	 * ou niveau du volume courant i
	 * stockage d'un brin de chaque face de celui-ci
untereiner's avatar
untereiner committed
412
	 */
413

414 415 416 417 418 419
	DartMarkerStore mf(map);		// Lock a face marker

	//Store faces that are traversed and start with the face of d
	std::vector<Dart> visitedFaces;
	visitedFaces.reserve(20);
	visitedFaces.push_back(old);
420

421 422 423
	//Store the edges before the cutEdge
	std::vector<Dart> oldEdges;
	oldEdges.reserve(20);
424

425 426
	mf.markOrbit(FACE_ORBIT, old) ;
	for(std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
427
	{
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442
		Dart e = *face ;

		do	// add all face neighbours to the table
		{

			oldEdges.push_back(e);

			Dart ee = map.phi2(e) ;
			if(!mf.isMarked(ee)) // not already marked
			{
				visitedFaces.push_back(ee) ;
				mf.markOrbit(FACE_ORBIT, ee) ;
			}
			e = map.phi1(e) ;
		} while(e != *face) ;
443 444
	}

445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545
	//First step : subdivide edges and faces
	//creates a i+1 edge level and i+1 face level
	for (std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
	{
		Dart d = *face;

		std::cout << "subdiv de toutes les faces" << std::endl;
		if(!map.faceIsSubdivided(d))
		{
			Dart cf = Algo::IHM::subdivideFaceTri<PFP>(map, *face, position);
		}
	}


//	map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision
//
//	std::vector<Dart> newEdges;	//save darts from inner edges
//	newEdges.reserve(100);
//
//	//Second step : deconnect each corner and close the hole
//	for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
//	{
//		Dart e = *edge;
//
//		do
//		{
//			map.unsewFaces(map.phi1(map.phi1(e)));
//			map.unsewFaces(map.phi1(e));
//			e = map.phi2(map.phi_1(e));
//		}while(e != *edge);
//
//		map.closeHole(map.phi1(e));
//
//		Dart old = map.phi1(map.phi2(map.phi1(e)));
//		Dart dd = map.phi1(old) ;
//		map.splitFace(dd, map.phi1(map.phi1(dd))) ;
//
//		Dart ne = map.phi2(map.phi_1(dd));
//		Dart ne2 = map.phi2(ne);
//
//		map.cutEdge(ne) ;
//		position[map.phi2(ne)] = volCenter;
//
//		newEdges.push_back(ne);
//		newEdges.push_back(map.phi1(ne));
//
//		dd = map.phi1(map.phi1(map.phi1(map.phi1(ne)))) ;
//		while(dd != ne)
//		{
//			Dart next = map.phi1(map.phi1(dd)) ;
//			map.splitFace(map.phi1(ne), dd) ;
//
//			newEdges.push_back(ne);
//			newEdges.push_back(map.phi1(ne));
//
//			dd = next ;
//		}
//	}
//
//	for (std::vector<Dart>::iterator nvol = quadfaces.begin(); nvol != quadfaces.end(); nvol = nvol + 2)
//	{
//		if(map.phi3(map.phi2(*nvol)) == map.phi2(*nvol) && map.phi3(map.phi2(*(nvol+1))) == map.phi2(*(nvol+1)))
//		{
//			//id pour toutes les faces interieures
//			map.sewVolumes(map.phi2(*nvol), map.phi2(*(nvol+1)));
//			unsigned int idface = map.getNewFaceId();
//			map.setFaceId(map.phi2(*nvol),idface, FACE_ORBIT);
//		}
//
//		//id pour toutes les aretes exterieurs des faces quadrangulees
//		unsigned int idedge = map.getEdgeId(*nvol);
//		map.setEdgeId(map.phi2(*nvol), idedge, DART_ORBIT);
//		map.setEdgeId( map.phi2(*(nvol+1)), idedge, DART_ORBIT);
//	}
//
//	//manque id pour les aretes interieurs : 6
//	DartMarker mne(map);
//	for(std::vector<Dart>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it)
//	{
//		if(!mne.isMarked(*it))
//		{
//			unsigned int idedge = map.getNewEdgeId();
//			map.setEdgeId(*it, idedge, EDGE_ORBIT);
//			mne.markOrbit(EDGE_ORBIT,*it);
//			std::cout << "plop" << std::endl;
//		}
//	}
//	std::cout << std::endl;

	map.setCurrentLevel(cur) ;
}


} //namespace IHM

} //namespace Algo

} //namespace CGoGN

/*
 	//Second step : deconnect the corners
546 547
	for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
	{
548 549 550
		Dart e = *edge;

		do
551
		{
552 553 554 555 556 557 558 559 560 561 562 563 564 565
			map.unsewFaces(map.phi1(map.phi1(e)));
			map.unsewFaces(map.phi1(e));
			e = map.phi2(map.phi_1(e));
		}while(e != *edge);

//		map.unsewFaces(map.phi1(*edge));
//		moe.markOrbit(DART_ORBIT,map.phi1(*edge));
//	}
//
//	//Third step : close the hole
//	for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
//	{
//		if(moe.isMarked(map.phi1(*edge)))
//		{
566 567 568
			map.closeHole(map.phi1(*edge));

			Dart h = map.phi2(map.phi1(*edge));
569
			unsigned int degree=0;
570 571
			do
			{
572
				map.setEdgeId(h,map.getEdgeId(map.phi2(h)), DART_ORBIT) ;
573
				h = map.phi1(h);
574
				++degree;
575 576 577
			}while(h != map.phi2(map.phi1(*edge)));


578 579 580
//			moe.unmark(map.phi1(*edge));
//			moe.unmark(map.phi1(map.phi2(map.phi_1(*edge))));
//			moe.unmark(map.phi1(map.phi1(map.phi2(*edge))));
581 582


583 584 585 586 587 588
			//si nouvelle diff de tri alors subdiviser la nouvelle face
			if(degree != 3)
			{
				Dart old = map.phi1(map.phi2(map.phi1(*edge)));
				Dart dd = map.phi1(old) ;
				map.splitFace(dd, map.phi1(map.phi1(dd))) ;
589

590 591
				Dart ne = map.phi2(map.phi_1(dd));
				Dart ne2 = map.phi2(ne);
592

593 594 595
				map.cutEdge(ne) ;
				//unsigned int id = map.getNewEdgeId() ;
				//map.setEdgeId(ne, id, EDGE_ORBIT) ;
596

597 598 599 600
				//id = map.getNewEdgeId() ;
				//map.setEdgeId(ne2, id, EDGE_ORBIT) ;

				position[map.phi2(ne)] = volCenter;
601

602 603 604 605 606 607
				dd = map.phi1(map.phi1(map.phi1(map.phi1(ne)))) ;
				while(dd != ne)
				{
					Dart next = map.phi1(map.phi1(dd)) ;
					map.splitFace(map.phi1(ne), dd) ;
					Dart nne = map.phi2(map.phi_1(dd)) ;
608

609 610 611 612 613
					//id = map.getNewEdgeId() ;
					//map.setEdgeId(nne, id, EDGE_ORBIT) ;

					dd = next ;
				}
614
			}
615
		//}
616 617
	}

618
//	moe.unmarkAll();
619 620 621 622 623 624

	for (std::vector<Dart>::iterator nvol = quadfaces.begin(); nvol != quadfaces.end(); nvol = nvol + 2)
	{
		if(map.phi3(map.phi2(*nvol)) == map.phi2(*nvol) && map.phi3(map.phi2(*(nvol+1))) == map.phi2(*(nvol+1)))
		{
			map.sewVolumes(map.phi2(*nvol), map.phi2(*(nvol+1)));
625 626 627 628 629
			unsigned int idface = map.getNewFaceId();
			map.setFaceId(map.phi2(*nvol),idface, FACE_ORBIT);

			unsigned int idedge = map.getNewEdgeId();
			map.setEdgeId(map.phi1(map.phi2(*nvol)), idedge, EDGE_ORBIT);
630 631 632
		}
	}

633 634 635 636 637 638 639 640
//	Dart ec = quadfaces.front();
//	unsigned int idedge = map.getNewEdgeId();
//	map.setEdgeId(map.phi1(map.phi1(ec)), idedge, EDGE_ORBIT);
//
//	ec = map.phi_1(map.phi2(map.phi3(map.phi2(map.phi1(ec)))));
//	idedge = map.getNewEdgeId();
//	map.setEdgeId(ec, idedge, EDGE_ORBIT);
 */
641