Addition of plugin DuplicateBoundary that allow one to modelize inter granular...

Addition of plugin DuplicateBoundary that allow one to modelize inter granular cracks with cohesive zones. 
Addition of the line command ComputeBestSeeds that allow one to optimize the placement of seeds in a microstructure to remove small edges from the Voronoï representation.
Small additions to cross3D and direction3D to compute the smoothness fo fields.

Common/CommandLine.cpp

@@ -6,6 +6,7 @@
 #include <string>
 #include <string.h>
 #include <stdlib.h>
+#include <iostream>
 #include <fstream>
 #include "GmshConfig.h"
 #include "GmshDefines.h"
@@ -527,6 +528,307 @@ void GetOptions(int argc, char *argv[])
           vm2.correspondance(0.00001,xMax,yMax,zMax);
         }
       }
+      else if(!strcmp(argv[i] + 1, "computeBestSeeds")) {
+              i++;
+              int j;
+              int radical;
+              double max;
+              double xMax;
+              double yMax;
+              double zMax;
+              std::vector<double> properties;
+              std::cout<<"entree dans computeBestSeeds"<<std::endl;
+              if(argv[i]){
+                std::ifstream file(argv[i++]);
+                file >> max;
+                file >> radical;
+                file >> xMax;
+                file >> yMax;
+                file >> zMax;
+                properties.clear();
+                properties.resize(4*max);
+                for(j=0;j<max;j++){
+                  file >> properties[4*j];
+                  file >> properties[4*j+1];
+                  file >> properties[4*j+2];
+                  file >> properties[4*j+3];
+                }
+                std::cout<<"Before count"<<std::endl;
+                std::vector<double> listDistances;
+                listDistances.clear();
+                int nbOfCount = 17;
+                listDistances.resize(nbOfCount);
+                for (unsigned int Count = 0; Count < nbOfCount; Count++){
+                	std::cout<<"Count"<<Count<<std::endl;
+                	double distMinGlobal = 0.0;
+                	int jMinGlobal = 0;
+                	int xORyORz = 0;
+                	int posORneg = 0;
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j] += 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 1;
+                			posORneg = 1;
+                		}
+                		delete m;
+                	}
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j + 1] += 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 2;
+                			posORneg = 1;
+                		}
+                		delete m;
+                	}
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j + 2] += 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 3;
+                			posORneg = 1;
+                		}
+                		delete m;
+                	}
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j] -= 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 1;
+                			posORneg = 2;
+                		}
+                		delete m;
+                	}
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j + 1] -= 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 2;
+                			posORneg = 2;
+                		}
+                		delete m;
+                	}
+                	for(j=0;j<max;j++){
+                		std::cout<<"j "<<j<<std::endl;
+                		std::vector<double> propertiesModified;
+                		propertiesModified.clear();
+                		propertiesModified.resize(4*max);
+                		std::cout<<"before assign propModif"<<std::endl;
+                		for(unsigned int k=0;k < properties.size();k++){
+                			propertiesModified[k] = properties[k];
+                		}
+                		std::cout<<"after assign propModif"<<std::endl;
+                		propertiesModified[4*j + 2] -= 0.01;
+                		voroMetal3D vm1;
+                		std::cout<<"before execute"<<std::endl;
+                		//std::remove("MicrostructurePolycrystal3D.geo");
+                		vm1.execute(propertiesModified,radical,0.1,xMax,yMax,zMax);
+                		//GModel::current()->destroy();
+                		GModel *m = new GModel();
+                		//GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                		m->load("MicrostructurePolycrystal3D.geo");
+                		double distMinTmp = 1000.0;
+                		//GModel *m = GModel::current();
+                		for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
+                			GEdge* eTmp = (*ite);
+                			GVertex* vTmp1 = eTmp->getBeginVertex();
+                			GVertex* vTmp2 = eTmp->getEndVertex();
+                			double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+                			if (distTmp < distMinTmp){
+                				distMinTmp = distTmp;
+                			}
+                		}
+                		if (distMinTmp > distMinGlobal){
+                			distMinGlobal = distMinTmp;
+                			jMinGlobal = j;
+                			xORyORz = 3;
+                			posORneg = 2;
+                		}
+                		delete m;
+                	}
+                	std::cout<<"distance minimale de "<<distMinGlobal<<std::endl;
+                	listDistances[Count] = distMinGlobal;
+                	if (xORyORz == 1){
+                		if (posORneg == 1){
+                			properties[4*jMinGlobal] += 0.01;
+                		} else if (posORneg == 2){
+                			properties[4*jMinGlobal] -= 0.01;
+                		}
+                	} else if (xORyORz == 2){
+                		if (posORneg == 1){
+                			properties[4*jMinGlobal + 1] += 0.01;
+                		} else if (posORneg == 2){
+                			properties[4*jMinGlobal + 1] -= 0.01;
+                		}
+                	} else if (xORyORz == 3){
+                		if (posORneg == 1){
+                			properties[4*jMinGlobal + 2] += 0.01;
+                		} else if (posORneg == 2){
+                			properties[4*jMinGlobal + 2] -= 0.01;
+                		}
+                	}
+                }
+                voroMetal3D vm1;
+                vm1.execute(properties,radical,0.1,xMax,yMax,zMax);
+                GModel::current()->load("MicrostructurePolycrystal3D.geo");
+                voroMetal3D vm2;
+                vm2.correspondance(0.00001,xMax,yMax,zMax);
+                for (unsigned int iTmp = 0; iTmp < listDistances.size();iTmp++){
+                	std::cout<<"distMinGlobal "<<iTmp<<" egale a "<<listDistances[iTmp]<<std::endl;
+                }
+                std::cout<<"liste des nouveaux seeds :"<<std::endl;
+                for(unsigned int iTmp = 0; iTmp < max;iTmp++){
+                	std::cout<<properties[4*iTmp]<<" "<<properties[4*iTmp + 1]<<" "<<properties[4*iTmp + 2]<<" "<<properties[4*iTmp + 3]<<std::endl;
+
+                }
+              }
+            }
 #endif
       else if(!strcmp(argv[i] + 1, "nopopup")) {
         CTX::instance()->noPopup = 1;

Geo/GFace.cpp

@@ -107,6 +107,22 @@ void GFace::delFreeEdge(GEdge *e)
   }
 }
 
