diff --git a/Common/VertexArray.cpp b/Common/VertexArray.cpp index d227bbae86407eb4124b1ece07fdb7a5c38c5065..987f5263fb950726556abb3faa922b17646db7f7 100644 --- a/Common/VertexArray.cpp +++ b/Common/VertexArray.cpp @@ -90,6 +90,7 @@ void VertexArray::add(double *x, double *y, double *z, SVector3 *n, unsigned cha return; } +#if 0 // removed this for now if(unique){ Barycenter pc(0.0F, 0.0F, 0.0F); for(int i = 0; i < npe; i++) @@ -99,6 +100,7 @@ void VertexArray::add(double *x, double *y, double *z, SVector3 *n, unsigned cha return; _barycenters.insert(pc); } +#endif for(int i = 0; i < npe; i++){ _addVertex((float)x[i], (float)y[i], (float)z[i]); diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index ad3e570db064bc392f74708ebfbb155242efc205..ec9fc2202f87c62d52594e4e963c591b5236cb55 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -1858,7 +1858,8 @@ MElement *MElementFactory::create(int num, int type, const std::vector<int> &dat return element; } -double MElement::skewness() { +double MElement::skewness() +{ double minsk = 1.0; for (int i=0;i<getNumFaces();i++){ MFace f = getFace(i); diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp index 26c57d70a1e174bbcedff69817fcb71b453b549b..aadfe0eee8c57d51dede13e4f6784228e9b12cee 100644 --- a/Geo/MTetrahedron.cpp +++ b/Geo/MTetrahedron.cpp @@ -9,7 +9,6 @@ #include "Context.h" #include "BasisFactory.h" - #if defined(HAVE_MESH) #include "qualityMeasures.h" #include "meshGFaceDelaunayInsertion.h" @@ -18,6 +17,36 @@ #define SQU(a) ((a)*(a)) +void MTetrahedron::getEdgeRep(bool curved, int num, double *x, double *y, double *z, + SVector3 *n) +{ + // don't use MElement::_getEdgeRep: it's slow due to the creation of MFaces + MVertex *v0 = getVertex(edges_tetra(num, 0)); + MVertex *v1 = getVertex(edges_tetra(num, 1)); + x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z(); + x[1] = v1->x(); y[1] = v1->y(); z[1] = v1->z(); +#if 0 // compute normals so that we can light the edges + static const int vv[6] = {2, 0, 1, 1, 0, 2}; + MVertex *v2 = getVertex(vv[num]); + double nn[3]; + normal3points(v0->x(), v0->y(), v0->z(), + v1->x(), v1->y(), v1->z(), + v2->x(), v2->y(), v2->z(), nn); + n[0] = n[1] = SVector3(nn[0], nn[1], nn[2]); +#else + n[0] = n[1] = SVector3(); +#endif +} + +void MTetrahedron::getFaceRep(bool curved, int num, double *x, double *y, double *z, + SVector3 *n) +{ + // don't use general MElement::_getFaceRep: it's slow due to the creation of + // MFaces + MFace f(getFace(num)); + _getFaceRep(f.getVertex(0), f.getVertex(1), f.getVertex(2), x, y, z, n); +} + SPoint3 MTetrahedron::circumcenter() { #if defined(HAVE_MESH) diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h index 1b25040f3f0f7a5024f87542194a9b6cad05c6e2..2adaf8aa354fd5eb83a357c9359a4186adba5bce 100644 --- a/Geo/MTetrahedron.h +++ b/Geo/MTetrahedron.h @@ -68,12 +68,7 @@ class MTetrahedron : public MElement { return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]); } virtual int getNumEdgesRep(bool curved){ return 6; } - virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) - { - static const int f[6] = {0, 0, 0, 1, 2, 3}; - MEdge e(getEdge(num)); - _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, f[num]); - } + virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n); virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const { v.resize(2); @@ -88,11 +83,7 @@ class MTetrahedron : public MElement { } virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign, int &rot) const; virtual int getNumFacesRep(bool curved){ return 4; } - virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) - { - MFace f(getFace(num)); - _getFaceRep(f.getVertex(0), f.getVertex(1), f.getVertex(2), x, y, z, n); - } + virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n); virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const { v.resize(3); diff --git a/Geo/SPoint3.h b/Geo/SPoint3.h index ea44797c98650e8d431e3353865b3805a7a3d2a2..dd69e00acd1438397e88098f00ad4edbc20889db 100644 --- a/Geo/SPoint3.h +++ b/Geo/SPoint3.h @@ -18,7 +18,12 @@ class SPoint3 { SPoint3(const SPoint3 &pt) { P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; } virtual ~SPoint3() {} void setPosition(double xx, double yy, double zz); - void setPosition(const SPoint3 &pt,const SPoint3 &dir,const double dist_) {P[0]=pt.P[0]; P[1]=pt.P[1]; P[2]=pt.P[2]; SPoint3 a(dir); a*=dist_; P[0]+=a[0]; P[1]+=a[1]; P[2]+=a[2];} + void setPosition(const SPoint3 &pt,const SPoint3 &dir,const double dist_) + { + P[0]=pt.P[0]; P[1]=pt.P[1]; P[2]=pt.P[2]; + SPoint3 a(dir); a*=dist_; + P[0]+=a[0]; P[1]+=a[1]; P[2]+=a[2]; + } void getPosition(double *xx, double *yy, double *zz) const; void position(double *) const; inline double x(void) const; diff --git a/Geo/SVector3.h b/Geo/SVector3.h index a8cc7dd84d6acc068027dce684a988e6d6f29831..967e92c84f0c3184e2be024440e00944c4fc081a 100644 --- a/Geo/SVector3.h +++ b/Geo/SVector3.h @@ -33,7 +33,8 @@ class SVector3 { // Beware that " w = v.normalize() " produces the vector // w = (v.norm(), v.norm(), v.norm()), which is NOT a unit vector! // Use " w = v.unit() " to affect to "w" the unit vector parallel to "v". - double normalize(){ + double normalize() + { double n = norm(); if(n){ P[0] /= n; P[1] /= n; P[2] /= n; } return n; } @@ -71,11 +72,14 @@ class SVector3 { } operator double *() { return P; } void print(std::string name="") const - { printf("Vector \'%s\': %f %f %f\n",name.c_str(),P[0],P[1],P[2]); } - - // Needed to allow the initialization of a SPoint3 from a SPoint3, a distance and a direction - SPoint3 point() const{return P;} - int getMaxValue(double& val) const{ + { + Msg::Info("Vector \'%s\': %f %f %f", name.c_str(), P[0], P[1], P[2]); + } + // Needed to allow the initialization of a SPoint3 from a SPoint3, a distance + // and a direction + SPoint3 point() const{ return P; } + int getMaxValue(double& val) const + { if ((P[0] >=P[1]) && (P[0]>=P[2])){ val = P[0]; return 0; @@ -166,8 +170,7 @@ inline void buildOrthoBasis_naive(SVector3 &dir, SVector3 &dir1, SVector3 &dir2) dir2 = SVector3(0.0, 0.0, 1.0); } else{ - Msg::Error("Problem with computing orthoBasis"); - Msg::Exit(1); + Msg::Error("Problem with computing orthoBasis"); } // printf("XYZ =%g %g %g r=%g dir0 = %g %g %g dir1 = %g %g %g dir2 =%g %g %g\n", // x,y,z,d, dir[0], dir[1], dir[2], dir1[0], dir1[1], dir1[2], dir2[0], dir2[1], dir2[2] ); diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp index 34ff4484b531147e7bcad97febea99d2e624ae36..709a286ffe2047e2e0126804006866504d08725d 100644 --- a/Mesh/Field.cpp +++ b/Mesh/Field.cpp @@ -1198,22 +1198,22 @@ class Popen2 { HANDLE hChildStd_OUT_Wr, hChildStd_IN_Rd; PROCESS_INFORMATION piProcInfo; SECURITY_ATTRIBUTES saAttr; - saAttr.nLength = sizeof(SECURITY_ATTRIBUTES); - saAttr.bInheritHandle = TRUE; - saAttr.lpSecurityDescriptor = NULL; - if (!CreatePipe(&_hIn, &hChildStd_OUT_Wr, &saAttr, 0)) - Msg::Error("StdoutRd CreatePipe"); - if (!CreatePipe(&hChildStd_IN_Rd, &_hOut, &saAttr, 0)) + saAttr.nLength = sizeof(SECURITY_ATTRIBUTES); + saAttr.bInheritHandle = TRUE; + saAttr.lpSecurityDescriptor = NULL; + if (!CreatePipe(&_hIn, &hChildStd_OUT_Wr, &saAttr, 0)) + Msg::Error("StdoutRd CreatePipe"); + if (!CreatePipe(&hChildStd_IN_Rd, &_hOut, &saAttr, 0)) Msg::Error("Stdin CreatePipe"); - if (!CreatePipe(&_hIn, &hChildStd_OUT_Wr, &saAttr, 0)) - Msg::Error("StdoutRd CreatePipe"); + if (!CreatePipe(&_hIn, &hChildStd_OUT_Wr, &saAttr, 0)) + Msg::Error("StdoutRd CreatePipe"); if (!SetHandleInformation(_hIn, HANDLE_FLAG_INHERIT, 0)) Msg::Error("Stdout SetHandleInformation"); STARTUPINFO siStartInfo; - BOOL bSuccess = FALSE; + BOOL bSuccess = FALSE; ZeroMemory( &piProcInfo, sizeof(PROCESS_INFORMATION) ); ZeroMemory( &siStartInfo, sizeof(STARTUPINFO) ); - siStartInfo.cb = sizeof(STARTUPINFO); + siStartInfo.cb = sizeof(STARTUPINFO); siStartInfo.hStdError = GetStdHandle(STD_ERROR_HANDLE); siStartInfo.hStdOutput = hChildStd_OUT_Wr; siStartInfo.hStdInput = hChildStd_IN_Rd; @@ -1301,7 +1301,7 @@ class ExternalProcessField : public Field Popen2 _pipes; void closePipes() { if (_pipes.started()) { - double xyz[3] = { + double xyz[3] = { std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN() }; @@ -2243,7 +2243,7 @@ void BoundaryLayerField::removeAttractors() } void BoundaryLayerField::setupFor1d(int iE) { - if (edges_id_saved.empty() ){ + if (edges_id_saved.empty() ){ edges_id_saved = edges_id; nodes_id_saved = nodes_id; } @@ -2253,7 +2253,7 @@ void BoundaryLayerField::setupFor1d(int iE) { bool found = std::find(edges_id_saved.begin(), edges_id_saved.end(), iE) != edges_id_saved.end(); - + if (!found) { GEdge *ge = GModel::current()->getEdgeByTag(iE); GVertex *gv0 = ge->getBeginVertex(); @@ -2277,21 +2277,21 @@ void BoundaryLayerField::setupFor2d(int iF) remove GFaces from the attarctors (only used in 2D) for edges and vertices */ - if (edges_id_saved.empty()){ + if (edges_id_saved.empty()){ edges_id_saved = edges_id; nodes_id_saved = nodes_id; } nodes_id.clear(); edges_id.clear(); - + // printf("have %d %d\n",faces_id_saved.size(),edges_id_saved.size()); - + /// FIXME : /// NOT REALLY A NICE WAY TO DO IT (VERY AD HOC) /// THIS COULD BE PART OF THE INPUT /// OR (better) CHANGE THE PHILOSOPHY - + GFace *gf = GModel::current()->getFaceByTag(iF); std::list<GEdge*> ed = gf->edges(); std::list<GEdge*> embedded_edges = gf->embeddedEdges(); @@ -2311,7 +2311,8 @@ void BoundaryLayerField::setupFor2d(int iF) // one only face --> 2D --> BL if (fc.size() <= 1) isIn = true; else { - Msg::Error ("Only 2D Boundary Layers are supported (edge %d is adjacet to %d faces\n",iE,fc.size()); + Msg::Error ("Only 2D Boundary Layers are supported (edge %d is adjacet to %d faces", + iE, fc.size()); } } if (isIn){ @@ -2358,7 +2359,7 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge double lc = dist*(ratio-1) + hwall_n; // double lc = hwall_n; - // double lc = hwall_n * pow (ratio, dist / hwall_t); + // double lc = hwall_n * pow (ratio, dist / hwall_t); return std::min (hfar,lc); } @@ -2403,7 +2404,7 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist, { // printf("WHAT THE FUCK\n"); - + // dist = hwall -> lc = hwall * ratio // dist = hwall (1+ratio) -> lc = hwall ratio ^ 2 // dist = hwall (1+ratio+ratio^2) -> lc = hwall ratio ^ 3 diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp index e0145cea9df6bdcc91a3d80489c80334e347d361..7b4d61e33b5c6f90aca1e81ae9dac5c938625cfd 100644 --- a/Numeric/Numeric.cpp +++ b/Numeric/Numeric.cpp @@ -34,7 +34,9 @@ double myacos(double a) else return acos(a); } -double norm2(double a[3][3]) { + +double norm2(double a[3][3]) +{ double norm2sq = gmsh_SQU(a[0][0])+ gmsh_SQU(a[0][1])+ @@ -96,7 +98,6 @@ void normal2points(double x0, double y0, double z0, norme(n); } - int sys2x2(double mat[2][2], double b[2], double res[2]) { double det, ud, norm; @@ -448,19 +449,18 @@ double computeInnerRadiusForQuad(double *x, double *y, int i) char float2char(float f) { // float normalized in [-1, 1], char in [-127, 127] - f *= 127.; - if(f > 127.) return 127; - else if(f < -127.) return -127; - else return (char)f; + float c = f * 127.; + if(c > 127.) return 127; + else if(c < -127.) return -127; + else return (int)c; } float char2float(char c) { float f = c; - f /= 127.; - if(f > 1.) return 1.; - else if(f < -1) return -1.; - else return f; + if(f > 127.) return 1.; + else if(f < -127) return -1.; + else return f / 127.; } void gradSimplex(double *x, double *y, double *z, double *v, double *grad) diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h index fe3d1de873eb7a86e9e4292dfc79ce4e9531fa52..fbace5405d09e293004621eb6bf8b74ad4178ba2 100644 --- a/Numeric/Numeric.h +++ b/Numeric/Numeric.h @@ -60,6 +60,7 @@ inline void prodve(double a[3], double b[3], double c[3]) c[1] = -a[0] * b[2] + a[2] * b[0]; c[0] = a[1] * b[2] - a[2] * b[1]; } + inline void prosca(double a[3], double b[3], double *c) { *c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt index 2018c483a8b9d20c0176154bf3ec34e74a715afb..0059dc46cd14b4ad86233bf5cab6a3642d9584ba 100644 --- a/doc/VERSIONS.txt +++ b/doc/VERSIONS.txt @@ -1,3 +1,6 @@ +2.14.2: fixed regression in multi-file partitioned grid export; fixed regression +for Mesh.SubdivisionAlgorithm=1; simplified 2D boundary layer field. + 2.14.1 (October 30, 2016): fixed regression in periodic meshes; small bug fixes and code cleanups.