diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index 3fda9d52e4d72be2aa4dca8998a1eae1e0ec4ee3..2ff748e8a085dc876046601776decd20d4a77d14 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -331,14 +331,19 @@ void GFace::computeMeanPlane() if(pts.size() < 3){ Msg::Info("Adding edge points to compute mean plane of face %d", tag()); std::list<GEdge*> edg = edges(); - std::list<GEdge*>::const_iterator ite = edg.begin(); - for(; ite != edg.end(); ite++){ + for(std::list<GEdge*>::const_iterator ite = edg.begin(); ite != edg.end(); ite++){ const GEdge *e = *ite; - Range<double> b = e->parBounds(0); - GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low())); - pts.push_back(SPoint3(p1.x(), p1.y(), p1.z())); - GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low())); - pts.push_back(SPoint3(p2.x(), p2.y(), p2.z())); + if(e->mesh_vertices.size() > 1){ + for(unsigned int i = 0; i < e->mesh_vertices.size(); i++) + pts.push_back(e->mesh_vertices[i]->point()); + } + else{ + Range<double> b = e->parBounds(0); + GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low())); + pts.push_back(SPoint3(p1.x(), p1.y(), p1.z())); + GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low())); + pts.push_back(SPoint3(p2.x(), p2.y(), p2.z())); + } } } @@ -499,10 +504,10 @@ end: double d = meanPlane.a * v->x() + meanPlane.b * v->y() + meanPlane.c * v->z() - meanPlane.d; if(fabs(d) > lc * 1.e-3) { - Msg::Error("Plane surface %d (%gx+%gy+%gz=%g) is not plane!", - tag(), meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d); - Msg::Error("Control point %d = (%g,%g,%g), val=%g", - v->tag(), v->x(), v->y(), v->z(), d); + Msg::Warning("Plane surface %d (%gx+%gy+%gz=%g) is not plane!", + tag(), meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d); + Msg::Warning("Control point %d = (%g,%g,%g), val=%g", + v->tag(), v->x(), v->y(), v->z(), d); break; } } diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp index d33009beee0d3d0c4def05864019cabf066cd772..2863ab8a7be3993f1d1973aae84699d5fbb989d2 100644 --- a/Mesh/BoundaryLayers.cpp +++ b/Mesh/BoundaryLayers.cpp @@ -304,6 +304,13 @@ int Mesh2DWithBoundaryLayers(GModel *m) } } + // recompute mean plane for plane surfaces just in case + for(std::set<GFace*>::iterator it = otherFaces.begin(); it != otherFaces.end(); it++){ + GFace *gf = *it; + if(gf->geomType() == GEntity::Plane) + gf->computeMeanPlane(); + } + // mesh non-source surfaces std::for_each(otherFaces.begin(), otherFaces.end(), meshGFace()); diff --git a/benchmarks/extrude/aorta_extrude.geo b/benchmarks/extrude/aorta_extrude.geo index f210adc3a829d03e5da6db997ca1ce34331551ba..8661a35bbb5e9f7e9ca98bcb27d91c4e04816faf 100644 --- a/benchmarks/extrude/aorta_extrude.geo +++ b/benchmarks/extrude/aorta_extrude.geo @@ -3,5 +3,17 @@ CreateTopology; Merge "aortaRADIUS2.bgm"; +out1[] = Extrude{Surface{1}; Layers{4, 0.8}; Using Index[0]; Using View[0]; }; +out2[] = Extrude{Surface{1}; Layers{4, -0.5}; Using Index[1]; Using View[0]; }; + +Line Loop(60) = {9}; Plane Surface(61) = {60}; +Line Loop(62) = {7}; Plane Surface(63) = {62}; +Line Loop(64) = {10}; Plane Surface(65) = {64}; +Line Loop(66) = {11}; Plane Surface(67) = {66}; +Line Loop(68) = {8}; Plane Surface(69) = {68}; + +Mesh.Algorithm3D = 4; + +Surface Loop(100) = {32, 61:69:2}; +Volume(100) = 100; -out[] = Extrude{Surface{1}; Layers{4, 1.5}; Using Index[0]; Using View[0]; }; diff --git a/benchmarks/extrude/t1_boundary_layer.geo b/benchmarks/extrude/t1_boundary_layer.geo index 067ef892e74eab900d9be2e586c42bec8f737501..e8664b500023fda9fe57e07c3f9e0b222b96d768 100644 --- a/benchmarks/extrude/t1_boundary_layer.geo +++ b/benchmarks/extrude/t1_boundary_layer.geo @@ -25,4 +25,7 @@ Line Loop(7) = {5,6,7,1}; Plane Surface(8) = {7}; // the minus sign inverts the orientation of surface 8 -Extrude { Surface{6, -8}; Layers{5, 0.01}; Recombine; } +Extrude { + Surface{6, -8}; Layers{5, 0.01}; Recombine; + // Using View[-3]; // Hack to force "box-type" boundary layer along x,y,z axes +}