+void GFace::replaceEdge(GEdge *e1, GEdge *e2){
+	std::list<GEdge*>::iterator ite = l_edges.begin();
+	std::list<GEdge*> newlist;
+	newlist.clear();
+	while(ite != l_edges.end()){
+		if(e1 == *ite){
+			newlist.push_back(e2);
+		}
+		else{
+			newlist.push_back((*ite));
+		}
+		ite++;
+	}
+	l_edges = newlist;
+}
+
 void GFace::deleteMesh()
 {
   for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i];

Geo/GFace.h

@@ -96,6 +96,10 @@ class GFace : public GEntity
   // edge in the face, not part of any edge loops--use with caution!)
   void delFreeEdge(GEdge *e);
 
+  //find the edge 1 from the list of edges and replace it by edge 2
+  //dont change the edge loops, and is susceptible to break some functionalities
+  void replaceEdge(GEdge *e1, GEdge *e2);
+
   // edge orientations
   virtual std::list<int> orientations() const { return l_dirs; }
 

Mesh/cross3D.h

@@ -39,6 +39,7 @@ public:
   double re() const{ return v[3]; }
   SVector3 im() const{ return SVector3(v[0], v[1], v[2]); }
   Qtn conj() { v[0] *= -1.; v[1] *= -1.; v[2] *= -1.; return *this; }
+  Qtn opp() { v[0] *= -1.; v[1] *= -1.; v[2] *= -1.; v[3] *= -1.; return *this; }
   double operator [](const unsigned int i) const { return v[i]; }
   void storeProduct(const Qtn &x, const Qtn &y);
   Qtn operator *(const Qtn &x) const;
@@ -137,6 +138,8 @@ public:
     scnd = im(R*scnd*conj(R));
     return *this;
   }
+  Qtn correspQuat();
+  Qtn rotationDirectTo(const cross3D &y) const;
   Qtn rotationTo(const cross3D &y) const;
   Qtn rotationToOnSurf(const cross3D &y) const;
 };
@@ -182,6 +185,71 @@ cross3D cross3D::get(int k) const{
   return cross3D(a,b);
 }
 
