diff --git a/CMakeLists.txt b/CMakeLists.txt
index bd8fe3846246e11fc18d6cfac1a3343f02ffcb9a..5b94ce6910f4efb94ec05b03de844ad10805f83a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -58,6 +58,7 @@ option(ENABLE_PARSER "Build the GEO file parser" ${DEFAULT})
option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ${DEFAULT})
option(ENABLE_PLUGINS "Build the post-processing plugins" ${DEFAULT})
option(ENABLE_POST "Build the post-processing module" ${DEFAULT})
+option(ENABLE_POPPLER "Enable POPPLER package for displaying PDF documents" ON)
option(ENABLE_QT "Build proof-of-concept QT GUI" OFF)
option(ENABLE_RTREE "Enable RTREE" ${DEFAULT})
option(ENABLE_SALOME "Enable Salome routines for CAD healing" ${DEFAULT})
@@ -558,6 +559,17 @@ if(ENABLE_MPI)
endif(MPI_FOUND)
endif(ENABLE_MPI)
+if(ENABLE_POPPLER)
+ find_library(POPPLER_LIB poppler)
+ find_library(POPPLER_CPP_LIB poppler-cpp)
+ if(POPPLER_LIB)
+ set_config_option(HAVE_POPPLER "Poppler")
+ list(APPEND EXTERNAL_LIBRARIES ${POPPLER_LIB})
+ list(APPEND EXTERNAL_LIBRARIES ${POPPLER_CPP_LIB})
+ endif(POPPLER_LIB)
+endif(ENABLE_POPPLER)
+
+
if(HAVE_MESH OR HAVE_SOLVER)
if(ENABLE_METIS)
add_subdirectory(contrib/Metis)
diff --git a/Common/CMakeLists.txt b/Common/CMakeLists.txt
index efddc80f63afbe4621ab75a90f65ca3878ee07c3..889c8b550378457dabc2f767010a66f943b2aefc 100644
--- a/Common/CMakeLists.txt
+++ b/Common/CMakeLists.txt
@@ -3,10 +3,11 @@
# See the LICENSE.txt file for license information. Please report all
# bugs and problems to <gmsh@geuz.org>.
-set(SRC
+set(SRC
Gmsh.cpp
GmshRemote.cpp
GmshMessage.cpp
+ gmshPopplerWrapper.cpp
Context.cpp
Options.cpp
CommandLine.cpp
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 35201fb8967b53a453a1e50446c9e22fdfa2f769..588890a30df5f3a66427375261512ccf2216e322 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -616,8 +616,10 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.algo2d = ALGO_2D_MESHADAPT_OLD;
else if(!strncmp(argv[i], "del2d", 5) || !strncmp(argv[i], "tri", 3))
CTX::instance()->mesh.algo2d = ALGO_2D_DELAUNAY;
- else if(!strncmp(argv[i], "delquad", 5) || !strncmp(argv[i], "tri", 3))
+ else if(!strncmp(argv[i], "delquad", 5))
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL_QUAD;
+ else if(!strncmp(argv[i], "pack", 4))
+ CTX::instance()->mesh.algo2d = ALGO_2D_PACK_PRLGRMS;
else if(!strncmp(argv[i], "front2d", 7) || !strncmp(argv[i], "frontal", 7))
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL;
else if(!strncmp(argv[i], "bamg",4))
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 79b05ce051379ff6ea23c028cb55aa13c4fac558..b29c6dab91aee92580331f49a5a9dbd65b3fd2be 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -50,6 +50,7 @@
#cmakedefine HAVE_PETSC
#cmakedefine HAVE_PLUGINS
#cmakedefine HAVE_POST
+#cmakedefine HAVE_POPPLER
#cmakedefine HAVE_QT
#cmakedefine HAVE_RTREE
#cmakedefine HAVE_SALOME
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 89f1a4e77622aa1e2ec66651d81dd4596ee0f607..7512cf32d62fee5e5e862188e5c5039561a678ef 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -220,6 +220,7 @@
#define ALGO_2D_FRONTAL 6
#define ALGO_2D_BAMG 7
#define ALGO_2D_FRONTAL_QUAD 8
+#define ALGO_2D_PACK_PRLGRMS 9
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 12420c090698a17b99fe911d1eeaea4311fa9ebe..1ea67bf5278cd231dd423c9120042fbd2e8a83fa 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -20,6 +20,10 @@
#include "GeomMeshMatcher.h"
#include "Field.h"
+#if defined(HAVE_POPPLER)
+#include "gmshPopplerWrapper.h"
+#endif
+
#if defined(HAVE_PARSER)
#include "Parser.h"
#endif
@@ -319,6 +323,14 @@ int MergeFile(const std::string &fileName, bool warnIfMissing)
else if(ext == ".diff" || ext == ".DIFF"){
status = GModel::current()->readDIFF(fileName);
}
+ else if(ext == ".pdf" || ext == ".PDF"){
+#if defined(HAVE_POPPLER)
+ status = gmshPopplerWrapper::instance()->load_from_file(fileName);
+#else
+ Msg::Error("Gmsh has to be compiled with POPPLER for displaying PDF documents");
+ status = 0;
+#endif
+ }
else if(ext == ".med" || ext == ".MED" || ext == ".mmed" || ext == ".MMED" ||
ext == ".rmed" || ext == ".RMED"){
status = GModel::readMED(fileName);
diff --git a/Common/Options.cpp b/Common/Options.cpp
index a84b043e336c629c54f32197b5537d38b3ca673a..6140ea6f31c7fd5be208d872324cbbd1f3251d28 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5070,6 +5070,9 @@ double opt_mesh_algo2d(OPT_ARGS_NUM)
case ALGO_2D_FRONTAL_QUAD:
FlGui::instance()->options->mesh.choice[2]->value(4);
break;
+ case ALGO_2D_PACK_PRLGRMS:
+ FlGui::instance()->options->mesh.choice[2]->value(5);
+ break;
case ALGO_2D_AUTO:
default:
FlGui::instance()->options->mesh.choice[2]->value(0);
diff --git a/Common/gmshPopplerWrapper.cpp b/Common/gmshPopplerWrapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..862152ce06f91f8a4d549803b67742bd1b4dbed6
--- /dev/null
+++ b/Common/gmshPopplerWrapper.cpp
@@ -0,0 +1,63 @@
+#include "gmshPopplerWrapper.h"
+
+#if defined(HAVE_POPPLER)
+
+#include "GmshMessage.h"
+#include <poppler/cpp/poppler-image.h>
+#include <poppler/cpp/poppler-page-renderer.h>
+
+gmshPopplerWrapper *gmshPopplerWrapper::_instance = 0;
+poppler::document *gmshPopplerWrapper::_current_doc = 0;
+#if defined(HAVE_OPENGL)
+std::map<int,GLuint> gmshPopplerWrapper::_pages2textures;
+int gmshPopplerWrapper::_w = -1;
+int gmshPopplerWrapper::_h = -1;
+int gmshPopplerWrapper::_current_page = 0;
+#endif
+
+gmshPopplerWrapper *gmshPopplerWrapper::instance() {
+ if (!_instance)_instance = new gmshPopplerWrapper;
+ return _instance;
+}
+
+int gmshPopplerWrapper::load_from_file (const std::string &file_name,
+ const std::string &owner_password,
+ const std::string &user_password){
+ if (_current_doc) delete _current_doc;
+ _current_doc = poppler::document::load_from_file (file_name,owner_password,
+ user_password);
+ if (!_current_doc) return 0;
+ Msg::Info("PDF File has been loaded in the Wrapper");
+ // createBitmap(1,72.,72.,-1,-1,-1,-1);
+ return 1;
+}
+
+#if defined(HAVE_OPENGL)
+GLuint gmshPopplerWrapper::getTextureForPage(double xres,
+ double yres) {
+ int iPage = _current_page;
+ std::map<int,GLuint>::iterator it = _pages2textures.find(iPage);
+ if (it != _pages2textures.end())return it->second;
+ if (!_current_doc)return 0;
+ poppler::page *_current_page = _current_doc->create_page (iPage);
+ poppler::page_renderer pr;
+ poppler::image im = pr.render_page (_current_page,xres,yres,-1,-1,-1);
+ _w = im.width();
+ _h = im.height();
+ // im.save("page.png","png");
+ GLuint texture;
+ glGenTextures(1, &texture);
+ glBindTexture(GL_TEXTURE_2D, texture);
+ _pages2textures[iPage] = texture;
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, im.width(),im.height(), 0,
+ GL_RGBA, GL_UNSIGNED_BYTE,im.const_data());
+ return texture;
+}
+#endif
+
+
+#endif
+
diff --git a/Common/gmshPopplerWrapper.h b/Common/gmshPopplerWrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..1dff680737e5408d211bf2ed0d1caf4f3edeec5a
--- /dev/null
+++ b/Common/gmshPopplerWrapper.h
@@ -0,0 +1,42 @@
+#ifndef _GMSHPOPPLERWRAPPER_PDF_H_
+#define _GMSHPOPPLERWRAPPER_PDF_H_
+
+#include "GmshConfig.h"
+
+#if defined(HAVE_POPPLER)
+#include <string>
+#include <map>
+#include <poppler/cpp/poppler-document.h>
+#include <poppler/cpp/poppler-page.h>
+#include <poppler/cpp/poppler-image.h>
+#if defined(HAVE_OPENGL)
+#include<OpenGL/gl.h>
+#endif
+
+class gmshPopplerWrapper {
+private:
+ static int _current_page;
+ static poppler::document *_current_doc;
+ static gmshPopplerWrapper *_instance;
+ static int _w,_h;
+#if defined(HAVE_OPENGL)
+ static std::map<int,GLuint> _pages2textures; // map pages to textures
+#endif
+
+public:
+ static gmshPopplerWrapper *instance();
+ static int load_from_file (const std::string &file_name,
+ const std::string &owner_password=std::string(),
+ const std::string &user_password=std::string());
+ static int width() {return _w;}
+ static int height() {return _h;}
+ static void setCurrentPageUp () {_current_page++;}
+ static void setCurrentPageDown () {if(_current_page > 0) _current_page--;}
+#if defined(HAVE_OPENGL)
+ static GLuint getTextureForPage(double xres,
+ double yres) ;
+#endif
+};
+
+#endif
+#endif
diff --git a/Fltk/FlGui.cpp b/Fltk/FlGui.cpp
index 38e850a07bf1af56bb99b2977d5460ab43cdf2d7..f31347fa576ef3936532df399b82185153893193 100644
--- a/Fltk/FlGui.cpp
+++ b/Fltk/FlGui.cpp
@@ -48,6 +48,7 @@ typedef unsigned long intptr_t;
#include "StringUtils.h"
#include "Generator.h"
#include "gl2ps.h"
+#include "gmshPopplerWrapper.h"
#if defined(HAVE_3M)
#include "3M.h"
#endif
@@ -686,6 +687,18 @@ int FlGui::testArrowShortcuts()
status_play_manual(0, CTX::instance()->post.animStep);
return 1;
}
+ else if(Fl::test_shortcut(FL_ALT + FL_Up)) {
+ gmshPopplerWrapper::setCurrentPageDown();
+ drawContext::global()->draw();
+ return 1;
+ }
+ else if(Fl::test_shortcut(FL_ALT + FL_Down)) {
+ gmshPopplerWrapper::setCurrentPageUp();
+ drawContext::global()->draw();
+ return 1;
+ }
+
+
return 0;
}
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 92ea57fff7aadbf10ffd2f4353750c6eac0ed569..cb9546ae3557905c595580ebca19f521a4e4cc66 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -511,6 +511,7 @@ static void mesh_options_ok_cb(Fl_Widget *w, void *data)
(o->mesh.choice[2]->value() == 2) ? ALGO_2D_DELAUNAY :
(o->mesh.choice[2]->value() == 3) ? ALGO_2D_FRONTAL :
(o->mesh.choice[2]->value() == 4) ? ALGO_2D_FRONTAL_QUAD :
+ (o->mesh.choice[2]->value() == 5) ? ALGO_2D_PACK_PRLGRMS :
ALGO_2D_AUTO);
opt_mesh_algo3d(0, GMSH_SET,
(o->mesh.choice[3]->value() == 0) ? ALGO_3D_DELAUNAY :
@@ -2134,6 +2135,7 @@ optionWindow::optionWindow(int deltaFontSize)
{"Delaunay", 0, 0, 0},
{"Frontal", 0, 0, 0},
{"Delaunay for quads", 0, 0, 0},
+ {"Packing Of Parallelograms", 0, 0, 0},
{0}
};
static Fl_Menu_Item menu_3d_algo[] = {
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 5e9b7786a52cd6fe77512d0aa16c2ddc05d19b82..337a8768d4bb88a5c3948084312a006edd84f753 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -2094,6 +2094,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
if(!lt && _mapping != RBF){
//printf("POINT no success %d tris %d quad \n", triangles.size(), quadrangles.size());
GPoint gp = pointInRemeshedOctree(par1,par2);
+ gp.setNoSuccess();
return gp;
}
else if (!lt && _mapping == RBF){
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 77e3415daa12e9d18872c14e0efe7b62a8d91eda..969d7bb1d29214f618cfc95ca51c123ac5920842 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -636,6 +636,7 @@ class GModel
// CEA triangulation
int writeMAIL(const std::string &name, bool saveAll, double scalingFactor);
+
};
#endif
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 06e7c5a9f65e82463d98dc910008c0bd6a20b482..a948b4b9b169f1069caa2e8c02b4acf38f9a99eb 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -243,7 +243,10 @@ double MQuadrangle::etaShapeMeasure()
SVector3 c = crossprod(v23,v30);
SVector3 d = crossprod(v30,v01);
- if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )return 0.0;
+ double sign = 1.0;
+ if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )sign = -1;
+ // FIXME ...
+ // if (a.z() > 0 || b.z() > 0 || c.z() > 0 || d.z() > 0) sign = -1;
double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[3]) / M_PI;
@@ -259,9 +262,61 @@ double MQuadrangle::etaShapeMeasure()
angle = std::max(fabs(90. - a3),angle);
angle = std::max(fabs(90. - a4),angle);
- return 1.-angle/90;
+ return sign*(1.-angle/90);
}
+/// a shape measure for quadrangles
+/// assume (for now) 2D elements --
+/// sf = (1 \pm xi) (1 \pm eta) / 4
+/// dsf_xi = \pm (1 \pm eta) / 4
+/// 1 + eta , -(1+eta) , -(1-eta), 1-eta
+/// dsf_eta = \pm (1 \pm xi) / 4
+/// 1 + xi , 1 - xi , -(1-xi), -(1+xi)
+double MQuadrangle::gammaShapeMeasure(){
+ return etaShapeMeasure();
+ /*
+ // xi = -1 eta = -1
+ const double dsf_corner1 [2][4] = {{0,0,-.5,.5},{0,0.5,-.5,0}};
+ // xi = 1 eta = -1
+ const double dsf_corner2 [2][4] = {{0,0,-.5,.5},{0.5,0,0,-.5}};
+ // xi = 1 eta = 1
+ const double dsf_corner3 [2][4] = {{.5,-.5,0,0},{0.5,0,0,-.5}};
+ // xi = -1 eta = 1
+ const double dsf_corner4 [2][4] = {{0,0,-.5,.5},{0.5,0,0,-.5}};
+ */
+ double QT[4] = {qmTriangle(_v[0],_v[1],_v[2],QMTRI_RHO),
+ qmTriangle(_v[1],_v[2],_v[3],QMTRI_RHO),
+ qmTriangle(_v[2],_v[3],_v[0],QMTRI_RHO),
+ qmTriangle(_v[3],_v[0],_v[1],QMTRI_RHO)} ;
+ std::sort(QT,QT+4);
+ double quality = QT[0]*QT[1] / (QT[2] * QT[3]);
+
+ SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
+ SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
+ SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
+ SVector3 v30 (_v[0]->x()-_v[3]->x(),_v[0]->y()-_v[3]->y(),_v[0]->z()-_v[3]->z());
+
+ SVector3 a = crossprod(v01,v12);
+ SVector3 b = crossprod(v12,v23);
+ SVector3 c = crossprod(v23,v30);
+ SVector3 d = crossprod(v30,v01);
+
+ if (a.z() < 0 || b.z() < 0 || c.z() < 0 || d.z() < 0) return -quality;
+
+ if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )return -quality;
+
+ return quality;
+ /*
+ double J[3][3];
+ double detJ = getJacobian(-1,-1, J);
+ double C[2][2] = {{0,0}{0,0}};
+ for (int i=0;i<2;i++){
+
+ }*/
+
+}
+
+
double MQuadrangle::distoShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 3c1d9a7939e42128cf498a568a2f40028afdb39e..5eb78ee5049b14cc7a8fade4d920bbcdddc37a98 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -52,6 +52,7 @@ class MQuadrangle : public MElement {
}
~MQuadrangle(){}
virtual double etaShapeMeasure();
+ virtual double gammaShapeMeasure();
virtual int getDim() const { return 2; }
virtual int getNumVertices() const { return 4; }
virtual MVertex *getVertex(int num){ return _v[num]; }
diff --git a/Graphics/drawContext.cpp b/Graphics/drawContext.cpp
index 03592a7d90aa7826439596cbf1e0c67bd892c030..771499495028bd3753492b3b39d8a791b0c99591 100644
--- a/Graphics/drawContext.cpp
+++ b/Graphics/drawContext.cpp
@@ -25,6 +25,11 @@
#include <FL/gl.h>
#endif
+#if defined(HAVE_POPPLER)
+#include "gmshPopplerWrapper.h"
+#endif
+
+
drawContextGlobal *drawContext::_global = 0;
drawContext::drawContext(drawTransform *transform)
@@ -334,7 +339,7 @@ void drawContext::drawBackgroundGradient()
glVertex2i(viewport[0], viewport[1]);
glEnd();
}
- else if(CTX::instance()->bgGradient == 3){ // radial
+ else if(CTX::instance()->bgGradient == 4){ // radial
double cx = 0.5 * (viewport[0] + viewport[2]);
double cy = 0.5 * (viewport[1] + viewport[3]);
double r = 0.5 * std::max(viewport[2] - viewport[0],
@@ -351,6 +356,32 @@ void drawContext::drawBackgroundGradient()
}
glEnd();
}
+#if defined(HAVE_POPPLER)
+ else if(CTX::instance()->bgGradient == 3){ // PDF @ background
+ GLuint texture=gmshPopplerWrapper::getTextureForPage(800,600);
+ glEnable( GL_TEXTURE_2D );
+ glBindTexture(GL_TEXTURE_2D,texture);
+ glBegin(GL_QUADS);
+ glColor4ubv((GLubyte *) & CTX::instance()->color.bg);
+
+ int dw =viewport[2] - viewport[0];
+ int dh =viewport[3] - viewport[1];
+
+ int dw_im = gmshPopplerWrapper::width();
+ int dh_im = gmshPopplerWrapper::height();
+
+ // consterve aspect ratio : dw / dh = dw_im / dh_im
+ // dw = dh * dw_im / dh_im;
+ dh = dw * dh_im / dw_im;
+
+ glTexCoord2f(1.0f,1.0f);glVertex2i(viewport[2], viewport[3]-dh);
+ glTexCoord2f(1.0f,0.0f);glVertex2i(viewport[2], viewport[3]);
+ glTexCoord2f(0.0f,0.0f);glVertex2i(viewport[2]-dw, viewport[3]);
+ glTexCoord2f(0.0f,1.0f);glVertex2i(viewport[2]-dw, viewport[3]-dh);
+ glEnd();
+ }
+#endif
+
}
void drawContext::drawBackgroundImage()
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index c89e474075662de5496e2426a7dfbf745c5e1e5d..a1f6f2fa652ae534dd0c23a05827208fb0fc5407 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -906,6 +906,11 @@ void backgroundMesh::updateSizes(GFace *_gf)
}
}
+bool backgroundMesh::inDomain (double u, double v, double w) const
+{
+ return _octree->find(u, v, w, 2, true) != 0;
+}
+
double backgroundMesh::operator() (double u, double v, double w) const
{
double uv[3] = {u, v, w};
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 5e2495f87484c589edcc0c19d96dd56a19814d70..ad73401b2a37ee33e1f771b4e99852ff56148e8e 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -71,6 +71,7 @@ class backgroundMesh : public simpleFunction<double>
void propagateCrossFieldByDistance(GFace *);
void updateSizes(GFace *);
double operator () (double u, double v, double w) const; // returns mesh size
+ bool inDomain (double u, double v, double w) const; // returns true if in domain
double getAngle(double u, double v, double w) const ;
void print(const std::string &filename, GFace *gf,
const std::map<MVertex*, double>&) const;
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 8ed4c732db65dddaf160614272287774fbc639fd..4eca4a6a8034c9ddd2298f62dc933f0013c5adbc 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -38,6 +38,7 @@ set(SRC
yamakawa.cpp
Field.cpp
CenterlineField.cpp
+ surfaceFiller.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 8c9d998450c26142865f0bbe79e446e6b3362fdc..a7333333b68aa4d8041482933675a36f327b8346 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1483,7 +1483,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
}
- highOrderTools hot(m);
+ // highOrderTools hot(m);
// now we smooth mesh the internal vertices of the faces
// we do that model face by model face
@@ -1511,7 +1511,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
std::vector<MElement*> v;
v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- hot.applySmoothingTo(v, (*it));
+ //hot.applySmoothingTo(v, (*it));
// hot.applySmoothingTo(v, .1,0);
}
// hot.ensureMinimumDistorsion(0.1);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 373cbf9c341fc72f355625829bae570d53b9783e..fb3bc0d54fac90781e408473676b039df4157cf6 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -390,7 +390,8 @@ void meshGEdge::operator() (GEdge *ge)
// force odd number of points for if blossom is used for recombination
if(ge->meshAttributes.Method != MESH_TRANSFINITE &&
CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){
- if(CTX::instance()->mesh.recombineAll){
+
+ if(1 || CTX::instance()->mesh.recombineAll){
N++;
}
else{
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index da7f0bebc26f5bd205a0461367853d92e4ac9ce8..7c1987f84aca5713078ebab6e2ec93d583c1b4f7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -307,6 +307,7 @@ static bool algoDelaunay2D(GFace *gf)
gf->getMeshingAlgo() == ALGO_2D_BAMG ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
gf->getMeshingAlgo() == ALGO_2D_BAMG)
return true;
@@ -1121,11 +1122,13 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// fill the small gmsh structures
BDS2GMSH(m, gf, recoverMap);
- // BOUNDARY LAYER
+ bool infty = false;
+ if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD || gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
+ infty = true;
if (!onlyInitialMesh) {
- if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
+ if (infty)
buildBackGroundMesh (gf);
- //backgroundMesh::setCrossFieldsByDistance(gf);
+ // BOUNDARY LAYER
modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf);
}
@@ -1138,6 +1141,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD){
bowyerWatsonFrontalLayers(gf,true);
}
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS){
+ bowyerWatsonParallelograms(gf);
+ }
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf);
@@ -1145,7 +1151,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
bowyerWatson(gf,15000);
meshGFaceBamg(gf);
}
- laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
+ if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
+ laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
}
if(debug){
@@ -1884,6 +1891,7 @@ void meshGFace::operator() (GFace *gf, bool print)
case ALGO_2D_DELAUNAY : algo = "Delaunay"; break;
case ALGO_2D_MESHADAPT_OLD : algo = "MeshAdapt (old)"; break;
case ALGO_2D_BAMG : algo = "Bamg"; break;
+ case ALGO_2D_PACK_PRLGRMS : algo = "Square Packing"; break;
case ALGO_2D_AUTO :
algo = (gf->geomType() == GEntity::Plane) ? "Delaunay" : "MeshAdapt";
break;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 2fd6d76494dad109d40fdf2b3bafbc8709be604c..b1208cdcfcddc6544275169c24a00a7d8ccdb2f5 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -22,6 +22,7 @@
#include "GModel.h"
#include "GFaceCompound.h"
#include "intersectCurveSurface.h"
+#include "surfaceFiller.h"
double LIMIT_ = 0.5 * sqrt(2.0) * 1;
@@ -38,6 +39,51 @@ static bool isBoundary(MTri3 *t, double limit_, int &active)
return false;
}
*/
+template <class ITERATOR>
+void _printTris(char *name, ITERATOR it, ITERATOR end,
+ std::vector<double> &Us, std::vector<double> &Vs, bool param=true)
+{
+ FILE *ff = fopen (name,"w");
+ fprintf(ff,"View\"test\"{\n");
+ while ( it != end ){
+ MTri3 *worst = *it;
+ if (!worst->isDeleted()){
+ if (param)
+ fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ Us[(worst)->tri()->getVertex(0)->getIndex()],
+ Vs[(worst)->tri()->getVertex(0)->getIndex()],
+ 0.0,
+ Us[(worst)->tri()->getVertex(1)->getIndex()],
+ Vs[(worst)->tri()->getVertex(1)->getIndex()],
+ 0.0,
+ Us[(worst)->tri()->getVertex(2)->getIndex()],
+ Vs[(worst)->tri()->getVertex(2)->getIndex()],
+ 0.0,
+ (worst)->getRadius(),
+ (worst)->getRadius(),
+ (worst)->getRadius());
+ else
+ fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ (worst)->tri()->getVertex(0)->x(),
+ (worst)->tri()->getVertex(0)->y(),
+ (worst)->tri()->getVertex(0)->z(),
+ (worst)->tri()->getVertex(1)->x(),
+ (worst)->tri()->getVertex(1)->y(),
+ (worst)->tri()->getVertex(1)->z(),
+ (worst)->tri()->getVertex(2)->x(),
+ (worst)->tri()->getVertex(2)->y(),
+ (worst)->tri()->getVertex(2)->z(),
+ (worst)->getRadius(),
+ (worst)->getRadius(),
+ (worst)->getRadius()
+ );
+ }
+ ++it;
+ }
+ fprintf(ff,"};\n");
+ fclose (ff);
+}
+
static bool isActive(MTri3 *t, double limit_, int &active)
{
@@ -560,7 +606,7 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
// avoid angles that are too obtuse
double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
- if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -.9) && !force) {
+ if ((d1 < LL * .0025 || d2 < LL * .0025 || cosv < -.9999999) && !force) {
onePointIsTooClose = true;
}
@@ -577,7 +623,9 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
newVolume += ss;
++it;
}
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() >= 3 &&
+
+
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() > 3 &&
!onePointIsTooClose){
connectTris(new_cavity.begin(), new_cavity.end());
allTets.insert(newTris, newTris + shell.size());
@@ -596,68 +644,34 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
// The cavity is NOT star shaped
else{
- for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
- delete [] newTris;
+
+ // printf("(%g %g) %22.15E %22.15E %d %d %d %d %d\n",v->x(),v->y(),oldVolume, newVolume, newVolume,shell.size(), onePointIsTooClose, cavity.size(), new_cavity.size(),allTets.size());
+ // printf("12.5E 12.5E 12.5E 12.5E 12.5E\n",d1,d2,LL,cosv);
+
ittet = cavity.begin();
ittete = cavity.end();
while(ittet != ittete){
(*ittet)->setDeleted(false);
++ittet;
}
+ // _printTris("cavity.pos", cavity.begin(), cavity.end(), Us, Vs, false);
+ // _printTris("new_cavity.pos", new_cavity.begin(), new_cavity.end(), Us, Vs, false);
+ // _printTris("newTris.pos", &newTris[0], newTris+shell.size(), Us, Vs, false);
+ // _printTris("allTris.pos", allTets.begin(),allTets.end(), Us, Vs, false);
+ for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
+ delete [] newTris;
+ // throw;
return false;
}
}
-void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
- std::vector<double> &Us, std::vector<double> &Vs, bool param=true)
-{
- FILE *ff = fopen (name,"w");
- fprintf(ff,"View\"test\"{\n");
- std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
- while (it != AllTris.end() ){
- MTri3 *worst = *it;
- if (!worst->isDeleted()){
- if (param)
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- Us[(worst)->tri()->getVertex(0)->getIndex()],
- Vs[(worst)->tri()->getVertex(0)->getIndex()],
- 0.0,
- Us[(worst)->tri()->getVertex(1)->getIndex()],
- Vs[(worst)->tri()->getVertex(1)->getIndex()],
- 0.0,
- Us[(worst)->tri()->getVertex(2)->getIndex()],
- Vs[(worst)->tri()->getVertex(2)->getIndex()],
- 0.0,
- (worst)->getRadius(),
- (worst)->getRadius(),
- (worst)->getRadius());
- else
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- (worst)->tri()->getVertex(0)->x(),
- (worst)->tri()->getVertex(0)->y(),
- (worst)->tri()->getVertex(0)->z(),
- (worst)->tri()->getVertex(1)->x(),
- (worst)->tri()->getVertex(1)->y(),
- (worst)->tri()->getVertex(1)->z(),
- (worst)->tri()->getVertex(2)->x(),
- (worst)->tri()->getVertex(2)->y(),
- (worst)->tri()->getVertex(2)->z(),
- (worst)->getRadius(),
- (worst)->getRadius(),
- (worst)->getRadius()
- );
- }
- ++it;
- }
- fprintf(ff,"};\n");
- fclose (ff);
-}
static MTri3* search4Triangle (MTri3 *t, double pt[2],
std::vector<double> &Us, std::vector<double> &Vs,
std::set<MTri3*,compareTri3Ptr> &AllTris, double uv[2]) {
bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
+ // printf("inside = %d (%g %g)\n",inside,pt[0],pt[1]);
if (inside) return t;
SPoint3 q1(pt[0],pt[1],0);
int ITER = 0;
@@ -682,10 +696,11 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
t = t->getNeigh(i);
if (!t) break;
bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
+ // printf("ITER %d %d\n",ITER,inside);
if (inside) return t;
if (ITER++ > (int)AllTris.size()) break;
}
- return 0;
+ // return 0;
for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin();
itx != AllTris.end();++itx){
if (!(*itx)->isDeleted()){
@@ -720,6 +735,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MTri3 *ptin = 0;
bool inside = false;
double uv[2];
+ // printf("inserting %g %g\n",center[0],center[1]);
ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
if (ptin) inside = true;
@@ -750,11 +766,13 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
Vs.push_back(center[1]);
if(!p.succeeded() || !insertVertex
- (false, gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
+ (false, gf, v, center, ptin, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
- Msg::Debug("Point %g %g cannot be inserted because %d",
- center[0], center[1], p.succeeded() );
+ // Msg::Debug("Point %g %g cannot be inserted because %d",
+ // center[0], center[1], p.succeeded() );
+ // printf("Point %g %g cannot be inserted because %d",
+ // center[0], center[1], p.succeeded() );
AllTris.erase(it);
worst->forceRadius(-1);
@@ -763,12 +781,14 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
return false;
}
else {
- //printf("ouf ! %d\n",AllTris.size());
+ // printf("done ! %d\n",AllTris.size());
+
gf->mesh_vertices.push_back(v);
return true;
}
}
else {
+ printf("NOT INSIDE\n");
/*
MTriangle *base = worst->tri();
Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
@@ -1100,8 +1120,16 @@ void bowyerWatsonFrontal(GFace *gf)
}
// insert points
+ int ITERATION = 0;
while (1){
-
+ ++ITERATION;
+ if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
+ char name[245];
+ sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ _printTris (name, AllTris.begin(), AllTris.end(), Us,Vs,true);
+ sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ _printTris (name, ActiveTris.begin(), ActiveTris.end(), Us,Vs,true);
+ }
/* if(ITER % 100== 0){
char name[245];
sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
@@ -1457,10 +1485,10 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
ITERATION ++;
if(ITERATION % 1== 0 && CTX::instance()->mesh.saveAll){
char name[245];
- sprintf(name,"denInfinite_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
- _printTris (name, AllTris, Us,Vs,true);
- sprintf(name,"denInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- _printTris (name, ActiveTris, Us,Vs,true);
+ sprintf(name,"delInfinite_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
+ _printTris (name, AllTris.begin(), AllTris.end(), Us,Vs,true);
+ sprintf(name,"delInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ _printTris (name, ActiveTris.begin(), ActiveTris.end(),Us,Vs,true);
}
std::set<MTri3*,compareTri3Ptr> ActiveTrisNotInFront;
@@ -1546,3 +1574,49 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
LIMIT_ = 0.5 * sqrt(2.0) * 1;
backgroundMesh::unset();
}
+
+void bowyerWatsonParallelograms(GFace *gf)
+{
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ std::vector<SMetric3> vMetricsBGM;
+ std::vector<MVertex*> packed;
+
+ printf("creating the points\n");
+ packingOfParallelograms(gf, packed);
+ printf("points created\n");
+
+ buildMeshGenerationDataStructures
+ (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
+
+ // delaunise the initial mesh
+ int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
+ Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
+
+ // printf("CARNAVAL !!!\n");
+
+ // std::sort(packed.begin(), packed.end(), lexicographicSort());
+
+ for (unsigned int i=0;i<packed.size();){
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
+ else{
+ double newPoint[2] ;
+ packed[i]->getParameter(0,newPoint[0]);
+ packed[i]->getParameter(1,newPoint[1]);
+ delete packed[i];
+ double metric[3];
+ buildMetric(gf, newPoint, metric);
+ bool success = insertAPoint(gf, AllTris.begin(), newPoint, metric, Us, Vs, vSizes,
+ vSizesBGM, vMetricsBGM, AllTris, 0, worst);
+ if (!success)printf("success %d %d\n",success,AllTris.size());
+ i++;
+ }
+ }
+ transferDataStructure(gf, AllTris, Us, Vs);
+ backgroundMesh::unset();
+}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 0c6952fe1727f914dcf04ed85b2e22afecaa02da..ab9878792d4d1ac5667fcba8059f5b73410657f5 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -98,6 +98,7 @@ void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
void bowyerWatson(GFace *gf, int MAXPNT= 1000000000);
void bowyerWatsonFrontal(GFace *gf);
void bowyerWatsonFrontalLayers(GFace *gf, bool quad);
+void bowyerWatsonParallelograms(GFace *gf);
void buildBackGroundMesh (GFace *gf);
struct edgeXface
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a18ea8b5b133f856b0c8d0d324a88d8bf6ea091a..8f549127a59c89f4dc474ad29f925592aec25e6c 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -630,7 +630,7 @@ static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
bool badQuality = false;
for (unsigned int j=0;j<e1.size();++j){
surface_new += surfaceFaceUV(e1[j],gf);
- qualityNew[j] = (e1[j]->getNumVertices() == 4) ? e1[j]->etaShapeMeasure() : e1[j]->gammaShapeMeasure();
+ qualityNew[j] = e1[j]->gammaShapeMeasure();
if (qualityNew[j] < 0.01) {
badQuality = true;
break;
@@ -648,6 +648,50 @@ static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
}
+static bool _isItAGoodIdeaToMoveThoseVertices (GFace *gf,
+ const std::vector<MElement*> &e1,
+ std::vector<MVertex *>v,
+ const std::vector<SPoint2> &before,
+ const std::vector<SPoint2> &after)
+{
+ double surface_old = 0;
+ double surface_new = 0;
+
+ std::vector<SPoint3> pafter(v.size()), pbefore(v.size());
+ for (unsigned int i=0; i<v.size();++i){
+ GPoint gp = gf->point(after[i]);
+ if (!gp.succeeded())return false;
+ pafter[i] = SPoint3 (gp.x(),gp.y(),gp.z());
+ pbefore[i] = SPoint3(v[i]->x(),v[i]->y(),v[i]->z());
+ }
+
+ for (unsigned int j=0;j<e1.size();++j)
+ surface_old += surfaceFaceUV(e1[j],gf);
+
+ for (unsigned int i=0; i<v.size();++i){
+ v[i]->setParameter(0,after[i].x());
+ v[i]->setParameter(1,after[i].y());
+ v[i]->setXYZ(pafter[i].x(),pafter[i].y(),pafter[i].z());
+ }
+
+ bool badQuality = false;
+ for (unsigned int j=0;j<e1.size();++j){
+ surface_new += surfaceFaceUV(e1[j],gf);
+ }
+
+ for (unsigned int i=0; i<v.size();++i){
+ v[i]->setParameter(0,before[i].x());
+ v[i]->setParameter(1,before[i].y());
+ v[i]->setXYZ(pbefore[i].x(),pbefore[i].y(),pbefore[i].z());
+ }
+
+ if (badQuality || (surface_new - surface_old) > 1.e-10 * surface_old) {
+ return false;
+ }
+ return true;
+}
+
+
static int _quadWithOneVertexOnBoundary (GFace *gf,
std::set<MVertex*> &touched,
std::set<MElement*> &diamonds,
@@ -1515,7 +1559,7 @@ static int _splitFlatQuads(GFace *gf, double minQual)
if(it->second.size()==1 && it->first->onWhat()->dim() == 1) {
const std::vector<MElement*> < = it->second;
MElement *e = lt[0];
- if (e->getNumVertices() == 4 && e->etaShapeMeasure() < minQual){
+ if (e->getNumVertices() == 4 && e->gammaShapeMeasure() < minQual){
int k=0;
while(1){
if (e->getVertex(k) == it->first){
@@ -1749,51 +1793,131 @@ struct p1p2p3 {
#if defined(HAVE_BFGS)
// Callback function for BFGS
+static void sort_edges (std::vector<MEdge> &eds){
+
+ std::list<MEdge> eds_sorted;
+ eds_sorted.push_back(*eds.begin());
+ eds.erase(eds.begin());
+
+ while(!eds.empty()){
+ MEdge first = (eds_sorted.front());
+ MEdge last = (eds_sorted.back());
+ for (unsigned int i=0;i<eds.size();i++){
+ MVertex *v1 = eds[i].getVertex(0);
+ MVertex *v2 = eds[i].getVertex(1);
+ bool found = true;
+ // v1 -- v2 -- first
+ if (first.getVertex(0) == v2) eds_sorted.push_front(MEdge(v1,v2));
+ // v2 -- v1 -- first
+ else if (first.getVertex(0) == v1) eds_sorted.push_front(MEdge(v2,v1));
+ // last -- v1 -- v2
+ else if (last.getVertex(1) == v1) eds_sorted.push_back(MEdge(v1,v2));
+ // last -- v2 -- v1
+ else if (last.getVertex(1) == v2) eds_sorted.push_back(MEdge(v2,v1));
+ else found = false;
+ if (found) {
+ eds.erase(eds.begin() + i);
+ break;
+ }
+ }
+ }
+ eds.insert(eds.begin(),eds_sorted.begin(),eds_sorted.end());
+}
+
+static int OPTI_NUMBER = 1;
struct opti_data_vertex_relocation {
+ int nv;
const std::vector<MElement*> & e;
- MVertex *v;
+ MVertex *v[4];
+ bool movable[4];
GFace *gf;
- double minq,eps;
opti_data_vertex_relocation (const std::vector<MElement*> & _e,
MVertex * _v, GFace *_gf)
- : e(_e), v(_v), gf(_gf)
+ : nv(1),e(_e), gf(_gf)
+ {
+ v[0] = _v;
+ movable[0] = true;
+ }
+
+ opti_data_vertex_relocation (const std::vector<MElement*> & _e,
+ MVertex * _v1, MVertex * _v2, MVertex * _v3, MVertex * _v4, GFace *_gf)
+ : nv(4),e(_e), gf(_gf)
{
- minq = -2;
- eps = minq / (1. - minq);
+ v[0] = _v1;
+ v[1] = _v2;
+ v[2] = _v3;
+ v[3] = _v4;
+ for (int i=0;i<4;i++){
+ movable[i] = v[0]->onWhat()->dim() == 2;
+ }
+ }
+
+ void print_cavity (char *name){
+ FILE *f = fopen(name,"w");
+ fprintf(f,"View \"\"{\n");
+ for (unsigned int i=0;i<e.size();++i){
+ MElement *el = e[i];
+ el->writePOS(f,false,false,true,false,false,false);
+ }
+ fprintf(f,"};");
+ fclose (f);
}
+
+ /// quality measure for a quad
+
double f() const
{
+ /*
+ double minq = 1.0;
+ for (unsigned int i=0;i<e.size();++i){
+ MElement *el = e[i];
+ minq = std::min(el->gammaShapeMeasure(), minq);
+ }
+ if (minq < 0)minq *= 1.1;
+ else minq *= 0.;
+ */
double val = 0.0;
- for (unsigned int i = 0; i < e.size(); i++){
- const double q = (e[i]->getNumVertices() == 4) ? e[i]->etaShapeMeasure() :
- e[i]->gammaShapeMeasure();
- const double dd = (1 + eps) * q - eps;
- const double l = (dd>0) ? log (dd) : 1.e22;
- const double m = (q-1);
- val += l * l + m*m;
+ for (unsigned int i=0;i<e.size();++i){
+ MElement *el = e[i];
+ // double m = log((el->gammaShapeMeasure()-minq)/(1.-minq));
+ // val += m*m;//el->gammaShapeMeasure();
+ val += el->gammaShapeMeasure();
}
+
return val;
}
- void df(const double &F, const double &U, const double &V, double &dfdu, double &dfdv)
+ void df(const double &F, const double U[], const double V[], double dfdu[], double dfdv[])
{
const double eps = 1.e-6;
- set_(U+eps,V);
- const double FU = f();
- set_(U,V+eps);
- const double FV = f();
- set_(U,V);
- dfdu = -(F-FU)/eps;
- dfdv = -(F-FV)/eps;
- }
- void set_(double U, double V)
+ for (int i=0;i<nv;i++){
+ if (!set_(U[i]+eps,V[i],i)){
+ for (int j=0;j<nv;j++)dfdu[j] = dfdv[j] = 0;
+ return;
+ }
+ // printf("coucou2 %22.15E, %22.15E\n", U[i]+eps,V[i]);
+ const double FU = f();
+ set_(U[i],V[i]+eps,i);
+ // printf("coucou3\n");
+ const double FV = f();
+ set_(U[i],V[i],i);
+ // printf("coucou4\n");
+ // printf("%12.5E %12.5E %12.5E\n",F,FU,FV);
+ dfdu[i] = movable[i] ? -(F-FU)/eps : 0 ;
+ dfdv[i] = movable[i] ? -(F-FV)/eps : 0;
+ }
+ }
+ bool set_(double U, double V, int iv)
{
+ // if (!movable[iv])return;
+ // printf("getting point of surface %d (%22.15E,%22.15E)\n",gf->tag(),U,V);
GPoint gp = gf->point(U,V);
- if (!gp.succeeded()) Msg::Error("point not OK \n");
- v->x() = gp.x();
- v->y() = gp.y();
- v->z() = gp.z();
- v->setParameter(0,U);
- v->setParameter(1,V);
+ if (!gp.succeeded()) return false;//Msg::Error("point not OK");
+ v[iv]->x() = gp.x();
+ v[iv]->y() = gp.y();
+ v[iv]->z() = gp.z();
+ v[iv]->setParameter(0,U);
+ v[iv]->setParameter(1,V);
+ return true;
}
};
@@ -1801,21 +1925,112 @@ void bfgs_callback_vertex_relocation (const alglib::real_1d_array& x,
double& func, alglib::real_1d_array& grad, void* ptr)
{
opti_data_vertex_relocation* w = static_cast<opti_data_vertex_relocation*>(ptr);
- w->set_(x[0],x[1]);
+
+ double u[4] ; for (int i=0;i<w->nv;i++)u[i] = x[2*i];
+ double v[4] ; for (int i=0;i<w->nv;i++)v[i] = x[2*i+1];
+ // printf("hoplaboum !!!\n");
+ for (int i=0;i<w->nv;i++) w->set_(u[i],v[i],i);
func = w->f();
// printf("F(%3.2f,%3.2f) = %12.5E\n",x[0],x[1],func);
- w->df(func,x[0],x[1],grad[0],grad[1]);
+ double dfdu[4],dfdv[4];
+ w->df(func,u,v,dfdu,dfdv);
+ for (int i=0;i<w->nv;i++) {
+ grad[2*i] = dfdu[i];
+ grad[2*i+1] = dfdv[i];
+ }
+}
+
+
+// use optimization for untangling one single quad
+// take all neighboring quads and move all vertices
+// when possible
+static int _untangleQuad (GFace *gf, MQuadrangle *q,v2t_cont & adj)
+{
+ std::set<MElement*> all;
+ for (int i=0;i<4;i++){
+ std::vector<MElement*> &adji = adj[q->getVertex(i)];
+ all.insert(adji.begin(),adji.end());
+ // FIXME
+ if (q->getVertex(i)->onWhat()->dim() != 2)return 0;
+ }
+ std::vector<MElement*> lt;
+ lt.insert(lt.begin(),all.begin(),all.end());
+
+ double minq_old = 100.;
+ for (unsigned int i = 0; i < lt.size(); i++){
+ const double q = lt[i]->gammaShapeMeasure();
+ minq_old = std::min(q,minq_old);
+ }
+ // printf("-------x--------x------------x-------------------\n");
+ // printf ("optimizing quadrangle (minq %12.5E) with BFGS %3d\n",minq_old,OPTI_NUMBER);
+
+
+ alglib::ae_int_t dim = 8;
+ alglib::ae_int_t corr = 2;
+ alglib::minlbfgsstate state;
+ alglib::real_1d_array x;
+ std::vector<double> initial_conditions(8);
+ opti_data_vertex_relocation data(lt,q->getVertex(0),q->getVertex(1),q->getVertex(2),q->getVertex(3),gf);
+
+ // char NNN[256];
+ // sprintf(NNN,"UNTANGLE_cavity_%d_before.pos",OPTI_NUMBER);
+ // data.print_cavity (NNN);
+
+ double U[4],V[4];
+ for (int i=0;i<4;i++){
+ q->getVertex(i)->getParameter(0,U[i]);
+ q->getVertex(i)->getParameter(1,V[i]);
+ initial_conditions[2*i] = U[i];
+ initial_conditions[2*i+1] = V[i];
+ }
+ x.setcontent(dim, &initial_conditions[0]);
+ minlbfgscreate(8, corr, x, state);
+ double epsg = 0.0;
+ double epsf = 0.0;
+ double epsx = 0.0;
+ alglib::ae_int_t maxits = 12;
+ minlbfgssetcond(state,epsg,epsf,epsx,maxits);
+ minlbfgsoptimize(state, bfgs_callback_vertex_relocation,NULL,&data);
+ alglib::minlbfgsreport rep;
+ minlbfgsresults(state,x,rep);
+
+ double minq = 100.;
+ for (unsigned int i = 0; i < data.e.size(); i++){
+ const double q = data.e[i]->gammaShapeMeasure();
+ minq = std::min(q,minq);
+ }
+ // printf("minq = %12.5E\n",minq);
+
+ std::vector<SPoint2> before;for(int i=0;i<4;i++)before.push_back(SPoint2(U[i],V[i]));
+ std::vector<SPoint2> after;for(int i=0;i<4;i++)after.push_back(SPoint2(x[2*i],x[2*i+1]));
+ std::vector<MVertex*> vs;for(int i=0;i<4;i++)vs.push_back(q->getVertex(i));
+ bool success = _isItAGoodIdeaToMoveThoseVertices (gf,lt,vs,before,after);
+
+ if (!success || minq < 0.0 || minq <= minq_old/2) for (int i=0;i<4;i++)data.set_(U[i],V[i],i);
+ else return 1;
+ return 0;
+ // else printf("OKIDOKI-UNTANGLE\n");
+ // sprintf(NNN,"UNTANGLE_cavity_%d_after.pos",OPTI_NUMBER++);
+ // data.print_cavity (NNN);
}
+
// use optimization
static void _relocateVertexOpti(GFace *gf, MVertex *ver,
const std::vector<MElement*> <)
{
- // DOES NOT WORK AT ALL !!!
- return;
+ // return;
if(ver->onWhat()->dim() != 2)return;
- printf("heyhey\n");
- // printf ("optimizing vertex position with BFGS\n");
+
+ double minq_old = 100.;
+ // printf("------------------------------------------------\n");
+ for (unsigned int i = 0; i < lt.size(); i++){
+ const double q = lt[i]->gammaShapeMeasure();
+ minq_old = std::min(q,minq_old);
+ // printf("Q(%d) = %12.5E\n",lt[i]->getNumVertices(),q);
+ }
+
+ // printf ("optimizing vertex position (minq %12.5E) with BFGS\n",minq_old);
// we optimize the local coordinates of the node
alglib::ae_int_t dim = 2;
alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7]
@@ -1823,6 +2038,12 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
alglib::real_1d_array x; // = "[0.5,0.5]";
std::vector<double> initial_conditions(dim);
opti_data_vertex_relocation data(lt,ver,gf);
+
+ // char NNN[256];
+ // sprintf(NNN,"cavity_%d_before.pos",OPTI_NUMBER);
+ // data.print_cavity (NNN);
+
+
double U, V;
ver->getParameter(0,U);
ver->getParameter(1,V);
@@ -1830,7 +2051,7 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
initial_conditions[1] = V;
x.setcontent(dim, &initial_conditions[0]);
minlbfgscreate(2, corr, x, state);
- double epsg = 1.e-12;
+ double epsg = 0.0;
double epsf = 0.0;
double epsx = 0.0;
alglib::ae_int_t maxits = 10;
@@ -1839,17 +2060,24 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
alglib::minlbfgsreport rep;
minlbfgsresults(state,x,rep);
- double minq = 0;
+ double minq = 100.;
for (unsigned int i = 0; i < data.e.size(); i++){
- const double q = (data.e[i]->getNumVertices() == 4) ? data.e[i]->etaShapeMeasure() :
- data.e[i]->gammaShapeMeasure();
+ const double q = data.e[i]->gammaShapeMeasure();
minq = std::min(q,minq);
}
- if (minq < 0.01) data.set_(U,V);
- else printf("OKIDOKI\n");
+ // printf("minq = %12.5E\n",minq);
+
+ bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,SPoint2(U,V),SPoint2(x[0],x[1]));
+
+ if (!success || minq < 0 || minq <= minq_old/2) data.set_(U,V,0);
+ else {
+//printf("OKIDOKI\n");
+ }
// if (rep.terminationtype != 4){
// data.set_(U,V);
// }
+ // sprintf(NNN,"cavity_%d_after.pos",OPTI_NUMBER++);
+ // data.print_cavity (NNN);
}
#endif
@@ -1915,7 +2143,7 @@ void _relocateVertex(GFace *gf, MVertex *ver,
}
else{
#if defined(HAVE_BFGS)
- _relocateVertexOpti(gf, ver, lt);
+ _relocateVertexOpti(gf, ver, lt);
#endif
}
}
@@ -1931,7 +2159,7 @@ void quadBlob::smooth (int niter)
}
}
-void laplaceSmoothing(GFace *gf, int niter)
+void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
{
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
@@ -1945,6 +2173,25 @@ void laplaceSmoothing(GFace *gf, int niter)
}
}
+
+int untangleInvalidQuads(GFace *gf, int niter)
+{
+ // return;
+ int N = 0;
+ v2t_cont adj;
+ buildVertexToElement(gf->triangles, adj);
+ buildVertexToElement(gf->quadrangles, adj);
+ for(int i = 0; i < niter; i++){
+ for (unsigned int j=0;j<gf->quadrangles.size();j++){
+ if (gf->quadrangles[j]->gammaShapeMeasure() < 0.1){
+ N += _untangleQuad (gf, gf->quadrangles[j], adj);
+ }
+ }
+ }
+ return N;
+}
+
+
int _edgeSwapQuadsForBetterQuality(GFace *gf)
{
e2t_cont adj;
@@ -1964,7 +2211,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
MVertex *v2 = it->first.getVertex(1);
MElement *e1 = it->second.first;
MElement *e2 = it->second.second;
- double worst_quality_old = std::min(e1->etaShapeMeasure(),e2->etaShapeMeasure());
+ double worst_quality_old = std::min(e1->gammaShapeMeasure(),e2->gammaShapeMeasure());
if (worst_quality_old < .3 && (
deleted.find(e1) == deleted.end() ||
@@ -2006,8 +2253,8 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
q1B = new MQuadrangle (v11,v12,v21,v1);
q2B = new MQuadrangle (v21,v12,v2,v22);
}
- double worst_quality_A = std::min(q1A->etaShapeMeasure(),q2A->etaShapeMeasure());
- double worst_quality_B = std::min(q1B->etaShapeMeasure(),q2B->etaShapeMeasure());
+ double worst_quality_A = std::min(q1A->gammaShapeMeasure(),q2A->gammaShapeMeasure());
+ double worst_quality_B = std::min(q1B->gammaShapeMeasure(),q2B->gammaShapeMeasure());
bool ca1,ca2,cb1,cb2;
double old_surface = surfaceFaceUV(e1,gf) + surfaceFaceUV(e2,gf) ;
@@ -2723,6 +2970,9 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
if(CTX::instance()->mesh.algoRecombine == 1){
#if defined(HAVE_BLOSSOM)
int ncount = gf->triangles.size();
+ if (ncount % 2 != 0) {
+ printf("strange\n");
+ }
if (ncount % 2 == 0) {
int ecount = cubicGraph ? pairs.size() + makeGraphPeriodic.size() : pairs.size();
Msg::Info("Blossom: %d internal %d closed",(int)pairs.size(), (int)makeGraphPeriodic.size());
@@ -3007,18 +3257,20 @@ void recombineIntoQuads(GFace *gf,
}
double t4 = Cpu();
- laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
+ untangleInvalidQuads(gf,CTX::instance()->mesh.nbSmoothing);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
double t5 = Cpu();
optistatus[3] = edgeSwapQuadsForBetterQuality(gf);
if(optistatus[3] && saveAll){
sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME);
}
+ double t5b = Cpu();
tFQ += (t2-t1);
tTQ += (t3-t2);
tDI += (t4-t3);
tLA += (t5-t4);
- tSW += (t6-t5);
+ tSW += (t5b-t5);
tBU += (t7-t6);
if (!(optistatus[0]+optistatus[1]+optistatus[2]+optistatus[3]+optistatus[4]))
break;
@@ -3027,15 +3279,14 @@ void recombineIntoQuads(GFace *gf,
if (optistatus[i] != oldoptistatus[i])theSame = false;
oldoptistatus[i] = optistatus[i];
}
- if(theSame)break;
+ if(theSame)break;
if (ITER++ >= 20)break;
}
Msg::Debug("Opti times FQ(%4.3f) tQ(%4.3f) DI(%4.3f) LA(%4.3f) "
"SW(%4.3f) BUN(%4.3f)",tFQ,tTQ,tDI,tLA,tSW,tBU);
}
- if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
-
+ if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing, true);
}
double t2 = Cpu();
Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
@@ -3046,9 +3297,9 @@ void recombineIntoQuads(GFace *gf,
// simple recombination algo
_recombineIntoQuads(gf, 0);
- if(haveParam) laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
+ if(haveParam) laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, true);
_recombineIntoQuads(gf, 0);
- if(haveParam) laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
+ if(haveParam) laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, true);
if (saveAll) gf->model()->writeMSH("after.msh");
@@ -3061,7 +3312,7 @@ void quadsToTriangles(GFace *gf, double minqual)
std::vector<MQuadrangle*> qds;
for (unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
- if (q->etaShapeMeasure() < minqual){
+ if (q->gammaShapeMeasure() < minqual){
MTriangle *t11 = new MTriangle (q->getVertex(0),q->getVertex(1),q->getVertex(2));
MTriangle *t12 = new MTriangle (q->getVertex(2),q->getVertex(3),q->getVertex(0));
MTriangle *t21 = new MTriangle (q->getVertex(1),q->getVertex(2),q->getVertex(3));
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index b41a9a39b77e3518f74c31187e808a80c305aba9..6c8942a162c0f013a07cbe7744a4fa45bf7bbbc4 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -56,7 +56,7 @@ void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
std::vector<edge_angle> &edges_lonly);
void buildEdgeToElements(std::vector<MElement*> &tris, e2t_cont &adj);
-void laplaceSmoothing(GFace *gf, int niter=1);
+void laplaceSmoothing(GFace *gf, int niter=1, bool infinity_norm = false);
void _relocateVertex(GFace *gf, MVertex *ver,
const std::vector<MElement*> <);
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 0c732af43a357b4b52a64a6f1d8666fdfa7f38c2..56a3c2dd8cf21fa67ebf60b4ea06e54de34eae96 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -64,6 +64,9 @@ bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2,
ring.push_back(lastinring);
cavity.push_back(t);
+ // a change here for hybrid meshes
+
+ int ITER = 0;
while (1){
MVertex *ov1 = t->tet()->getVertex(edges[5 - iLocalEdge][0]);
MVertex *ov2 = t->tet()->getVertex(edges[5 - iLocalEdge][1]);
@@ -100,6 +103,9 @@ bool buildEdgeCavity(MTet4 *t, int iLocalEdge, MVertex **v1, MVertex **v2,
Msg::Error("loc = %d", iLocalEdge);
return false;
}
+ // FIXME when hybrid mesh, this loops for ever
+ // if (cavity.size() > 10)return false;
+ // printf("%d %d\n",ITER++, cavity.size());
}
computeNeighboringTetsOfACavity (cavity,outside);
return true;
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index fdef2b6416b9ffddcd61659ebca8240e2d16dbcb..38527c783dce169b3ef72a44ba01bceec0b30123 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -691,3 +691,5 @@ void Filler::print_node(Node* node,std::ofstream& file){
/*********static declarations*********/
std::vector<MVertex*> Filler::new_vertices;
+
+
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d143564e11a05d99d54d514fa64a02698e04be4
--- /dev/null
+++ b/Mesh/surfaceFiller.cpp
@@ -0,0 +1,342 @@
+#include "GmshConfig.h"
+#include "surfaceFiller.h"
+#include <queue>
+#include <stack>
+
+/// Here, we aim at producing a set of points that
+/// enables to generate a nice quad mesh
+
+#if defined(HAVE_RTREE)
+#include "rtree.h"
+#endif
+
+#include "MVertex.h"
+//#include "directions3D.h"
+#include "BackgroundMesh.h"
+#include "intersectCurveSurface.h"
+
+static const double FACTOR = .81;
+static const int NUMDIR = 3;
+static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
+
+/// a rectangle in the tangent plane is transformed
+/// into a parallelogram. We define an exclusion zone
+/// that is centered around a vertex and that is used
+/// in a r-tree structure for generating points with the
+/// right spacing in the tangent plane
+
+
+struct surfacePointWithExclusionRegion {
+ MVertex *_v;
+ SPoint2 _center;
+ SPoint2 _p[4][NUMDIR];
+ SPoint2 _q[4];
+ /*
+ + p3
+ p4 |
+ +----c-----+ p2
+ |
+ + p1
+
+*/
+
+ surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR]){
+ _v = v;
+ _center = (p[0][0]+p[1][0]+p[2][0]+p[3][0])*.25;
+ for (int i=0;i<4;i++)_q[i] = _center + (p[i][0]+p[(i+1)%4][0]-_center*2)*FACTOR;
+ for (int i=0;i<4;i++)for (int j=0;j<NUMDIR;j++)_p[i][j] = p[i][j];
+ }
+ bool inExclusionZone (const SPoint2 &p){
+ double mat[2][2];
+ double b[2] , uv[2];
+ mat[0][0]= _q[1].x()-_q[0].x();
+ mat[0][1]= _q[2].x()-_q[0].x();
+ mat[1][0]= _q[1].y()-_q[0].y();
+ mat[1][1]= _q[2].y()-_q[0].y();
+ b[0] = p.x() - _q[0].x();
+ b[1] = p.y() - _q[0].y();
+ sys2x2(mat, b, uv);
+ // printf("inversion 1 : %g %g \n",uv[0],uv[1]);
+ if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
+ mat[0][0]= _q[3].x()-_q[2].x();
+ mat[0][1]= _q[0].x()-_q[2].x();
+ mat[1][0]= _q[3].y()-_q[2].y();
+ mat[1][1]= _q[0].y()-_q[2].y();
+ b[0] = p.x() - _q[2].x();
+ b[1] = p.y() - _q[2].y();
+ sys2x2(mat, b, uv);
+ // printf("inversion 2 : %g %g \n",uv[0],uv[1]);
+ if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
+ return false;
+ }
+ void minmax (double _min[2], double _max[2]) const{
+ _min[0] = std::min(std::min(std::min(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
+ _min[1] = std::min(std::min(std::min(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
+ _max[0] = std::max(std::max(std::max(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
+ _max[1] = std::max(std::max(std::max(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
+ }
+};
+
+struct my_wrapper {
+ bool _tooclose;
+ SPoint2 _p;
+ my_wrapper (SPoint2 sp) : _tooclose (false), _p(sp) {}
+};
+
+
+bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
+ my_wrapper *w = static_cast<my_wrapper*>(point);
+
+ if (neighbour->inExclusionZone(w->_p)){
+ w->_tooclose = true;
+ return false;
+ }
+
+ return true;
+}
+
+bool inExclusionZone (SPoint2 &p,
+ RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
+ std::vector<surfacePointWithExclusionRegion*> & all ){
+ // should assert that the point is inside the domain
+ if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true;
+
+ my_wrapper w (p);
+ double _min[2] = {p.x()-1.e-8, p.y()-1.e-8},_max[2] = {p.x()+1.e-8,p.y()+1.e-8};
+ rtree.Search(_min,_max,rtree_callback,&w);
+
+ return w._tooclose;
+
+ for (unsigned int i=0;i<all.size();++i){
+ if (all[i]->inExclusionZone(p)){
+ // printf("%g %g is in exclusion zone of %g %g\n",p.x(),p.y(),all[i]._center.x(),all[i]._center.y());
+ return true;
+ }
+ }
+ return false;
+}
+
+
+
+
+// assume a point on the surface, compute the 4 possible neighbors.
+//
+// ^ t2
+// |
+// |
+// v2
+// |
+// |
+// v1-----+------v3 -------> t1
+// |
+// |
+// v4
+//
+// we aim at generating a rectangle with sizes size_1 and size_2 along t1 and t2
+
+bool compute4neighbors (GFace *gf, // the surface
+ MVertex *v_center, // the wertex for which we wnt to generate 4 neighbors
+ bool goNonLinear, // do we compute the position in the real surface which is nonlinear
+ SPoint2 newP[4][NUMDIR]) // look into other directions
+{
+ // we assume that v is on surface gf
+
+ // get the parameter of the point on the surface
+ SPoint2 midpoint;
+ reparamMeshVertexOnFace(v_center, gf, midpoint);
+
+ double size_1 = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0.0);
+ double size_2 = size_1;
+
+ // get the unit normal at that point
+ Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(midpoint[0],midpoint[1]));
+ SVector3 s1 = der.first();
+ SVector3 s2 = der.second();
+ SVector3 n = crossprod(s1,s2);
+ n.normalize();
+
+ for (int DIR = 0 ; DIR < NUMDIR ; DIR ++){
+ double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) + DIRS[DIR];
+
+ // normalize vector t1 that is tangent to gf at midpoint
+ SVector3 t1 = s1 * cos (quadAngle) + s2 * sin (quadAngle);
+ t1.normalize();
+
+ // compute the second direction t2 and normalize (t1,t2,n) is the tangent frame
+ SVector3 t2 = crossprod(t1,n);
+ t2.normalize();
+
+ // compute the first fundamental form i.e. the metric tensor at the point
+ // M_{ij} = s_i \cdot s_j
+ double M = dot(s1,s1);
+ double N = dot(s2,s2);
+ double E = dot(s1,s2);
+ double metric[2][2] = {{M,E},{E,N}};
+
+
+ // compute covariant coordinates of t1 and t2
+ // t1 = a s1 + b s2 -->
+ // t1 . s1 = a M + b E
+ // t1 . s2 = a E + b N --> solve the 2 x 2 system
+ // and get covariant coordinates a and b
+ double rhs1[2] = {dot(t1,s1),dot(t1,s2)}, covar1[2];
+ sys2x2(metric,rhs1,covar1);
+ double rhs2[2] = {dot(t2,s1),dot(t2,s2)}, covar2[2];
+ sys2x2(metric,rhs2,covar2);
+
+ // transform the sizes with respect to the metric
+ // consider a vector v of size 1 in the parameter plane
+ // its length is sqrt (v^T M v) --> if I want a real size
+ // of size1 in direction v, it should be sqrt(v^T M v) * size1
+ double l1 = sqrt(covar1[0]*covar1[0]+covar1[1]*covar1[1]);
+ double l2 = sqrt(covar2[0]*covar2[0]+covar2[1]*covar2[1]);
+ covar1[0] /= l1;covar1[1] /= l1;
+ covar2[0] /= l2;covar2[1] /= l2;
+
+
+ const double size_param_1 = size_1 / sqrt ( M*covar1[0]*covar1[0]+
+ 2*E*covar1[1]*covar1[0]+
+ N*covar1[1]*covar1[1]);
+ const double size_param_2 = size_2 / sqrt ( M*covar2[0]*covar2[0]+
+ 2*E*covar2[1]*covar2[0]+
+ N*covar2[1]*covar2[1]);
+
+ // printf("%12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());
+ // this is the rectangle in the parameter plane.
+ double r1 = 1.e-8*(double)rand() / RAND_MAX;
+ double r2 = 1.e-8*(double)rand() / RAND_MAX;
+ double r3 = 1.e-8*(double)rand() / RAND_MAX;
+ double r4 = 1.e-8*(double)rand() / RAND_MAX;
+ double r5 = 1.e-8*(double)rand() / RAND_MAX;
+ double r6 = 1.e-8*(double)rand() / RAND_MAX;
+ double r7 = 1.e-8* (double)rand() / RAND_MAX;
+ double r8 = 1.e-8*(double)rand() / RAND_MAX;
+ double newPoint[4][2] = {{midpoint[0] - covar1[0] * size_param_1 +r1,
+ midpoint[1] - covar1[1] * size_param_1 +r2},
+ {midpoint[0] - covar2[0] * size_param_2 +r3,
+ midpoint[1] - covar2[1] * size_param_2 +r4},
+ {midpoint[0] + covar1[0] * size_param_1 +r5,
+ midpoint[1] + covar1[1] * size_param_1 +r6},
+ {midpoint[0] + covar2[0] * size_param_2 +r7,
+ midpoint[1] + covar2[1] * size_param_2 +r8 }};
+ // We could stop here. Yet, if the metric varies a lot, we can solve
+ // a nonlinear problem in order to find a better approximation in the real
+ // surface
+ if (goNonLinear){//---------------------------------------------------//
+ surfaceFunctorGFace ss (gf); //
+ SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1,t2}; //
+ for (int i=0;i<4;i++){ //
+ double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0}; //
+ curveFunctorCircle cf (n,dirs[i],
+ SVector3(v_center->x(),v_center->y(),v_center->z()),
+ (i%2==1 ?size_param_1:size_param_2)*.5); //
+ if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-8)){ //
+ newPoint[i][0] = uvt[0]; //
+ newPoint[i][1] = uvt[1]; //
+ } //
+ } //
+ } /// end non linear -------------------------------------------------//
+
+ // return the four new vertices
+ for (int i=0;i<4;i++){
+ newP[i][DIR] = SPoint2(newPoint[i][0],newPoint[i][1]);
+ }
+ }
+ return true;
+}
+
+// fills a surface with points in order to build a nice
+// quad mesh ------------
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed){
+
+ // get all the boundary vertices
+ std::set<MVertex*> bnd_vertices;
+ for(unsigned int i=0;i<gf->getNumMeshElements();i++){
+ MElement* element = gf->getMeshElement(i);
+ for(int j=0;j<element->getNumVertices();j++){
+ MVertex *vertex = element->getVertex(j);
+ if (vertex->onWhat()->dim() < 2)bnd_vertices.insert(vertex);
+ }
+ }
+
+ // put boundary vertices in a fifo queue
+
+ // let us start without r-tree for now
+
+ std::queue<surfacePointWithExclusionRegion*> fifo;
+ std::vector<surfacePointWithExclusionRegion*> vertices;
+ // put the RTREE
+ RTree<surfacePointWithExclusionRegion*,double,2,double> rtree;
+
+ SVector3 t1;
+ SPoint2 newp[4][NUMDIR];
+ std::set<MVertex*>::iterator it = bnd_vertices.begin() ;
+ for (; it != bnd_vertices.end() ; ++it){
+ compute4neighbors (gf, *it, false, newp);
+ surfacePointWithExclusionRegion *sp =
+ new surfacePointWithExclusionRegion (*it, newp);
+ fifo.push(sp);
+ vertices.push_back(sp);
+ double _min[2],_max[2];
+ sp->minmax(_min,_max);
+ rtree.Insert(_min,_max,sp);
+ }
+
+ // printf("initially : %d vertices in the domain\n",vertices.size());
+
+
+ while(!fifo.empty()){
+ //surfacePointWithExclusionRegion & parent = fifo.top();
+ surfacePointWithExclusionRegion * parent = fifo.front();
+ fifo.pop();
+ for (int dir=0;dir<NUMDIR;dir++){
+ // printf("dir = %d\n",dir);
+ int countOK = 0;
+ for (int i=0;i<4;i++){
+ // printf("i = %d %12.5E %12.5E \n",i,parent._p[i][dir].x(),parent._p[i][dir].y());
+
+ // if (!w._tooclose){
+ if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){
+ countOK++;
+ GPoint gp = gf->point(parent->_p[i][dir]);
+ MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
+ // printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
+ compute4neighbors (gf, v, false, newp);
+ surfacePointWithExclusionRegion *sp =
+ new surfacePointWithExclusionRegion (v, newp);
+ fifo.push(sp);
+ vertices.push_back(sp);
+ double _min[2],_max[2];
+ sp->minmax(_min,_max);
+ rtree.Insert(_min,_max,sp);
+ }
+ }
+ if (countOK)break;
+ }
+ // printf("%d\n",vertices.size());
+ }
+
+ // printf("done\n");
+
+ // add the vertices as additional vertices in the
+ // surface mesh
+ FILE *f = fopen("points.pos","w");
+ fprintf(f,"View \"\"{\n");
+ for (int i=0;i<vertices.size();i++){
+ if(vertices[i]->_v->onWhat() == gf) {
+ packed.push_back(vertices[i]->_v);
+ SPoint2 midpoint;
+ reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
+ // fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",vertices[i]._v->x(),vertices[i]._v->y(),vertices[i]._v->z());
+ fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",midpoint.x(),midpoint.y(),0.0);
+ }
+ delete vertices[i];
+ }
+ fprintf(f,"};");
+ fclose(f);
+ // printf("packed.size = %d\n",packed.size());
+ // delete rtree;
+}
+
+
+
diff --git a/Mesh/surfaceFiller.h b/Mesh/surfaceFiller.h
new file mode 100644
index 0000000000000000000000000000000000000000..c5b215be935fcb336f9e6123ae6ef0a1c4247599
--- /dev/null
+++ b/Mesh/surfaceFiller.h
@@ -0,0 +1,10 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+
+#include <vector>
+class GFace;
+class MVertex;
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed);
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index ff5cc82aa5debc5b6f574eb4d8c9696889c3161f..5078de84475cc5a8289d4b45b8452c9c22c3e090 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -1,12 +1,13 @@
fact = 100;
lc = .1 * fact;
-Point(1) = {0.0,0.0,0,lc*.000001};
-Point(2) = {1* fact,0.0,0,lc};
-Point(3) = {1* fact,1* fact,0,lc*.00000001};
-Point(4) = {0,1* fact,0,lc};
+Point(1) = {0.0,0.0,0,lc*.3};
+Point(2) = {1* fact,0.0,0,lc*1};
+Point(3) = {1* fact,1* fact,0,lc};
+Point(4) = {0,1* fact,0,lc*.1};
Line(1) = {3,2};
Line(2) = {2,1};
Line(3) = {1,4};
Line(4) = {4,3};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
+Recombine Surface {6};
diff --git a/benchmarks/2d/Square-02.geo b/benchmarks/2d/Square-02.geo
index f781bd1dba09cdf7c8da3e94322b901a4a46667d..70c16d58ad5b2437d44faef176bdefb25e685eaa 100644
--- a/benchmarks/2d/Square-02.geo
+++ b/benchmarks/2d/Square-02.geo
@@ -3,8 +3,8 @@ Square non uniformly
meshed
******************************/
lc = .1;
-Point(1) = {0.0,0.0,0,lc*.01};
-Point(2) = {1,0.0,0,lc};
+Point(1) = {0.0,0.0,0,lc*.1};
+Point(2) = {1,0.0,0,lc*.1};
Point(3) = {1,1,0,lc};
Point(4) = {0,1,0,lc};
Line(1) = {3,2};
@@ -13,3 +13,4 @@ Line(3) = {1,4};
Line(4) = {4,3};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
+Recombine Surface {6};
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 05eb0a27cf98e14610763b2e8da6da63cd499566..21f9ae08d4a67f97ea461030a0d8116932e926d3 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -1,4 +1,4 @@
-unit = 1.0e-02 ;
+unit = 1.0e0 ;
e1 = 4.5 * unit ;
e2 = 6.0 * unit / 2.0 ;
@@ -14,8 +14,8 @@ r = 1.0 * unit ;
ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
ssin = ( 1.0 - ccos*ccos )^0.5 ;
-Lc1 = 0.003 ;
-Lc2 = 0.003 ;
+Lc1 = 0.3 * unit;
+Lc2 = 0.1 * unit;
Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};
Point(2) = { -e1-e2, h1 , 0.0 , Lc1};
diff --git a/benchmarks/3d/Torus2.geo b/benchmarks/3d/Torus2.geo
index 227ea82289155fc8d5e5914ce14299966c0d8900..667eb7e656a049e7b7995fc8011dbe8a4ee9e7ae 100644
--- a/benchmarks/3d/Torus2.geo
+++ b/benchmarks/3d/Torus2.geo
@@ -1,5 +1,5 @@
-Mesh.CharacteristicLengthFactor = 0.1;
+Mesh.CharacteristicLengthFactor = 0.05;
lc = 0.7;
Point(1) = {0,0,0,lc};
@@ -21,3 +21,4 @@ Surface Loop(72) = {45,23,67,71,49,27,15,59,37,41,19,63};
//Compound Surface(100)={45,23,67,71,49,27,15,59,37,41,19,63};
+Recombine Surface {50, 71, 67, 63, 59, 6, 45, 41, 23, 19, 37, 49, 28, 15, 27};