Commit a7f94810 by Amaury Johnen

Merge branch 'master' into validity

parents eda1c6e5 b3cbb0a7
Pipeline #173 passed with stage
in 8 minutes 48 seconds
......@@ -67,7 +67,7 @@ opt(NATIVE_FILE_CHOOSER "Enable native file chooser in GUI" ${DEFAULT})
opt(NETGEN "Enable Netgen 3D frontal mesh generator" ${DEFAULT})
opt(NUMPY "Enable conversion between fullMatrix and numpy array (requires SWIG)" OFF)
opt(PETSC4PY "Enable petsc4py wrappers for petsc matrices" ON)
opt(OCC "Enable Open CASCADE geometrical models" ${DEFAULT})
opt(OCC "Enable OpenCASCADE geometrical models" ${DEFAULT})
opt(ONELAB "Enable ONELAB solver interface" ${DEFAULT})
opt(ONELAB_METAMODEL "Enable ONELAB metamodels (experimental)" ${DEFAULT})
opt(OPENMP "Enable OpenMP (experimental)" OFF)
......@@ -107,6 +107,7 @@ set(GMSH_API
Common/VertexArray.h Common/Octree.h Common/OctreeInternals.h
Common/OS.h Common/StringUtils.h Common/OpenFile.h Common/Hash.h
Common/onelab.h Common/GmshSocket.h Common/onelabUtils.h Common/Options.h
Common/picojson.h
Numeric/Numeric.h Numeric/GaussIntegration.h Numeric/polynomialBasis.h
Numeric/JacobianBasis.h Numeric/bezierBasis.h Numeric/fullMatrix.h
Numeric/FuncSpaceData.h Numeric/pointsGenerators.h
......@@ -1089,6 +1090,7 @@ if(HAVE_SOLVER)
endif(HAVE_SOLVER)
if(ENABLE_OCC)
set(OCC_MINIMAL_VERSION "6.9.1")
if(WIN32 OR CYGWIN)
if(HAVE_64BIT_SIZE_T)
set(OCC_SYS_NAME win64)
......@@ -1098,53 +1100,72 @@ if(ENABLE_OCC)
else(WIN32 OR CYGWIN)
set(OCC_SYS_NAME ${CMAKE_SYSTEM_NAME})
endif(WIN32 OR CYGWIN)
set(OCC_LIBS_REQUIRED
# subset of DataExchange
TKSTEP TKSTEP209 TKSTEPAttr TKSTEPBase TKIGES TKXSBase
# ModelingAlgorithms
TKOffset TKFeat TKFillet TKBool TKShHealing TKMesh TKHLR TKBO TKPrim
TKTopAlgo TKGeomAlgo
# ModelingData
TKBRep TKGeomBase TKG3d TKG2d
# FoundationClasses
# TKAdvTools -- not necessary? (and removed from OCC 6.8)
TKMath TKernel)
list(LENGTH OCC_LIBS_REQUIRED NUM_OCC_LIBS_REQUIRED)
if(OCC_LIBS)
message(STATUS "OCC libraries specified explicitly: " ${OCC_LIBS})
list(LENGTH OCC_LIBS_REQUIRED NUM_OCC_LIBS)
else(OCC_LIBS)
set(OCC_LIBS)
foreach(OCC ${OCC_LIBS_REQUIRED})
find_library(OCC_LIB ${OCC} HINTS ENV CASROOT PATH_SUFFIXES lib
${OCC_SYS_NAME}/lib ${OCC_SYS_NAME}/vc8/lib)
if(OCC_LIB)
list(APPEND OCC_LIBS ${OCC_LIB})
else(OCC_LIB)
message(STATUS "OCC lib " ${OCC} " not Found")
endif(OCC_LIB)
unset(OCC_LIB CACHE)
endforeach(OCC)
list(LENGTH OCC_LIBS NUM_OCC_LIBS)
endif(OCC_LIBS)
if(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
find_path(OCC_INC "BRep_Tool.hxx" PATHS ENV CASROOT PATH_SUFFIXES inc
include include/oce opencascade)
if(OCC_INC)
set_config_option(HAVE_OCC "OpenCascade")
find_path(OCC_INC "Standard_Version.hxx" HINTS ENV CASROOT PATH_SUFFIXES inc
include include/oce opencascade include/opencascade)
if(OCC_INC)
file(STRINGS ${OCC_INC}/Standard_Version.hxx
OCC_MAJOR REGEX "#define OCC_VERSION_MAJOR.*")
file(STRINGS ${OCC_INC}/Standard_Version.hxx
OCC_MINOR REGEX "#define OCC_VERSION_MINOR.*")
file(STRINGS ${OCC_INC}/Standard_Version.hxx
OCC_MAINT REGEX "#define OCC_VERSION_MAINTENANCE.*")
if(OCC_MAJOR AND OCC_MINOR AND OCC_MAINT)
string(REGEX MATCH "[0-9]+" OCC_MAJOR "${OCC_MAJOR}")
string(REGEX MATCH "[0-9]+" OCC_MINOR "${OCC_MINOR}")
string(REGEX MATCH "[0-9]+" OCC_MAINT "${OCC_MAINT}")
set(OCC_VERSION "${OCC_MAJOR}.${OCC_MINOR}.${OCC_MAINT}")
message(STATUS "Found OpenCASCADE version ${OCC_VERSION} in ${OCC_INC}")
endif(OCC_MAJOR AND OCC_MINOR AND OCC_MAINT)
endif(OCC_INC)
if(OCC_VERSION AND OCC_VERSION STRLESS ${OCC_MINIMAL_VERSION})
message(WARNING "Gmsh requires OpenCASCADE >= ${OCC_MINIMAL_VERSION}")
message(WARNING "Use CMAKE_PREFIX_PATH or the CASROOT environment variable "
"to explicitely specify the installation path of OpenCASCADE")
else(OCC_VERSION AND OCC_VERSION STRLESS ${OCC_MINIMAL_VERSION})
set(OCC_LIBS_REQUIRED
# subset of DataExchange
TKSTEP TKSTEP209 TKSTEPAttr TKSTEPBase TKIGES TKXSBase
# ModelingAlgorithms
TKOffset TKFeat TKFillet TKBool TKShHealing TKMesh TKHLR TKBO TKPrim
TKTopAlgo TKGeomAlgo
# ModelingData
TKBRep TKGeomBase TKG3d TKG2d
# FoundationClasses
# TKAdvTools -- not necessary? (and removed from OCC 6.8)
TKMath TKernel)
list(LENGTH OCC_LIBS_REQUIRED NUM_OCC_LIBS_REQUIRED)
if(OCC_LIBS)
message(STATUS "OCC libraries specified explicitly: " ${OCC_LIBS})
list(LENGTH OCC_LIBS_REQUIRED NUM_OCC_LIBS)
else(OCC_LIBS)
set(OCC_LIBS)
foreach(OCC ${OCC_LIBS_REQUIRED})
find_library(OCC_LIB ${OCC} HINTS ENV CASROOT PATH_SUFFIXES
lib ${OCC_SYS_NAME}/lib ${OCC_SYS_NAME}/vc8/lib)
if(OCC_LIB)
list(APPEND OCC_LIBS ${OCC_LIB})
else(OCC_LIB)
message(STATUS "OCC lib " ${OCC} " not Found")
endif(OCC_LIB)
unset(OCC_LIB CACHE)
endforeach(OCC)
list(LENGTH OCC_LIBS NUM_OCC_LIBS)
endif(OCC_LIBS)
if(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
set_config_option(HAVE_OCC "OpenCASCADE")
list(APPEND EXTERNAL_LIBRARIES ${OCC_LIBS})
list(APPEND EXTERNAL_INCLUDES ${OCC_INC})
if(HAVE_64BIT_SIZE_T)
add_definitions(-D_OCC64)
endif(HAVE_64BIT_SIZE_T)
if(CYGWIN)
list(APPEND EXTERNAL_LIBRARIES "winspool")
list(APPEND EXTERNAL_LIBRARIES "winspool")
add_definitions(-DOCC_CONVERT_SIGNALS)
elseif(MSVC)
add_definitions(-DWNT)
add_definitions(-DWNT)
endif(CYGWIN)
endif(OCC_INC)
endif(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
endif(NUM_OCC_LIBS EQUAL NUM_OCC_LIBS_REQUIRED)
endif(OCC_VERSION AND OCC_VERSION STRLESS ${OCC_MINIMAL_VERSION})
endif(ENABLE_OCC)
if(ENABLE_ACIS)
......@@ -1455,9 +1476,12 @@ if(ENABLE_BUILD_ANDROID)
set(LIBRARY_OUTPUT_PATH_ROOT ${CMAKE_CURRENT_BINARY_DIR})
set(LIBRARY_OUTPUT_PATH ${CMAKE_CURRENT_BINARY_DIR}/libs/)
add_definitions(-DBUILD_ANDROID)
add_definitions(-DPICOJSON_USE_LOCALE=0)
add_library(androidGmsh SHARED ${GMSH_SRC})
set_target_properties(androidGmsh PROPERTIES OUTPUT_NAME Gmsh)
target_link_libraries(androidGmsh ${EXTERNAL_LIBRARIES} ${LAPACK_LIBRARIES})
add_custom_command(TARGET androidGmsh POST_BUILD COMMAND
${CMAKE_STRIP} ${LIBRARY_OUTPUT_PATH}/libGmsh.so)
endif(ENABLE_BUILD_ANDROID)
# shared library target
......
......@@ -913,7 +913,7 @@ StringXNumber GeometryOptions_Number[] = {
"Display geometry surfaces?" },
{ F|O, "SurfaceNumbers" , opt_geometry_surfaces_num , 0. ,
"Display surface numbers?" },
{ F|O, "SurfaceType" , opt_geometry_surface_type , 2. ,
{ F|O, "SurfaceType" , opt_geometry_surface_type , 0. ,
"Surface display type (0=cross, 1=wireframe, 2=solid)" },
{ F|O, "Tangents" , opt_geometry_tangents , 0. ,
......
......@@ -413,6 +413,10 @@ int MergeFile(const std::string &fileName, bool warnIfMissing, bool setBoundingB
else if(ext == ".dat" || ext == ".DAT"){
if(!strncmp(header, "BEGIN ACTRAN", 12))
status = GModel::current()->readACTRAN(fileName);
else if(!strncmp(header, "!", 1) ||
!strncmp(header, ";ECHO", 5) ||
!strncmp(header, ".NOE", 4))
status = GModel::current()->readSAMCEF(fileName);
else
status = GModel::current()->readBDF(fileName);
}
......
......@@ -9,10 +9,11 @@
#include "avl.h"
#include "ListUtils.h"
typedef struct {
class Tree_T {
public:
int size;
avl_tree *root;
} Tree_T;
};
Tree_T *Tree_Create(int size, int (*fcmp)(const void *a, const void *b));
void Tree_Delete(Tree_T *Tree);
......
......@@ -56,10 +56,10 @@ static void draw_stl(std::vector<SPoint3> &vertices, std::vector<SVector3> &norm
{
GLint mode[2];
glGetIntegerv(GL_POLYGON_MODE, mode);
if(CTX::instance()->geom.surfaceType > 1)
glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
else
if(CTX::instance()->geom.surfaceType == 1)
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
else
glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
glEnable(GL_LIGHTING);
glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
glColor4ubv((GLubyte *) & CTX::instance()->color.geom.highlight[0]);
......
......@@ -2187,6 +2187,13 @@ static void mesh_smooth_cb(Fl_Widget *w, void *data)
drawContext::global()->draw();
}
static void mesh_recombine_cb(Fl_Widget *w, void *data)
{
RecombineMesh(GModel::current());
drawContext::global()->draw();
}
static void mesh_optimize_netgen_cb(Fl_Widget *w, void *data)
{
if(CTX::instance()->lock) {
......@@ -4284,6 +4291,8 @@ static menuItem static_modules[] = {
#endif
{"0Modules/Mesh/Smooth 2D",
(Fl_Callback *)mesh_smooth_cb} ,
{"0Modules/Mesh/Recombine 2D",
(Fl_Callback *)mesh_recombine_cb} ,
{"0Modules/Mesh/Reclassify 2D",
(Fl_Callback *)mesh_classify_cb} ,
#if defined(HAVE_FOURIER_MODEL)
......
......@@ -31,7 +31,7 @@ set(SRC
GModelIO_PLY.cpp GModelIO_VRML.cpp GModelIO_UNV.cpp GModelIO_BDF.cpp
GModelIO_IR3.cpp GModelIO_DIFF.cpp GModelIO_GEOM.cpp GModelIO_INP.cpp
GModelIO_MAIL.cpp GModelIO_P3D.cpp GModelIO_CELUM.cpp
GModelIO_ACTRAN.cpp GModelIO_SU2.cpp
GModelIO_ACTRAN.cpp GModelIO_SU2.cpp GModelIO_SAMCEF.cpp
ExtrudeParams.cpp
Geo.cpp
GeoStringInterface.cpp GeoInterpolation.cpp
......
......@@ -669,6 +669,9 @@ class GModel {
// Actran mesh
int readACTRAN(const std::string &name);
// Sameced mesh
int readSAMCEF(const std::string &name);
// Plot3D structured mesh format
int readP3D(const std::string &name);
int writeP3D(const std::string &name,
......
......@@ -6,15 +6,15 @@
#ifndef _GMODELIO_GEO_H_
#define _GMODELIO_GEO_H_
#include "ListUtils.h"
#include "TreeUtils.h"
class GModel;
class ExtrudeParams;
class gmshSurface;
class List_T;
class Tree_T;
class GEO_Internals{
public:
// this will become private
Tree_T *Points, *Curves, *EdgeLoops, *Surfaces, *SurfaceLoops, *Volumes;
Tree_T *DelPoints, *DelCurves, *DelSurfaces, *DelVolumes;
List_T *PhysicalGroups, *DelPhysicalGroups;
......
......@@ -905,8 +905,8 @@ bool OCC_Internals::_addSpline(int &tag, const std::vector<int> &vertexTags, int
Msg::Error("OpenCASCADE edge with tag %d already exists", tag);
return false;
}
if(vertexTags.size() < 2 || vertexTags.size() > 20){
Msg::Error("Number of control points should be in [2,20]");
if(vertexTags.size() < 2){
Msg::Error("Number of control points should be at least 2");
return false;
}
......
// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
#include <stdlib.h>
#include <string.h>
#include <sstream>
#include "GModel.h"
#include "OS.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "Context.h"
static bool getVertices(GModel *m, int num, int *n, std::vector<MVertex*> &vec)
{
for(int i = 0; i < num; i++){
MVertex *v = m->getMeshVertexByTag(n[i]);
if(!v){
Msg::Error("Wrong vertex number %d", n[i]);
return false;
}
else
vec.push_back(v);
}
return true;
}
int GModel::readSAMCEF(const std::string &name)
{
FILE *fp = Fopen(name.c_str(), "r");
if(!fp){
Msg::Error("Unable to open file '%s'", name.c_str());
return 0;
}
_vertexMapCache.clear();
std::map<int, std::vector<MElement*> > elements[2];
char buffer[256], dummy[256];
while(!feof(fp)) {
if(!fgets(buffer, 256, fp)) break;
if(!strncmp(buffer, ".NOE", 4)){
while(!feof(fp)) {
if(!fgets(buffer, 256, fp)) break;
if(buffer[0] == '!') continue;
if(buffer[0] == ';') break;
int num;
double x, y, z;
if(sscanf(buffer, "%s %d %s %lf %s %lf %s %lf", dummy, &num,
dummy, &x, dummy, &y, dummy, &z) != 8)
return 0;
_vertexMapCache[num] = new MVertex(x, y, z, 0, num);
}
Msg::Info("Read %d mesh vertices", (int)_vertexMapCache.size());
}
else if(!strncmp(buffer, ".MAI", 4)){
while(!feof(fp)) {
if(!fgets(buffer, 256, fp)) break;
if(buffer[0] == '!') continue;
if(buffer[0] == ';') break;
std::stringstream s(buffer);
std::string word;
int nn = 0;
while(s >> word) nn++;
int num, reg, n[4];
std::vector<MVertex*> vertices;
if(nn == 8){ // TRIA3
if(sscanf(buffer, "%s %d %s %d %s %d %d %d", dummy, &num,
dummy, &reg, dummy, &n[0], &n[1], &n[2]) != 8)
return 0;
if(!getVertices(this, 3, n, vertices)) break;
elements[0][reg].push_back(new MTriangle(vertices, num));
}
else if(nn == 9){ // QUAD4
if(sscanf(buffer, "%s %d %s %d %s %d %d %d %d", dummy, &num,
dummy, &reg, dummy, &n[0], &n[1], &n[2], &n[3]) != 9)
return 0;
if(!getVertices(this, 4, n, vertices)) break;
elements[1][reg].push_back(new MQuadrangle(vertices, num));
}
}
}
}
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(_vertexMapCache);
fclose(fp);
return 1;
}
......@@ -181,7 +181,7 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
bool haveCells = true;
bool haveLines = false;
if( !strcmp(buffer, "CELLS") && numElements > 0)
if(!strcmp(buffer, "CELLS") && numElements > 0)
Msg::Info("Reading %d cells", numElements);
else if (!strcmp(buffer, "POLYGONS") && numElements > 0)
Msg::Info("Reading %d polygons", numElements);
......@@ -249,11 +249,13 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
case 13: elements[6][1].push_back(new MPrism(cells[i])); break;
case 14: elements[7][1].push_back(new MPyramid(cells[i])); break;
// second order elements
case 21: elements[1][1].push_back(new MLine(cells[i])); break;
case 22: elements[2][1].push_back(new MTriangle(cells[i])); break;
case 23: elements[3][1].push_back(new MQuadrangle(cells[i])); break;
case 24: elements[4][1].push_back(new MTetrahedron(cells[i])); break;
case 25: elements[5][1].push_back(new MHexahedron(cells[i])); break;
case 21: elements[1][1].push_back(new MLine3(cells[i])); break;
case 22: elements[2][1].push_back(new MTriangle6(cells[i])); break;
case 23: elements[3][1].push_back(new MQuadrangle8(cells[i])); break;
case 28: elements[3][1].push_back(new MQuadrangle9(cells[i])); break;
case 24: elements[4][1].push_back(new MTetrahedron10(cells[i])); break;
case 25: elements[5][1].push_back(new MHexahedron20(cells[i])); break;
case 29: elements[5][1].push_back(new MHexahedron27(cells[i])); break;
default:
Msg::Error("Unknown type of cell %d", type);
break;
......@@ -297,7 +299,7 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
}
}
else{
Msg::Error("TODO: implement reading lines for binary files \n");
Msg::Error("Line import not done for binary VTK files");
}
}
......
......@@ -27,6 +27,25 @@ gmshFace::gmshFace(GModel *m, Surface *face)
resetMeshAttributes();
}
// a face is degenerate if
bool gmshFace::degenerate(int dim) const {
std::list<GEdge*> eds = edges();
int numNonDegenerate = 0;
std::set<GEdge*> t;
for(std::list<GEdge*>::iterator it = eds.begin(); it != eds.end(); ++it){
GEdge *e = *it;
GVertex *start = e->getBeginVertex();
GVertex *next = e->getEndVertex();
if (start != next && t.find(e) == t.end()){
numNonDegenerate++;
}
t.insert(e);
}
// printf("%d \n",numNonDegenerate);
return numNonDegenerate <= 1;
}
void gmshFace::resetNativePtr(Surface *face)
{
s = face;
......
......@@ -36,6 +36,7 @@ class gmshFace : public GFace {
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
virtual void resetMeshAttributes();
void resetNativePtr(Surface *_s);
bool degenerate(int dim) const;
};
#endif
......@@ -982,7 +982,7 @@ void RecombineMesh(GModel *m)
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GFace *gf = *it;
recombineIntoQuads(gf);
recombineIntoQuads(gf,true,true,.01);
}
CTX::instance()->mesh.changed = ENT_ALL;
......
......@@ -20,13 +20,17 @@
#endif
typedef unsigned char CHECKTYPE;
struct Tet;
struct Vert {
private :
double _x[3];
double _lc;
unsigned int _num;
Tet *_t;
public :
inline void setT(Tet*t) { _t = t;}
inline Tet* getT() const { return _t;}
inline unsigned int getNum () const {return _num;}
inline void setNum (unsigned int n) {_num=n;}
unsigned char _thread;
......@@ -40,7 +44,7 @@ public :
inline double &lc() {return _lc;}
inline operator double *() { return _x; }
Vert (double X=0, double Y=0, double Z=0, double lc=0, int num = 0)
: _num(num), _thread(0)
: _num(num),_t(NULL), _thread(0)
{
_x[0] = X; _x[1] = Y; _x[2] = Z; _lc = lc;
}
......@@ -223,6 +227,7 @@ struct Tet {
{
_modified=true;
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
......@@ -232,6 +237,7 @@ struct Tet {
double val = robustPredicates::orient3d((double*)v0, (double*)v1,
(double*)v2, (double*)v3);
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
if (val > 0){
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
......@@ -397,6 +403,7 @@ void delaunayTrgl(const unsigned int numThreads,
unsigned int Npts,
std::vector<Vert*> assignTo[],
tetContainer &allocator, edgeContainer *embedded = 0);
bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread);
void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embeddedEdges);
void vertexRelocationPass (int numThreads, std::vector<Vert*> &v);
#endif
......@@ -142,20 +142,10 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
bgm = fields->get(fields->getBackgroundField());
}
if (0){
const double length = p1.distance(p2);
const double lc_avg = (e.first->lc() + e.second->lc())*.5;
const double l_adim_approx = length / lc_avg;
if (l_adim_approx > 10.){
SPoint3 pmid = (p1+p2)*0.5;
const double lc = GMSHSIZE(pmid, bgm, lc_avg);
S.push_back(new Vert(pmid.x(),pmid.y(),pmid.z(),lc));
return;
}
}
// double _a = Cpu();
// printf("%g %g \n",e.first->lc(), e.second->lc());
const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(),
_result, dl, temp, bgm);
// _C1 += Cpu() - _a;
......@@ -449,15 +439,25 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){
updateSize ((*it)->triangles[i]->getVertex((j+1)%3), s, d);
}
}
}
}
// there may be steiner points
for (unsigned int i=0;i<gr->tetrahedra.size();i++){
for (unsigned int j = 0; j < 4; j++){
MVertex *v = gr->tetrahedra[i]->getVertex(j);
if (s.find(v) == s.end()){
s[v] = 1.0;
}
}
}
}
extern double tetQuality (Tet* t, double *volume);
void edgeBasedRefinement(const int numThreads,
const int nptsatonce,
GRegion *gr)
{
// fill up old Datastructures
edgeContainer embeddedEdges (10000);
std::map<MVertex*, double> sizes;
......@@ -525,7 +525,6 @@ void edgeBasedRefinement(const int numThreads,
}
}
// do not allow to saturate boundary edges
{
for (unsigned int i=0;i< allocator.size(0);i++) {
......@@ -550,7 +549,6 @@ void edgeBasedRefinement(const int numThreads,
for (unsigned int i=0;i<_vertices.size();i++){
_filter.insert( _vertices[i] );
}
int iter = 1;
Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
......@@ -561,10 +559,10 @@ void edgeBasedRefinement(const int numThreads,
double t1 = Cpu();
// C_COUNT = 0;
saturateEdges (ec, allocator, numThreads, add);
// printf("%d calls %d points added",C_COUNT,add.size());
// printf("%d points added",add.size());
double t2 = Cpu();
filterVertices (numThreads, _filter, add);
// printf(" %d points remain\n",add.size());
// printf(" %d points remain\n",add.size());
double t3 = Cpu();
if (add.empty())break;
// randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
......@@ -572,6 +570,7 @@ void edgeBasedRefinement(const int numThreads,
// sort them using BRIO
std::vector<int> indices;
SortHilbert(add, indices);
// printf("%d indices\n",indices.size());
double t4 = Cpu();
delaunayTrgl (1,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
double t5 = Cpu();
......@@ -588,12 +587,45 @@ void edgeBasedRefinement(const int numThreads,
}
}
std::vector<Vert*> vv;
for (int myThread=0; myThread < numThreads; myThread++)
for (unsigned int i=0;i<allocator.size(myThread);i++)
for (unsigned int j=0;j<4;j++)
if (allocator(myThread,i)->V[j])
allocator(myThread,i)->V[j]->setNum(0);
for (int myThread=0; myThread < numThreads; myThread++)
for (unsigned int i=0;i<allocator.size(myThread);i++)
for (unsigned int j=0;j<4;j++)
if (allocator(myThread,i)->V[j] && allocator(myThread,i)->V[j]->getNum() == 0){
allocator(myThread,i)->V[j]->setNum(1);
std::map<Vert*,MVertex*>::iterator it = _ma.find(allocator(myThread,i)->V[j]);
if (it == _ma.end())vv.push_back(allocator(myThread,i)->V[j]);
}
double t6 = Cpu();
Msg::Info("Optimizing");
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
edgeSwapPass (numThreads, allocator, embeddedEdges);
vertexRelocationPass (numThreads,vv);
double t7 = Cpu();
Msg::Info("Optimization done (%g seconds)",t7-t6);
int cat [10] = {0,0,0,0,0,0,0,0,0,0};
double MIN = 1.0;
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
MVertex *mvs[4];
if (tt->V[0]){
double vol;
double q = tetQuality (tt,&vol);
cat [(int)(q*10)] ++;
MIN = std::min(MIN,q);
for (int j=0;j<4;j++){
Vert *v = tt->V[j];
std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
......@@ -609,6 +641,11 @@ void edgeBasedRefinement(const int numThreads,
}
}
Msg::Info ("Min Tet Quality %22.15E",MIN);
for (int i=0;i<10;i++){
Msg::Info ("Tet Quality [%4.3f,%4.3f] %8d",0.1*i,0.1*(i+1),cat[i]);
}
if (Msg::GetVerbosity() == 99) {
std::map<Edge,double> _sizes;
for (unsigned int i=0; i< allocator.size(0);i++){
......@@ -638,4 +675,6 @@ void edgeBasedRefinement(const int numThreads,
Msg::Info("MESH EFFICIENCY : %22.15E %6d edges among %d are out of range",
exp (sum / _sizes.size()),nbBad,_sizes.size());
}
}
......@@ -1006,6 +1006,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if (isMeshed) return true;
}
// if (gf->degenerate(0))return 0;
// build a set with all points of the boundaries
std::set<MVertex*, MVertexLessThanNum> all_vertices, boundary;
std::list<GEdge*>::iterator ite = edges.begin();
......
......@@ -510,8 +510,8 @@ static int _removeTwoQuadsNodes(GFace *gf)
return 0;
}
MQuadrangle *q = new MQuadrangle(v1,v2,v3,v4);
double s1 = surfaceFaceUV(q,gf);
double s2 = surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);;
double s1 = 0;//surfaceFaceUV(q,gf);
double s2 = 1;//surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);;
if (s1 > s2){
delete q;
}
......@@ -564,6 +564,66 @@ int removeTwoQuadsNodes(GFace *gf)
return nbRemove;
}
static bool _tryToCollapseThatVertex2 (GFace *gf,
std::vector<MElement*> &e1,
std::vector<MElement*> &e2,
MElement *q,
MVertex *v1,
MVertex *v2)
{
std::vector<MElement*> e = e1;
e.insert(e.end(), e2.begin(), e2.end());
double x1 = v1->x();
double y1 = v1->y();
double z1 = v1->z();
double x2 = v2->x();
double y2 = v2->y();
double z2 = v2->z();
// new position of v1 && v2
double initialGuess[2]={0,0};
GPoint pp = gf->closestPoint(SPoint3(0.5*(x1+x2),0.5*(y1+y2),0.5*(z1+z2)),initialGuess);
// double surface_old = 0;
// double surface_new = 0;
double worst_quality_old = 1.0;
double worst_quality_new = 1.0;
// surface_old = surfaceFaceUV(q,gf,false);