Commit e0abd072 authored by Christophe Geuzaine's avatar Christophe Geuzaine

added option to specify random perturbation factor for initial 3D delaunay + cleanup

parent 44f4bc8c
Pipeline #1197 passed with stage
in 35 minutes and 32 seconds
......@@ -18,7 +18,7 @@ class GamePad;
struct contextMeshOptions {
// mesh algorithms
int optimize, optimizeNetgen, smoothCrossField, refineSteps;
double lcFactor, randFactor, lcIntegrationPrecision;
double lcFactor, randFactor, randFactor3d, lcIntegrationPrecision;
double optimizeThreshold, normals, tangents, explode, angleSmoothNormals;
double allowSwapEdgeAngle;
double qualityInf, qualitySup, radiusInf, radiusSup;
......
......@@ -1181,6 +1181,8 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "RandomFactor" , opt_mesh_rand_factor , 1.e-9 ,
"Random factor used in the 2D meshing algorithm (should be increased if "
"RandomFactor * size(triangle)/size(model) approaches machine accuracy)" },
{ F|O, "RandomFactor3D" , opt_mesh_rand_factor3d , 1.e-12 ,
"Random factor used in the 3D meshing algorithm" },
{ F|O, "PreserveNumberingMsh2" , opt_mesh_preserve_numbering_msh2 , 0. ,
"Preserve element numbering in MSH2 format (will break meshes with multiple "
"physical groups for a single elementary entity)"},
......
......@@ -5191,6 +5191,16 @@ double opt_mesh_rand_factor(OPT_ARGS_NUM)
return CTX::instance()->mesh.randFactor;
}
double opt_mesh_rand_factor3d(OPT_ARGS_NUM)
{
if(action & GMSH_SET){
if(!(action & GMSH_SET_DEFAULT) && val != CTX::instance()->mesh.randFactor3d)
Msg::SetOnelabChanged(2);
CTX::instance()->mesh.randFactor3d = val;
}
return CTX::instance()->mesh.randFactor3d;
}
double opt_mesh_quality_type(OPT_ARGS_NUM)
{
if(action & GMSH_SET) {
......
......@@ -433,6 +433,7 @@ double opt_mesh_lc_from_points(OPT_ARGS_NUM);
double opt_mesh_lc_extend_from_boundary(OPT_ARGS_NUM);
double opt_mesh_lc_integration_precision(OPT_ARGS_NUM);
double opt_mesh_rand_factor(OPT_ARGS_NUM);
double opt_mesh_rand_factor3d(OPT_ARGS_NUM);
double opt_mesh_quality_inf(OPT_ARGS_NUM);
double opt_mesh_quality_sup(OPT_ARGS_NUM);
double opt_mesh_quality_type(OPT_ARGS_NUM);
......
This diff is collapsed.
This diff is collapsed.
......@@ -404,7 +404,8 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
mid->lcBGM() = BGM_MeshSize (gf,U,V, mid->X,mid->Y,mid->Z);
mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
// if (isSphere && !canWeSplitAnEdgeOnASphere (e, mid, center,radius))m.del_point(mid);
// if(isSphere && !canWeSplitAnEdgeOnASphere(e, mid, center,radius))
// m.del_point(mid);
if(!m.split_edge(e, mid)) m.del_point(mid);
else nb_split++;
}
......@@ -441,7 +442,8 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point
return maxLc;
}
static bool revertTriangleSphere (SPoint3 &center, BDS_Point *p, BDS_Point *o){
static bool revertTriangleSphere(SPoint3 &center, BDS_Point *p, BDS_Point *o)
{
std::list<BDS_Face*> t;
p->getTriangles(t);
std::list<BDS_Face*>::iterator it = t.begin();
......@@ -468,7 +470,6 @@ static bool revertTriangleSphere (SPoint3 &center, BDS_Point *p, BDS_Point *o){
return false;
}
void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
......@@ -532,7 +533,6 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
}
}
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
// FIXME SUPER HACK
......@@ -684,7 +684,8 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
Msg::Debug(" ---------------------------------------");
}
void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
void invalidEdgesPeriodic(BDS_Mesh &m,
std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
std::set<BDS_Edge*, EdgeLessThan> &toSplit)
{
// first look for degenerated vertices
......@@ -801,9 +802,12 @@ void TRYTOFIXSPHERES(GFace *gf, BDS_Mesh &m,
BDS_Point *n[4];
f->getNodes(n);
MVertex *v1 = (recoverMap->find(n[0])==recoverMap->end()) ? NULL : (*recoverMap)[n[0]];
MVertex *v2 = (recoverMap->find(n[1])==recoverMap->end()) ? NULL : (*recoverMap)[n[1]];
MVertex *v3 = (recoverMap->find(n[2])==recoverMap->end()) ? NULL : (*recoverMap)[n[2]];
MVertex *v1 = (recoverMap->find(n[0])==recoverMap->end()) ? NULL
: (*recoverMap)[n[0]];
MVertex *v2 = (recoverMap->find(n[1])==recoverMap->end()) ? NULL
: (*recoverMap)[n[1]];
MVertex *v3 = (recoverMap->find(n[2])==recoverMap->end()) ? NULL
: (*recoverMap)[n[2]];
if ((!v1 || (v1 != v2 && v1 != v3)) && (!v2 || v2 != v3)){
normal_triangle(n[0], n[1], n[2], norm);
......
......@@ -61,7 +61,6 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
J.mult(M,W);
W.mult(JT,R);
bamg::Metric M1(R(0,0),R(1,0),R(1,1));
// printf("%12.5E %12.5E %12.5E vs %12.5E %12.5E %12.5E \n",m(0,0),m(1,0),m(1,1),M1.a11,M1.a21,M1.a22);
mm11[i] = M1.a11;
mm12[i] = M1.a21;
mm22[i] = M1.a22;
......@@ -202,10 +201,12 @@ void meshGFaceBamg(GFace *gf)
Vertex2 &v = refinedBamgMesh->vertices[i];
if (i >= nbFixedVertices){
GPoint gp = gf->point(SPoint2(v[0], v[1]));
// if (gp.x() > 2.){ printf("wrong vertex index=%d %g %g %g (%g %g)\n", i, gp.x(), gp.y(), gp.z(), v[0], v[1]);
// if (gp.x() > 2.){
// printf("wrong vertex index=%d %g %g %g (%g %g)\n",
// i, gp.x(), gp.y(), gp.z(), v[0], v[1]);
// }
//If point not found because compound edges have been remeshed and boundary triangles have changed
//then we call our new octree
// If point not found because compound edges have been remeshed and
//boundary triangles have changed then we call our new octree
MFaceVertex *x = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, v[0], v[1]);
yetAnother[i] = x;
gf->mesh_vertices.push_back(x);
......
This diff is collapsed.
......@@ -38,7 +38,6 @@ struct bidimMeshData
if (it != parametricCoordinates->end()){
u = it->second.x();
v = it->second.y();
// printf("%g %g\n",u,v);
}
}
Us.push_back(u);
......@@ -46,11 +45,13 @@ struct bidimMeshData
vSizes.push_back(size);
vSizesBGM.push_back(sizeBGM);
}
inline int getIndex (MVertex *mv) {
inline int getIndex (MVertex *mv)
{
if (mv->onWhat()->dim() == 2)return mv->getIndex();
return indices[mv];
}
inline MVertex * equivalent (MVertex *v1) const {
inline MVertex * equivalent (MVertex *v1) const
{
if (equivalence){
std::map<MVertex* , MVertex*>::iterator it = equivalence->find(v1);
if (it == equivalence->end())return 0;
......@@ -58,12 +59,13 @@ struct bidimMeshData
}
return 0;
}
bidimMeshData (std::map<MVertex* , MVertex*>* e = 0, std::map<MVertex*, SPoint2> *p = 0) : equivalence(e), parametricCoordinates(p)
bidimMeshData(std::map<MVertex* , MVertex*>* e = 0,
std::map<MVertex*, SPoint2> *p = 0)
: equivalence(e), parametricCoordinates(p)
{
}
};
void buildMetric(GFace *gf, double *uv, double *metric);
int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
double *p4, double *metric);
......@@ -73,7 +75,8 @@ void circumCenterMetric(double *pa, double *pb, double *pc, const double *metric
double *x, double &Radius2);
void circumCenterMetric(MTriangle *base, const double *metric, bidimMeshData & data,
double *x, double &Radius2);
bool circumCenterMetricInTriangle(MTriangle *base, const double *metric, bidimMeshData &data);
bool circumCenterMetricInTriangle(MTriangle *base, const double *metric,
bidimMeshData &data);
bool invMapUV(MTriangle *t, double *p, bidimMeshData &data,
double *uv, double tol);
......@@ -96,10 +99,12 @@ class MTri3
MVertex *v1 = base->getVertex((i+2)%3);
MVertex *v2 = base->getVertex(i);
for (int j=0;j<3;j++)
if (n->tri()->getVertex(j) != v1 && n->tri()->getVertex(j) != v2)return n->tri()->getVertex(j);
if (n->tri()->getVertex(j) != v1 && n->tri()->getVertex(j) != v2)
return n->tri()->getVertex(j);
return 0;
}
MTri3(MTriangle *t, double lc, SMetric3 *m = 0, bidimMeshData * data = 0, GFace *gf = 0);
MTri3(MTriangle *t, double lc, SMetric3 *m = 0, bidimMeshData * data = 0,
GFace *gf = 0);
inline void setTri(MTriangle *t) { base = t; }
inline MTriangle *tri() const { return base; }
inline void setNeigh(int iN , MTri3 *n) { neigh[iN] = n; }
......@@ -144,26 +149,27 @@ class compareTri3Ptr
void connectTriangles(std::list<MTri3*> &);
void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
void connectTriangles(std::set<MTri3*, compareTri3Ptr> &AllTris);
void bowyerWatson(GFace *gf, int MAXPNT= 1000000000,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::map<MVertex*, MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatsonFrontal(GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::map<MVertex*, MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::map<MVertex*, MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatsonParallelograms(GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::map<MVertex*, MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatsonParallelogramsConstrained(GFace *gf,
std::set<MVertex*> constr_vertices,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::set<MVertex*> constr_vertices,
std::map<MVertex*, MVertex*> *equivalence = 0,
std::map<MVertex*, SPoint2>
*parametricCoordinates = 0);
void buildBackGroundMesh (GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
std::map<MVertex*, MVertex*> *equivalence = 0,
std::map<MVertex*, SPoint2> *parametricCoordinates = 0);
void delaunayMeshIn2D(std::vector<MVertex*> &,
std::vector<MTriangle*> &,
......@@ -174,7 +180,7 @@ void delaunayMeshIn2D(std::vector<MVertex*> &,
struct edgeXface
{
MVertex *v[2];
MTri3 * t1;
MTri3 *t1;
int i1;
edgeXface(MTri3 *_t, int iFac) : t1(_t), i1(iFac)
{
......@@ -196,7 +202,9 @@ struct edgeXface
}
inline bool operator == ( const edgeXface &other) const
{
if(v[0]->getNum() == other.v[0]->getNum() && v[1]->getNum() == other.v[1]->getNum()) return true;
if(v[0]->getNum() == other.v[0]->getNum() &&
v[1]->getNum() == other.v[1]->getNum())
return true;
return false;
}
};
......
......@@ -682,7 +682,8 @@ static bool _tryToCollapseThatVertex (GFace *gf,
}
}
// printf("%d %g %g %g %g\n", count, surface_old, surface_new, worst_quality_old , worst_quality_new);
// printf("%d %g %g %g %g\n", count, surface_old, surface_new,
// worst_quality_old , worst_quality_new);
if (worst_quality_new > worst_quality_old ) {
v1->x() = pp.x();v1->y() = pp.y();v1->z() = pp.z();
......
......@@ -109,14 +109,14 @@ class splitQuadRecovery {
MVertex *v = it2->second.first;
v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
v->onWhat()->mesh_vertices.end(), v));
std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0),
f.getVertex(1),v));
std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
f.getVertex(2),v));
std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
f.getVertex(3),v));
std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3),
f.getVertex(0),v));
std::set<MFace, Less_Face>::iterator itf0 =
allFaces.find(MFace(f.getVertex(0), f.getVertex(1),v));
std::set<MFace, Less_Face>::iterator itf1 =
allFaces.find(MFace(f.getVertex(1), f.getVertex(2),v));
std::set<MFace, Less_Face>::iterator itf2 =
allFaces.find(MFace(f.getVertex(2), f.getVertex(3),v));
std::set<MFace, Less_Face>::iterator itf3 =
allFaces.find(MFace(f.getVertex(3), f.getVertex(0),v));
if (itf0 != allFaces.end() && itf1 != allFaces.end() &&
itf2 != allFaces.end() && itf3 != allFaces.end()){
(*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1),
......@@ -136,7 +136,8 @@ class splitQuadRecovery {
// A quad face connected to an hex or a primsm --> leave the quad face as is
if (_toDelete.find(f) == _toDelete.end()){
(*it)->getRegion(iReg)->pyramids.push_back
(new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v));
(new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2),
f.getVertex(3), v));
(*it)->getRegion(iReg)->mesh_vertices.push_back(v);
NBPY++;
}
......@@ -259,18 +260,19 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
getBoundingInfoAndSplitQuads(gr, allBoundingFaces, allBoundingVertices, sqr);
//// TEST
// TEST
{
// std::vector<MVertex*>ALL;
// std::vector<MTetrahedron*> MESH;
// ALL.insert(ALL.begin(),allBoundingVertices.begin(),allBoundingVertices.end());
// delaunayMeshIn3D (ALL,MESH);
// exit(1);
// std::vector<MVertex*>ALL;
// std::vector<MTetrahedron*> MESH;
// ALL.insert(ALL.begin(),allBoundingVertices.begin(),allBoundingVertices.end());
// delaunayMeshIn3D (ALL,MESH);
// exit(1);
}
in.mesh_dim = 3;
in.firstnumber = 1;
int nbvertices_filler = (old_algo_hexa()) ? Filler::get_nbr_new_vertices() : Filler3D::get_nbr_new_vertices();
int nbvertices_filler = (old_algo_hexa()) ? Filler::get_nbr_new_vertices() :
Filler3D::get_nbr_new_vertices();
in.numberofpoints = allBoundingVertices.size() + nbvertices_filler +
LpSmoother::get_nbr_interior_vertices();
in.pointlist = new REAL[in.numberofpoints * 3];
......@@ -400,7 +402,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// re-create the triangular meshes FIXME: this can lead to hanging nodes for
// non manifold geometries (single surface connected to volume)
for(int i = 0; i < out.numberoftrifaces; i++){
// printf("%d %d %d\n",i,out.numberoftrifaces,needParam);
// printf("%d %d %d\n",i,out.numberoftrifaces,needParam);
if (out.trifacemarkerlist[i] > 0) {
MVertex *v[3];
......@@ -482,13 +484,13 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[j] = v1b;
}
}
// printf("new triangle %d %d %d\n",v[0]->onWhat()->dim(), v[1]->onWhat()->dim(), v[2]->onWhat()->dim());
// printf("new triangle %d %d %d\n", v[0]->onWhat()->dim(),
// v[1]->onWhat()->dim(), v[2]->onWhat()->dim());
MTriangle *t = new MTriangle(v[0], v[1], v[2]);
gf->triangles.push_back(t);
}// do not do anything with prismatic faces
} // do not do anything with prismatic faces
}
for(int i = 0; i < out.numberoftetrahedra; i++){
MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] - 1];
MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] - 1];
......@@ -565,7 +567,8 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> &regions)
allEmbVerticesSet.insert(e.begin(), e.end());
}
std::list<GVertex*> allEmbVertices;
allEmbVertices.insert(allEmbVertices.end(), allEmbVerticesSet.begin(), allEmbVerticesSet.end());
allEmbVertices.insert(allEmbVertices.end(), allEmbVerticesSet.begin(),
allEmbVerticesSet.end());
std::list<GVertex*> oldEmbVertices = gr->embeddedVertices();
gr->embeddedVertices() = allEmbVertices;
......@@ -697,7 +700,8 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
allEmbVerticesSet.insert(e.begin(), e.end());
}
std::list<GVertex*> allEmbVertices;
allEmbVertices.insert(allEmbVertices.end(), allEmbVerticesSet.begin(), allEmbVerticesSet.end());
allEmbVertices.insert(allEmbVertices.end(), allEmbVerticesSet.begin(),
allEmbVerticesSet.end());
std::list<GVertex*> oldEmbVertices = gr->embeddedVertices();
gr->embeddedVertices() = allEmbVertices;
......@@ -741,7 +745,8 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
// just to remove tets that are not to be meshed
insertVerticesInRegion(gr, 0);
for(unsigned int i = 0; i < regions.size(); i++){
Msg::Info("Refining volume %d with %d threads",regions[i]->tag(),Msg::GetMaxThreads());
Msg::Info("Refining volume %d with %d threads", regions[i]->tag(),
Msg::GetMaxThreads());
edgeBasedRefinement(Msg::GetMaxThreads(), 1, regions[i]);
}
// RelocateVertices (regions,-1);
......@@ -774,7 +779,8 @@ namespace nglib {
using namespace nglib;
static void getAllBoundingVertices(GRegion *gr,
std::set<MVertex*, MVertexLessThanNum> &allBoundingVertices)
std::set<MVertex*, MVertexLessThanNum>
&allBoundingVertices)
{
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
......@@ -1201,7 +1207,6 @@ GFace *findInFaceSearchStructure(const MFace &ff,
return it->second;
}
GEdge *findInEdgeSearchStructure(MVertex *p1, MVertex *p2, const es_cont &search)
{
MVertex *p = std::min(p1, p2);
......
......@@ -380,9 +380,10 @@ int makeCavityStarShaped (std::list<faceXtet> & shell,
}
}
if (wrong.empty()) return 0;
// printf ("cavity %p (shell size %d cavity size %d)is not star shaped (%d faces not visible), correcting it\n",
// v,shell.size(),cavity.size(),wrong.size());
// bool doneNothing = true;
// printf("cavity %p (shell size %d cavity size %d)is not star shaped "
// "(%d faces not visible), correcting it\n",
// v,shell.size(),cavity.size(),wrong.size());
// bool doneNothing = true;
while (!wrong.empty()){
faceXtet &fxt = *(wrong.begin());
std::list<faceXtet>::iterator its = std::find(shell.begin(),shell.end(),fxt);
......@@ -401,7 +402,7 @@ int makeCavityStarShaped (std::list<faceXtet> & shell,
}
wrong.erase(wrong.begin());
}
// printf("after : shell size %d cavity size %d\n",shell.size(),cavity.size());
// printf("after : shell size %d cavity size %d\n",shell.size(),cavity.size());
return 1;
}
......@@ -922,7 +923,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
double worstA = 0.0;
while (1){
// printf("coucou\n");
std::vector<MTet4*> newTets;
/* for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
......@@ -938,7 +938,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
}
}
*/
// printf("coucou\n");
illegals.clear();
for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
......@@ -974,7 +973,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
break;
}
// printf("coucou\n");
// add all the new tets in the container
for(unsigned int i = 0; i < newTets.size(); i++){
if(!newTets[i]->isDeleted()){
......@@ -985,7 +983,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
delete newTets[i];
}
}
// printf("coucou\n");
// relocate vertices
if (!gr->hexahedra.size() &&
......@@ -1002,7 +999,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
}
}
// printf("coucou\n");
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
......@@ -1017,7 +1013,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
totalVolumeb += vol;
}
}
// printf("coucou\n");
double t2 = Cpu();
Msg::Info("Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",
nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
......@@ -1105,14 +1101,6 @@ double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
circumcenter[1] = ycirca + a[1];
circumcenter[2] = zcirca + a[2];
/*
printf(" %g %g %g %g\n",
sqrt((a[0]-xcirca)*(a[0]-xcirca)+(a[1]-ycirca)*(a[1]-ycirca)+(a[2]-zcirca)*(a[2]-zcirca)),
sqrt((b[0]-xcirca)*(b[0]-xcirca)+(b[1]-ycirca)*(b[1]-ycirca)+(b[2]-zcirca)*(b[2]-zcirca)),
sqrt((c[0]-xcirca)*(c[0]-xcirca)+(c[1]-ycirca)*(c[1]-ycirca)+(c[2]-zcirca)*(c[2]-zcirca)),
sqrt((d[0]-xcirca)*(d[0]-xcirca)+(d[1]-ycirca)*(d[1]-ycirca)+(d[2]-zcirca)*(d[2]-zcirca)) );
*/
if (xi != (double *) NULL) {
/* To interpolate a linear function at the circumcenter, define a */
/* coordinate system with a xi-axis directed from `a' to `b', */
......
......@@ -29,7 +29,7 @@ class MTet4Factory;
// Memory usage for 1 million tets:
//
// * sizeof(MTet4) = 36 Bytes and sizeof(MTetrahedron) = 28 Bytes
// * sizeof(MTet4) = 36 Bytes and sizeof(MTetrahedron) = 28 Bytes
// -> 64 MB
// * rb tree containing all pointers sorted with respect to tet
// radius: each bucket of the tree contains 4 pointers (16 Bytes)
......@@ -37,10 +37,10 @@ class MTet4Factory;
// * sizeof(MVertex) = 44 Bytes and there are about 200000 verts per
// million tet -> 9MB
// * vector of char lengths per vertex -> 1.6Mb
// * vectors in GEntities to store the element and vertex pointers
// * vectors in GEntities to store the element and vertex pointers
// -> 5Mb
//
// Grand total should thus be about 100 MB.
// Grand total should thus be about 100 MB.
//
// The observed mem usage with "demos/cube.geo -clscale 0.61" is
// 157MB. Where do the extra 57 MB come from?
......@@ -60,14 +60,14 @@ class MTet4
MTet4 *neigh[4];
GRegion *gr;
public :
static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm
static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm
~MTet4(){}
MTet4()
MTet4()
: deleted(false), circum_radius(0.0), base(0), gr(0)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
}
MTet4(MTetrahedron *t, double qual)
MTet4(MTetrahedron *t, double qual)
: deleted(false), circum_radius(qual), base(t), gr(0)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
......@@ -114,7 +114,7 @@ class MTet4
base->getVertex(3)->getIndex());
}
*/
*/
double lc1 = 0.25*(sizes[base->getVertex(0)->getIndex()]+
sizes[base->getVertex(1)->getIndex()]+
sizes[base->getVertex(2)->getIndex()]+
......@@ -126,9 +126,10 @@ class MTet4
double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
circum_radius /= lc;
deleted = false;
}
}
void setup(MTetrahedron *t, std::vector<double> &sizes, std::vector<double> &sizesBGM, double lcA, double lcB)
void setup(MTetrahedron *t, std::vector<double> &sizes,
std::vector<double> &sizesBGM, double lcA, double lcB)
{
base = t;
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
......@@ -150,7 +151,7 @@ class MTet4
base->getVertex(3)->getIndex());
}
*/
*/
double lc1 = 0.25*(sizes[base->getVertex(0)->getIndex()]+
sizes[base->getVertex(1)->getIndex()]+
sizes[base->getVertex(2)->getIndex()]+
......@@ -162,22 +163,22 @@ class MTet4
double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
circum_radius /= lc;
deleted = false;
}
}
inline GRegion *onWhat() const { return gr; }
inline void setOnWhat(GRegion *g) { gr = g; }
inline bool isDeleted() const { return deleted; }
inline void forceRadius(double r){ circum_radius = r; }
inline double getRadius() const { return circum_radius; }
inline double getQuality() const { return circum_radius; }
inline void setQuality(const double &q){ circum_radius = q; }
inline double getQuality() const { return circum_radius; }
inline void setQuality(const double &q){ circum_radius = q; }
inline MTetrahedron *tet() const { return base; }
inline MTetrahedron *&tet() { return base; }
inline void setTet(MTetrahedron *t) { base=t; }
inline void setNeigh(int iN, MTet4 *n) { neigh[iN] = n; }
inline MTet4 *getNeigh(int iN) const { return neigh[iN]; }
int inCircumSphere(const double *p) const;
inline int inCircumSphere(double x, double y, double z) const
int inCircumSphere(const double *p) const;
inline int inCircumSphere(double x, double y, double z) const
{
const double p[3] = {x, y, z};
return inCircumSphere(p);
......@@ -186,27 +187,27 @@ class MTet4
{
return inCircumSphere(v->x(), v->y(), v->z());
}
inline double getVolume() const {
double pa[3] = {base->getVertex(0)->x(),
base->getVertex(0)->y(),
inline double getVolume() const {
double pa[3] = {base->getVertex(0)->x(),
base->getVertex(0)->y(),
base->getVertex(0)->z()};
double pb[3] = {base->getVertex(1)->x(),
double pb[3] = {base->getVertex(1)->x(),
base->getVertex(1)->y(),
base->getVertex(1)->z()};
double pc[3] = {base->getVertex(2)->x(),
base->getVertex(2)->y(),
base->getVertex(2)->y(),
base->getVertex(2)->z()};
double pd[3] = {base->getVertex(3)->x(),
base->getVertex(3)->y(),
base->getVertex(3)->z()};
return fabs(robustPredicates::orient3d(pa, pb, pc, pd))/6.0;
return fabs(robustPredicates::orient3d(pa, pb, pc, pd))/6.0;
}
inline void setDeleted(bool d)
{
deleted = d;
}
inline bool assertNeigh() const
inline bool assertNeigh() const
{
if (deleted) return true;
for (int i = 0; i < 4; i++)
......@@ -224,17 +225,19 @@ class MTet4
void connectTets(std::list<MTet4*> &);
void connectTets(std::vector<MTet4*> &);
// IN --> Vertices ---- OUT --> Tets
void delaunayMeshIn3D(std::vector<MVertex*> &, std::vector<MTetrahedron*> &, bool removeBox = true);
void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000, bool _classify = true);
void delaunayMeshIn3D(std::vector<MVertex*> &, std::vector<MTetrahedron*> &,
bool removeBox = true);
void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000,
bool _classify = true);
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
GRegion *getRegionFromBoundingFaces(GModel *model,
std::set<GFace *> &faces_bound);
class compareTet4Ptr
class compareTet4Ptr
{
public:
inline bool operator () (const MTet4 *a, const MTet4 *b) const
{
{
if (a->getRadius() > b->getRadius()) return true;
if (a->getRadius() < b->getRadius()) return false;
return a->tet()->getNum() < b->tet()->getNum();
......@@ -267,7 +270,7 @@ class MTet4Factory
return t;
}
return getANewSlot();
};
};
#endif
public :
MTet4Factory(int _size = 1000000)
......@@ -277,13 +280,13 @@ class MTet4Factory
allSlots = new MTet4[s_alloc];
#endif
}
~MTet4Factory()
~MTet4Factory()
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
delete [] allSlots;
#endif
}
MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
std::vector<double> &sizesBGM)
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
......@@ -294,7 +297,7 @@ class MTet4Factory
t4->setup(t, sizes, sizesBGM);
return t4;
}
MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
std::vector<double> &sizesBGM, double lc1, double lc2)
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
......
......@@ -11,9 +11,10 @@
#include "meshGFaceOptimize.h"
#include "meshGRegionRelocateVertex.h"
static double objective_function (double xi, MVertex *ver,
double xTarget, double yTarget, double zTarget,
const std::vector<MElement*> &lt){
static double objective_function(double xi, MVertex *ver,
double xTarget, double yTarget, double zTarget,
const std::vector<MElement*> &lt)
{
double x = ver->x();