+Qtn cross3D::correspQuat(){
+	/*
+	   returns the quaternion representing the rotation from the canonical basis to this
+	 */
+	SVector3 frst = SVector3(1,0,0);
+	SVector3 scnd = SVector3(0,1,0);
+	cross3D xx = cross3D(frst,scnd);
+	cross3D yy = *this;
+	Qtn R = xx.rotationDirectTo(yy);
+	return R;
+
+}
+
+Qtn cross3D::rotationDirectTo(const cross3D &y) const{
+	  double d, dmin, jmin, kmin, th1, th2;
+	  SVector3 axis;
+	  Qtn Rxy1, Rxy2;
+
+	  cross3D xx = *this;
+	  cross3D yy = y;
+
+	  //ind1 = 0; jmin=0; dmin = angle(xx.get(kmin).frst, vect[jmin]);
+	  dmin = angle(xx.frst, yy.frst);
+
+	  th1 = dmin;
+
+	  if(th1 > 1e-8)
+	    axis = crossprod(xx.frst, yy.frst).unit();
+	  else {
+	    axis = SVector3(1,0,0);
+	    th1 = 0.;
+	  }
+
+	  Rxy1 = Qtn(axis, th1);
+	  xx.rotate(Rxy1);
+
+	  dmin = angle(xx.scnd, yy.scnd);
+	  // xx.scnd and yy.scnd now form the smallest angle among the 4^2 possible.
+
+	  th2 = dmin;
+
+	  if(th2 > 1e-8)
+	    axis = crossprod(xx.scnd, yy.scnd).unit();
+	  else {
+	    axis = SVector3(1,0,0);
+	    th2 = 0.;
+	  }
+
+	  Rxy2 = Qtn(axis, th2);
+	  Qtn R = Rxy2*Rxy1;
+
+	  // test
+	  double theta = eulerAngleFromQtn(R);
+	  if(theta > 1.07 /*0.988*/){ //
+	    std::cout << "Ouch! th1 = " << th1 << " th2 = " << th2 << std::endl;
+	    std::cout << "x = " << *this << std::endl;
+	    std::cout << "y = " << y << std::endl;
+	    std::cout << "R = " << R << std::endl;
+	    std::cout << "u = " << eulerAngleFromQtn(R) << std::endl;
+	    std::cout << "axis = " << eulerAxisFromQtn(R) << std::endl;
+	  }
+
+	  return R;
+}
+
 Qtn cross3D::rotationTo(const cross3D &y) const{
   /* x.rotationTo(y) returns a quaternion R representing the rotation 
      such that y = R x, x = conj(R) y
@@ -326,6 +394,93 @@ STensor3 convert(const cross3D x){
   return m;
 }
 
+double computeSetSmoothness(std::vector<cross3D > S){
+	  /* this->frst and y.frst are assumed to be the normal to the face
+	     R1 is the rotation such that R1(this->frst) = y.frst.
+	     R2 is then the rotation such that R2 o R1(this->scnd) = y.scnd
+	  */
+	double result = 1.0;
+	double qmean[4];
+	qmean[0] = 0.0;
+	qmean[1] = 0.0;
+	qmean[2] = 0.0;
+	qmean[3] = 0.0;
+	std::vector<cross3D>::iterator it1 = S.begin();
+	cross3D cInit = (*it1);
+	Qtn qInit = cInit.correspQuat();
+	for (it1 = S.begin();it1 != S.end();it1++){
+		//pour chaque element du set
+		cross3D cTmp = (*it1);
+		Qtn qTmp = cTmp.correspQuat();
+		double prodVecMin = 0.0;
+		for (unsigned int i = 1;i < 24;i++){
+			//on trouve la cross appropriee
+			Qtn qTmpi = cTmp.get(i).correspQuat();
+			double prodVeci = 0.0;
+			for (unsigned int j = 0;j < 4;j++){
+				prodVeci += qInit.v[j] * qTmpi.v[j];
+			}
+			if (prodVeci >= 0.0){
+				if (prodVeci > prodVecMin){
+					prodVecMin = prodVeci;
+					qTmp = qTmpi;
+				}
+			}
+			else{
+				prodVeci = - prodVeci;
+				qTmpi = qTmpi.opp();
+				if (prodVeci > prodVecMin){
+					prodVecMin = prodVeci;
+					qTmp = qTmpi;
+				}
+			}
+		}
+		//on a trouve le quat approprie
+		for (unsigned int j = 0;j < 4;j++){
+			qmean[j] += qTmp[j];
+		}
+	}
+	double normQt = sqrt(qmean[0] * qmean[0] + qmean[1] * qmean[1] + qmean[2] * qmean[2] + qmean[3] * qmean[3]);
+	if (normQt != 0.0){
+		for (unsigned int j = 0;j < 4;j++){
+			qmean[j] = qmean[j] / normQt;
+		}
+	}
+	for (it1 = S.begin();it1 != S.end();it1++){
+		//pour chaque element du set
+		cross3D cTmp = (*it1);
+		Qtn qTmp = cTmp.correspQuat();
+		double prodVecMin = 0.0;
+		for (unsigned int i = 0;i < 24;i++){
+			//on trouve la cross appropriee
+			Qtn qTmpi = cTmp.get(i).correspQuat();
+			double prodVeci = 0.0;
+			for (unsigned int j = 0;j < 4;j++){
+				prodVeci += qmean[j] * qTmpi.v[j];
+			}
+			if (prodVeci >= 0.0){
+				if (prodVeci > prodVecMin){
+					prodVecMin = prodVeci;
+					qTmp = qTmpi;
+				}
+			}
+			else{
+				prodVeci = - prodVeci;
+				qTmpi = qTmpi.opp();
+				if (prodVeci > prodVecMin){
+					prodVecMin = prodVeci;
+					qTmp = qTmpi;
+				}
+			}
+		}
+		//on a trouve le quat approprie
+		if (prodVecMin < result){
+			result = prodVecMin;
+		}
+	}
+	return result;
+}
+
 
 
 #endif

