Commit 5b0e1e17 authored by untereiner's avatar untereiner

Merge cgogn:~cgogn/CGoGN into HEAD

Conflicts:
	include/Algo/Modelisation/subdivision.h
	include/Algo/Modelisation/subdivision.hpp
parents 32045a8a 74ad6e7b
......@@ -200,7 +200,7 @@ void MyQT::cb_initGL()
// create the renders
m_topo_render = new Algo::Render::GL2::Topo3Render();
m_explode_render = new Algo::Render::GL2::ExplodeVolumeRender(true,true);
m_explode_render = new Algo::Render::GL2::ExplodeVolumeRender(true,true,true);
SelectorDartNoBoundary<PFP::MAP> nb(myMap);
m_topo_render->updateData<PFP>(myMap, position, 0.8f, 0.8f, 0.8f, nb);
......
......@@ -49,7 +49,7 @@ void filterAverageAttribute_OneRing(
const FunctorSelect& select = allDarts)
{
FunctorAverage<T, VERTEX> fa(attIn) ;
Algo::Selection::Collector_OneRing<PFP> col(map) ;
Algo::Surface::Selection::Collector_OneRing<PFP> col(map) ;
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
......@@ -98,7 +98,7 @@ void filterAverageVertexAttribute_WithinSphere(
{
FunctorAverage<T, VERTEX> faInside(attIn) ;
FunctorAverageOnSphereBorder<PFP, T> faBorder(map, attIn, position) ;
Algo::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
Algo::Surface::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
......@@ -139,7 +139,7 @@ void filterAverageEdgeAttribute_WithinSphere(
const FunctorSelect& select = allDarts)
{
FunctorAverage<T, EDGE> fa(attIn) ;
Algo::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
Algo::Surface::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
TraversorE<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
......@@ -167,7 +167,7 @@ void filterAverageFaceAttribute_WithinSphere(
const FunctorSelect& select = allDarts)
{
FunctorAverage<T, FACE> fa(attIn) ;
Algo::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
Algo::Surface::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
TraversorF<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
......
......@@ -87,9 +87,9 @@ void filterAverageNormals(typename PFP::MAP& map, const VertexAttribute<typename
FaceAutoAttribute<VEC3> faceNormal(map, "faceNormal") ;
FaceAutoAttribute<VEC3> faceCentroid(map, "faceCentroid") ;
Algo::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
FaceAutoAttribute<VEC3> faceNewNormal(map, "faceNewNormal") ;
......@@ -131,9 +131,9 @@ void filterMMSE(typename PFP::MAP& map, float sigmaN2, const VertexAttribute<typ
FaceAutoAttribute<VEC3> faceNormal(map, "faceNormal") ;
FaceAutoAttribute<VEC3> faceCentroid(map, "faceCentroid") ;
Algo::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
FaceAutoAttribute<VEC3> faceNewNormal(map, "faceNewNormal") ;
......@@ -216,9 +216,9 @@ void filterTNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
FaceAutoAttribute<VEC3> faceNormal(map, "faceNormal") ;
FaceAutoAttribute<VEC3> faceCentroid(map, "faceCentroid") ;
Algo::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
FaceAutoAttribute<VEC3> faceNewNormal(map, "faceNewNormal") ;
......@@ -335,9 +335,9 @@ void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
FaceAutoAttribute<VEC3> faceNormal(map, "faceNormal") ;
FaceAutoAttribute<VEC3> faceCentroid(map, "faceCentroid") ;
Algo::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea, select) ;
Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal, select) ;
Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid, select) ;
VertexAutoAttribute<REAL> vertexArea(map, "vertexArea") ;
FaceAutoAttribute<VEC3> faceNewNormal(map, "faceNewNormal") ;
......@@ -369,7 +369,7 @@ void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
float angle = Geom::angle(normV, neighborNormal) ;
if( angle <= SUSANthreshold )
{
REAL umbArea = Algo::Geometry::vertexOneRingArea<PFP>(map, it, position) ;
REAL umbArea = Algo::Surface::Geometry::vertexOneRingArea<PFP>(map, it, position) ;
vertexArea[it] = umbArea ;
sumArea += umbArea ;
......
......@@ -49,7 +49,7 @@ void sigmaBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP::
TraversorE<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
sumLengths += Algo::Geometry::edgeLength<PFP>(map, d, position) ;
sumLengths += Algo::Surface::Geometry::edgeLength<PFP>(map, d, position) ;
sumAngles += Geom::angle(normal[d], normal[map.phi1(d)]) ;
++nbEdges ;
}
......@@ -80,7 +80,7 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP:
Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next())
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) ) ;
......@@ -125,7 +125,7 @@ void filterSUSAN(typename PFP::MAP& map, float SUSANthreshold, const VertexAttri
float angle = Geom::angle(normal_d, neighborNormal) ;
if( angle <= SUSANthreshold )
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) );
......
......@@ -42,7 +42,7 @@ void filterTaubin(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& p
{
typedef typename PFP::VEC3 VEC3 ;
Algo::Selection::Collector_OneRing<PFP> c(map) ;
Algo::Surface::Selection::Collector_OneRing<PFP> c(map) ;
const float lambda = 0.6307 ;
const float mu = -0.6732 ;
......@@ -109,7 +109,7 @@ void filterTaubin_modified(typename PFP::MAP& map, VertexAttribute<typename PFP:
CellMarkerNoUnmark<VERTEX> mv(map) ;
FunctorAverageOnSphereBorder<PFP, VEC3> fa1(map, position, position) ;
Algo::Selection::Collector_WithinSphere<PFP> c1(map, position, radius) ;
Algo::Surface::Selection::Collector_WithinSphere<PFP> c1(map, position, radius) ;
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !mv.isMarked(d))
......@@ -133,7 +133,7 @@ void filterTaubin_modified(typename PFP::MAP& map, VertexAttribute<typename PFP:
// unshrinking step
FunctorAverageOnSphereBorder<PFP, VEC3> fa2(map, position2, position2) ;
Algo::Selection::Collector_WithinSphere<PFP> c2(map, position2, radius) ;
Algo::Surface::Selection::Collector_WithinSphere<PFP> c2(map, position2, radius) ;
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && mv.isMarked(d))
......
......@@ -27,6 +27,15 @@
#include "Topology/generic/traversorCell.h"
#include "Topology/generic/traversor2.h"
extern "C"
{
#include "C_BLAS_LAPACK/INCLUDE/f2c.h"
#include "C_BLAS_LAPACK/INCLUDE/clapack.h"
}
#undef max
#undef min
namespace CGoGN
{
......
......@@ -14,6 +14,10 @@ namespace CGoGN
namespace Algo
{
namespace Surface
{
namespace Geometry
{
......@@ -111,6 +115,7 @@ protected :
}// end namespace Geometry
}// end namespace Surface
}// end namespace Algo
}// end namespace CGoGN
......
......@@ -4,6 +4,9 @@ namespace CGoGN
namespace Algo
{
namespace Surface
{
namespace Geometry
{
......@@ -543,5 +546,6 @@ unsigned int CentroidalVoronoiDiagram<PFP>::moveSeed(unsigned int numSeed){
}
*/
}// end namespace Geometry
} // Surface
}// end namespace Algo
}// end namespace CGoGN
......@@ -116,7 +116,8 @@ template <typename PFP>
void DooSabin(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position);
///**
// * Sqrt(3) subdivision scheme
// * Reverse the orientation of the map
// * NOW IN THE MAP
// */
//template <typename PFP>
//void Sqrt3Subdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected = allDarts) ;
......@@ -124,14 +125,10 @@ void DooSabin(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
///**
// * Dual mesh computation
// * Sqrt(3) subdivision scheme
// */
//template <typename PFP>
//void computeDual(typename PFP::MAP& map, const FunctorSelect& selected = allDarts) ;
//
//template <typename PFP>
//void computeDualV2(typename PFP::MAP& map, const FunctorSelect& selected = allDarts) ;
//
//void Sqrt3Subdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected = allDarts) ;
......
......@@ -698,108 +698,6 @@ inline double sqrt3_K(unsigned int n)
//}
//template <typename PFP>
//void computeDual(typename PFP::MAP& map, const FunctorSelect& selected)
//{
// DartAttribute<Dart> phi1 = map.template getAttribute<Dart, DART>("phi1") ;
// DartAttribute<Dart> phi_1 = map.template getAttribute<Dart, DART>("phi_1") ;
// DartAttribute<Dart> new_phi1 = map.template addAttribute<Dart, DART>("new_phi1") ;
// DartAttribute<Dart> new_phi_1 = map.template addAttribute<Dart, DART>("new_phi_1") ;
//
// for(Dart d = map.begin(); d != map.end(); map.next(d))
// {
// Dart dd = map.phi1(map.phi2(d));
//
// new_phi1[d] = dd ;
// new_phi_1[dd] = d ;
// }
//
// map.template swapAttributes<Dart>(phi1, new_phi1) ;
// map.template swapAttributes<Dart>(phi_1, new_phi_1) ;
//
// map.removeAttribute(new_phi1) ;
// map.removeAttribute(new_phi_1) ;
//
// map.swapEmbeddingContainers(VERTEX, FACE) ;
//
// reverseOrientation<PFP>(map) ;
//
//// //boundary management
//// for(Dart d = map.begin(); d != map.end(); map.next(d))
//// {
//// if(map.isBoundaryMarked(d))
//// {
//// map.template boundaryMarkOrbit<FACE>(map.deleteVertex(map.phi2(d))); //map.deleteCycle();
//// }
//// }
//}
//
//template <typename PFP>
//void computeDualV2(typename PFP::MAP& map, const FunctorSelect& selected)
//{
// DartAttribute<Dart> phi1 = map.template getAttribute<Dart, DART>("phi1") ;
// DartAttribute<Dart> phi_1 = map.template getAttribute<Dart, DART>("phi_1") ;
// DartAttribute<Dart> new_phi1 = map.template addAttribute<Dart, DART>("new_phi1") ;
// DartAttribute<Dart> new_phi_1 = map.template addAttribute<Dart, DART>("new_phi_1") ;
//
// for(Dart d = map.begin(); d != map.end(); map.next(d))
// {
// Dart dd = map.phi1(map.phi2(d));
//
// new_phi1[d] = dd ;
// new_phi_1[dd] = d ;
// }
//
// map.template swapAttributes<Dart>(phi1, new_phi1) ;
// map.template swapAttributes<Dart>(phi_1, new_phi_1) ;
//
// map.removeAttribute(new_phi1) ;
// map.removeAttribute(new_phi_1) ;
//
// //boundary management
// for(Dart d = map.begin(); d != map.end(); map.next(d))
// {
// if(map.isBoundaryMarked(d))
// {
// Dart d1 = map.phi1(d);
// Dart dd = map.phi_1(map.phi2(d));
//
// //marquer le brin d1 et le bin dd
//
// phi1[dd] = d1 ;
// phi_1[d1] = dd ;
//
// phi1[d] = map.phi2(d);
// phi_1[map.phi2(d)] = d;
// }
// }
//
// map.swapEmbeddingContainers(VERTEX, FACE) ;
//
// reverseOrientation<PFP>(map) ;
//
// for(Dart d = map.begin(); d != map.end(); map.next(d))
// {
// if(map.isBoundaryMarked(d))
// {
// map.deleteCycle(map.phi2(d));
// }
// }
//
// //transformer le marker des brins d1 et dd en boundaryMarker
//}
} // namespace Modelisation
}
......
......@@ -93,7 +93,8 @@ public:
/**
* Constructor
* @param withColorPerFace affect a color per face
* @param withExplodeFace shrinj each face
* @param withExplodeFace shrink each face
* @param withSmoothFaces use a smooth gouraud interpolation between triangles of a faces
*/
ExplodeVolumeRender(bool withColorPerFace = false, bool withExplodeFace = false, bool withSmoothFaces = false) ;
......
......@@ -25,6 +25,7 @@
#ifndef __INCLUSION__
#define __INCLUSION__
#include "Geometry/basic.h"
#include "Geometry/orientation.h"
......
......@@ -439,8 +439,6 @@ public:
* @param fonct functor obj ref
*/
bool foreach_dart_of_edge1(Dart d, FunctorType& fonct, unsigned int thread = 0);
//@}
/*! @name Close map after import or creation
......@@ -469,6 +467,17 @@ public:
*/
unsigned int closeMap();
//@}
/*! @name Compute dual
* These functions compute the dual mesh
*************************************************************************/
//@{
//! Dual mesh computation
/*!
*/
void computeDual();
//@}
};
} // namespace CGoGN
......
......@@ -378,6 +378,19 @@ public:
* @return the number of closed holes
*/
unsigned int closeMap();
//@}
/*! @name Compute dual
* These functions compute the dual mesh
*************************************************************************/
//@{
//! Dual mesh computation
/*!
*/
void computeDual();
//@}
};
} // namespace CGoGN
......
......@@ -432,7 +432,6 @@ public:
* These functions must be used with care, generally only by import algorithms
*************************************************************************/
//@{
/**
* create a face of map1 marked as boundary
......
......@@ -27,6 +27,7 @@
#include <assert.h>
#include <vector>
#include <list>
namespace CGoGN {
......
......@@ -947,4 +947,20 @@ unsigned int GMap2::closeMap()
return nb ;
}
/*! @name Compute dual
* These functions compute the dual mesh
*************************************************************************/
void GMap2::computeDual()
{
// DartAttribute<Dart> old_beta0 = this->getAttribute<Dart, DART>("beta0");
// DartAttribute<Dart> old_beta2 = this->getAttribute<Dart, DART>("beta2") ;
//
// swapAttributes<Dart>(old_beta0, old_beta2) ;
//
// swapEmbeddingContainers(VERTEX, FACE) ;
//
// //boundary management ?
}
} // namespace CGoGN
......@@ -996,4 +996,24 @@ unsigned int GMap3::closeMap()
return nb ;
}
/*! @name Compute dual
* These functions compute the dual mesh
*************************************************************************/
void GMap3::computeDual()
{
// DartAttribute<Dart> old_beta0 = getAttribute<Dart, DART>("beta0");
// DartAttribute<Dart> old_beta1 = getAttribute<Dart, DART>("beta1");
// DartAttribute<Dart> old_beta2 = getAttribute<Dart, DART>("beta2");
// DartAttribute<Dart> old_beta3 = getAttribute<Dart, DART>("beta3") ;
//
// swapAttributes<Dart>(old_beta0, old_beta3) ;
// swapAttributes<Dart>(old_beta1, old_beta2) ;
//
// swapEmbeddingContainers(VERTEX, FACE) ;
//
// //boundary management ?
}
} // namespace CGoGN
......@@ -278,8 +278,8 @@ Dart Map3::splitVertex(std::vector<Dart>& vd)
Dart Map3::deleteVertex(Dart d)
{
if(isBoundaryVertex(d))
return NIL ;
//if(isBoundaryVertex(d))
// return NIL ;
// Save the darts around the vertex
// (one dart per face should be enough)
......@@ -301,27 +301,55 @@ Dart Map3::deleteVertex(Dart d)
}
}
std::cout << "nb faces " << fstore.size() << std::endl;
Dart res = NIL ;
for(std::vector<Dart>::iterator it = fstore.begin() ; it != fstore.end() ; ++it)
{
Dart fit = *it ;
Dart end = phi_1(fit) ;
fit = phi1(fit) ;
while(fit != end)
if(fit == end)
{
Dart d2 = phi2(fit) ;
Dart d3 = phi3(fit) ;
Dart d32 = phi2(d3) ;
std::cout << " mmmmmmmmmmmmmmmmmmmmmerrrrrrrrrrrrrrrrrde !!!!!!!!!!!! " << std::endl;
if(res == NIL)
res = d2 ;
// Dart d2 = phi2(fit) ;
// Dart d23 = phi3(d2) ;
// Dart d3 = phi3(fit) ;
// Dart d32 = phi2(d3) ;
//
// //phi3unsew()
// phi3sew(d3,23);
//
// fit = phi_1(fit);
//
// d2 = phi2(fit) ;
// d23 = phi3(d2) ;
// d3 = phi3(fit) ;
// d32 = phi2(d3) ;
// phi3sew(d3,23);
phi2unsew(d2) ;
phi2unsew(d32) ;
phi2sew(d2, d32) ;
phi2sew(fit, d3) ;
// Map2::deleteCC(fit);
}
else
{
while(fit != end)
{
Dart d2 = phi2(fit) ;
Dart d3 = phi3(fit) ;
Dart d32 = phi2(d3) ;
fit = phi1(fit) ;
if(res == NIL)
res = d2 ;
phi2unsew(d2) ;
phi2unsew(d32) ;
phi2sew(d2, d32) ;
phi2sew(fit, d3) ;
fit = phi1(fit) ;
}
}
}
......@@ -1270,60 +1298,86 @@ void Map3::computeDual()
unsigned int count = 0;
std::vector<Dart> vbound;
for(Dart d = begin(); d != end(); next(d))
{
if(isBoundaryMarked3(d))
boundaryUnmark<3>(d);
if(d == 1569)
{
std::cout << "d " << std::endl;
Traversor3WE<Map3> t(*this,d);
for(Dart dit = t.begin() ; dit != t.end() ; dit = t.next())
{
Dart temp = dit;
do
{
if(isBoundaryMarked3(d))
std::cout << "d boundary " << std::endl;
temp = alpha2(temp);
}while(temp != dit);
}
if(isBoundaryMarked3(d))
std::cout << "d boundary " << std::endl;
if(isBoundaryMarked3(phi1(d)))
std::cout << "phi1(d) boundary " << std::endl;
if(isBoundaryMarked3(phi_1(d)))
std::cout << "phi_1(d) boundary " << std::endl;
if(isBoundaryMarked3(phi2(d)))
std::cout << "phi2(d) boundary " << std::endl;
if(isBoundaryMarked3(phi3(d)))
std::cout << "phi3(d) boundary " << std::endl;
if(isBoundaryMarked3(phi2(phi3(d))))
std::cout << "phi2(phi3(d)) boundary " << std::endl;
if(isBoundaryMarked3(d) && !isBoundaryMarked3(phi3(d)))
{
vbound.push_back(d);
}
}
if(isBoundaryMarked3(phi3(phi2(d))))
std::cout << "phi3(phi2(d)) boundary " << std::endl;
std::cout << "vbound size = " << vbound.size() << std::endl;
if(isBoundaryMarked3(phi1(phi3(d))))
std::cout << "phi1(phi3(d)) boundary " << std::endl;
for(std::vector<Dart>::iterator it = vbound.begin() ; it != vbound.end() ; ++it)
{
Dart d = *it;
//Dart d3 = phi3(d);
phi3unsew(d);
//phi3unsew(d3);
}
if(isBoundaryMarked3(phi3(phi1(d))))
std::cout << "phi3(phi1(d)) boundary " << std::endl;
}
if(isBoundaryMarked3(d))
{
//std::cout << "nb faces : " << closeMap() << std::endl;
// if(d == 14208)
// {
// std::cout << "yeahhhhhhhh" << std::endl;
// std::cout << "isBoundaryMarked ? " << isBoundaryMarked3(phi3(phi2(14208))) << std::endl;
//
// }
//
// //boundaryUnmark<3>(d);
//
// }
// if(d == 1569)
// {
// std::cout << "d " << std::endl;
//
// Traversor3WE<Map3> t(*this,d);
// for(Dart dit = t.begin() ; dit != t.end() ; dit = t.next())
// {
// Dart temp = dit;
// do
// {
// if(isBoundaryMarked3(d))
// std::cout << "d boundary " << std::endl;
//
// temp = alpha2(temp);
// }while(temp != dit);
// }
//
// if(isBoundaryMarked3(d))
// std::cout << "d boundary " << std::endl;
//
// if(isBoundaryMarked3(phi1(d)))
// std::cout << "phi1(d) boundary " << std::endl;
//
// if(isBoundaryMarked3(phi_1(d)))
// std::cout << "phi_1(d) boundary " << std::endl;
//
// if(isBoundaryMarked3(phi2(d)))
// std::cout << "phi2(d) boundary " << std::endl;
//
// if(isBoundaryMarked3(phi3(d)))
// std::cout << "phi3(d) boundary " << std::endl;
//
// if(isBoundaryMarked3(phi2(phi3(d))))
// std::cout << "phi2(phi3(d)) boundary " << std::endl;
//