From 02b11a027bce92f1619b68c32ff2c08d5fdb2202 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Fri, 14 Nov 2008 15:13:29 +0000 Subject: [PATCH] *** empty log message *** --- Geo/GEdge.h | 5 +++-- Geo/GModelIO_Mesh.cpp | 11 +++++++++-- Mesh/meshGEdge.cpp | 37 ++++++++++++++++++------------------- 3 files changed, 30 insertions(+), 23 deletions(-) diff --git a/Geo/GEdge.h b/Geo/GEdge.h index d9cf681516..834283d28b 100644 --- a/Geo/GEdge.h +++ b/Geo/GEdge.h @@ -103,8 +103,9 @@ class GEdge : public GEntity { virtual double prescribedMeshSizeAtVertex() const { return meshAttributes.meshSize; } // true if start == end and no more than 2 segments - void setTooSmall ( bool b) {_tooSmall = b;} - bool isMeshDegenerated() const{ return _tooSmall || (v0 == v1 && mesh_vertices.size() < 2); } + void setTooSmall(bool b) { _tooSmall = b; } + bool isMeshDegenerated() const + { return _tooSmall || (v0 == v1 && mesh_vertices.size() < 2); } // number of types of elements int getNumElementTypes() const { return 1; } diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp index a2a4f68e9d..3cef450ce8 100644 --- a/Geo/GModelIO_Mesh.cpp +++ b/Geo/GModelIO_Mesh.cpp @@ -2426,13 +2426,20 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll, std::vector<GEntity*> entities; getEntities(entities); + // find max dimension of mesh elements we need to save + int dim = 0; + for(unsigned int i = 0; i < entities.size(); i++) + if(entities[i]->physicals.size() || saveAll) + for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++) + dim = std::max(dim, entities[i]->getMeshElement(j)->getDim()); + // loop over all elements we need to save int numElements = 0, maxNumNodesPerElement = 0; for(unsigned int i = 0; i < entities.size(); i++){ if(entities[i]->physicals.size() || saveAll){ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){ MElement *e = entities[i]->getMeshElement(j); - if(e->getStringForDIFF()){ + if(e->getStringForDIFF() && e->getDim() == dim){ numElements++; maxNumNodesPerElement = std::max(maxNumNodesPerElement, e->getNumVertices()); } @@ -2494,7 +2501,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll, if(entities[i]->physicals.size() || saveAll){ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){ MElement *e = entities[i]->getMeshElement(j); - if(e->getStringForDIFF()) + if(e->getStringForDIFF() && e->getDim() == dim) e->writeDIFF(fp, ++num, binary, entities[i]->tag()); } } diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index 3a55d93127..1bed85769b 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -127,15 +127,20 @@ static double F_Lc(GEdge *ge, double t) static double F_Transfinite(GEdge *ge, double t) { - double val, r; + double length = ge->length(); + if(length == 0.0){ + Msg::Error("Zero-length curve in transfinite mesh"); + return 1.; + } SVector3 der = ge->firstDer(t) ; double d = norm(der); - double coef = ge->meshAttributes.coeffTransfinite; int type = ge->meshAttributes.typeTransfinite; int nbpt = ge->meshAttributes.nbPointsTransfinite; + double val; + if(coef <= 0.0 || coef == 1.0) { // coef < 0 should never happen val = d * coef / ge->length(); @@ -145,12 +150,9 @@ static double F_Transfinite(GEdge *ge, double t) case 1: // Geometric progression ar^i; Sum of n terms = length = a (r^n-1)/(r-1) { - if(sign(type) >= 0) - r = coef; - else - r = 1. / coef; - double a = ge->length() * (r - 1.) / (pow(r, nbpt - 1.) - 1.); - int i = (int)(log(t * ge->length() / a * (r - 1.) + 1.) / log(r)); + double r = (sign(type) >= 0) ? coef : 1. / coef; + double a = length * (r - 1.) / (pow(r, nbpt - 1.) - 1.); + int i = (int)(log(t * length / a * (r - 1.) + 1.) / log(r)); val = d / (a * pow(r, (double)i)); } break; @@ -159,18 +161,16 @@ static double F_Transfinite(GEdge *ge, double t) { double a; if(coef > 1.0) { - a = -4. * sqrt(coef - 1.) * - atan2(1., sqrt(coef - 1.)) / - ((double)nbpt * ge->length()); + a = -4. * sqrt(coef - 1.) * atan2(1., sqrt(coef - 1.)) / + ((double)nbpt * length); } else { - a = 2. * sqrt(1. - coef) * - log(fabs((1. + 1. / sqrt(1. - coef)) - / (1. - 1. / sqrt(1. - coef)))) - / ((double)nbpt * ge->length()); + a = 2. * sqrt(1. - coef) * log(fabs((1. + 1. / sqrt(1. - coef)) / + (1. - 1. / sqrt(1. - coef)))) + / ((double)nbpt * length); } - double b = -a * ge->length() * ge->length() / (4. * (coef - 1.)); - val = d / (-a * SQU(t * ge->length() - (ge->length()) * 0.5) + b); + double b = -a * length * length / (4. * (coef - 1.)); + val = d / (-a * SQU(t * length - (length) * 0.5) + b); } break; @@ -288,8 +288,7 @@ void meshGEdge::operator() (GEdge *ge) Msg::Debug("Curve %d has a zero length", ge->tag()); // TEST - if (length < CTX.mesh.tolerance_edge_length)ge->setTooSmall(true); - + if (length < CTX.mesh.tolerance_edge_length) ge->setTooSmall(true); // Integrate detJ/lc du double a; -- GitLab