diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 7a95867f..9d2c0ddf 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -17,10 +17,34 @@ namespace netgen #define IGNORECURVELENGTH 1e-4 - extern OCCParameters occparam; - bool merge_solids = 1; + double Line :: Dist (Line l) + { + Vec<3> n = p1-p0; + Vec<3> q = l.p1-l.p0; + double nq = n*q; + + Point<3> p = p0 + 0.5*n; + double lambda = (p-l.p0)*n / nq; + + if (lambda >= 0 && lambda <= 1) + { + double d = (p-l.p0-lambda*q).Length(); + // if (d < 1e-3) d = 1e99; + return d; + } + else + return 1e99; + } + + + + double Line :: Length () + { + return (p1-p0).Length(); + } + inline Point<3> occ2ng (const gp_Pnt & p) @@ -30,6 +54,179 @@ namespace netgen + double ComputeH (double kappa) + { + double hret; + kappa *= mparam.curvaturesafety; + + if (mparam.maxh * kappa < 1) + hret = mparam.maxh; + else + hret = 1 / kappa; + + if (mparam.maxh < hret) + hret = mparam.maxh; + + return (hret); + } + + + + + void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, + BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0) + { + int ls = -1; + + gp_Pnt pnt0,pnt1,pnt2; + + prop->SetParameters (par0.X(), par0.Y()); + pnt0 = prop->Value(); + + prop->SetParameters (par1.X(), par1.Y()); + pnt1 = prop->Value(); + + prop->SetParameters (par2.X(), par2.Y()); + pnt2 = prop->Value(); + + double aux; + double maxside = pnt0.Distance(pnt1); + ls = 2; + aux = pnt1.Distance(pnt2); + if(aux > maxside) + { + maxside = aux; + ls = 0; + } + aux = pnt2.Distance(pnt0); + if(aux > maxside) + { + maxside = aux; + ls = 1; + } + + + + gp_Pnt2d parmid; + + parmid.SetX(0.3*(par0.X()+par1.X()+par2.X())); + parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y())); + + if (depth%3 == 0) + { + double curvature = 0; + + prop->SetParameters (parmid.X(), parmid.Y()); + if (!prop->IsCurvatureDefined()) + { + (*testout) << "curvature not defined!" << endl; + return; + } + curvature = max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature())); + + prop->SetParameters (par0.X(), par0.Y()); + if (!prop->IsCurvatureDefined()) + { + (*testout) << "curvature not defined!" << endl; + return; + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); + + prop->SetParameters (par1.X(), par1.Y()); + if (!prop->IsCurvatureDefined()) + { + (*testout) << "curvature not defined!" << endl; + return; + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); + + prop->SetParameters (par2.X(), par2.Y()); + if (!prop->IsCurvatureDefined()) + { + (*testout) << "curvature not defined!" << endl; + return; + } + curvature = max(curvature,max(fabs(prop->MinCurvature()), + fabs(prop->MaxCurvature()))); + + //(*testout) << "curvature " << curvature << endl; + + if (curvature < 1e-3) + { + //(*testout) << "curvature too small (" << curvature << ")!" << endl; + return; + // return war bis 10.2.05 auskommentiert + } + + + + h = ComputeH (curvature+1e-10); + + if(h < 1e-4*maxside) + return; + + + if (h > 30) return; + } + + if (h < maxside && depth < 10) + { + //cout << "\r h " << h << flush; + gp_Pnt2d pm; + + //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; + //cout << "par0 " << par0.X() << " " << par0.Y() + //<< " par1 " << par1.X() << " " << par1.Y() + // << " par2 " << par2.X() << " " << par2.Y()<< endl; + + if(ls == 0) + { + pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); + RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); + RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); + } + else if(ls == 1) + { + pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); + RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); + RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); + } + else if(ls == 2) + { + pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); + RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); + RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); + } + + } + else + { + gp_Pnt pnt; + Point3d p3d; + + prop->SetParameters (parmid.X(), parmid.Y()); + pnt = prop->Value(); + p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); + mesh.RestrictLocalH (p3d, h); + + p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); + mesh.RestrictLocalH (p3d, h); + + p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); + mesh.RestrictLocalH (p3d, h); + + p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); + mesh.RestrictLocalH (p3d, h); + + //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; + + } + } + + void DivideEdge (TopoDS_Edge & edge, Array & ps, Array & params, Mesh & mesh) @@ -53,7 +250,9 @@ namespace netgen double olddist = 0; double dist = 0; - for (int i = 1; i <= DIVIDEEDGESECTIONS; i++) + int tmpVal = (int)(DIVIDEEDGESECTIONS); + + for (int i = 1; i <= tmpVal; i++) { oldpnt = pnt; pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); @@ -64,7 +263,6 @@ namespace netgen //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; - olddist = dist; dist = pnt.Distance(oldpnt); } @@ -111,7 +309,7 @@ namespace netgen - static void FindEdges (OCCGeometry & geom, Mesh & mesh) + void OCCFindEdges (OCCGeometry & geom, Mesh & mesh) { const char * savetask = multithread.task; multithread.task = "Edge meshing"; @@ -180,8 +378,6 @@ namespace netgen total++; - - int facenr = 0; int edgenr = 0; @@ -272,8 +468,7 @@ namespace netgen Array params; DivideEdge (edge, mp, params, mesh); - - + Array pnums; pnums.SetSize (mp.Size()+2); @@ -404,7 +599,7 @@ namespace netgen - static void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend) + void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend) { int i, j, k; int changed; @@ -774,339 +969,300 @@ namespace netgen multithread.task = savetask; } - double ComputeH (double kappa) + + + void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh) { - double hret; - kappa *= mparam.curvaturesafety; + mesh.SetGlobalH (mparam.maxh); + mesh.SetMinimalH (mparam.minh); - if (mparam.maxh * kappa < 1) - hret = mparam.maxh; - else - hret = 1 / kappa; + Array maxhdom; + maxhdom.SetSize (geom.NrSolids()); + maxhdom = mparam.maxh; - if (mparam.maxh < hret) - hret = mparam.maxh; + mesh.SetMaxHDomain (maxhdom); - return (hret); - } + Box<3> bb = geom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); + mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); - class Line - { - public: - Point<3> p0, p1; - double Dist (Line l) + if (mparam.uselocalh) { - Vec<3> n = p1-p0; - Vec<3> q = l.p1-l.p0; - double nq = n*q; + const char * savetask = multithread.task; + multithread.percent = 0; - Point<3> p = p0 + 0.5*n; - double lambda = (p-l.p0)*n / nq; + mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); - if (lambda >= 0 && lambda <= 1) + int nedges = geom.emap.Extent(); + + double maxedgelen = 0; + double minedgelen = 1e99; + + multithread.task = "Setting local mesh size (elements per edge)"; + + // setting elements per edge + + for (int i = 1; i <= nedges && !multithread.terminate; i++) { - double d = (p-l.p0-lambda*q).Length(); - // if (d < 1e-3) d = 1e99; - return d; + TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); + multithread.percent = 100 * (i-1)/double(nedges); + if (BRep_Tool::Degenerated(e)) continue; + + GProp_GProps system; + BRepGProp::LinearProperties(e, system); + double len = system.Mass(); + + if (len < IGNORECURVELENGTH) + { + (*testout) << "ignored" << endl; + continue; + } + + double localh = len/mparam.segmentsperedge; + double s0, s1; + + // Philippose - 23/01/2009 + // Find all the parent faces of a given edge + // and limit the mesh size of the edge based on the + // mesh size limit of the face + TopTools_IndexedDataMapOfShapeListOfShape edge_face_map; + edge_face_map.Clear(); + + TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); + const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); + + TopTools_ListIteratorOfListOfShape parent_face_list; + + for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) + { + TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); + + int face_index = geom.fmap.FindIndex(parent_face); + + if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); + } + + Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); + + maxedgelen = max (maxedgelen, len); + minedgelen = min (minedgelen, len); + + // Philippose - 23/01/2009 + // Modified the calculation of maxj, because the + // method used so far always results in maxj = 2, + // which causes the localh to be set only at the + // starting, mid and end of the edge. + // Old Algorithm: + // int maxj = 2 * (int) ceil (localh/len); + int maxj = max((int) ceil(len/localh), 2); + + for (int j = 0; j <= maxj; j++) + { + gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); + mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); + } } - else - return 1e99; + + multithread.task = "Setting local mesh size (edge curvature)"; + + // setting edge curvature + + int nsections = 20; + + for (int i = 1; i <= nedges && !multithread.terminate; i++) + { + double maxcur = 0; + multithread.percent = 100 * (i-1)/double(nedges); + TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); + if (BRep_Tool::Degenerated(edge)) continue; + double s0, s1; + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + BRepAdaptor_Curve brepc(edge); + BRepLProp_CLProps prop(brepc, 2, 1e-5); + + for (int j = 1; j <= nsections; j++) + { + double s = s0 + j/(double) nsections * (s1-s0); + prop.SetParameter (s); + double curvature = prop.Curvature(); + if(curvature> maxcur) maxcur = curvature; + + if (curvature >= 1e99) + continue; + + gp_Pnt pnt = c->Value (s); + + mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature))); + } + // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; + } + + multithread.task = "Setting local mesh size (face curvature)"; + + // setting face curvature + + int nfaces = geom.fmap.Extent(); + + for (int i = 1; i <= nfaces && !multithread.terminate; i++) + { + multithread.percent = 100 * (i-1)/double(nfaces); + TopoDS_Face face = TopoDS::Face(geom.fmap(i)); + TopLoc_Location loc; + Handle(Geom_Surface) surf = BRep_Tool::Surface (face); + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); + + if (triangulation.IsNull()) continue; + + BRepAdaptor_Surface sf(face, Standard_True); + BRepLProp_SLProps prop(sf, 2, 1e-5); + + int ntriangles = triangulation -> NbTriangles(); + for (int j = 1; j <= ntriangles; j++) + { + gp_Pnt p[3]; + gp_Pnt2d par[3]; + + for (int k = 1; k <=3; k++) + { + int n = triangulation->Triangles()(j)(k); + p[k-1] = triangulation->Nodes()(n).Transformed(loc); + par[k-1] = triangulation->UVNodes()(n); + } + + //double maxside = 0; + //maxside = max (maxside, p[0].Distance(p[1])); + //maxside = max (maxside, p[0].Distance(p[2])); + //maxside = max (maxside, p[1].Distance(p[2])); + //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; + + RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0); + //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; + } + } + + // setting close edges + + if (occparam.resthcloseedgeenable) + { + multithread.task = "Setting local mesh size (close edges)"; + + int sections = 100; + + Array lines(sections*nedges); + + Box3dTree* searchtree = + new Box3dTree (bb.PMin(), bb.PMax()); + + int nlines = 0; + for (int i = 1; i <= nedges && !multithread.terminate; i++) + { + TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); + if (BRep_Tool::Degenerated(edge)) continue; + + double s0, s1; + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + BRepAdaptor_Curve brepc(edge); + BRepLProp_CLProps prop(brepc, 1, 1e-5); + prop.SetParameter (s0); + + gp_Vec d0 = prop.D1().Normalized(); + double s_start = s0; + int count = 0; + for (int j = 1; j <= sections; j++) + { + double s = s0 + (s1-s0)*(double)j/(double)sections; + prop.SetParameter (s); + gp_Vec d1 = prop.D1().Normalized(); + double cosalpha = fabs(d0*d1); + if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) + { + count++; + gp_Pnt p0 = c->Value (s_start); + gp_Pnt p1 = c->Value (s); + lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); + lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); + + Box3d box; + box.SetPoint (Point3d(lines[nlines].p0)); + box.AddPoint (Point3d(lines[nlines].p1)); + + searchtree->Insert (box.PMin(), box.PMax(), nlines+1); + nlines++; + + s_start = s; + d0 = d1; + } + } + } + + Array linenums; + + for (int i = 0; i < nlines; i++) + { + multithread.percent = (100*i)/double(nlines); + Line & line = lines[i]; + + Box3d box; + box.SetPoint (Point3d(line.p0)); + box.AddPoint (Point3d(line.p1)); + double maxhline = max (mesh.GetH(box.PMin()), + mesh.GetH(box.PMax())); + box.Increase(maxhline); + + double mindist = 1e99; + linenums.SetSize(0); + searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); + + for (int j = 0; j < linenums.Size(); j++) + { + int num = linenums[j]-1; + if (i == num) continue; + if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; + if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; + if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; + if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; + mindist = min (mindist, line.Dist(lines[num])); + } + + mindist *= occparam.resthcloseedgefac; + + if (mindist < 1e-3) + { + (*testout) << "extremely small local h: " << mindist + << " --> setting to 1e-3" << endl; + (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; + mindist = 1e-3; + } + + mesh.RestrictLocalHLine(line.p0, line.p1, mindist); + } + } + + multithread.task = savetask; + } - double Length () - { - return (p1-p0).Length(); - }; - }; - - /* - void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, - BRepLProp_SLProps * prop, Mesh & mesh, const double maxside, int depth, double h = 0) - { - - - gp_Pnt2d parmid; - - parmid.SetX(0.3*(par0.X()+par1.X()+par2.X())); - parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y())); - - if (depth == 0) - { - double curvature = 0; - - prop->SetParameters (parmid.X(), parmid.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature())); - - prop->SetParameters (par0.X(), par0.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - prop->SetParameters (par1.X(), par1.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - prop->SetParameters (par2.X(), par2.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - //(*testout) << "curvature " << curvature << endl; - - if (curvature < 1e-3) - { - //(*testout) << "curvature too small (" << curvature << ")!" << endl; - return; - // return war bis 10.2.05 auskommentiert + // Philippose - 09/03/2009 + // Added the capability to load the mesh size from a + // file also for OpenCascade Geometry + // Note: + // ** If the "uselocalh" option is ticked in + // the "mesh options...insider" menu, the mesh + // size will be further modified by the topology + // analysis routines. + // ** To use the mesh size file as the sole source + // for defining the mesh size, uncheck the "uselocalh" + // option. + mesh.LoadLocalMeshSize (mparam.meshsizefilename); } - h = ComputeH (curvature+1e-10); - - if(h < 1e-4*maxside) - return; - - - if (h > 30) return; - } - - if (h < maxside) // && depth < 5) + int OCCGenerateMesh (OCCGeometry & geom, Mesh *& mesh, + int perfstepsstart, int perfstepsend, + char * optstr) { - //cout << "\r h " << h << flush; - gp_Pnt2d pm0; - gp_Pnt2d pm1; - gp_Pnt2d pm2; - - //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; - //cout << "par0 " << par0.X() << " " << par0.Y() - //<< " par1 " << par1.X() << " " << par1.Y() - // << " par2 " << par2.X() << " " << par2.Y()<< endl; - - - pm0.SetX(0.5*(par1.X()+par2.X())); pm0.SetY(0.5*(par1.Y()+par2.Y())); - pm1.SetX(0.5*(par0.X()+par2.X())); pm1.SetY(0.5*(par0.Y()+par2.Y())); - pm2.SetX(0.5*(par1.X()+par0.X())); pm2.SetY(0.5*(par1.Y()+par0.Y())); - - RestrictHTriangle (pm0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h); - RestrictHTriangle (par0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h); - RestrictHTriangle (par1, pm0, pm2, prop, mesh, 0.5*maxside, depth+1, h); - RestrictHTriangle (par2, pm1, pm0, prop, mesh, 0.5*maxside, depth+1, h); - } - else - { - gp_Pnt pnt; - Point3d p3d; - - prop->SetParameters (parmid.X(), parmid.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); - - - prop->SetParameters (par0.X(), par0.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); - - prop->SetParameters (par1.X(), par1.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); - - prop->SetParameters (par2.X(), par2.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); - - } - } - */ - - - void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, - BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0) - { - int ls = -1; - - gp_Pnt pnt0,pnt1,pnt2; - - prop->SetParameters (par0.X(), par0.Y()); - pnt0 = prop->Value(); - - prop->SetParameters (par1.X(), par1.Y()); - pnt1 = prop->Value(); - - prop->SetParameters (par2.X(), par2.Y()); - pnt2 = prop->Value(); - - double aux; - double maxside = pnt0.Distance(pnt1); - ls = 2; - aux = pnt1.Distance(pnt2); - if(aux > maxside) - { - maxside = aux; - ls = 0; - } - aux = pnt2.Distance(pnt0); - if(aux > maxside) - { - maxside = aux; - ls = 1; - } - - - - gp_Pnt2d parmid; - - parmid.SetX(0.3*(par0.X()+par1.X()+par2.X())); - parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y())); - - if (depth%3 == 0) - { - double curvature = 0; - - prop->SetParameters (parmid.X(), parmid.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature())); - - prop->SetParameters (par0.X(), par0.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - prop->SetParameters (par1.X(), par1.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - prop->SetParameters (par2.X(), par2.Y()); - if (!prop->IsCurvatureDefined()) - { - (*testout) << "curvature not defined!" << endl; - return; - } - curvature = max(curvature,max(fabs(prop->MinCurvature()), - fabs(prop->MaxCurvature()))); - - //(*testout) << "curvature " << curvature << endl; - - if (curvature < 1e-3) - { - //(*testout) << "curvature too small (" << curvature << ")!" << endl; - return; - // return war bis 10.2.05 auskommentiert - } - - - - h = ComputeH (curvature+1e-10); - - if(h < 1e-4*maxside) - return; - - - if (h > 30) return; - } - - if (h < maxside && depth < 10) - { - //cout << "\r h " << h << flush; - gp_Pnt2d pm; - - //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; - //cout << "par0 " << par0.X() << " " << par0.Y() - //<< " par1 " << par1.X() << " " << par1.Y() - // << " par2 " << par2.X() << " " << par2.Y()<< endl; - - if(ls == 0) - { - pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); - RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); - } - else if(ls == 1) - { - pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); - RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h); - } - else if(ls == 2) - { - pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); - RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); - RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); - } - - } - else - { - gp_Pnt pnt; - Point3d p3d; - - prop->SetParameters (parmid.X(), parmid.Y()); - pnt = prop->Value(); - p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); - mesh.RestrictLocalH (p3d, h); - - p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); - mesh.RestrictLocalH (p3d, h); - - p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); - mesh.RestrictLocalH (p3d, h); - - p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); - mesh.RestrictLocalH (p3d, h); - - //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; - - } - } - - - - - int OCCGenerateMesh (OCCGeometry & geom, - Mesh *& mesh, - int perfstepsstart, int perfstepsend, - char * optstr) - { - // int i, j; - multithread.percent = 0; if (perfstepsstart <= MESHCONST_ANALYSE) @@ -1114,289 +1270,8 @@ namespace netgen delete mesh; mesh = new Mesh(); mesh->geomtype = Mesh::GEOM_OCC; - mesh->SetGlobalH (mparam.maxh); - mesh->SetMinimalH (mparam.minh); - Array maxhdom; - maxhdom.SetSize (geom.NrSolids()); - maxhdom = mparam.maxh; - - mesh->SetMaxHDomain (maxhdom); - - Box<3> bb = geom.GetBoundingBox(); - bb.Increase (bb.Diam()/10); - - mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5); - - if (mparam.uselocalh) - { - - const char * savetask = multithread.task; - multithread.percent = 0; - - mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); - - int nedges = geom.emap.Extent(); - - double maxedgelen = 0; - double minedgelen = 1e99; - - multithread.task = "Setting local mesh size (elements per edge)"; - - // setting elements per edge - - for (int i = 1; i <= nedges && !multithread.terminate; i++) - { - TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); - multithread.percent = 100 * (i-1)/double(nedges); - if (BRep_Tool::Degenerated(e)) continue; - - GProp_GProps system; - BRepGProp::LinearProperties(e, system); - double len = system.Mass(); - - if (len < IGNORECURVELENGTH) - { - (*testout) << "ignored" << endl; - continue; - } - - double localh = len/mparam.segmentsperedge; - double s0, s1; - - // Philippose - 23/01/2009 - // Find all the parent faces of a given edge - // and limit the mesh size of the edge based on the - // mesh size limit of the face - TopTools_IndexedDataMapOfShapeListOfShape edge_face_map; - edge_face_map.Clear(); - - TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); - const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); - - TopTools_ListIteratorOfListOfShape parent_face_list; - - for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) - { - TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); - - int face_index = geom.fmap.FindIndex(parent_face); - - if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); - } - - Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); - - maxedgelen = max (maxedgelen, len); - minedgelen = min (minedgelen, len); - - // Philippose - 23/01/2009 - // Modified the calculation of maxj, because the - // method used so far always results in maxj = 2, - // which causes the localh to be set only at the - // starting, mid and end of the edge. - // Old Algorithm: - // int maxj = 2 * (int) ceil (localh/len); - int maxj = max((int) ceil(len/localh), 2); - - for (int j = 0; j <= maxj; j++) - { - gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); - mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); - } - } - - multithread.task = "Setting local mesh size (edge curvature)"; - - // setting edge curvature - - int nsections = 20; - - for (int i = 1; i <= nedges && !multithread.terminate; i++) - { - double maxcur = 0; - multithread.percent = 100 * (i-1)/double(nedges); - TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); - if (BRep_Tool::Degenerated(edge)) continue; - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - BRepAdaptor_Curve brepc(edge); - BRepLProp_CLProps prop(brepc, 2, 1e-5); - - for (int j = 1; j <= nsections; j++) - { - double s = s0 + j/(double) nsections * (s1-s0); - prop.SetParameter (s); - double curvature = prop.Curvature(); - if(curvature> maxcur) maxcur = curvature; - - if (curvature >= 1e99) - continue; - - gp_Pnt pnt = c->Value (s); - - mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), - ComputeH (fabs(curvature))); - } - // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; - } - - multithread.task = "Setting local mesh size (face curvature)"; - - // setting face curvature - - int nfaces = geom.fmap.Extent(); - - for (int i = 1; i <= nfaces && !multithread.terminate; i++) - { - multithread.percent = 100 * (i-1)/double(nfaces); - TopoDS_Face face = TopoDS::Face(geom.fmap(i)); - TopLoc_Location loc; - Handle(Geom_Surface) surf = BRep_Tool::Surface (face); - Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); - if (triangulation.IsNull()) continue; - - BRepAdaptor_Surface sf(face, Standard_True); - BRepLProp_SLProps prop(sf, 2, 1e-5); - - int ntriangles = triangulation -> NbTriangles(); - for (int j = 1; j <= ntriangles; j++) - { - gp_Pnt p[3]; - gp_Pnt2d par[3]; - - for (int k = 1; k <=3; k++) - { - int n = triangulation->Triangles()(j)(k); - p[k-1] = triangulation->Nodes()(n).Transformed(loc); - par[k-1] = triangulation->UVNodes()(n); - } - - //double maxside = 0; - //maxside = max (maxside, p[0].Distance(p[1])); - //maxside = max (maxside, p[0].Distance(p[2])); - //maxside = max (maxside, p[1].Distance(p[2])); - //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; - - RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, 0); - //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; - } - } - - // setting close edges - - if (occparam.resthcloseedgeenable) - { - multithread.task = "Setting local mesh size (close edges)"; - - int sections = 100; - - Array lines(sections*nedges); - - Box3dTree* searchtree = - new Box3dTree (bb.PMin(), bb.PMax()); - - int nlines = 0; - for (int i = 1; i <= nedges && !multithread.terminate; i++) - { - TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); - if (BRep_Tool::Degenerated(edge)) continue; - - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - BRepAdaptor_Curve brepc(edge); - BRepLProp_CLProps prop(brepc, 1, 1e-5); - prop.SetParameter (s0); - - gp_Vec d0 = prop.D1().Normalized(); - double s_start = s0; - int count = 0; - for (int j = 1; j <= sections; j++) - { - double s = s0 + (s1-s0)*(double)j/(double)sections; - prop.SetParameter (s); - gp_Vec d1 = prop.D1().Normalized(); - double cosalpha = fabs(d0*d1); - if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) - { - count++; - gp_Pnt p0 = c->Value (s_start); - gp_Pnt p1 = c->Value (s); - lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); - lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); - - Box3d box; - box.SetPoint (Point3d(lines[nlines].p0)); - box.AddPoint (Point3d(lines[nlines].p1)); - - searchtree->Insert (box.PMin(), box.PMax(), nlines+1); - nlines++; - - s_start = s; - d0 = d1; - } - } - } - - Array linenums; - - for (int i = 0; i < nlines; i++) - { - multithread.percent = (100*i)/double(nlines); - Line & line = lines[i]; - - Box3d box; - box.SetPoint (Point3d(line.p0)); - box.AddPoint (Point3d(line.p1)); - double maxhline = max (mesh->GetH(box.PMin()), - mesh->GetH(box.PMax())); - box.Increase(maxhline); - - double mindist = 1e99; - linenums.SetSize(0); - searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); - - for (int j = 0; j < linenums.Size(); j++) - { - int num = linenums[j]-1; - if (i == num) continue; - if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; - if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; - if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; - if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; - mindist = min (mindist, line.Dist(lines[num])); - } - - mindist *= occparam.resthcloseedgefac; - - if (mindist < 1e-3) - { - (*testout) << "extremely small local h: " << mindist - << " --> setting to 1e-3" << endl; - (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; - mindist = 1e-3; - } - - mesh->RestrictLocalHLine(line.p0, line.p1, mindist); - } - } - - multithread.task = savetask; - - } - - // Philippose - 09/03/2009 - // Added the capability to load the mesh size from a - // file also for OpenCascade Geometry - // Note: - // ** If the "uselocalh" option is ticked in - // the "mesh options...insider" menu, the mesh - // size will be further modified by the topology - // analysis routines. - // ** To use the mesh size file as the sole source - // for defining the mesh size, uncheck the "uselocalh" - // option. - mesh->LoadLocalMeshSize (mparam.meshsizefilename); + OCCSetLocalMeshSize(geom,*mesh); } if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE) @@ -1404,7 +1279,7 @@ namespace netgen if (perfstepsstart <= MESHCONST_MESHEDGES) { - FindEdges (geom, *mesh); + OCCFindEdges (geom, *mesh); /* cout << "Removing redundant points" << endl; @@ -1503,8 +1378,7 @@ namespace netgen { multithread.task = "Volume meshing"; - MESHING3_RESULT res = - MeshVolume (mparam, *mesh); + MESHING3_RESULT res = MeshVolume (mparam, *mesh); ofstream problemfile("occmesh.rep",ios_base::app); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 6eb89933..d35782e8 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1603,8 +1603,6 @@ namespace netgen return false; } - extern int OCCGenerateMesh (OCCGeometry & occgeometry, Mesh*& mesh, - int perfstepsstart, int perfstepsend, char* optstring); int OCCGeometry :: GenerateMesh (Mesh*& mesh, int perfstepsstart, int perfstepsend, char* optstring) diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 2aa91f14..b43e80e0 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -80,6 +80,7 @@ #include "ShapeAnalysis.hxx" #include "ShapeBuild_ReShape.hxx" + // Philippose - 29/01/2009 // OpenCascade XDE Support // Include support for OpenCascade XDE Features @@ -164,6 +165,17 @@ namespace netgen + class Line + { + public: + Point<3> p0, p1; + + double Dist (Line l); + + double Length (); + }; + + inline double Det3 (double a00, double a01, double a02, double a10, double a11, double a12, @@ -386,6 +398,20 @@ namespace netgen OCCGeometry * LoadOCC_BREP (const char * filename); extern OCCParameters occparam; + + + // Philippose - 31.09.2009 + // External access to the mesh generation functions within the OCC + // subsystem (Not sure if this is the best way to implement this....!!) + extern int OCCGenerateMesh (OCCGeometry & occgeometry, Mesh*& mesh, + int perfstepsstart, int perfstepsend, + char* optstring); + + extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); + + extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend); + + extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh); } #endif diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 223f5c4d..9ed170b7 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -27,8 +27,8 @@ namespace netgen { extern void MeshFromSpline2D (SplineGeometry2d & geometry, - Mesh *& mesh, - MeshingParameters & mp); + Mesh *& mesh, + MeshingParameters & mp); } @@ -131,7 +131,7 @@ namespace nglib // Manually add a surface element of a given type to an existing mesh object DLL_HEADER void Ng_AddSurfaceElement (Ng_Mesh * mesh, Ng_Surface_Element_Type et, - int * pi) + int * pi) { Mesh * m = (Mesh*)mesh; Element2d el (3); @@ -146,7 +146,7 @@ namespace nglib // Manually add a volume element of a given type to an existing mesh object DLL_HEADER void Ng_AddVolumeElement (Ng_Mesh * mesh, Ng_Volume_Element_Type et, - int * pi) + int * pi) { Mesh * m = (Mesh*)mesh; Element el (4); @@ -270,8 +270,10 @@ namespace nglib { Mesh * m = (Mesh*)mesh; - - MeshingParameters mparam; + // Philippose - 30/08/2009 + // Do not locally re-define "mparam" here... "mparam" is a global + // object + //MeshingParameters mparam; mparam.maxh = mp->maxh; mparam.meshsizefilename = mp->meshsize_filename; @@ -369,9 +371,12 @@ namespace nglib return (Ng_Geometry_2D *)geom; } + + + DLL_HEADER Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom, - Ng_Mesh ** mesh, - Ng_Meshing_Parameters * mp) + Ng_Mesh ** mesh, + Ng_Meshing_Parameters * mp) { // use global variable mparam // MeshingParameters mparam; @@ -388,6 +393,9 @@ namespace nglib return NG_OK; } + + + DLL_HEADER void Ng_HP_Refinement (Ng_Geometry_2D * geom, Ng_Mesh * mesh, int levels) @@ -396,6 +404,9 @@ namespace nglib HPRefinement (*(Mesh*)mesh, &ref, levels); } + + + DLL_HEADER void Ng_HP_Refinement (Ng_Geometry_2D * geom, Ng_Mesh * mesh, int levels, double parameter) @@ -461,12 +472,16 @@ namespace nglib return geo2; } + + // generate new STL Geometry DLL_HEADER Ng_STL_Geometry * Ng_STL_NewGeometry () { return (Ng_STL_Geometry*)(void*)new STLGeometry; } + + // after adding triangles (and edges) initialize DLL_HEADER Ng_Result Ng_STL_InitSTLGeometry (Ng_STL_Geometry * geom) { @@ -489,10 +504,12 @@ namespace nglib return NG_SURFACE_INPUT_ERROR; } + + // automatically generates edges: DLL_HEADER Ng_Result Ng_STL_MakeEdges (Ng_STL_Geometry * geom, - Ng_Mesh* mesh, - Ng_Meshing_Parameters * mp) + Ng_Mesh* mesh, + Ng_Meshing_Parameters * mp) { STLGeometry* stlgeometry = (STLGeometry*)geom; Mesh* me = (Mesh*)mesh; @@ -507,8 +524,8 @@ namespace nglib me -> SetGlobalH (mparam.maxh); me -> SetLocalH (stlgeometry->GetBoundingBox().PMin() - Vec3d(10, 10, 10), - stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10), - 0.3); + stlgeometry->GetBoundingBox().PMax() + Vec3d(10, 10, 10), + 0.3); me -> LoadLocalMeshSize (mp->meshsize_filename); /* @@ -533,8 +550,8 @@ namespace nglib // generates mesh, empty mesh be already created. DLL_HEADER Ng_Result Ng_STL_GenerateSurfaceMesh (Ng_STL_Geometry * geom, - Ng_Mesh* mesh, - Ng_Meshing_Parameters * mp) + Ng_Mesh* mesh, + Ng_Meshing_Parameters * mp) { STLGeometry* stlgeometry = (STLGeometry*)geom; Mesh* me = (Mesh*)mesh; @@ -593,7 +610,8 @@ namespace nglib // positive orientation // normal vector may be null-pointer DLL_HEADER void Ng_STL_AddTriangle (Ng_STL_Geometry * geom, - double * p1, double * p2, double * p3, double * nv) + double * p1, double * p2, double * p3, + double * nv) { Point<3> apts[3]; apts[0] = Point<3>(p1[0],p1[1],p1[2]); @@ -621,7 +639,6 @@ namespace nglib - #ifdef OCCGEOMETRY // --------------------- OCC Geometry / Meshing Utility Functions ------------------- @@ -632,46 +649,166 @@ namespace nglib } + // Loads geometry from STEP File DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_STEP (const char * filename) { Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry(); - geo = (Ng_OCC_Geometry *)LoadOCC_STEP(filename); + OCCGeometry * occgeo = LoadOCC_STEP(filename); + + // Create the initial triangulation for the OCC + BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1); + + geo = (Ng_OCC_Geometry *)occgeo; return (geo); } + // Loads geometry from IGES File DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_IGES (const char * filename) { Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry(); - geo = (Ng_OCC_Geometry *)LoadOCC_IGES(filename); + OCCGeometry * occgeo = LoadOCC_IGES(filename); + + // Create the initial triangulation for the OCC + BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1); + + geo = (Ng_OCC_Geometry *)occgeo; return (geo); } + // Loads geometry from BREP File DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_BREP (const char * filename) { Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry(); - geo = (Ng_OCC_Geometry *)LoadOCC_BREP(filename); + OCCGeometry * occgeo = LoadOCC_BREP(filename); + + // Create the initial triangulation for the OCC + BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1); + + geo = (Ng_OCC_Geometry *)occgeo; return (geo); } - // Extract a + + // + DLL_HEADER Ng_Result Ng_OCC_SetLocalMeshSize (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp) + { + OCCGeometry * occgeom = (OCCGeometry*)geom; + Mesh * me = (Mesh*)mesh; + + mparam.uselocalh = mp->uselocalh; + + mparam.maxh = mp->maxh; + mparam.minh = mp->minh; + + mparam.segmentsperedge = mp->elementsperedge; + mparam.curvaturesafety = mp->elementspercurve; + + mparam.grading = mp->grading; + mparam.meshsizefilename = mp->meshsize_filename; + + occparam.resthcloseedgeenable = mp->closeedgeenable; + occparam.resthcloseedgefac = mp->closeedgefact; + + delete me; + me = new Mesh; + + OCCSetLocalMeshSize(*occgeom, *me); + + return(NG_OK); + } + + + + // Mesh the edges and add Face descriptors to prepare for surface meshing + DLL_HEADER Ng_Result Ng_OCC_GenerateEdgeMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp) + { + OCCGeometry * occgeom = (OCCGeometry*)geom; + Mesh * me = (Mesh*)mesh; + + mparam.uselocalh = mp->uselocalh; + + OCCFindEdges(*occgeom, *me); + + if((me->GetNP()) && (me->GetNFD())) + { + return NG_OK; + } + else + { + return NG_ERROR; + } + } + + + + + // Mesh the edges and add Face descriptors to prepare for surface meshing + DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp) + { + int numpoints = 0; + + OCCGeometry * occgeom = (OCCGeometry*)geom; + Mesh * me = (Mesh*)mesh; + + mparam.uselocalh = mp->uselocalh; + + // Only go into surface meshing if the face descriptors have already been added + if(!me->GetNFD()) + return NG_ERROR; + + numpoints = me->GetNP(); + + // Initially set up only for surface meshing without any optimisation + int perfstepsend = MESHCONST_MESHSURFACE; + + // Check and if required, enable surface mesh optimisation step + if(mp->optsurfmeshenable) + { + perfstepsend = MESHCONST_OPTSURFACE; + } + + OCCMeshSurface(*occgeom, *me, perfstepsend); + + me->CalcSurfacesOfNode(); + + if(me->GetNP() <= numpoints) + return NG_ERROR; + + if(me->GetNSE() <= 0) + return NG_ERROR; + + return NG_OK; + } + + + // Extract the face map from the OCC geometry + // The face map basically gives an index to each face in the geometry, + // which can be used to access a specific face DLL_HEADER Ng_Result Ng_OCC_GetFMap(Ng_OCC_Geometry * geom, - Ng_OCC_TopTools_IndexedMapOfShape * FMap) + Ng_OCC_TopTools_IndexedMapOfShape * FMap) { OCCGeometry* occgeom = (OCCGeometry*)geom; TopTools_IndexedMapOfShape *occfmap = (TopTools_IndexedMapOfShape *)FMap; + // Copy the face map from the geometry to the given variable occfmap->Assign(occgeom->fmap); if(occfmap->Extent()) @@ -694,13 +831,28 @@ namespace nglib DLL_HEADER Ng_Meshing_Parameters :: Ng_Meshing_Parameters() { + uselocalh = 1; + maxh = 1000; + minh = 0.0; + fineness = 0.5; grading = 0.3; + + elementsperedge = 2.0; + elementspercurve = 2.0; + + closeedgeenable = 0; + closeedgefact = 2.0; + secondorder = 0; - meshsize_filename = 0; quad_dominated = 0; + meshsize_filename = 0; + + optsurfmeshenable = 1; + optvolmeshenable = 1; + optsteps_2d = 3; optsteps_3d = 3; } @@ -709,13 +861,28 @@ namespace nglib DLL_HEADER void Ng_Meshing_Parameters :: Reset_Parameters() { + uselocalh = 1; + maxh = 1000; + minh = 0; + fineness = 0.5; grading = 0.3; + + elementsperedge = 2.0; + elementspercurve = 2.0; + + closeedgeenable = 0; + closeedgefact = 2.0; + secondorder = 0; - meshsize_filename = 0; quad_dominated = 0; + meshsize_filename = 0; + + optsurfmeshenable = 1; + optvolmeshenable = 1; + optsteps_2d = 3; optsteps_3d = 3; } @@ -723,6 +890,8 @@ namespace nglib } // End of namespace nglib + + // compatibility functions: namespace netgen @@ -762,3 +931,4 @@ namespace netgen } } // End of namespace netgen + diff --git a/nglib/nglib.h b/nglib/nglib.h index 2c851a49..289ea595 100644 --- a/nglib/nglib.h +++ b/nglib/nglib.h @@ -92,12 +92,28 @@ enum Ng_Result class Ng_Meshing_Parameters { public: - double maxh; //!< Maximum global mesh size limit + int uselocalh; //!< Switch to enable / disable usage of local mesh size modifiers + + double maxh; //!< Maximum global mesh size allowed + double minh; //!< Minimum global mesh size allowed + double fineness; //!< Mesh density: 0...1 (0 => coarse; 1 => fine) double grading; //!< Mesh grading: 0...1 (0 => uniform mesh; 1 => aggressive local grading) + + double elementsperedge; //!< Number of elements to generate per edge of the geometry + double elementspercurve; //!< Elements to generate per curvature radius + + int closeedgeenable; //!< Enable / Disable mesh refinement at close edges + double closeedgefact; //!< Factor to use for refinement at close edges (STL: larger => finer ; OCC: larger => coarser) + int secondorder; //!< Generate second-order surface and volume elements - char * meshsize_filename; //!< Optional external mesh size file int quad_dominated; //!< Creates a Quad-dominated mesh + + char * meshsize_filename; //!< Optional external mesh size file + + int optsurfmeshenable; //!< Enable / Disable automatic surface mesh optimization + int optvolmeshenable; //!< Enable / Disable automatic volume mesh optimization + int optsteps_3d; //!< Number of optimize steps to use for 3-D mesh optimization int optsteps_2d; //!< Number of optimize steps to use for 2-D mesh optimization @@ -107,12 +123,19 @@ public: Note: This constructor initialises the variables in the class with the following default values + - #uselocalh: 1 - #maxh: 1000.0 - #fineness: 0.5 - #grading: 0.3 + - #elementsperedge: 2.0 + - #elementspercurve: 2.0 + - #closeedgeenable: 0 + - #closeedgefact: 2.0 - #secondorder: 0.0 - #meshsize_filename: null - #quad_dominated: 0 + - #optsurfmeshenable: 1 + - #optvolmeshenable: 1 - #optsteps_2d: 3 - #optsteps_3d: 3 */ @@ -543,8 +566,8 @@ DLL_HEADER Ng_Result Ng_STL_MakeEdges (Ng_STL_Geometry * geom, // generates mesh, empty mesh must be already created. DLL_HEADER Ng_Result Ng_STL_GenerateSurfaceMesh (Ng_STL_Geometry * geom, - Ng_Mesh * mesh, - Ng_Meshing_Parameters * mp); + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); #ifdef ACIS @@ -587,10 +610,25 @@ DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_IGES (const char * filename); // Loads geometry from BREP file DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_BREP (const char * filename); +// Set the local mesh size based on geometry / topology +DLL_HEADER Ng_Result Ng_OCC_SetLocalMeshSize (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + +// Mesh the edges and add Face descriptors to prepare for surface meshing +DLL_HEADER Ng_Result Ng_OCC_GenerateEdgeMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + +// Mesh the surfaces of an OCC geometry +DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom, + Ng_Mesh * mesh, + Ng_Meshing_Parameters * mp); + // Get the face map of an already loaded OCC geometry DLL_HEADER Ng_Result Ng_OCC_GetFMap(Ng_OCC_Geometry * geom, Ng_OCC_TopTools_IndexedMapOfShape * FMap); -#endif +#endif // OCCGEOMETRY -#endif +#endif // NGLIB