diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index 46bcf30383efdc03631de7d73339fd851d8f1110..0cbbb498d80e812e373a666cd72baea16d56c2d1 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -519,8 +519,8 @@ bool GFaceCompound::parametrize() const else if (_mapping == CONFORMAL){ Msg::Debug("Parametrizing surface %d with 'conformal map'", tag()); fillNeumannBCS(); - bool withoutFolding = parametrize_conformal_spectral() ; - //bool withoutFolding = parametrize_conformal(); + //bool withoutFolding = parametrize_conformal_spectral() ; + bool withoutFolding = parametrize_conformal(); if ( withoutFolding == false ){ //printStuff(); exit(1); Msg::Warning("$$$ Parametrization switched to harmonic map"); @@ -1822,7 +1822,10 @@ bool GFaceCompound::checkTopology() const // FIXME!!! I think those things are wrong with cross-patch reparametrization //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true; + + //TODO: smthg to exit here for lloyd remeshing + bool correctTopo = true; if(allNodes.empty()) buildAllNodes(); @@ -1913,7 +1916,7 @@ double GFaceCompound::checkAspectRatio() const } double AR = 2*3.14*area3D/(tot_length*tot_length); - if (areaMin < limit && nb > 2) { + if (areaMin > 0 && areaMin < limit && nb > 2) { Msg::Warning("Too small triangles in mapping (a_2D=%g)", areaMin); } else { diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h index eeab87d002a610459cf2a44f661814b0bbdb2424..b92b76d8e8233cd3b7685ebd545094738d5e5b1c 100644 --- a/Geo/GFaceCompound.h +++ b/Geo/GFaceCompound.h @@ -121,6 +121,7 @@ class GFaceCompound : public GFace { void partitionFaceCM(); virtual std::list<GFace*> getCompounds() const {return _compound;}; mutable int nbSplit; + int getNbSplit() const {return nbSplit;} int allowPartition() const{return _allowPartition;}; private: typeOfIsomorphism _type; diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp index 6aa6d5d697f7aff28ed3080892c53058643811ca..f9759546a97bc735205897666c4cc4061134c924 100644 --- a/Mesh/Generator.cpp +++ b/Mesh/Generator.cpp @@ -435,12 +435,27 @@ static void Mesh2D(GModel *m) // lloyd optimization if (CTX::instance()->mesh.optimizeLloyd){ + Msg::Info("------------------------------"); for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){ GFace *gf = *it; + if(gf->geomType() == GEntity::DiscreteSurface) continue; + if(gf->geomType() == GEntity::CompoundSurface ) { + GFaceCompound *gfc = (GFaceCompound*) gf; + if (gfc->getNbSplit() != 0) continue; + } int recombine = gf->meshAttributes.recombine; - Msg::Info("Lloyd optimization for face %d", gf->tag()); - gf->lloyd(40, recombine); - if(recombine) recombineIntoQuads(gf); + + Msg::StatusBar(2, true, "Lloyd optimization for face %d", gf->tag()); + gf->lloyd(25,recombine); + + if(recombine) recombineIntoQuads(gf); + + // computeElementShapes(gf, gf->meshStatistics.worst_element_shape, +// gf->meshStatistics.average_element_shape, +// gf->meshStatistics.best_element_shape, +// gf->meshStatistics.nbTriangle, +// gf->meshStatistics.nbGoodQuality); + } } diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp index d7e254695595697f8c40d81b676f4c1a91c23017..7f4b30d93376a5a664260e99fb4155010900162d 100644 --- a/Mesh/meshGFaceLloyd.cpp +++ b/Mesh/meshGFaceLloyd.cpp @@ -68,6 +68,7 @@ void lloydAlgorithm::operator () ( GFace * gf) { fullMatrix<double> cgs(triangulator.numPoints,2); // now iterate on internal vertices double ENERGY = 0.0; + double criteria = 0.0; for (int i=0; i<triangulator.numPoints;i++){ // get the ith vertex PointRecord &pt = triangulator.points[i]; @@ -76,26 +77,23 @@ void lloydAlgorithm::operator () ( GFace * gf) { // get the voronoi corners std::vector<SPoint2> pts; triangulator.voronoiCell (i,pts); - double E; + double E, A; SPoint2 p(pt.where.h,pt.where.v); if (!infiniteNorm){ - centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E,backgroundMesh::current()); + centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E, A, NULL); //backgroundMesh::current()); } else { - centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E); + centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A); } ENERGY += E; + double d = sqrt((p.x()-cgs(i,0))*(p.x()-cgs(i,0))+ + (p.y()-cgs(i,1))*(p.y()-cgs(i,1))); + criteria = d/A; }// if (v->onWhat() == gf) else { } }// for all points - if (ITER % 10 == 0){ - char name[234]; - sprintf(name,"LloydIter%d.pos",ITER); - triangulator.makePosView(name); - } - for(PointNumero i = 0; i < triangulator.numPoints; i++) { MVertex *v = (MVertex*)triangulator.points[i].data; if (v->onWhat() == gf && !triangulator.onHull(i)){ @@ -104,11 +102,18 @@ void lloydAlgorithm::operator () ( GFace * gf) { } } - Msg::Debug("GFace %d Lloyd iteration %d energy %g",gf->tag(),ITER++,ENERGY); + Msg::Debug("GFace %d Lloyd-iter %d Inertia=%g Convergence=%g ",gf->tag(),ITER++,ENERGY, criteria); if (ITER > ITER_MAX)break; // compute the Voronoi diagram triangulator.Voronoi(); + + if (ITER % 10 == 0){ + char name[234]; + sprintf(name,"LloydIter%d.pos",ITER); + triangulator.makePosView(name); + } + } // now create the vertices @@ -132,6 +137,7 @@ void lloydAlgorithm::operator () ( GFace * gf) { // vertices gf->_additional_vertices = mesh_vertices; // remesh the face with all those vertices in + Msg::Info("Lloyd remeshing of face %d ", gf->tag()); meshGFace mesher; mesher(gf); // assign those vertices to the face internal vertices diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp index ddc3f1c8b9794d26bec8422a1cc3e3a9f53a3e02..6bd532988ef5d9562b6b24d0e97dd0d9af96434f 100644 --- a/Numeric/DivideAndConquer.cpp +++ b/Numeric/DivideAndConquer.cpp @@ -592,7 +592,7 @@ void DocRecord::makePosView(std::string fileName) } void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle, - double &xc, double &yc, double &inertia) + double &xc, double &yc, double &inertia, double &area) { const int N = pts.size(); @@ -614,13 +614,15 @@ void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle, xc = XC*cosa - YC *sina; yc = XC*sina + YC *cosa; inertia = std::max(xmax-xmin,ymax-ymin); + area = (xmax-xmin) * (ymax-ymin); } void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts, - double &xc, double &yc, double &inertia, + double &xc, double &yc, double &inertia, double &areaCell, simpleFunction<double> *bgm) { double area_tot = 0; + areaCell = 0.0; SPoint2 center(0,0); for (unsigned int j = 0; j < pts.size(); j++){ SPoint2 &pa = pts[j]; @@ -629,7 +631,8 @@ void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts, const double lc = bgm ? (*bgm)((pa.x()+pb.x()+pc.x())/3.0, (pa.y()+pb.y()+pc.y())/3.0, 0.0) : 1.0; - const double fact = 1./(lc*lc); + const double fact = 1./(lc*lc*lc*lc);//rho = 1/lc^4 (emi) + areaCell += area; area_tot += area*fact; center += ((pa+pb+pc) * (area*fact/3.0)); } @@ -664,14 +667,14 @@ double DocRecord::Lloyd(int type) PointRecord &pt = points[i]; std::vector<SPoint2> pts; voronoiCell (i,pts); - double E; + double E, A; if (!points[i].data){ SPoint2 p (pt.where.h,pt.where.v); if (type == 0) - centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E); + centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E, A); else - centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E); + centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A); } inertia_tot += E; } diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h index e44784dcf98c21226b24ad5da2c596c055e3ce51..b8093c27875ac5e5a26c9765129a07a7a416a267 100644 --- a/Numeric/DivideAndConquer.h +++ b/Numeric/DivideAndConquer.h @@ -104,12 +104,14 @@ void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle, double &xc, double &yc, - double &inertia); + double &inertia, + double &area); void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts, double &xc, double &yc, double &inertia, + double &areaCell, simpleFunction<double> *bgm = 0); diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp index 353665900f35516fa90abd3a519271a4df29dd13..4e58cbad6c56eef03eb04cc24c363a06a094e3be 100644 --- a/Solver/eigenSolver.cpp +++ b/Solver/eigenSolver.cpp @@ -53,7 +53,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which) // set some default options _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); - _try(EPSSetTolerances(eps,1.e-7,50)); + _try(EPSSetTolerances(eps,1.e-6,20));//1.e-7 50 //_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default _try(EPSSetType(eps,EPSARNOLDI)); //_try(EPSSetType(eps,EPSARPACK));