From a029b66711d752ff039b4fb4976fd19fc1e48eec Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 1 Feb 2012 11:32:23 +0000
Subject: [PATCH] Worked a lot oncurvilinear meshing : projection of point on
surafces work, matching meshes andgeometry works better, added some test
examples...
---
Common/OpenFile.cpp | 2 +-
Fltk/statisticsWindow.cpp | 329 ++++++++++++++--------------
Geo/GFace.cpp | 74 +++++--
Geo/Geo.cpp | 2 +-
Geo/GeoInterpolation.cpp | 79 ++++++-
Geo/GeomMeshMatcher.cpp | 37 ++--
Geo/MElement.cpp | 11 +-
Geo/gmshEdge.cpp | 4 +-
Geo/gmshFace.cpp | 5 +-
Geo/gmshFace.h | 2 -
Geo/gmshVertex.cpp | 37 ++++
Mesh/BackgroundMesh.cpp | 20 ++
Mesh/HighOrder.cpp | 80 ++++++-
Mesh/meshGEdge.cpp | 1 +
Mesh/meshGFaceDelaunayInsertion.cpp | 64 ++++++
Mesh/meshGFaceOptimize.cpp | 11 +-
Mesh/qualityMeasures.cpp | 21 +-
Solver/elasticitySolver.cpp | 2 +-
Solver/elasticitySolver.h | 2 +-
benchmarks/2d/recombine.geo | 2 +-
benchmarks/2d/transfinite_b.geo | 2 +-
21 files changed, 543 insertions(+), 244 deletions(-)
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index a1881b2695..eba6ade46b 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -366,7 +366,7 @@ int MergeFile(std::string fileName, bool warnIfMissing)
// mesh matcher
if(CTX::instance()->geom.matchGeomAndMesh && !GModel::current()->empty()) {
GModel* tmp2 = GModel::current();
- GModel* tmp = new GModel();
+ GModel* tmp = new GModel();
tmp->readMSH(fileName);
status = GeomMeshMatcher::instance()->match(tmp2, tmp);
delete tmp;
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 89cf92139d..789230bfdf 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -199,183 +199,188 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
void statisticsWindow::compute(bool elementQuality)
{
-#if 0
- // MINIMUM ANGLES
- double minAngle = 120.0;
- double meanAngle = 0.0;
- int count = 0;
- std::vector<GEntity*> entities;
- GModel::current()->getEntities(entities);
- std::map<int, std::vector<double> > d;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- double angle = e->angleShapeMeasure();
- minAngle = std::min(minAngle, angle);
- meanAngle += angle;
- count++;
+#if 1
+ {
+ // MINIMUM ANGLES
+ double minAngle = 120.0;
+ double meanAngle = 0.0;
+ int count = 0;
+ std::vector<GEntity*> entities;
+ GModel::current()->getEntities(entities);
+ std::map<int, std::vector<double> > d;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ if(entities[i]->dim() < 2) continue;
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ double angle = e->angleShapeMeasure();
+ minAngle = std::min(minAngle, angle);
+ meanAngle += angle;
+ count++;
+ }
}
+ meanAngle = meanAngle / count;
+ printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
}
- meanAngle = meanAngle / count;
- printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
- // MESH DEGREE VERTICES
- std::set<MEdge, Less_Edge> edges;
- GModel::current()->getEntities(entities);
- std::map<MVertex*, int > vert2Deg;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2 ) continue;
- if(entities[i]->tag() < 100) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- for(unsigned int k = 0; k < e->getNumEdges(); k++){
- edges.insert(e->getEdge(k));
- }
- for(unsigned int k = 0; k < e->getNumVertices(); k++){
- MVertex *v = e->getVertex(k);
- if (v->onWhat()->dim() < 2) continue;
- std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
- if (it == vert2Deg.end()) {
- vert2Deg.insert(std::make_pair(v,1));
- }
- else{
- int nbE = it->second+1;
- it->second = nbE;
- }
+ {
+ // MESH DEGREE VERTICES
+ std::vector<GEntity*> entities;
+ std::set<MEdge, Less_Edge> edges;
+ GModel::current()->getEntities(entities);
+ std::map<MVertex*, int > vert2Deg;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ if(entities[i]->dim() < 2 ) continue;
+ // if(entities[i]->tag() < 100) continue;
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ edges.insert(e->getEdge(k));
+ }
+ for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ MVertex *v = e->getVertex(k);
+ if (v->onWhat()->dim() < 2) continue;
+ std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ if (it == vert2Deg.end()) {
+ vert2Deg.insert(std::make_pair(v,1));
+ }
+ else{
+ int nbE = it->second+1;
+ it->second = nbE;
+ }
+ }
}
}
- }
- int dMin = 10;
- int dMax = 0;
- int d4 = 0;
- int nbElems = vert2Deg.size();
- std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
- for(; itmap !=vert2Deg.end(); itmap++){
- MVertex *v = itmap->first;
- int nbE = itmap->second;
- dMin = std::min(nbE, dMin);
- dMax = std::max(nbE, dMax);
- if (nbE == 4) d4 += 1;
- }
- if (nbElems > 0)
- printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
- dMin, dMax, (double)d4/nbElems);
-
- FieldManager *fields = GModel::current()->getFields();
- Field *f = fields->get(fields->background_field);
- int nbEdges = edges.size();
- printf("nb edges =%d \n", nbEdges);
- system("rm qualEdges.txt");
- FILE *fp = fopen("qualEdges.txt", "w");
- std::vector<int> qualE;
- int nbS = 50;
- qualE.resize(nbS);
- if(fields->background_field > 0){
- printf("found field \n");
- std::set<MEdge, Less_Edge>::iterator it = edges.begin();
- double sum = 0;
- for (; it !=edges.end();++it){
- MVertex *v0 = it->getVertex(0);
- MVertex *v1 = it->getVertex(1);
- double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
- (v0->y()-v1->y())*(v0->y()-v1->y())+
- (v0->z()-v1->z())*(v0->z()-v1->z()));
- double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
- 0.5*(v0->z()+v1->z()),v0->onWhat());
- double el = l/lf;
- int index = (int) ceil(el*nbS*0.5);
- qualE[index]+= 1;
- double e = (l>lf) ? lf/l : l/lf;
- sum += e - 1.0;
+ int dMin = 10;
+ int dMax = 0;
+ int d4 = 0;
+ int nbElems = vert2Deg.size();
+ std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ for(; itmap !=vert2Deg.end(); itmap++){
+ MVertex *v = itmap->first;
+ int nbE = itmap->second;
+ dMin = std::min(nbE, dMin);
+ dMax = std::max(nbE, dMax);
+ if (nbE == 4) d4 += 1;
}
- double tau = exp ((1./edges.size()) * sum);
- printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
-
- double ibegin = 2./(2*nbS);
- double inext = 2./nbS;
- for (int i= 0; i< qualE.size(); i++){
- fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ if (nbElems > 0)
+ printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
+ dMin, dMax, (double)d4/nbElems);
+ FieldManager *fields = GModel::current()->getFields();
+ Field *f = fields->get(fields->background_field);
+ int nbEdges = edges.size();
+ printf("nb edges =%d \n", nbEdges);
+ system("rm qualEdges.txt");
+ FILE *fp = fopen("qualEdges.txt", "w");
+ std::vector<int> qualE;
+ int nbS = 50;
+ qualE.resize(nbS);
+ if(fields->background_field > 0){
+ printf("found field \n");
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ double sum = 0;
+ for (; it !=edges.end();++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
+ (v0->y()-v1->y())*(v0->y()-v1->y())+
+ (v0->z()-v1->z())*(v0->z()-v1->z()));
+ double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
+ 0.5*(v0->z()+v1->z()),v0->onWhat());
+ double el = l/lf;
+ int index = (int) ceil(el*nbS*0.5);
+ qualE[index]+= 1;
+ double e = (l>lf) ? lf/l : l/lf;
+ sum += e - 1.0;
+ }
+ double tau = exp ((1./edges.size()) * sum);
+ printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
+
+ double ibegin = 2./(2*nbS);
+ double inext = 2./nbS;
+ for (int i= 0; i< qualE.size(); i++){
+ fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ }
+
}
-
+ fclose(fp);
}
- fclose(fp);
-
- std::vector<GEntity*> entities;
- std::set<MEdge, Less_Edge> edges;
- GModel::current()->getEntities(entities);
- std::map<MVertex*, int > vert2Deg;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2 ) continue;
- if(entities[i]->tag() != 10) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- for(unsigned int k = 0; k < e->getNumEdges(); k++){
- edges.insert(e->getEdge(k));
+ {
+ std::vector<GEntity*> entities;
+ std::set<MEdge, Less_Edge> edges;
+ GModel::current()->getEntities(entities);
+ std::map<MVertex*, int > vert2Deg;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ if(entities[i]->dim() < 2 ) continue;
+ if(entities[i]->tag() != 10) continue;
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ edges.insert(e->getEdge(k));
+ }
+ for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ MVertex *v = e->getVertex(k);
+ if (v->onWhat()->dim() < 2) continue;
+ std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ if (it == vert2Deg.end()){
+ vert2Deg.insert(std::make_pair(v,1));
}
- for(unsigned int k = 0; k < e->getNumVertices(); k++){
- MVertex *v = e->getVertex(k);
- if (v->onWhat()->dim() < 2) continue;
- std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
- if (it == vert2Deg.end()){
- vert2Deg.insert(std::make_pair(v,1));
- }
- else{
- int nbE = it->second+1;
- it->second = nbE;
- }
+ else{
+ int nbE = it->second+1;
+ it->second = nbE;
}
}
- }
- int dMin = 10;
- int dMax = 0;
- int d4 = 0;
- int nbElems = vert2Deg.size();
- std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
- for(; itmap !=vert2Deg.end(); itmap++){
- MVertex *v = itmap->first;
- int nbE = itmap->second;
- dMin = std::min(nbE, dMin);
- dMax = std::max(nbE, dMax);
- if (nbE == 4) d4 += 1;
- }
- if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
- dMin, dMax, (double)d4/nbElems);
- FieldManager *fields = GModel::current()->getFields();
- Field *f = fields->get(fields->background_field);
- int nbEdges = edges.size();
- system("rm qualEdges.txt");
- FILE *fp = fopen("qualEdges.txt", "w");
- std::vector<int> qualE;
- int nbS = 50;
- qualE.resize(nbS);
- if(fields->background_field > 0){
- std::set<MEdge, Less_Edge>::iterator it = edges.begin();
- double sum = 0;
- for (; it !=edges.end();++it){
- MVertex *v0 = it->getVertex(0);
- MVertex *v1 = it->getVertex(1);
- double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
- (v0->y()-v1->y())*(v0->y()-v1->y())+
- (v0->z()-v1->z())*(v0->z()-v1->z()));
- double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
- 0.5*(v0->z()+v1->z()),v0->onWhat());
- double el = l/lf;
- int index = (int) ceil(el*nbS*0.5);
- qualE[index]+= 1;
- double e = (l>lf) ? lf/l : l/lf;
- sum += e - 1.0;
+ }
+ }
+ int dMin = 10;
+ int dMax = 0;
+ int d4 = 0;
+ int nbElems = vert2Deg.size();
+ std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ for(; itmap !=vert2Deg.end(); itmap++){
+ MVertex *v = itmap->first;
+ int nbE = itmap->second;
+ dMin = std::min(nbE, dMin);
+ dMax = std::max(nbE, dMax);
+ if (nbE == 4) d4 += 1;
}
- double tau = exp ((1./edges.size()) * sum);
- double ibegin = 2./(2*nbS);
- double inext = 2./nbS;
- for (int i= 0; i< qualE.size(); i++){
- fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
+ dMin, dMax, (double)d4/nbElems);
+ FieldManager *fields = GModel::current()->getFields();
+ Field *f = fields->get(fields->background_field);
+ int nbEdges = edges.size();
+ system("rm qualEdges.txt");
+ FILE *fp = fopen("qualEdges.txt", "w");
+ std::vector<int> qualE;
+ int nbS = 50;
+ qualE.resize(nbS);
+ if(fields->background_field > 0){
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ double sum = 0;
+ for (; it !=edges.end();++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
+ (v0->y()-v1->y())*(v0->y()-v1->y())+
+ (v0->z()-v1->z())*(v0->z()-v1->z()));
+ double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
+ 0.5*(v0->z()+v1->z()),v0->onWhat());
+ double el = l/lf;
+ int index = (int) ceil(el*nbS*0.5);
+ qualE[index]+= 1;
+ double e = (l>lf) ? lf/l : l/lf;
+ sum += e - 1.0;
+ }
+ double tau = exp ((1./edges.size()) * sum);
+ double ibegin = 2./(2*nbS);
+ double inext = 2./nbS;
+ for (int i= 0; i< qualE.size(); i++){
+ fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ }
}
+ fclose(fp);
}
- fclose(fp);
#endif
-
+
int num = 0;
static double s[50];
static char label[50][256];
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 455a53e182..53cefb61b4 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -908,7 +908,7 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p, bool onSurface) const
return SPoint2(U, V);
}
-#if defined(QUASINEWTON)
+#if defined(HAVE_BFGS)
// Callback function for BFGS
void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr) {
data_wrapper* w = static_cast<data_wrapper*>(ptr);
@@ -917,26 +917,29 @@ void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d
// Value of the objective
GPoint pnt = gf->point(x[0], x[1]);
- SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
- func = p.distance(spnt);
+ func = 0.5 *
+ ( p.x() - pnt.x() ) * ( p.x() - pnt.x() ) +
+ ( p.y() - pnt.y() ) * ( p.y() - pnt.y() ) +
+ ( p.z() - pnt.z() ) * ( p.z() - pnt.z() ) ;
//printf("func : %f\n", func);
// Value of the gradient
- std::vector<double> values(2);
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(x[0], x[1]));
- grad[0] = -2./(2.0*func) * (der.left().x() + der.left().y() + der.left().z());
- grad[1] = -2./(2.0*func) * (der.right().x() + der.right().y() + der.right().z());
- //printf("Gradients %f %f\n", grad[0], grad[1]);
+ grad[0] = -(p.x() - pnt.x()) * der.left().x() - (p.y() - pnt.y()) * der.left().y() - (p.z() - pnt.z()) * der.left().z();
+ grad[1] = -(p.x() - pnt.x()) * der.right().x() - (p.y() - pnt.y()) * der.right().y() - (p.z() - pnt.z()) * der.right().z();
+ // printf("func %22.15E Gradients %22.15E %22.15E der %g %g %g\n",func, grad[0], grad[1],der.left().x(),der.left().y(),der.left().z());
}
#endif
GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
{
-#if defined(QUASINEWTON)
-
+#if defined(HAVE_BFGS)
//
// Creating the optimisation problem
//
+
+ // printf("STARTING OPTIMIZATION\n");
+
alglib::ae_int_t dim = 2;
alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7]
alglib::minlbfgsstate state;
@@ -950,16 +953,21 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
Range<double> uu = parBounds(0);
Range<double> vv = parBounds(1);
- double initial_guesses = 1000.0;
- for(double u = uu.low(); u<=uu.high(); u+=(uu.high()-uu.low())/initial_guesses) {
- printf("%f\n", u);
- for(double v = vv.low(); v<=vv.high(); v+=(vv.high()-vv.low())/initial_guesses) {
+ // FILE *F = fopen ("hop.pos","w");
+ // fprintf(F,"View \" \" {\n");
+ // fprintf(F,"SP(%g,%g,%g) {%g};\n",queryPoint.x(),queryPoint.y(),queryPoint.z(),0.0);
+ double initial_guesses = 10.0;
+ for(double u = uu.low(); u<=uu.high()+1.e-5; u+=(uu.high()-uu.low())/initial_guesses) {
+ // printf("%f\n", u);
+ for(double v = vv.low(); v<=vv.high()+1.e-5; v+=(vv.high()-vv.low())/initial_guesses) {
GPoint pnt = point(u, v);
- printf("(point) : %f %f %f\n", pnt.x(), pnt.y(), pnt.z());
SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
double dist = queryPoint.distance(spnt);
+ // fprintf(F,"SP(%g,%g,%g) {%g};\n",pnt.x(), pnt.y(), pnt.z(),dist);
+ // printf("lmocal (%12.5E %12.5E) (point) : %12.5E %12.5E %12.5E (query) : %12.5E %12.5E %12.5E DSIT %12.5E\n",u,v, pnt.x(), pnt.y(), pnt.z(),queryPoint.x(),queryPoint.y(),queryPoint.z(),
+ // dist);
if (dist < min_dist) {
- printf("min_dist %f\n", dist);
+ // printf("min_dist %f\n", dist);
min_dist = dist;
min_u = u;
min_v = v;
@@ -967,25 +975,27 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
}
}
}
-
+ // fprintf(F,"};\n");
+ // fclose(F);
+ // getchar();
initial_conditions[0] = min_u;
initial_conditions[1] = min_v;
- printf("Initial conditions : %f %f\n", min_u, min_v);
+ // printf("Initial conditions : %f %f %12.5E\n", min_u, min_v,min_dist);
GPoint pnt = point(min_u, min_v);
- printf("Initial conditions (point) : %f %f %f\n", pnt.x(), pnt.y(), pnt.z());
+ // printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f DIST = %12.5E\n",
+ // pnt.x(), pnt.y(), pnt.z(),min_u,min_v,
+ // queryPoint.x(),queryPoint.y(),queryPoint.z(),
+ // min_dist);
x.setcontent(dim, &initial_conditions[0]);
-
-
-
minlbfgscreate(2, corr, x, state);
//
// Defining the stopping criterion
//
- double epsg = 0.0;
+ double epsg = 1.e-12;
double epsf = 0.0;
double epsx = 0.0;
alglib::ae_int_t maxits = 500;
@@ -1009,7 +1019,25 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
minlbfgsresults(state,x,rep);
- return point(x[0],x[1]);
+ GPoint pntF = point(x[0], x[1]);
+ if (rep.terminationtype != 4){
+ // printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f DIST = %12.5E at the end (%f %f) %f %f %f\n",
+ // pnt.x(), pnt.y(), pnt.z(),min_u,min_v,
+ // queryPoint.x(),queryPoint.y(),queryPoint.z(),
+ // min_dist,x[0],x[1],pntF.x(), pntF.y(), pntF.z());
+ // double DDD =
+ // ( queryPoint.x() - pntF.x()) * ( queryPoint.x() - pntF.x()) +
+ // ( queryPoint.y() - pntF.y()) * ( queryPoint.y() - pntF.y()) +
+ // ( queryPoint.z() - pntF.z()) * ( queryPoint.z() - pntF.z());
+ // if (sqrt(DDD) > 1.e-4)
+ /*
+ printf("Initial conditions (point) : %f %f %f local (%g %g) Looking for %f %f %f DIST = %12.5E at the end (%f %f) %f %f %f termination %d\n",
+ pnt.x(), pnt.y(), pnt.z(),min_u,min_v,
+ queryPoint.x(),queryPoint.y(),queryPoint.z(),
+ min_dist,x[0],x[1],pntF.x(), pntF.y(), pntF.z(),rep.terminationtype);
+ */
+ }
+ return pntF;
#else
Msg::Error("Closest point not implemented for this type of surface");
SPoint2 p = parFromPoint(queryPoint, false);
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index c17f99e359..607b603716 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3259,7 +3259,7 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
p = InterpolateSurface(s, x(0), x(1), 0, 0);
uv[0] = x(0);
uv[1] = x(1);
- if (ITER > 0)
+ if (ITER >= 0)
Msg::Info("ProjectPoint (%g,%g,%g) On Surface %d converged after %d iterations",
p.Pos.X, p.Pos.Y, p.Pos.Z, s->Num, ITER);
return true;
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 464dcd318f..ed1e23b88a 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -295,8 +295,8 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
switch (c->Typ) {
case MSH_SEGM_LINE:
-#if defined(QUASINEWTON)
- printf("MSH_SEGM_LINE\n");
+#if defined(HAVE_BFGS)
+ // printf("MSH_SEGM_LINE\n");
#endif
N = List_Nbr(c->Control_Points);
i = (int)((double)(N - 1) * u);
@@ -418,8 +418,8 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
break;
case MSH_SEGM_DISCRETE:
-#if defined(QUASINEWTON)
- printf("MSH_SEGM_DISCRETE\n");
+#if defined(HAVE_BFGS)
+ // printf("MSH_SEGM_DISCRETE\n");
#endif
Msg::Debug("Cannot interpolate discrete curve");
@@ -463,9 +463,41 @@ static Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
// Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4
// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
+//
+// s3(1,1)
+// +
+// / |
+// c3 / |
+// / | c2
+// / |
+// / c1 |
+// +----------+
+// s1(0,0) s2(1,0)
+
+// u = v = 0 -----> x = c1(0) = s1 --> OK
+// u = 1 ; v = 0 -----> x = c2(0) + c1(1) - s2 = s2 --> OK
+// u = 1 ; v = 1 -----> x = c2(1) + c3(1) - s3 = s3 --> OK
+// v = 0 --> u s2 + c1(u) - u s2 --> x = c1(u) --> OK
+// u = 1 --> c2(v) + (1-v) s2 + v s3 -(1-v) s2 - v s3 --> x = c2(v) --> OK
+// u = 0 --> (1-v) s1 + v s1 = s1 !!! not terrible (singular)
+
+// Transfinite approximation on a triangle
+
+// f(u,v) = (1-u) (c1(u-v) + c3(1-v) - s1) +
+// (u-v) (c2(v) + c1(u) - s2) +
+// v (c3(1-u) + c2(1-u+v) - s3)
+//
+// u = v = 0 --> f(0,0) = s1 + s1 - s1 = s1
+// u = v = 1 --> f(1,1) = s3 + s3 - s3 = s3
+// u = 1 ; v = 0 --> f(1,1) = s2 + s2 - s2 = s2
+// v = 0 --> (1-u)c1(u) + u c1(u) = c1(u) --> COOL
+
#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
+#define TRAN_TRIB(c1,c1b,c2,c2b,c3,c3b,s1,s2,s3,u,v) (1.-u)*(c1+c3b-s1)+(u-v)*(c2+c1b-s2)+v*(c3+c2b-s3);
+
+
static Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v)
{
@@ -480,6 +512,23 @@ static Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, u, v);
return V;
}
+static Vertex TransfiniteTriB(Vertex c1, Vertex c1b, Vertex c2,
+ Vertex c2b, Vertex c3, Vertex c3b,
+ Vertex s1, Vertex s2, Vertex s3,
+ double u, double v)
+{
+ Vertex V;
+ V.lc = TRAN_TRIB(c1.lc,c1b.lc, c2.lc,c2b.lc, c3.lc, c3b.lc, s1.lc, s2.lc, s3.lc, u, v);
+ V.w = TRAN_TRIB(c1.w, c1b.w, c2.w,c2b.w, c3.w,c3b.w, s1.w, s2.w, s3.w, u, v);
+ V.Pos.X = TRAN_TRIB(c1.Pos.X, c1b.Pos.X,c2.Pos.X, c2b.Pos.X,c3.Pos.X,c3b.Pos.X,
+ s1.Pos.X, s2.Pos.X, s3.Pos.X, u, v);
+ V.Pos.Y = TRAN_TRIB(c1.Pos.Y, c1b.Pos.Y,c2.Pos.Y, c2b.Pos.Y,c3.Pos.Y,c3b.Pos.Y,
+ s1.Pos.Y, s2.Pos.Y, s3.Pos.Y, u, v);
+ V.Pos.Z = TRAN_TRIB(c1.Pos.Z, c1b.Pos.Z,c2.Pos.Z, c2b.Pos.Z,c3.Pos.Z,c3b.Pos.Z,
+ s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, u, v);
+ return V;
+}
+
static void TransfiniteSph(Vertex S, Vertex center, Vertex *T)
{
@@ -592,7 +641,7 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
}
}
- Vertex *S[4], V[4], T;
+ Vertex *S[4], V[4],VB[3], T;
if(s->Typ == MSH_SURF_REGL && List_Nbr(s->Generatrices) >= 4){
S[0] = C[0]->beg;
S[1] = C[1]->beg;
@@ -602,8 +651,9 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
-#if defined(QUASINEWTON)
- printf("----- u %f v %f\n", u, v);
+ /*
+#if defined(HAVE_BFGS)
+ printf("----- u %f v %f %g %g\n", u, v,C[1]->ubeg,C[1]->uend);
printf("S[0] %f %f %f\n", S[0]->Pos.X, S[0]->Pos.Y,S[0]->Pos.Z);
printf("S[1] %f %f %f\n", S[1]->Pos.X, S[1]->Pos.Y,S[1]->Pos.Z);
printf("S[2] %f %f %f\n", S[2]->Pos.X, S[2]->Pos.Y,S[2]->Pos.Z);
@@ -613,6 +663,7 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
printf("V[2] %f %f %f\n", V[2].Pos.X, V[2].Pos.Y,V[2].Pos.Z);
printf("V[3] %f %f %f\n", V[3].Pos.X, V[3].Pos.Y,V[3].Pos.Z);
#endif
+ */
T = TransfiniteQua(V[0], V[1], V[2], V[3], *S[0], *S[1], *S[2], *S[3], u, v);
if(isSphere) TransfiniteSph(*S[0], *O, &T);
}
@@ -620,11 +671,23 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
S[0] = C[0]->beg;
S[1] = C[1]->beg;
S[2] = C[2]->beg;
+ /*
V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
T = TransfiniteTri(V[0], V[1], V[2], *S[0], *S[1], *S[2], u, v);
- if(isSphere) TransfiniteSph(*S[0], *O, &T);
+ */
+ V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * (u-v), 0);
+ V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
+ V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
+ VB[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
+ VB[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * (1-u+v), 0);
+ VB[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - v), 0);
+ T = TransfiniteTriB(V[0],VB[0], V[1],VB[1], V[2],VB[2], *S[0], *S[1], *S[2], u, v);
+
+ if(isSphere) {
+ TransfiniteSph(*S[0], *O, &T);
+ }
}
return T;
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index e498931d95..76f9ec1030 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -132,7 +132,7 @@ GeomMeshMatcher::matchVertices(GModel* m1, GModel *m2, bool& ok)
}
if (best_score != DBL_MAX) {
- Msg::Info("Vertices %i (geom) and %i (mesh) match.\n",
+ Msg::Info("Model Vertex %i (geom) and %i (mesh) match",
(*entity1)->tag(),
best_candidate_ge->tag());
@@ -154,7 +154,7 @@ GeomMeshMatcher::matchVertices(GModel* m1, GModel *m2, bool& ok)
//vert++)
//m2->add(*vert);
- Msg::Info("Vertices matched : %i / %i (%i)\n", num_matched_vertices, num_total_vertices, num_mesh_vertices);
+ Msg::Info("Vertices matched : %i / %i (%i)", num_matched_vertices, num_total_vertices, num_mesh_vertices);
//if(num_matched_vertices != num_total_vertices) ok = false;
return (coresp_v);
}
@@ -570,7 +570,7 @@ int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
}
}
}
- printf("creating faces\n");
+ // printf("creating faces\n");
for (GModel::fiter it = mesh->firstFace(); it != mesh->lastFace(); ++it){
for (int i=0;i<(*it)->triangles.size();i++){
MVertex *v1 = geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(0)->getNum());
@@ -608,6 +608,7 @@ static void copy_vertices (GVertex *to, GVertex *from, std::map<MVertex*,MVertex
}
}
static void copy_vertices (GRegion *to, GRegion *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+ to->deleteMesh();
for (int i=0;i<from->mesh_vertices.size();i++){
MVertex *v_from = from->mesh_vertices[i];
MVertex *v_to = new MVertex (v_from->x(),v_from->y(),v_from->z(), to);
@@ -621,23 +622,30 @@ static void copy_vertices (GEdge *to, GEdge *from, std::map<MVertex*,MVertex*> &
MVertex *v_from = from->mesh_vertices[i];
double t;
GPoint gp = to->closestPoint(SPoint3(v_from->x(),v_from->y(),v_from->z()), t );
- MEdgeVertex *v_to = new MEdgeVertex (v_from->x(),v_from->y(),v_from->z(), to, t );
+ MEdgeVertex *v_to = new MEdgeVertex (gp.x(),gp.y(),gp.z(), to, gp.u() );
to->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
}
+ // printf("Ending Edge %d %d vertices to match\n",from->tag(),from->mesh_vertices.size());
}
-static void copy_vertices (GFace *to, GFace *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
- printf("Starting Face %d, with %d vertices\n", to->tag(), from->mesh_vertices.size());
- for (int i=0;i<from->mesh_vertices.size();i++){
- MVertex *v_from = from->mesh_vertices[i];
+static void copy_vertices (GFace *geom, GFace *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+ // printf("Starting Face %d, with %d vertices\n", geom->tag(), mesh->mesh_vertices.size());
+ for (int i=0;i<mesh->mesh_vertices.size();i++){
+ MVertex *v_from = mesh->mesh_vertices[i];
double uv[2];
- GPoint gp = to->closestPoint ( SPoint3(v_from->x(),v_from->y(),v_from->z()), uv );
- printf("Original point %f %f %f\n", v_from->x(), v_from->y(), v_from->z());
- printf("New point %f %f %f\n", gp.x(), gp.y(), gp.z());
- MFaceVertex *v_to = new MFaceVertex (v_from->x(),v_from->y(),v_from->z(), to, uv[0],uv[1] );
- to->mesh_vertices.push_back(v_to);
+ GPoint gp = geom->closestPoint ( SPoint3(v_from->x(),v_from->y(),v_from->z()), uv );
+ // printf("Original point %f %f %f\n", v_from->x(), v_from->y(), v_from->z());
+ // printf("New point %f %f %f\n", gp.x(), gp.y(), gp.z());
+ double DDD = ( v_from->x() - gp.x()) * ( v_from->x() - gp.x()) +
+ ( v_from->y() - gp.y()) * ( v_from->y() - gp.y()) +
+ ( v_from->z() - gp.z()) * ( v_from->z() - gp.z()) ;
+ if (sqrt(DDD) > 1.e-3)Msg::Error("Impossible to match one point Original point %f %f %f New point %f %f %f",
+ v_from->x(), v_from->y(), v_from->z(),gp.x(), gp.y(), gp.z());
+ MFaceVertex *v_to = new MFaceVertex (v_from->x(),v_from->y(),v_from->z(), geom, gp.u(),gp.v() );
+ geom->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
}
+ // printf("Ending Face %d %d vertices to match\n",geom->tag(),geom->mesh_vertices.size());
}
template <class ELEMENT>
@@ -669,7 +677,7 @@ void copy_vertices (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_me
void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom){
// copy all elements
for (GModel::viter vit = geom->firstVertex() ; vit != geom->lastVertex(); ++vit)
- copy_elements<MPoint>((*vit)->points,mesh->getVertexByTag((*vit)->tag())->points,_mesh_to_geom);
+ if (mesh->getVertexByTag((*vit)->tag())) copy_elements<MPoint>((*vit)->points,mesh->getVertexByTag((*vit)->tag())->points,_mesh_to_geom);
for (GModel::eiter eit = geom->firstEdge() ; eit != geom->lastEdge(); ++eit)
copy_elements<MLine>((*eit)->lines,mesh->getEdgeByTag((*eit)->tag())->lines,_mesh_to_geom);
for (GModel::fiter fit = geom->firstFace() ; fit != geom->lastFace(); ++fit) {
@@ -688,6 +696,7 @@ void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_me
int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
{
mesh->createTopologyFromMesh();
+ GModel::setCurrent(geom);
bool ok = true;
// This will match VERTICES
std::vector<Pair<GVertex*, GVertex*> > *coresp_v = matchVertices(geom, mesh,ok);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index a47a61cc9e..5881e05096 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -271,13 +271,14 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
double gsf[1256][3];
getGradShapeFunctions(u, v, w, gsf);
for (int i = 0; i < getNumShapeFunctions(); i++) {
- const MVertex *v = getShapeFunctionNode(i);
+ const MVertex *ver = getShapeFunctionNode(i);
double* gg = gsf[i];
- for (int j = 0; j < 3; j++) {
- jac[j][0] += v->x() * gg[j];
- jac[j][1] += v->y() * gg[j];
- jac[j][2] += v->z() * gg[j];
+ for (int j = 0; j < getDim(); j++) {
+ jac[j][0] += ver->x() * gg[j];
+ jac[j][1] += ver->y() * gg[j];
+ jac[j][2] += ver->z() * gg[j];
}
+ // printf("GSF (%d,%g %g) = %g %g \n",i,u,v,gg[0],gg[1]);
}
return _computeDeterminantAndRegularize(this, jac);
}
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 41a1e62572..956763aae2 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -270,11 +270,11 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const
}
else if (C[2]->Num == c->Num) {
U = 1-(epar - C[2]->ubeg) / (C[2]->uend - C[2]->ubeg) ;
- V = 1;
+ V = U;
}
else if (C[2]->Num == -c->Num) {
U = 1-(C[2]->uend - epar - C[2]->ubeg) / (C[2]->uend - C[2]->ubeg) ;
- V = 1;
+ V = U;
}
else{
Msg::Info("Reparameterizing edge %d on face %d", c->Num, s->Num);
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index d2ae3d350d..4057afbdbd 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -217,9 +217,11 @@ GPoint gmshFace::point(double par1, double par2) const
}
}
-#if not defined(QUASINEWTON)
GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
{
+#if defined(HAVE_BFGS)
+ return GFace::closestPoint(qp, initialGuess);
+#endif
if (s->Typ == MSH_SURF_PLAN && !s->geometry){
double XP = qp.x();
double YP = qp.y();
@@ -258,7 +260,6 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
return GPoint(-1.e22, -1.e22, -1.e22, 0, u);
return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
}
-#endif
SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
{
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index ecd59489b9..b70242b3ef 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -22,10 +22,8 @@ class gmshFace : public GFace {
Range<double> parBounds(int i) const;
void setModelEdges(std::list<GEdge*> &);
virtual GPoint point(double par1, double par2) const;
-#if not defined(QUASINEWTON)
virtual GPoint closestPoint(const SPoint3 &queryPoint,
const double initialGuess[2]) const;
-#endif
virtual bool containsPoint(const SPoint3 &pt) const;
virtual double getMetricEigenvalue(const SPoint2 &);
virtual SVector3 normal(const SPoint2 ¶m) const;
diff --git a/Geo/gmshVertex.cpp b/Geo/gmshVertex.cpp
index b4bf4ac2cd..985580e0ee 100644
--- a/Geo/gmshVertex.cpp
+++ b/Geo/gmshVertex.cpp
@@ -39,8 +39,11 @@ GEntity::GeomType gmshVertex::geomType() const
SPoint2 gmshVertex::reparamOnFace(const GFace *face, int dir) const
{
+
Surface *s = (Surface*)face->getNativePtr();
+ // printf("coucou HERE I AM %d\n",s->Typ);
+
if(s->geometry){
// It is not always right if it is periodic.
if(l_edges.size() == 1 &&
@@ -91,6 +94,40 @@ SPoint2 gmshVertex::reparamOnFace(const GFace *face, int dir) const
}
return SPoint2(U, V);
}
+ else if(s->Typ == MSH_SURF_TRIC){
+ // printf("coucou HERE I AM TRIC\n");
+ Curve *C[3];
+ for(int i = 0; i < 3; i++)
+ List_Read(s->Generatrices, i, &C[i]);
+
+ double U, V;
+ if ((C[0]->beg == v && C[2]->beg == v) ||
+ (C[0]->end == v && C[2]->beg == v) ||
+ (C[0]->beg == v && C[2]->end == v) ||
+ (C[0]->end == v && C[2]->end == v)){
+ U = V = 0;
+ }
+ else if ((C[0]->beg == v && C[1]->beg == v) ||
+ (C[0]->end == v && C[1]->beg == v) ||
+ (C[0]->beg == v && C[1]->end == v) ||
+ (C[0]->end == v && C[1]->end == v)){
+ U = 1;
+ V = 0;
+ }
+ else if ((C[2]->beg == v && C[1]->beg == v) ||
+ (C[2]->end == v && C[1]->beg == v) ||
+ (C[2]->beg == v && C[1]->end == v) ||
+ (C[2]->end == v && C[1]->end == v)){
+ U = 1;
+ V = 1;
+ }
+ else{
+ Msg::Info("Reparameterizing point %d on face %d", v->Num, s->Num);
+ return GVertex::reparamOnFace(face, dir);
+ }
+ // printf("coucou1 %g %g\n",U,V);
+ return SPoint2(U, V);
+ }
else{
return GVertex::reparamOnFace(face, dir);
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 605db6ccf9..4787d8071e 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -735,6 +735,17 @@ double backgroundMesh::operator() (double u, double v, double w) const
double backgroundMesh::getAngle(double u, double v, double w) const
{
+
+ // HACK FOR LEWIS
+ // h = 1+30(y-x^2)^2 + (1-x)^2
+ // double x = u;
+ // double y = v;
+ // double dhdx = 30 * 2 * (y-x*x) * (-2*x) - 2 * (1-x);
+ // double dhdy = 30 * 2 * (y-x*x);
+ // double angles = atan2(y,x*x);
+ // crossField2d::normalizeAngle (angles);
+ // return angles;
+
double uv[3] = {u, v, w};
double uv2[3];
MElement *e = _octree->find(u, v, w, 2, true);
@@ -776,6 +787,15 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
v1->x(),v1->y(),v1->z(),
v2->x(),v2->y(),v2->z(),
v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
+ /*
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ v1->x(),v1->y(),v1->z(),
+ v2->x(),v2->y(),v2->z(),
+ v3->x(),v3->y(),v3->z(),
+ getAngle(v1->x(),v1->y(),v1->z()),
+ getAngle(v2->x(),v2->y(),v2->z()),
+ getAngle(v3->x(),v3->y(),v3->z()));
+ */
}
else {
/*
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 50461ab573..c7e17705f0 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -41,7 +41,28 @@ static void myresid(int N, GEdge *ge, double *u, fullVector<double> &r)
for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
+static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N,
+ double *u, double underRelax)
+{
+ GPoint p0 = ge->point(u0);
+ GPoint p1 = ge->point(uN);
+ double du = 1. / (N - 1);
+ u[0] = u0;
+ // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
+ // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
+ for (int i = 1; i < N; i++){
+ SPoint3 pi (p0.x() + i * du * (p1.x()-p0.x()),
+ p0.y() + i * du * (p1.y()-p0.y()),
+ p0.z() + i * du * (p1.z()-p0.z()));
+ double t;
+ GPoint gp = ge->closestPoint(pi, t);
+ u[i] = gp.u();
+ // printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
+ }
+ return true;
+}
+
+static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N,
double *u, double underRelax)
{
const double PRECISION = 1.e-6;
@@ -100,6 +121,15 @@ static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
}
return false;
}
+// 1 = geodesics
+static int method_for_computing_intermediary_points = 0;
+static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
+ double *u, double underRelax){
+ if (method_for_computing_intermediary_points == 0) // use linear abscissa
+ return computeEquidistantParameters0(ge,u0,uN,N,u,underRelax);
+ else if (method_for_computing_intermediary_points == 1) // use projection
+ return computeEquidistantParameters1(ge,u0,uN,N,u,underRelax);
+}
static double mylength(GFace *gf, int i, double *u, double *v)
{
@@ -113,7 +143,31 @@ static void myresid(int N, GFace *gf, double *u, double *v, fullVector<double> &
for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
+static bool computeEquidistantParameters1(GFace *gf, double u0, double uN,
+ double v0, double vN, int N,
+ double *u, double *v){
+ // printf("coucou\n");
+ GPoint p0 = gf->point(u0,v0);
+ GPoint p1 = gf->point(uN,vN);
+ double du = 1. / (N - 1);
+ u[0] = u0;
+ u[0] = u0;
+ v[0] = v0;
+ // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
+ // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
+ for (int i = 1; i < N; i++){
+ SPoint3 pi (p0.x() + i * du * (p1.x()-p0.x()),
+ p0.y() + i * du * (p1.y()-p0.y()),
+ p0.z() + i * du * (p1.z()-p0.z()));
+ SPoint2 t;
+ GPoint gp = gf->closestPoint(pi, t);
+ u[i] = gp.u();
+ v[i] = gp.v();
+ }
+ return true;
+}
+
+static bool computeEquidistantParameters0(GFace *gf, double u0, double uN,
double v0, double vN, int N,
double *u, double *v)
{
@@ -191,6 +245,16 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
return false;
}
+static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
+ double v0, double vN, int N,
+ double *u, double *v){
+ if (method_for_computing_intermediary_points == 0) // use linear abscissa
+ return computeEquidistantParameters0(gf,u0,uN,v0,vN,N,u,v);
+ else if (method_for_computing_intermediary_points == 1) // use projection
+ return computeEquidistantParameters1(gf,u0,uN,v0,vN,N,u,v);
+}
+
+
static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
edgeContainer &edgeVertices, bool linear,
int nPts = 1, highOrderSmoother *displ2D = 0,
@@ -236,6 +300,9 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
"for edge %d-%d parametrized with %g %g on GEdge %d linear %d",
relax, US[1], v0->getNum(), v1->getNum(),u0,u1,ge->tag(), linear);
}
+ else{
+ Msg::Error("Cannot reparam a mesh Vertex in high order meshing");
+ }
}
for(int j = 0; j < nPts; j++){
const double t = (double)(j + 1)/(nPts + 1);
@@ -247,12 +314,13 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
v0->getNum(), v1->getNum());
// we don't have a (valid) parameter on the curve
SPoint3 pc = edge.interpolate(t);
- v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
+ v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
}
else {
- GPoint pc = ge->point(US[u0<u1? j + 1 : nPts - 1 - (j + 1)]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,
- US[u0 < u1 ? j + 1 : nPts - 1 - (j + 1)]);
+ int count = u0<u1? j + 1 : nPts + 1 - (j + 1);
+ GPoint pc = ge->point(US[count]);
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,pc.u());
+ // printf("Edge %d(%g) %d(%g) new vertex %g %g at %g\n",v0->getNum(),u0,v1->getNum(),u1,v->x(),v->y(), US[count]);
if (displ2D || displ3D){
SPoint3 pc2 = edge.interpolate(t);
if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 80d87b475c..088898a79b 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -397,6 +397,7 @@ void meshGEdge::operator() (GEdge *ge)
// FIXME JF : MAYBE WE SHOULD NOT ALWAYS SMOOTH THE 1D MESH SIZE FIELD ??
//printFandPrimitive(ge->tag(), Points);
// a = smoothPrimitive (ge,CTX::instance()->mesh.smoothRatio*CTX::instance()->mesh.smoothRatio,Points);
+ // printf(" coucou %g\n",a);
//printFandPrimitive(ge->tag()+10000, Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 83f1026d9d..5849d5c0f5 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -999,6 +999,70 @@ void optimalPointFrontal (GFace *gf,
}
+void optimalPointFrontalB (GFace *gf,
+ MTri3* worst,
+ int active_edge,
+ std::vector<double> &Us,
+ std::vector<double> &Vs,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM,
+ double newPoint[2],
+ double metric[3]){
+ double center[2],r2;
+ MTriangle *base = worst->tri();
+ circUV(base, Us, Vs, center, gf);
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
+ Us[base->getVertex(2)->getIndex()]) / 3.,
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
+ Vs[base->getVertex(2)->getIndex()]) / 3.};
+ buildMetric(gf, pa, metric);
+ circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+ // compute the middle point of the edge
+ int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
+ int ip2 = active_edge;
+
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+ Vs[base->getVertex(ip1)->getIndex()]};
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ Vs[base->getVertex(ip2)->getIndex()]};
+ double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
+
+ // now we have the edge center and the center of the circumcircle,
+ // we try to find a point that would produce a perfect triangle while
+ // connecting the 2 points of the active edge
+
+ double dir[2] = {center[0] - midpoint[0], center[1] - midpoint[1]};
+ double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
+ dir[0] /= norm;
+ dir[1] /= norm;
+ const double RATIO = sqrt(dir[0] * dir[0] * metric[0] +
+ 2 * dir[1] * dir[0] * metric[1] +
+ dir[1] * dir[1] * metric[2]);
+
+ const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
+ const double q = lengthMetric(center, midpoint, metric);
+ const double rhoM1 = 0.5 *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
+ vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
+ const double rhoM2 = 0.5 *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
+ const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
+
+ const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
+ const double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
+
+ // printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
+
+
+ newPoint[0] = midpoint[0] + d * dir[0];
+ newPoint[1] = midpoint[1] + d * dir[1];
+}
+
+
+
void bowyerWatsonFrontal(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e11c4d7c7e..59fae4b688 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -798,7 +798,6 @@ static int _removeDiamonds(GFace *gf)
touched.insert(gf->triangles[i]->getVertex(2));
}
-
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
MVertex *v1 = q->getVertex(0);
@@ -836,7 +835,7 @@ static int _removeDiamonds(GFace *gf)
deleted.insert(v1);
diamonds.insert(q);
}
- else if (v2->onWhat()->dim() == 2 &&
+ else if (0 && v2->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
v1->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
@@ -908,9 +907,11 @@ static int _removeDiamonds(GFace *gf)
mesh_vertices2.push_back(gf->mesh_vertices[i]);
}
else {
- delete gf->mesh_vertices[i];
+ // FIXME : GMSH SOMETIMES CRASHES IF DELETED ....
+ // delete gf->mesh_vertices[i];
}
}
+
gf->mesh_vertices = mesh_vertices2;
return diamonds.size();
@@ -1109,7 +1110,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
{
e2t_cont adj;
// buildEdgeToElement(gf->triangles, adj);
- buildEdgeToElement(gf->quadrangles, adj);
+ buildEdgeToElement(gf, adj);
std::vector<MQuadrangle*>created;
std::set<MElement*>deleted;
@@ -1196,7 +1197,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
}
static int edgeSwapQuadsForBetterQuality ( GFace *gf ) {
- // return 0;
+ // return 0;
int COUNT = 0;
while(1){
int k = _edgeSwapQuadsForBetterQuality (gf);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 8301a624e5..31ee8bfcd4 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -188,9 +188,7 @@ double mesh_functional_distorsion(MElement *t, double u, double v)
double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal1[3];
prodve(v1, v2, normal1);
- double nn = sqrt(normal1[0]*normal1[0] +
- normal1[1]*normal1[1] +
- normal1[2]*normal1[2]);
+ double nn = norm3(normal1);
// compute uncurved element jacobian d_u x and d_v x
@@ -200,14 +198,19 @@ double mesh_functional_distorsion(MElement *t, double u, double v)
double normal[3];
prodve(v1b, v2b, normal);
+ // printf("%g %g %g -- %g %g %g - %g\n",mat[0][0], mat[0][1], mat[0][2],mat[1][0], mat[1][1], mat[1][2],nn);
+
double sign = 1.0;
prosca(normal1, normal, &sign);
- double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
+ // double det = (norm3(normal) / nn) * (sign > 0 ? 1. : -1.);
// printf("%g %g : %g : %g n1 (%g,%g,%g)\n",u,v,sign,det, normal1[0], normal1[1], normal1[2]);
- // printf("n (%g,%g,%g)\n", normal[0], normal[1], normal[2]);
+ // for (int i=0;i<t->getNumVertices();i++){
+ // printf("COORD (%d) = %g %g %g\n",i,t->getVertex(i)->x(),t->getVertex(i)->y(),t->getVertex(i)->z());
+ // }
+ // printf("n (%g,%g,%g)\n", normal1[0], normal1[1], normal1[2]);
- return det;
+ return sign/(nn*nn);
}
static double MINQ (double a, double b, double c){
@@ -324,9 +327,9 @@ double mesh_functional_distorsion_pN(MElement *t)
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
jac->matrixLag2Bez.mult(Ji,Bi);
- // for (int i=0;i<jac->points.size1();i++){
- // printf("J(%d) = %12.5E\n",i,Bi(i));
- // }
+ // for (int i=0;i<jac->points.size1();i++){
+ // printf("J(%d) = %12.5E B(%d) = %12.5E\n",i,Ji(i),i,Bi(i));
+ // }
/*
if (t->getType() == TYPE_QUA){
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index eb9a2145ac..d9cf9cd6c6 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -380,7 +380,7 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
}
// Neumann conditions
GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
- std::cout << "Neumann BC"<< std::endl;
+
for (unsigned int i = 0; i < allNeumann.size(); i++)
{
LoadTerm<SVector3> Lterm(*LagSpace,*allNeumann[i]._f);
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 899dc9c68a..7b4f7e683e 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -57,7 +57,7 @@ struct neumannBC : public BoundaryCondition
// an elastic solver ...
class elasticitySolver
{
- protected:
+ public:
GModel *pModel;
int _dim, _tag;
dofManager<double> *pAssembler;
diff --git a/benchmarks/2d/recombine.geo b/benchmarks/2d/recombine.geo
index 7594af00aa..1ac40be9eb 100644
--- a/benchmarks/2d/recombine.geo
+++ b/benchmarks/2d/recombine.geo
@@ -9,5 +9,5 @@ Line(3) = {3,4} ;
Line(4) = {4,1} ;
Line Loop(5) = {4,1,-2,3} ;
Plane Surface(6) = {5} ;
-//Recombine Surface{6};
+Recombine Surface{6};
//Mesh.RecombinationAlgorithm=1;
diff --git a/benchmarks/2d/transfinite_b.geo b/benchmarks/2d/transfinite_b.geo
index 1a209a14a2..8efc8c4535 100644
--- a/benchmarks/2d/transfinite_b.geo
+++ b/benchmarks/2d/transfinite_b.geo
@@ -16,4 +16,4 @@ Transfinite Line {4} = 20 Using Progression 1./1.2;
Transfinite Line {3,2} = 10 Using Progression 1;
Transfinite Line {5} = 19 Using Progression 1;
Transfinite Surface {7} = {3,6,4,2} Alternated;
-//Recombine Surface {7};
+Recombine Surface {7};
--
GitLab