From d0b4ede439e2024716c77fdcc6174611df00501e Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Thu, 24 Jan 2013 19:17:46 +0000 Subject: [PATCH] tweaks --- Geo/GEdge.cpp | 2 +- Geo/GFace.cpp | 18 +++++----- Geo/GPoint.cpp | 5 --- Geo/GRegion.cpp | 88 ++++++++++++++++++++++------------------------ Geo/gmshEdge.cpp | 2 +- Mesh/meshGEdge.cpp | 6 ++-- 6 files changed, 56 insertions(+), 65 deletions(-) delete mode 100644 Geo/GPoint.cpp diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp index 940d6e59ca..76ae91e04c 100644 --- a/Geo/GEdge.cpp +++ b/Geo/GEdge.cpp @@ -172,7 +172,7 @@ void GEdge::setVisibility(char val, bool recursive) std::string GEdge::getAdditionalInfoString() { std::ostringstream sstream; - if(v0 && v1) sstream << "{" << v0->tag() << "," << v1->tag() << "}"; + if(v0 && v1) sstream << "{" << v0->tag() << " " << v1->tag() << "}"; return sstream.str(); } diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index 6dcb922de7..d90e9b0933 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -614,8 +614,6 @@ void GFace::computeMeshSizeFieldAccuracy(double &avg,double &max_e, double &min_ min_e = std::min(min_e, lone); } #endif - //printf("Emi efficiency tau (%g) =%g nE=%d nT=%d \n", avg, 100*exp(avg/(double)nE), nE, getNumMeshElements()); - } double GFace::curvatureDiv(const SPoint2 ¶m) const @@ -902,7 +900,8 @@ class data_wrapper{ }; // Callback function for BFGS -void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr) +void bfgs_callback(const alglib::real_1d_array& x, double& func, + alglib::real_1d_array& grad, void* ptr) { data_wrapper* w = static_cast<data_wrapper*>(ptr); SPoint3 p = w->get_point(); @@ -950,9 +949,11 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[ // fprintf(F,"View \" \" {\n"); // fprintf(F,"SP(%g,%g,%g) {%g};\n",queryPoint.x(),queryPoint.y(),queryPoint.z(),0.0); double initial_guesses = 10.0; - for(double u = uu.low(); u<=uu.high()+1.e-5; u+=(uu.high()-uu.low())/initial_guesses) { + for(double u = uu.low(); u <= uu.high() + 1.e-5; + u += (uu.high() - uu.low()) / initial_guesses) { // printf("%f\n", u); - for(double v = vv.low(); v<=vv.high()+1.e-5; v+=(vv.high()-vv.low())/initial_guesses) { + for(double v = vv.low(); v <= vv.high() + 1.e-5; + v += (vv.high() - vv.low()) / initial_guesses) { GPoint pnt = point(u, v); SPoint3 spnt(pnt.x(), pnt.y(), pnt.z()); double dist = queryPoint.distance(spnt); @@ -993,19 +994,19 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[ double epsx = 0.0; alglib::ae_int_t maxits = 500; - minlbfgssetcond(state,epsg,epsf,epsx,maxits); + minlbfgssetcond(state, epsg, epsf, epsx, maxits); // Solving the problem data_wrapper w; w.set_point(queryPoint); w.set_face(this); - minlbfgsoptimize(state,bfgs_callback,NULL,&w); + minlbfgsoptimize(state, bfgs_callback, NULL, &w); // Getting the results alglib::minlbfgsreport rep; - minlbfgsresults(state,x,rep); + minlbfgsresults(state, x, rep); GPoint pntF = point(x[0], x[1]); if (rep.terminationtype != 4){ @@ -1286,7 +1287,6 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint2> *uvpoints, std::vector<SVector3> *normals) { - if(!points) return false; if (buildSTLTriangulation()){ diff --git a/Geo/GPoint.cpp b/Geo/GPoint.cpp deleted file mode 100644 index 4774096e7f..0000000000 --- a/Geo/GPoint.cpp +++ /dev/null @@ -1,5 +0,0 @@ -// Gmsh - Copyright (C) 1997-2013 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@geuz.org>. - diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp index 8331eb8fbd..c77b8edfee 100644 --- a/Geo/GRegion.cpp +++ b/Geo/GRegion.cpp @@ -209,13 +209,10 @@ void GRegion::setVisibility(char val, bool recursive) std::string GRegion::getAdditionalInfoString() { std::ostringstream sstream; - if(l_faces.size() > 20){ - sstream << "{" << l_faces.front()->tag() << ",...," << l_faces.back()->tag() << "}"; - } - else if(l_faces.size()){ + if(l_faces.size()){ sstream << "{"; for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); ++it){ - if(it != l_faces.begin()) sstream << ","; + if(it != l_faces.begin()) sstream << " "; sstream << (*it)->tag(); } sstream << "}"; @@ -288,15 +285,15 @@ bool GRegion::edgeConnected(GRegion *r) const return false; } -void GRegion::replaceFaces (std::list<GFace*> &new_faces) +void GRegion::replaceFaces(std::list<GFace*> &new_faces) { - replaceFacesInternal (new_faces); + replaceFacesInternal(new_faces); if (l_faces.size() != new_faces.size()){ - Msg::Fatal("impossible to replace faces in region %d (%d vs %d)",tag(), - l_faces.size(),new_faces.size()); + Msg::Error("Impossible to replace faces in region %d (%d vs %d)", + tag(), l_faces.size(), new_faces.size()); } - std::list<GFace*>::iterator it = l_faces.begin(); + std::list<GFace*>::iterator it = l_faces.begin(); std::list<GFace*>::iterator it2 = new_faces.begin(); std::list<int>::iterator it3 = l_dirs.begin(); std::list<int> newdirs; @@ -312,18 +309,17 @@ void GRegion::replaceFaces (std::list<GFace*> &new_faces) l_dirs = newdirs; } -double GRegion::computeSolidProperties (std::vector<double> cg, - std::vector<double> inertia) +double GRegion::computeSolidProperties(std::vector<double> cg, + std::vector<double> inertia) { std::list<GFace*>::iterator it = l_faces.begin(); - std::list<int>::iterator itdir = l_dirs.begin(); + std::list<int>::iterator itdir = l_dirs.begin(); double volumex = 0; double volumey = 0; double volumez = 0; double surface = 0; cg[0] = cg[1] = cg[2] = 0.0; for ( ; it != l_faces.end(); ++it,++itdir){ - //printf("face %d dir %d %d elements\n",(*it)->tag(),*itdir,(int)(*it)->triangles.size()); for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){ MTriangle *e = (*it)->triangles[i]; int npt; @@ -332,65 +328,65 @@ double GRegion::computeSolidProperties (std::vector<double> cg, for (int j=0;j<npt;j++){ SPoint3 pt; // compute x,y,z of the integration point - e->pnt(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],pt); + e->pnt(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], pt); double jac[3][3]; // compute normal - double detJ = e->getJacobian(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],jac); - SVector3 n (jac[2][0],jac[2][1],jac[2][2]); + double detJ = e->getJacobian(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], jac); + SVector3 n(jac[2][0], jac[2][1], jac[2][2]); n.normalize(); n *= (double)*itdir; - surface += detJ*pts[j].weight; - volumex += detJ*n.x()*pt.x()*pts[j].weight; - volumey += detJ*n.y()*pt.y()*pts[j].weight; - volumez += detJ*n.z()*pt.z()*pts[j].weight; - cg[0] += detJ*n.x()*(pt.x()*pt.x())*pts[j].weight*0.5; - cg[1] += detJ*n.y()*(pt.y()*pt.y())*pts[j].weight*0.5; - cg[2] += detJ*n.z()*(pt.z()*pt.z())*pts[j].weight*0.5; + surface += detJ* pts[j].weight; + volumex += detJ * n.x() * pt.x() * pts[j].weight; + volumey += detJ * n.y() * pt.y() * pts[j].weight; + volumez += detJ * n.z() * pt.z() * pts[j].weight; + cg[0] += detJ * n.x() * (pt.x() * pt.x()) * pts[j].weight * 0.5; + cg[1] += detJ * n.y() * (pt.y() * pt.y()) * pts[j].weight * 0.5; + cg[2] += detJ * n.z() * (pt.z() * pt.z()) * pts[j].weight * 0.5; } } } - printf("%g -- %g %g %g\n",surface,volumex,volumey,volumez); + printf("%g -- %g %g %g\n", surface, volumex, volumey, volumez); double volume = volumex; - cg[0]/=volume; - cg[1]/=volume; - cg[2]/=volume; + cg[0] /= volume; + cg[1] /= volume; + cg[2] /= volume; it = l_faces.begin(); - itdir = l_dirs.begin(); - inertia[0] = - inertia[1] = - inertia[2] = - inertia[3] = - inertia[4] = - inertia[5] = 0.0; + itdir = l_dirs.begin(); + inertia[0] = inertia[1] = inertia[2] = inertia[3] = inertia[4] = inertia[5] = 0.0; for ( ; it != l_faces.end(); ++it,++itdir){ for (unsigned int i = 0; i < (*it)->getNumMeshElements(); ++i){ MElement *e = (*it)->getMeshElement(i); int npt; IntPt *pts; - e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts); + e->getIntegrationPoints(2 * (e->getPolynomialOrder() - 1) + 3, &npt, &pts); for (int j = 0; j < npt; j++){ SPoint3 pt; // compute x,y,z of the integration point - e->pnt(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],pt); + e->pnt(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], pt); double jac[3][3]; // compute normal - double detJ = e->getJacobian(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],jac); - SVector3 n (jac[2][0],jac[2][1],jac[2][2]); + double detJ = e->getJacobian(pts[j].pt[0], pts[j].pt[1], pts[j].pt[2], jac); + SVector3 n(jac[2][0], jac[2][1], jac[2][2]); n *= (double)*itdir; - inertia[0] += pts[j].weight*detJ*n.x()*(pt.x()-cg[0])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0; - inertia[1] += pts[j].weight*detJ*n.y()*(pt.y()-cg[1])*(pt.y()-cg[1])*(pt.y()-cg[1])/3.0; - inertia[2] += pts[j].weight*detJ*n.z()*(pt.z()-cg[2])*(pt.z()-cg[2])*(pt.z()-cg[2])/3.0; - inertia[3] += pts[j].weight*detJ*n.x()*(pt.y()-cg[1])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0; - inertia[4] += pts[j].weight*detJ*n.x()*(pt.z()-cg[2])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0; - inertia[5] += pts[j].weight*detJ*n.y()*(pt.z()-cg[2])*(pt.y()-cg[1])*(pt.y()-cg[1])/3.0; + inertia[0] += pts[j].weight * detJ * n.x() * + (pt.x() - cg[0]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0; + inertia[1] += pts[j].weight * detJ * n.y() * + (pt.y() - cg[1]) * (pt.y() - cg[1]) * (pt.y() - cg[1]) / 3.0; + inertia[2] += pts[j].weight * detJ * n.z() * + (pt.z() - cg[2]) * (pt.z() - cg[2]) * (pt.z() - cg[2]) / 3.0; + inertia[3] += pts[j].weight * detJ * n.x() * + (pt.y() - cg[1]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0; + inertia[4] += pts[j].weight * detJ * n.x() * + (pt.z() - cg[2]) * (pt.x() - cg[0]) * (pt.x() - cg[0]) / 3.0; + inertia[5] += pts[j].weight * detJ * n.y() * + (pt.z() - cg[2]) * (pt.y() - cg[1]) * (pt.y() - cg[1]) / 3.0; } } } return volume; } - diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp index 69b29c8e17..394d2d419c 100644 --- a/Geo/gmshEdge.cpp +++ b/Geo/gmshEdge.cpp @@ -77,7 +77,7 @@ std::string gmshEdge::getAdditionalInfoString() std::ostringstream sstream; sstream << "{"; for(int i = 0; i < List_Nbr(c->Control_Points); i++){ - if(i) sstream << ","; + if(i) sstream << " "; Vertex *v; List_Read(c->Control_Points, i, &v); sstream << v->Num; diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index 8ef3ca968e..21e4d12410 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -372,8 +372,8 @@ void meshGEdge::operator() (GEdge *ge) N = ge->meshAttributes.nbPointsTransfinite; } else{ - if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG - // || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS + if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG + // || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS || blf){ a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points, CTX::instance()->mesh.lcIntegrationPrecision); @@ -393,7 +393,7 @@ void meshGEdge::operator() (GEdge *ge) N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99)); } - // force odd number of points for if blossom is used for recombination + // force odd number of points if blossom is used for recombination if(ge->meshAttributes.Method != MESH_TRANSFINITE && CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){ if(/* 1 ||*/ CTX::instance()->mesh.recombineAll){ -- GitLab