diff --git a/Geo/MElement.h b/Geo/MElement.h
index 90de79149e84abe11f99ff9851e6d7f25381bc30..45f3ffdca3eda10045fe12e73e6bfab4573b0423 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -159,6 +159,7 @@ class MElement
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
virtual double distoShapeMeasure(){ return 1.0; }
+ virtual double angleShapeMeasure() { return 1.0; }
// compute the barycenter
virtual SPoint3 barycenter();
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 45b1b501c7e1fbea9de5a532a033a3612c6b1e04..04be90f28b0a152f7bdf67e49e63c8d1dfc4d7d4 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -32,6 +32,15 @@ double MTriangle::distoShapeMeasure()
#endif
}
+double MTriangle::angleShapeMeasure()
+{
+#if defined(HAVE_MESH)
+ return qmTriangleAngles(this);
+#else
+ return 0.;
+#endif
+}
+
double MTriangle::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 67172ff3581b9e3ca7b213f6d151579eb3b0554f..b8ec8fef1d87ff88312036c4e80698f8762da5ad 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -52,6 +52,7 @@ class MTriangle : public MElement {
virtual int getDim(){ return 2; }
virtual double gammaShapeMeasure();
virtual double distoShapeMeasure();
+ virtual double angleShapeMeasure();
virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index e6e5168e5325ad08b6436ab5f67ef4bccf5e996e..c593b323401a2ddd2919b4f67df43a06d6290c68 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -28,6 +28,7 @@
#include "dofManager.h"
#include "elasticityTerm.h"
#include "linearSystemCSR.h"
+#include "linearSystemFull.h"
#define SQU(a) ((a)*(a))
@@ -43,7 +44,8 @@ extern double angle3Points(MVertex *p1, MVertex *p2, MVertex *p3);
static double shapeMeasure(MElement *e)
{
- const double d1 = e->distoShapeMeasure();
+ //const double d1 = e->distoShapeMeasure();
+ const double d1 = e->angleShapeMeasure();
//const double d2 = e->gammaShapeMeasure();
return d1;
}
@@ -365,7 +367,7 @@ void highOrderSmoother::optimize(GFace * gf,
// }
// then try to swap for better configurations
- smooth(gf, true);
+ //smooth(gf, true);
//int nbSwap =
swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
// smooth(gf,true);
@@ -711,6 +713,128 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
delete lsys;
}
+void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
+ std::vector<MElement*>& old_elems,
+ GFace* gf) {
+
+ printf("Smoothing a cavity...");
+
+ linearSystemFull<double> *lsys = new linearSystemFull<double>;
+
+ dofManager<double> myAssembler(lsys);
+ elasticityTerm El(0,1.0,0.333, getTag());
+
+ std::vector<MElement*> layer,v;
+
+ v.insert(v.begin(), gf->triangles.begin(),gf->triangles.end());
+ v.insert(v.begin(), gf->quadrangles.begin(),gf->quadrangles.end());
+
+ addOneLayer(v,cavity,layer);
+
+ for (unsigned int i = 0; i < layer.size(); i++){
+ bool skip = false;
+ for (std::vector<MElement*>::iterator itold = old_elems.begin();
+ itold != old_elems.end();
+ itold++)
+ if (*itold == layer[i]) {
+ skip = true;
+ break;
+ }
+ if (skip) continue;
+ for (int j = 0; j < layer[i]->getNumVertices(); j++){
+ MVertex *vert = layer[i]->getVertex(j);
+ //printf("i = %d, j = %d \n",i,j);
+ myAssembler.fixVertex(vert, 0, getTag(), 0);
+ myAssembler.fixVertex(vert, 1, getTag(), 0);
+ }
+ }
+
+ std::map<MVertex*,SVector3>::iterator its;
+ std::set<MVertex*> verticesToMove;
+
+ for (unsigned int i = 0; i < cavity.size(); i++){
+ for (int j = 0; j < cavity[i]->getNumVertices(); j++){
+ MVertex *vert = cavity[i]->getVertex(j);
+ //printf("Vertex on dim %d with tag %d\n",vert->onWhat()->dim(),vert->onWhat()->tag());
+
+ // printf("%d %d %d v\n", i, j, v[i]->getNumVertices());
+ its = _straightSidedLocation.find(vert);
+ if (its == _straightSidedLocation.end()){
+ //printf("Ping\n");
+ _straightSidedLocation[vert] =
+ SVector3(vert->x(), vert->y(), vert->z());
+ _targetLocation[vert] =
+ SVector3(vert->x(), vert->y(), vert->z());
+ if (vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex(vert, 0, getTag(), 0);
+ myAssembler.fixVertex(vert, 1, getTag(), 0);
+ }
+ }
+ else{
+ //printf("Pong\n");
+ vert->x() = its->second.x();
+ vert->y() = its->second.y();
+ vert->z() = its->second.z();
+ if (vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex(vert, 0, getTag(), 0);
+ myAssembler.fixVertex(vert, 1, getTag(), 0);
+ }
+ }
+ }
+ }
+
+ // number the other DOFs
+ for (unsigned int i = 0; i < cavity.size(); i++){
+ for (int j = 0; j < cavity[i]->getNumVertices(); j++){
+ //printf("Numbering i = %d, j = %d \n",i,j);
+ //printf("%d vertices FIXED %d NUMBERED\n", myAssembler.sizeOfF()
+ // , myAssembler.sizeOfR());
+ MVertex *vert = cavity[i]->getVertex(j);
+ myAssembler.numberVertex(vert, 0, getTag());
+ myAssembler.numberVertex(vert, 1, getTag());
+ verticesToMove.insert(vert);
+ }
+ }
+
+ //Msg::Info("%d vertices FIXED %d NUMBERED\n", myAssembler.sizeOfF()
+ // , myAssembler.sizeOfR());
+
+ double dx0 = smooth_metric_(cavity, gf, myAssembler, verticesToMove, El);
+ double dx = dx0;
+ printf(" dx0 = %12.5E\n", dx0);
+ int iter = 0;
+ while(1){
+ double dx2 = smooth_metric_(cavity, gf, myAssembler, verticesToMove, El);
+ printf(" dx2 = %12.5E\n", dx2);
+ if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
+ if (iter++ > 2)break;
+ dx = dx2;
+ }
+
+ std::set<MVertex*>::iterator it;
+
+ for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
+ SPoint2 param;
+ reparamMeshVertexOnFace(*it, gf, param);
+ GPoint gp = gf->point(param);
+ if ((*it)->onWhat()->dim() == 2){
+ (*it)->x() = gp.x();
+ (*it)->y() = gp.y();
+ (*it)->z() = gp.z();
+ _targetLocation[*it] = SVector3(gp.x(), gp.y(), gp.z());
+ }
+ else{
+ SVector3 p = _targetLocation[(*it)];
+ (*it)->x() = p.x();
+ (*it)->y() = p.y();
+ (*it)->z() = p.z();
+ }
+ }
+
+ delete lsys;
+}
+
+
/*
n3 n23 n2
x-----x-------x
@@ -988,7 +1112,13 @@ struct swap_triangles_pN
!t1->getNumFaceVertices(),
t1->getPolynomialOrder()-1,s);
- optimalLocationPN_ (gf,me, t3, t4,s);
+ //optimalLocationPN_ (gf,me, t3, t4,s);
+ std::vector<MElement*> cavity, old_elems;
+ cavity.push_back(t3);
+ cavity.push_back(t4);
+ old_elems.push_back(t1);
+ old_elems.push_back(t2);
+ s->smooth_cavity(cavity, old_elems, gf);
const double qnew1 = shapeMeasure(t3);
const double qnew2 = shapeMeasure(t4);
@@ -1047,7 +1177,8 @@ static int optimalLocationP2_(GFace *gf,
n12->y() = gp12.y();
n12->z() = gp12.z();
n12->setParameter(0,pp(0));
- n12->setParameter(1,pp(1));
+ n12->setParameter(1,pp(1));
+ printf("Hum, order 2 here...");
return 1;
}
@@ -1151,11 +1282,10 @@ static int swapHighOrderTriangles(GFace *gf,
if (it->second.second){
MTriangle *t1 = dynamic_cast<MTriangle*>(it->second.first);
MTriangle *t2 = dynamic_cast<MTriangle*>(it->second.second);
+
const double qold1 = shapeMeasure(t1);
const double qold2 = shapeMeasure(t2);
- // printf("swap : %g %g\n",qold1,qold2);
-
if (qold1 < 0.6 || qold2 < 0.6)
pairs.insert(swap_triangles_pN(it->first,t1,t2,gf,
edgeVertices,faceVertices,s));
@@ -1174,29 +1304,101 @@ static int swapHighOrderTriangles(GFace *gf,
itp = pairs.begin();
while(itp != pairs.end()){
double diff = fabs(itp->s_before - itp->s_after);
+
+ std::vector<MVertex*> v1,v2,v3,v4;
+ int o1 = itp->t1->getPolynomialOrder();
+ int o2 = itp->t2->getPolynomialOrder();
+ int o3 = itp->t3->getPolynomialOrder();
+ int o4 = itp->t4->getPolynomialOrder();
+
+ itp->t1->getFaceVertices(0,v1);
+ itp->t2->getFaceVertices(0,v2);
+ itp->t3->getFaceVertices(0,v3);
+ itp->t4->getFaceVertices(0,v4);
+ std::vector<MVertex*> ve1(v1.begin()+3,v1.begin()+3*o1);
+ std::vector<MVertex*> ve2(v2.begin()+3,v2.begin()+3*o2);
+ std::vector<MVertex*> ve3(v3.begin()+3,v3.begin()+3*o3);
+ std::vector<MVertex*> ve4(v4.begin()+3,v4.begin()+3*o4);
+ std::vector<MVertex*> vf1(v1.begin()+3*o1,v1.end());
+ std::vector<MVertex*> vf2(v2.begin()+3*o2,v2.end());
+ std::vector<MVertex*> vf3(v3.begin()+3*o3,v3.end());
+ std::vector<MVertex*> vf4(v4.begin()+3*o4,v4.end());
+
+ // for (int ed = 0; ed < 3; ed++) {
+ // std::vector<MVertex*> v1,v2,v3,v4;
+ // itp->t1->getEdgeVertices(ed,v1);
+ // itp->t2->getEdgeVertices(ed,v2);
+ // itp->t3->getEdgeVertices(ed,v3);
+ // itp->t4->getEdgeVertices(ed,v4);
+ // ve1.insert(ve1.end(),v1.begin()+2,v1.end());
+ // ve2.insert(ve2.end(),v2.begin()+2,v2.end());
+ // ve3.insert(ve3.end(),v3.begin()+2,v3.end());
+ // ve4.insert(ve4.end(),v4.begin()+2,v4.end());
+ // }
+
+
if ( t_removed.find(itp->t1) == t_removed.end() &&
t_removed.find(itp->t2) == t_removed.end() &&
itp->quality_new > itp->quality_old &&
diff < 1.e-9){
// itp->print();
+ printf("quality was %f, is now %f\n",itp->quality_old,itp->quality_new);
t_removed.insert(itp->t1);
t_removed.insert(itp->t2);
triangles2.push_back(itp->t3);
triangles2.push_back(itp->t4);
+
+ v_removed.insert(vf1.begin(),vf1.end());
+ v_removed.insert(vf2.begin(),vf2.end());
+ mesh_vertices2.insert(mesh_vertices2.end(),vf3.begin(),vf3.end());
+ mesh_vertices2.insert(mesh_vertices2.end(),vf4.begin(),vf4.end());
+
+ for(std::vector<MVertex*>::iterator vit = ve1.begin(); vit != ve1.end(); vit++) {
+ if (find(ve2.begin(),ve2.begin(),*vit)!=ve2.end())
+ v_removed.insert(*vit);
+ }
+ for(std::vector<MVertex*>::iterator vit = ve3.begin(); vit != ve3.end(); vit++) {
+ if (find(ve4.begin(),ve4.end(),*vit)!=ve4.end())
+ mesh_vertices2.push_back(*vit);
+ }
+
// if (itp->n34 != itp->n12){
- // v_removed.insert(itp->n12);
- // mesh_vertices2.push_back(itp->n34);
+ // v_removed.insert(itp->n12);
+ // mesh_vertices2.push_back(itp->n34);
// }
nbSwap++;
}
else{
+ for(std::vector<MVertex*>::iterator vit = ve1.begin(); vit != ve1.end(); vit++) {
+ if (find(ve2.begin(),ve2.end(),*vit)!=ve2.end())
+ mesh_vertices2.push_back(*vit);
+ }
+ for(std::vector<MVertex*>::iterator vit = ve3.begin(); vit != ve3.end(); vit++) {
+ if (find(ve4.begin(),ve4.end(),*vit)!=ve4.end())
+ v_removed.insert(*vit);
+ }
delete itp->t3;
delete itp->t4;
+
// if (itp->n34 != itp->n12) delete itp->n34;
}
++itp;
}
-
+ int c1=0,c2=0;
+ for (unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
+ if (v_removed.find(gf->mesh_vertices[i]) == v_removed.end()){
+ //mesh_vertices2.push_back(gf->mesh_vertices[i]);
+ c1++;
+ }
+ else {
+ delete gf->mesh_vertices[i];
+ c2++;
+ }
+ }
+ gf->mesh_vertices = mesh_vertices2;
+ printf("Deleted %d vertices from %d\n",c2,gf->mesh_vertices.size());
+ printf("Added %d vertices\n",c1);
+
for (unsigned int i = 0; i < gf->triangles.size(); i++){
if (t_removed.find(gf->triangles[i]) == t_removed.end()){
triangles2.push_back(gf->triangles[i]);
@@ -1205,7 +1407,8 @@ static int swapHighOrderTriangles(GFace *gf,
delete gf->triangles[i];
}
}
- // printf("replacing %d by %d\n",gf->triangles.size(),triangles2.size());
+
+ printf("replacing %d by %d\n",gf->triangles.size(),triangles2.size());
gf->triangles = triangles2;
return nbSwap;
}
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index 3b4379fbc00a6e4678dc1174fb8c6e2c42c3b046..be8cf2f307e2b3236c8f22ae92c3324683a15bd1 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -52,6 +52,9 @@ public:
void smooth_p2point(GFace *);
void smooth_pNpoint(GFace *);
void smooth(GRegion*);
+ void smooth_cavity(std::vector<MElement*> &,
+ std::vector<MElement*> &,
+ GFace *gf);
int getTag() const { return _tag; }
void swap(GFace *,
edgeContainer &edgeVertices,
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 8cf1826f0fc394cbf8781c28951e48cd4388b622..6075ae69a90f94ee99151cbd99e9f86bfe8c0d84 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -102,8 +102,8 @@ static int getAspectRatio (std::vector<MElement *> &elements,
vs.insert(e->getVertex(j));
}
}
- SBoundingBox3d bb;
- SOrientedBoundingBox obbox ;
+ SBoundingBox3d bb;
+ SOrientedBoundingBox obbox ;
std::vector<SPoint3> vertices;
for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index a4bc0b5b97e4826fd920ede441b377ce7723d548..ec701c9c61a06c08a3bacd8b86ea12ba4056bb8e 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -11,6 +11,7 @@
#include "Numeric.h"
#include "polynomialBasis.h"
#include "GmshMessage.h"
+#include <limits>
double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
const qualityMeasure4Triangle &cr)
@@ -293,3 +294,59 @@ double qmDistorsionOfMapping(MTetrahedron *e)
return dmin;
}
+
+double qmTriangleAngles (MTriangle *e) {
+ double a = 100;
+ double worst_quality = std::numeric_limits<double>::max();
+ double mat[3][3];
+ double mat2[3][3];
+ double den = atan(a*(M_PI/9)) + atan(a*(2*M_PI/9 - (M_PI/9)));
+
+ // This matrix is used to "rotate" the triangle to get each vertex
+ // as the "origin" of the mapping in turn
+ double rot[3][3];
+ rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
+ rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
+ rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
+ double tmp[3][3];
+
+ for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
+ const double u = i == 1 ? 1 : 0;
+ const double v = i == 2 ? 1 : 0;
+ const double w = 0;
+ e->getJacobian(u, v, w, mat);
+ e->getPrimaryJacobian(u,v,w,mat2);
+ for (int j = 0; j < i; j++) {
+ matmat(rot,mat,tmp);
+ memcpy(mat, tmp, sizeof(mat));
+ }
+ //get angle
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
+ double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
+ double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
+ norme(v1);
+ norme(v2);
+ norme(v3);
+ norme(v4);
+ double v12[3], v34[3];
+ prodve(v1,v2,v12);
+ prodve(v3,v4,v34);
+ norme(v12);
+ norme(v34);
+ double orientation;
+ prosca(v12,v34,&orientation);
+
+ // If the if the triangle is "flipped" it's no good
+ if (orientation < 0)
+ return -std::numeric_limits<double>::max();
+
+ double c;
+ prosca(v1,v2,&c);
+ double x = acos(c)-M_PI/3;
+ double quality = (atan(a*(x+M_PI/9)) + atan(a*(2*M_PI/9 - (x+M_PI/9))))/den;
+ worst_quality = std::min(worst_quality, quality);
+ }
+
+ return worst_quality;
+}
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 0daf095b5021429e87f0b4eea2722e1297f52ccb..cacf55add0a5e8d777fc16daae34f6acafd961b2 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -19,6 +19,8 @@ enum qualityMeasure4Tet{QMTET_1, QMTET_2, QMTET_3, QMTET_ONE, QMTET_COND};
double qmDistorsionOfMapping (MTriangle *e);
double qmDistorsionOfMapping (MTetrahedron *e);
+double qmTriangleAngles (MTriangle *e);
+
double qmTriangle(MTriangle *f, const qualityMeasure4Triangle &cr);
double qmTriangle(BDS_Face *f, const qualityMeasure4Triangle &cr);
double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3,
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 3f110c426f955324b1557134b23fcde4a789d1a8..3de18b979b942f1ee731ba5a546db3511814e13c 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -43,6 +43,19 @@ void matvec(double mat[3][3], double vec[3], double res[3])
res[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2];
}
+void matmat(double mat1[3][3], double mat2[3][3], double res[3][3])
+{
+ res[0][0] = mat1[0][0]*mat2[0][0] + mat1[0][1]*mat2[1][0] + mat1[0][2]*mat2[2][0];
+ res[0][1] = mat1[0][0]*mat2[0][1] + mat1[0][1]*mat2[1][1] + mat1[0][2]*mat2[2][1];
+ res[0][2] = mat1[0][0]*mat2[0][2] + mat1[0][1]*mat2[1][2] + mat1[0][2]*mat2[2][2];
+ res[1][0] = mat1[1][0]*mat2[0][0] + mat1[1][1]*mat2[1][0] + mat1[1][2]*mat2[2][0];
+ res[1][1] = mat1[1][0]*mat2[0][1] + mat1[1][1]*mat2[1][1] + mat1[1][2]*mat2[2][1];
+ res[1][2] = mat1[1][0]*mat2[0][2] + mat1[1][1]*mat2[1][2] + mat1[1][2]*mat2[2][2];
+ res[2][0] = mat1[2][0]*mat2[0][0] + mat1[2][1]*mat2[1][0] + mat1[2][2]*mat2[2][0];
+ res[2][1] = mat1[2][0]*mat2[0][1] + mat1[2][1]*mat2[1][1] + mat1[2][2]*mat2[2][1];
+ res[2][2] = mat1[2][0]*mat2[0][2] + mat1[2][1]*mat2[1][2] + mat1[2][2]*mat2[2][2];
+}
+
void normal3points(double x0, double y0, double z0,
double x1, double y1, double z1,
double x2, double y2, double z2,
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 29db8d3d60a1da570bb6415f37759ec2370f862f..9678825457fccf2761f8591f74e6865389ffd83e 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -29,6 +29,7 @@ inline void prosca(double a[3], double b[3], double *c)
*c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
void matvec(double mat[3][3], double vec[3], double res[3]);
+void matmat(double mat1[3][3], double mat2[3][3], double res[3][3]);
inline double norm3(double a[3])
{
return sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);