Mesh/directions3D.cpp

@@ -674,6 +674,89 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
   return energy;
 }
 
+void Frame_field::buildSmoothness(){
+	GModel *m = GModel::current();
+	std::vector<GEntity*> entities;
+	m->getEntities(entities);
+	//pour commencer on va creer une map de connectique des Mvertex
+	std::map<MVertex*, std::vector<MVertex* > > Neighbours;
+	for(unsigned int i = 0; i < entities.size(); i++){
+		GEntity* eTmp = entities[i];
+		for (int j = 0; j < eTmp->getNumMeshElements();j++){
+			MElement* elem = eTmp->getMeshElement(j);
+			for (unsigned int k = 0;k < elem->getNumVertices();k++){
+				for (unsigned int l = k;l < elem->getNumVertices();l++){
+					if (k != l){
+						MVertex* v1 = elem->getVertex(k);
+						MVertex* v2 = elem->getVertex(l);
+						Neighbours[v1].push_back(v2);
+						Neighbours[v2].push_back(v1);
+					}
+				}
+			}
+		}
+	}
+	for(unsigned int i = 0; i < entities.size(); i++){
+		for (unsigned int j = 0;j < entities[i]->mesh_vertices.size();j++){
+			//on va traiter chaque noeud
+			std::set<MVertex*> V1;
+			std::set<MVertex*> V2;
+			std::set<MVertex*> V3;
+			MVertex* v0 = entities[i]->mesh_vertices[j];
+			V1.insert(v0);
+			std::vector<MVertex* > v0vec = Neighbours[v0];
+			for (unsigned int k = 0;k < v0vec.size();k++){
+				V1.insert(v0vec[k]);
+			}
+			for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
+				MVertex* vTmp = (*itSet);
+				V2.insert(vTmp);
+				v0vec = Neighbours[vTmp];
+				for (unsigned int k = 0;k < v0vec.size();k++){
+					V2.insert(v0vec[k]);
+				}
+			}
+			for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
+				MVertex* vTmp = (*itSet);
+				V3.insert(vTmp);
+				v0vec = Neighbours[vTmp];
+				for (unsigned int k = 0;k < v0vec.size();k++){
+					V3.insert(v0vec[k]);
+				}
+			}
+			//we have all three set here, time to compute the smoothnesses for each one
+			std::vector<cross3D> C1;
+			std::vector<cross3D> C2;
+			std::vector<cross3D> C3;
+			double S1 = 0.0;
+			double S2 = 0.0;
+			double S3 = 0.0;
+			for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
+				MVertex* vTmp = (*itSet);
+				STensor3 tTmp = crossField[vTmp];
+				cross3D cTmp = cross3D(tTmp);
+				C1.push_back(cTmp);
+			}
+			for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
+				MVertex* vTmp = (*itSet);
+				STensor3 tTmp = crossField[vTmp];
+				cross3D cTmp = cross3D(tTmp);
+				C2.push_back(cTmp);
+			}
+			for (std::set<MVertex*>::iterator itSet = V3.begin();itSet != V3.end();itSet++){
+				MVertex* vTmp = (*itSet);
+				STensor3 tTmp = crossField[vTmp];
+				cross3D cTmp = cross3D(tTmp);
+				C3.push_back(cTmp);
+			}
+			S1 = computeSetSmoothness(C1);
+			S2 = computeSetSmoothness(C2);
+			S3 = computeSetSmoothness(C3);
+			double finalSmoothness = (9.0 * S1 + 3.0 * S2 + S3) / 13.0 ;
+			crossFieldSmoothness[v0] = finalSmoothness;
+		}
+	}
+}
 
 double Frame_field::smoothFace(GFace *gf, int n){
   double energy=0;
@@ -1642,6 +1725,7 @@ void Nearest_point::clear(){
 std::vector<std::pair<SPoint3,STensor3> > Frame_field::field;
 std::vector<int> Frame_field::labels;
 std::map<MVertex*, STensor3> Frame_field::crossField;
+std::map<MVertex*, double> Frame_field::crossFieldSmoothness;
 std::map<MEdge, double, Less_Edge> Frame_field::crossDist;
 std::map<MVertex*,std::set<MVertex*> > Frame_field::vertex_to_vertices;
 std::map<MVertex*,std::set<MElement*> > Frame_field::vertex_to_elements;

Mesh/directions3D.h

@@ -28,6 +28,7 @@ class Frame_field{
   static std::vector<std::pair<SPoint3,STensor3> > field;
   static std::vector<int> labels;
   static std::map<MVertex*, STensor3> crossField;
+  static std::map<MVertex*, double> crossFieldSmoothness;
   static std::map<MEdge, double, Less_Edge> crossDist;
   static std::vector<MVertex*> listVertices;
 #if defined(HAVE_ANN)
@@ -54,6 +55,7 @@ class Frame_field{
   static STensor3 findCross(double x, double y, double z);
   static void initFace(GFace* gf);
   static void initRegion(GRegion* gr, int n);
+  static void buildSmoothness();
   static double smoothFace(GFace *gf, int n);
   static double smoothRegion(GRegion *gr, int n);
   static double smooth();

Plugin/CMakeLists.txt

@@ -32,7 +32,7 @@ set(SRC
   Scal2Tens.cpp Scal2Vec.cpp
   CutMesh.cpp
   NewView.cpp
-  SimplePartition.cpp Crack.cpp
+  SimplePartition.cpp Crack.cpp DuplicateBoundaries.cpp
   FaultZone.cpp
 )
 

Plugin/DuplicateBoundaries.h

+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _DUPLICATEBOUNDARIES_H_
+#define _DUPLICATEBOUNDARIES_H_
+
+#include "Plugin.h"
+
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterDuplicateBoundariesPlugin();
+}
+
+class GMSH_DuplicateBoundariesPlugin : public GMSH_PostPlugin
+{
+ public:
+  GMSH_DuplicateBoundariesPlugin(){}
+  std::string getName() const { return "DuplicateBoundaries"; }
+  std::string getShortHelp() const
+  {
+    return "Duplicate Boundaries";
+  }
+  std::string getHelp() const;
+  int getNbOptions() const;
+  StringXNumber* getOption(int iopt);
+  double abs(double x);
+  PView *execute(PView *);
+  PView *execute2D(PView *);
+  PView *execute2DWithBound(PView *);
+  PView *executeDuplicate(PView *);
+  PView *executeBis(PView *);
+  PView *executeTer(PView *);
+  PView *executeFourth(PView *);
+  PView *ComputeBestSeeds(PView *);
+};
+
+#endif

Plugin/PluginManager.cpp

@@ -25,6 +25,7 @@
 #include "ExtractElements.h"
 #include "SimplePartition.h"
 #include "Crack.h"
+#include "DuplicateBoundaries.h"
 #include "HarmonicToTime.h"
 #include "ModulusPhase.h"
 #include "Integrate.h"
@@ -257,6 +258,8 @@ void PluginManager::registerDefaultPlugins()
                       ("Crack", GMSH_RegisterCrackPlugin()));
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("FaultZone", GMSH_RegisterFaultZonePlugin()));
+    allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+                      ("DuplicateBoundaries", GMSH_RegisterDuplicateBoundariesPlugin()));
 #if defined(HAVE_TETGEN)
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("Tetrahedralize", GMSH_RegisterTetrahedralizePlugin()));