diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 2817ba18..ed63cf8c 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -27,7 +27,7 @@ namespace netgen { return Point<3> (p.X(), p.Y(), p.Z()); } - + void DivideEdge (TopoDS_Edge & edge, Array & ps, Array & params, @@ -35,7 +35,7 @@ namespace netgen { double s0, s1; double maxh = mparam.maxh; - int nsubedges = 1; + int nsubedges = 1; gp_Pnt pnt, oldpnt; double svalue[DIVIDEEDGESECTIONS]; @@ -60,7 +60,7 @@ namespace netgen 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* pnt.Distance(oldpnt); - //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) + //(*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; @@ -106,7 +106,7 @@ namespace netgen params[nsubedges] = s1; } } - + @@ -138,13 +138,13 @@ namespace netgen exists = 1; break; } - + if (!exists) mesh.AddPoint (mp); } (*testout) << "different vertices = " << mesh.GetNP() << endl; - + int first_ep = mesh.GetNP()+1; @@ -177,17 +177,17 @@ namespace netgen for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) total++; - - + + int facenr = 0; int edgenr = 0; - + (*testout) << "faces = " << geom.fmap.Extent() << endl; int curr = 0; - + for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) { TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); @@ -207,16 +207,16 @@ namespace netgen } */ - mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); + mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) Handle(Geom_Surface) occface = BRep_Tool::Surface(face); - + for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) { TopoDS_Shape wire = exp2.Current(); - + for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next()) { curr++; @@ -226,13 +226,13 @@ namespace netgen if (multithread.terminate) return; TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); - if (BRep_Tool::Degenerated(edge)) + if (BRep_Tool::Degenerated(edge)) { //(*testout) << "ignoring degenerated edge" << endl; continue; } - if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == + if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == geom.vmap.FindIndex(TopExp::LastVertex (edge))) { GProp_GProps system; @@ -257,11 +257,11 @@ namespace netgen Array params; DivideEdge (edge, mp, params, mesh); - - + + Array pnums; pnums.SetSize (mp.Size()+2); - + if (!merge_solids) { pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)); @@ -271,7 +271,7 @@ namespace netgen { Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); - + pnums[0] = -1; pnums.Last() = -1; for (PointIndex pi = 1; pi < first_ep; pi++) @@ -280,14 +280,14 @@ namespace netgen if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; } } - + for (int i = 1; i <= mp.Size(); i++) { bool exists = 0; int j; for (j = first_ep; j <= mesh.GetNP(); j++) - if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) + if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) { exists = 1; break; @@ -303,14 +303,14 @@ namespace netgen } } (*testout) << "NP = " << mesh.GetNP() << endl; - + //(*testout) << pnums[pnums.Size()-1] << endl; - + for (int i = 1; i <= mp.Size()+1; i++) { edgenr++; Segment seg; - + seg.p1 = pnums[i-1]; seg.p2 = pnums[i]; seg.edgenr = edgenr; @@ -319,7 +319,7 @@ namespace netgen seg.epgeominfo[1].dist = params[i]; seg.epgeominfo[0].edgenr = geomedgenr; seg.epgeominfo[1].edgenr = geomedgenr; - + gp_Pnt2d p2d; p2d = cof->Value(params[i-1]); // if (i == 1) p2d = cof->Value(s0); @@ -359,7 +359,7 @@ namespace netgen seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); } */ - + if (edge.Orientation() == TopAbs_REVERSED) { swap (seg.p1, seg.p2); @@ -378,13 +378,13 @@ namespace netgen } // for(i=1; i<=mesh.GetNSeg(); i++) - // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si + // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si // << " p1 " << mesh.LineSegment(i).p1 << " p2 " << mesh.LineSegment(i).p2 << endl; // exit(10); mesh.CalcSurfacesOfNode(); multithread.task = savetask; - } + } @@ -396,7 +396,7 @@ namespace netgen const char * savetask = multithread.task; multithread.task = "Surface meshing"; - + geom.facemeshstatus = 0; int noldp = mesh.GetNP(); @@ -406,7 +406,7 @@ namespace netgen Array glob2loc(noldp); //int projecttype = PARAMETERSPACE; - + int projecttype = PARAMETERSPACE; int notrys = 1; @@ -426,16 +426,16 @@ namespace netgen multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); geom.facemeshstatus[k-1] = -1; - - /* + + /* if (k != 42) { cout << "skipped" << endl; continue; } */ - - + + FaceDescriptor & fd = mesh.GetFaceDescriptor(k); int oldnf = mesh.GetNSE(); @@ -443,8 +443,8 @@ namespace netgen Box<3> bb = geom.GetBoundingBox(); // int projecttype = PLANESPACE; - - Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype); + + Meshing2OCCSurfaces meshing(TopoDS::Face(geom.fmap(k)), bb, projecttype); if (meshing.GetProjectionType() == PLANESPACE) PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); @@ -454,7 +454,7 @@ namespace netgen if (surfmesherror) cout << "Surface meshing error occured before (in " << surfmesherror << " faces)" << endl; - // Meshing2OCCSurfaces meshing(f2, bb); + // Meshing2OCCSurfaces meshing(f2, bb); meshing.SetStartTime (starttime); //(*testout) << "Face " << k << endl << endl; @@ -481,7 +481,7 @@ namespace netgen } } } - + for (i = 1; i <= mesh.GetNSeg(); i++) { Segment & seg = mesh.LineSegment(i); @@ -493,7 +493,7 @@ namespace netgen gi0.v = seg.epgeominfo[0].v; gi1.u = seg.epgeominfo[1].u; gi1.v = seg.epgeominfo[1].v; - + meshing.AddBoundaryElement (glob2loc.Get(seg.p1), glob2loc.Get(seg.p2), gi0, gi1); //(*testout) << gi0.u << " " << gi0.v << endl; //(*testout) << gi1.u << " " << gi1.v << endl; @@ -503,17 +503,17 @@ namespace netgen else { int cntp = 0; - + for (i = 1; i <= mesh.GetNSeg(); i++) if (mesh.LineSegment(i).si == k) cntp+=2; - - + + Array< PointGeomInfo > gis; - + gis.SetAllocSize (cntp); gis.SetSize (0); - + for (i = 1; i <= mesh.GetNSeg(); i++) { Segment & seg = mesh.LineSegment(i); @@ -525,13 +525,13 @@ namespace netgen gi0.v = seg.epgeominfo[0].v; gi1.u = seg.epgeominfo[1].u; gi1.v = seg.epgeominfo[1].v; - + int locpnum[2] = {0, 0}; - + for (j = 0; j < 2; j++) { PointGeomInfo gi = (j == 0) ? gi0 : gi1; - + int l; for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) { @@ -540,7 +540,7 @@ namespace netgen if (dist < 1e-10) locpnum[j] = l+1; } - + if (locpnum[j] == 0) { int pi = (j == 0) ? seg.p1 : seg.p2; @@ -551,11 +551,11 @@ namespace netgen locpnum[j] = l+1; } } - + meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); //(*testout) << gi0.u << " " << gi0.v << endl; //(*testout) << gi1.u << " " << gi1.v << endl; - + } } } @@ -564,8 +564,9 @@ namespace netgen - - double maxh = mparam.maxh; + // Philippose - 15/01/2009 + double maxh = geom.face_maxh[k-1]; + //double maxh = mparam.maxh; mparam.checkoverlap = 0; // int noldpoints = mesh->GetNP(); int noldsurfel = mesh.GetNSE(); @@ -600,7 +601,7 @@ namespace netgen { for (int i = noldsurfel+1; i <= mesh.GetNSE(); i++) mesh.DeleteSurfaceElement (i); - + mesh.Compress(); cout << "retry Surface " << k << endl; @@ -624,7 +625,7 @@ namespace netgen } notrys = 1; - + for (i = oldnf+1; i <= mesh.GetNSE(); i++) mesh.SurfaceElement(i).SetIndex (k); @@ -657,10 +658,10 @@ namespace netgen gp_Pnt point = BRep_Tool::Pnt(vertex); problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; } - } + } } problemfile << endl; - + } cout << endl << endl; cout << "for more information open IGES/STEP Topology Explorer" << endl; @@ -672,8 +673,8 @@ namespace netgen problemfile << "OK" << endl << endl; problemfile.close(); } - - + + if (multithread.terminate || perfstepsend < MESHCONST_OPTSURFACE) @@ -683,7 +684,7 @@ namespace netgen static int timer_opt2d = NgProfiler::CreateTimer ("Optimization 2D"); NgProfiler::StartTimer (timer_opt2d); - + for (k = 1; k <= mesh.GetNFD(); k++) { // if (k != 42) continue; @@ -691,26 +692,26 @@ namespace netgen // (*testout) << "optimize face " << k << endl; multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); - + FaceDescriptor & fd = mesh.GetFaceDescriptor(k); - + PrintMessage (1, "Optimize Surface ", k); for (i = 1; i <= mparam.optsteps2d; i++) { // (*testout) << "optstep " << i << endl; if (multithread.terminate) return; - + { MeshOptimize2dOCCSurfaces meshopt(geom); meshopt.SetFaceIndex (k); meshopt.SetImproveEdges (0); meshopt.SetMetricWeight (0.2); meshopt.SetWriteStatus (0); - + // (*testout) << "EdgeSwapping (mesh, (i > mparam.optsteps2d/2))" << endl; meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2)); } - + if (multithread.terminate) return; { MeshOptimize2dOCCSurfaces meshopt(geom); @@ -718,22 +719,22 @@ namespace netgen meshopt.SetImproveEdges (0); meshopt.SetMetricWeight (0.2); meshopt.SetWriteStatus (0); - + // (*testout) << "ImproveMesh (mesh)" << endl; meshopt.ImproveMesh (mesh); } - + { MeshOptimize2dOCCSurfaces meshopt(geom); meshopt.SetFaceIndex (k); meshopt.SetImproveEdges (0); meshopt.SetMetricWeight (0.2); meshopt.SetWriteStatus (0); - + // (*testout) << "CombineImprove (mesh)" << endl; meshopt.CombineImprove (mesh); } - + if (multithread.terminate) return; { MeshOptimize2dOCCSurfaces meshopt(geom); @@ -741,15 +742,15 @@ namespace netgen meshopt.SetImproveEdges (0); meshopt.SetMetricWeight (0.2); meshopt.SetWriteStatus (0); - + // (*testout) << "ImproveMesh (mesh)" << endl; meshopt.ImproveMesh (mesh); } } - + } - - + + mesh.CalcSurfacesOfNode(); mesh.Compress(); @@ -762,15 +763,15 @@ namespace netgen { 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); } @@ -784,7 +785,7 @@ namespace netgen 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; @@ -797,7 +798,7 @@ namespace netgen else return 1e99; } - + double Length () { return (p1-p0).Length(); @@ -814,7 +815,7 @@ namespace netgen 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; @@ -867,14 +868,14 @@ namespace netgen h = ComputeH (curvature+1e-10); - + if(h < 1e-4*maxside) return; if (h > 30) return; } - + if (h < maxside) // && depth < 5) { //cout << "\r h " << h << flush; @@ -887,7 +888,7 @@ namespace netgen //<< " 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())); @@ -907,7 +908,7 @@ namespace netgen 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()); @@ -927,23 +928,23 @@ namespace netgen } */ - + 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; @@ -959,14 +960,14 @@ namespace netgen 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; @@ -1019,14 +1020,14 @@ namespace netgen 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; @@ -1055,7 +1056,7 @@ namespace netgen RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h); RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h); } - + } else { @@ -1085,456 +1086,465 @@ namespace netgen int OCCGenerateMesh (OCCGeometry & geom, - Mesh *& mesh, - int perfstepsstart, int perfstepsend, - char * optstr) - { - // int i, j; + Mesh *& mesh, + int perfstepsstart, int perfstepsend, + char * optstr) + { + // int i, j; - multithread.percent = 0; + multithread.percent = 0; - if (perfstepsstart <= MESHCONST_ANALYSE) + if (perfstepsstart <= MESHCONST_ANALYSE) { - delete mesh; - mesh = new Mesh(); - mesh->geomtype = Mesh::GEOM_OCC; - mesh->SetGlobalH (mparam.maxh); - mesh->SetMinimalH (mparam.minh); + 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); - + Array maxhdom; + maxhdom.SetSize (geom.NrSolids()); + maxhdom = mparam.maxh; - if (mparam.uselocalh) - { + mesh->SetMaxHDomain (maxhdom); - const char * savetask = multithread.task; - multithread.percent = 0; + Box<3> bb = geom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); - mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); + mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5); - int nedges = geom.emap.Extent(); + if (mparam.uselocalh) + { - double maxedgelen = 0; - double minedgelen = 1e99; + const char * savetask = multithread.task; + multithread.percent = 0; - - multithread.task = "Setting local mesh size (elements per edge)"; + mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); - // 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; + int nedges = geom.emap.Extent(); - GProp_GProps system; - BRepGProp::LinearProperties(e, system); - double len = system.Mass(); + double maxedgelen = 0; + double minedgelen = 1e99; - if (len < IGNORECURVELENGTH) - { - (*testout) << "ignored" << endl; - continue; - } + multithread.task = "Setting local mesh size (elements per edge)"; + // setting elements per edge - double localh = len/mparam.segmentsperedge; - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); + 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; - maxedgelen = max (maxedgelen, len); - minedgelen = min (minedgelen, len); + GProp_GProps system; + BRepGProp::LinearProperties(e, system); + double len = system.Mass(); - int maxj = 2 * (int) ceil (localh/len); - 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); - } - } + 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(); - - multithread.task = "Setting local mesh size (edge curvature)"; - - - // setting edge curvature + TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); + const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); - int nsections = 20; + TopTools_ListIteratorOfListOfShape parent_face_list; - 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; + 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()); - if (curvature >= 1e99) - continue; - + int face_index = geom.fmap.FindIndex(parent_face); - gp_Pnt pnt = c->Value (s); + if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); + } - mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), - ComputeH (fabs(curvature))); - } - // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; - } + Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); - - - multithread.task = "Setting local mesh size (face curvature)"; + maxedgelen = max (maxedgelen, len); + minedgelen = min (minedgelen, len); - // 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; + // 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); - RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, 0); - //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; - } - } - + 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); + } + } - // setting close edges + multithread.task = "Setting local mesh size (edge curvature)"; - if (stlparam.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; + // setting edge curvature - 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); + int nsections = 20; - 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)); + 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); - searchtree->Insert (box.PMin(), box.PMax(), nlines+1); - nlines++; + 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; - 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); + if (curvature >= 1e99) + continue; - double mindist = 1e99; - linenums.SetSize(0); - searchtree->GetIntersecting(box.PMin(),box.PMax(),linenums); + gp_Pnt pnt = c->Value (s); - 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 *= stlparam.resthcloseedgefac; + mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), + ComputeH (fabs(curvature))); + } + // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl; + } - 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; - } + multithread.task = "Setting local mesh size (face curvature)"; - mesh->RestrictLocalHLine(line.p0, line.p1, mindist); - } - } - + // setting face curvature - multithread.task = savetask; + 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 (stlparam.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 *= stlparam.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; + + } } - - if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE) + if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE) return TCL_OK; - if (perfstepsstart <= MESHCONST_MESHEDGES) + if (perfstepsstart <= MESHCONST_MESHEDGES) { - FindEdges (geom, *mesh); + FindEdges (geom, *mesh); - /* - cout << "Removing redundant points" << endl; - - int i, j; - int np = mesh->GetNP(); - Array equalto; + /* + cout << "Removing redundant points" << endl; - equalto.SetSize (np); - equalto = 0; + int i, j; + int np = mesh->GetNP(); + Array equalto; - for (i = 1; i <= np; i++) - { - for (j = i+1; j <= np; j++) - { - if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12)) - equalto[j-1] = i; - } - } + equalto.SetSize (np); + equalto = 0; - for (i = 1; i <= np; i++) - if (equalto[i-1]) - { - cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl; - for (j = 1; j <= mesh->GetNSeg(); j++) - { - Segment & seg = mesh->LineSegment(j); - if (seg.p1 == i) seg.p1 = equalto[i-1]; - if (seg.p2 == i) seg.p2 = equalto[i-1]; - } - } + for (i = 1; i <= np; i++) + { + for (j = i+1; j <= np; j++) + { + if (!equalto[j-1] && (Dist2 (mesh->Point(i), mesh->Point(j)) < 1e-12)) + equalto[j-1] = i; + } + } - cout << "Removing degenerated segments" << endl; - for (j = 1; j <= mesh->GetNSeg(); j++) - { - Segment & seg = mesh->LineSegment(j); - if (seg.p1 == seg.p2) - { - mesh->DeleteSegment(j); - cout << "Deleting Segment " << j << endl; - } - } + for (i = 1; i <= np; i++) + if (equalto[i-1]) + { + cout << "Point " << i << " is equal to Point " << equalto[i-1] << endl; + for (j = 1; j <= mesh->GetNSeg(); j++) + { + Segment & seg = mesh->LineSegment(j); + if (seg.p1 == i) seg.p1 = equalto[i-1]; + if (seg.p2 == i) seg.p2 = equalto[i-1]; + } + } - mesh->Compress(); - */ + cout << "Removing degenerated segments" << endl; + for (j = 1; j <= mesh->GetNSeg(); j++) + { + Segment & seg = mesh->LineSegment(j); + if (seg.p1 == seg.p2) + { + mesh->DeleteSegment(j); + cout << "Deleting Segment " << j << endl; + } + } - /* - for (int i = 1; i <= geom.fmap.Extent(); i++) - { - Handle(Geom_Surface) hf1 = - BRep_Tool::Surface(TopoDS::Face(geom.fmap(i))); - for (int j = i+1; j <= geom.fmap.Extent(); j++) - { - Handle(Geom_Surface) hf2 = - BRep_Tool::Surface(TopoDS::Face(geom.fmap(j))); - if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl; - } - } - */ + mesh->Compress(); + */ + /* + for (int i = 1; i <= geom.fmap.Extent(); i++) + { + Handle(Geom_Surface) hf1 = + BRep_Tool::Surface(TopoDS::Face(geom.fmap(i))); + for (int j = i+1; j <= geom.fmap.Extent(); j++) + { + Handle(Geom_Surface) hf2 = + BRep_Tool::Surface(TopoDS::Face(geom.fmap(j))); + if (hf1 == hf2) cout << "face " << i << " and face " << j << " lie on same surface" << endl; + } + } + */ - -#ifdef LOG_STREAM - (*logout) << "Edges meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - } - - if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES) - return TCL_OK; - - if (perfstepsstart <= MESHCONST_MESHSURFACE) - { - OCCMeshSurface (geom, *mesh, perfstepsend); - if (multithread.terminate) return TCL_OK; - #ifdef LOG_STREAM - (*logout) << "Surfaces meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; -#endif - -#ifdef STAT_STREAM - (*statout) << mesh->GetNSeg() << " & " - << mesh->GetNSE() << " & - &" - << GetTime() << " & " << endl; -#endif - - // MeshQuality2d (*mesh); - mesh->CalcSurfacesOfNode(); + (*logout) << "Edges meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif } - if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE) + if (multithread.terminate || perfstepsend <= MESHCONST_MESHEDGES) return TCL_OK; - - - if (perfstepsstart <= MESHCONST_MESHVOLUME) + if (perfstepsstart <= MESHCONST_MESHSURFACE) { - multithread.task = "Volume meshing"; - - MESHING3_RESULT res = - MeshVolume (mparam, *mesh); + OCCMeshSurface (geom, *mesh, perfstepsend); + if (multithread.terminate) return TCL_OK; - ofstream problemfile("occmesh.rep",ios_base::app); - - problemfile << "VOLUMEMESHING" << endl << endl; - if(res != MESHING3_OK) - problemfile << "ERROR" << endl << endl; - else - problemfile << "OK" << endl - << mesh->GetNE() << " elements" << endl << endl; - - problemfile.close(); - - if (res != MESHING3_OK) return TCL_ERROR; - - if (multithread.terminate) return TCL_OK; - - RemoveIllegalElements (*mesh); - if (multithread.terminate) return TCL_OK; - - MeshQuality3d (*mesh); - -#ifdef STAT_STREAM - (*statout) << GetTime() << " & "; -#endif - #ifdef LOG_STREAM - (*logout) << "Volume meshed" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; + (*logout) << "Surfaces meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; #endif + +#ifdef STAT_STREAM + (*statout) << mesh->GetNSeg() << " & " + << mesh->GetNSE() << " & - &" + << GetTime() << " & " << endl; +#endif + + // MeshQuality2d (*mesh); + mesh->CalcSurfacesOfNode(); } - if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME) + if (multithread.terminate || perfstepsend <= MESHCONST_OPTSURFACE) return TCL_OK; - - if (perfstepsstart <= MESHCONST_OPTVOLUME) + if (perfstepsstart <= MESHCONST_MESHVOLUME) { - multithread.task = "Volume optimization"; - - OptimizeVolume (mparam, *mesh); - if (multithread.terminate) return TCL_OK; - -#ifdef STAT_STREAM - (*statout) << GetTime() << " & " - << mesh->GetNE() << " & " - << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl; -#endif + multithread.task = "Volume meshing"; -#ifdef LOG_STREAM - (*logout) << "Volume optimized" << endl - << "time = " << GetTime() << " sec" << endl - << "points: " << mesh->GetNP() << endl; + MESHING3_RESULT res = + MeshVolume (mparam, *mesh); + + ofstream problemfile("occmesh.rep",ios_base::app); + + problemfile << "VOLUMEMESHING" << endl << endl; + if(res != MESHING3_OK) + problemfile << "ERROR" << endl << endl; + else + problemfile << "OK" << endl + << mesh->GetNE() << " elements" << endl << endl; + + problemfile.close(); + + if (res != MESHING3_OK) return TCL_ERROR; + + if (multithread.terminate) return TCL_OK; + + RemoveIllegalElements (*mesh); + if (multithread.terminate) return TCL_OK; + + MeshQuality3d (*mesh); + +#ifdef STAT_STREAM + (*statout) << GetTime() << " & "; #endif - - // cout << "Optimization complete" << endl; - +#ifdef LOG_STREAM + (*logout) << "Volume meshed" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif } - (*testout) << "NP: " << mesh->GetNP() << endl; - for (int i = 1; i <= mesh->GetNP(); i++) + if (multithread.terminate || perfstepsend <= MESHCONST_MESHVOLUME) + return TCL_OK; + + if (perfstepsstart <= MESHCONST_OPTVOLUME) + { + multithread.task = "Volume optimization"; + + OptimizeVolume (mparam, *mesh); + if (multithread.terminate) return TCL_OK; + +#ifdef STAT_STREAM + (*statout) << GetTime() << " & " + << mesh->GetNE() << " & " + << mesh->GetNP() << " " << '\\' << '\\' << " \\" << "hline" << endl; +#endif + +#ifdef LOG_STREAM + (*logout) << "Volume optimized" << endl + << "time = " << GetTime() << " sec" << endl + << "points: " << mesh->GetNP() << endl; +#endif + + // cout << "Optimization complete" << endl; + + } + + (*testout) << "NP: " << mesh->GetNP() << endl; + for (int i = 1; i <= mesh->GetNP(); i++) (*testout) << mesh->Point(i) << endl; - - (*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl; - for (int i = 1; i <= mesh->GetNSeg(); i++) + + (*testout) << endl << "NSegments: " << mesh->GetNSeg() << endl; + for (int i = 1; i <= mesh->GetNSeg(); i++) (*testout) << mesh->LineSegment(i) << endl; - - - return TCL_OK; - } + return TCL_OK; + } } #endif diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 8b0399ee..770ce1b9 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -2,7 +2,7 @@ #ifdef OCCGEOMETRY #include -#include +#include #include "ShapeAnalysis_ShapeTolerance.hxx" #include "ShapeAnalysis_ShapeContents.hxx" #include "ShapeAnalysis_CheckSmallFace.hxx" @@ -62,9 +62,9 @@ namespace netgen for (e.Init(geom->shape, TopAbs_COMPSOLID); e.More(); e.Next()) count++; (*testout) << "CompSolids: " << count << endl; - + (*testout) << endl; - + cout << "Highest entry in topology hierarchy: " << endl; if (count) cout << count << " composite solid(s)" << endl; @@ -102,7 +102,7 @@ namespace netgen nrw = wmap.Extent(), nre = emap.Extent(), nrv = vmap.Extent(); - + TopExp_Explorer exp0; TopExp_Explorer exp1; @@ -111,10 +111,10 @@ namespace netgen for (exp0.Init(shape, TopAbs_COMPSOLID); exp0.More(); exp0.Next()) nrcs++; double surfacecont = 0; - - - + + + { Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); @@ -129,16 +129,16 @@ namespace netgen BuildFMap(); - + for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next()) { TopoDS_Face face = TopoDS::Face(exp0.Current()); - + GProp_GProps system; BRepGProp::SurfaceProperties(face, system); surfacecont += system.Mass(); } - + cout << "Starting geometry healing procedure (tolerance: " << tolerance << ")" << endl << "-----------------------------------" << endl; @@ -149,12 +149,12 @@ namespace netgen Handle(ShapeFix_Face) sff; Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); - + for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next()) { TopoDS_Face face = TopoDS::Face (exp0.Current()); - + sff = new ShapeFix_Face (face); sff->FixAddNaturalBoundMode() = Standard_True; sff->FixSmallAreaWireMode() = Standard_True; @@ -178,7 +178,7 @@ namespace netgen else if(sff->Status(ShapeExtend_DONE5)) cout << "(natural bounds added)" <Face(); - + rebuild->Replace(face, newface, Standard_False); } @@ -188,7 +188,7 @@ namespace netgen //delete rebuild; rebuild = NULL; } - + { Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; @@ -210,8 +210,8 @@ namespace netgen Handle(ShapeFix_Wire) sfw; Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); - - + + for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next()) { TopoDS_Face face = TopoDS::Face(exp0.Current()); @@ -225,12 +225,12 @@ namespace netgen sfw->ClosedWireMode() = Standard_True; bool replace = false; - + replace = sfw->FixReorder() || replace; replace = sfw->FixConnected() || replace; - + if (sfw->FixSmall (Standard_False, tolerance) && ! (sfw->StatusSmall(ShapeExtend_FAIL1) || sfw->StatusSmall(ShapeExtend_FAIL2) || @@ -238,7 +238,7 @@ namespace netgen { cout << "Fixed small edge in wire " << wmap.FindIndex (oldwire) << endl; replace = true; - + } else if (sfw->StatusSmall(ShapeExtend_FAIL1)) cerr << "Failed to fix small edge in wire " << wmap.FindIndex (oldwire) @@ -265,14 +265,14 @@ namespace netgen } //delete sfw; sfw = NULL; - + } } shape = rebuild->Apply(shape); - - + + { BuildFMap(); Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; @@ -281,14 +281,14 @@ namespace netgen for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next()) { TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); - if (vmap.FindIndex(TopExp::FirstVertex (edge)) == + if (vmap.FindIndex(TopExp::FirstVertex (edge)) == vmap.FindIndex(TopExp::LastVertex (edge))) { GProp_GProps system; BRepGProp::LinearProperties(edge, system); if (system.Mass() < tolerance) { - cout << "removing degenerated edge " << emap.FindIndex(edge) + cout << "removing degenerated edge " << emap.FindIndex(edge) << " from vertex " << vmap.FindIndex(TopExp::FirstVertex (edge)) << " to vertex " << vmap.FindIndex(TopExp::LastVertex (edge)) << endl; rebuild->Remove(edge, false); @@ -299,9 +299,9 @@ namespace netgen //delete rebuild; rebuild = NULL; } - - + + { Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); @@ -321,7 +321,7 @@ namespace netgen sfwf->SetPrecision(tolerance); sfwf->Load (shape); sfwf->ModeDropSmallEdges() = Standard_True; - + sfwf->SetPrecision(boundingbox.Diam()); if (sfwf->FixWireGaps()) @@ -333,10 +333,10 @@ namespace netgen if (sfwf->StatusWireGaps(ShapeExtend_FAIL1)) cout << "failed to fix some 2D gaps" << endl; if (sfwf->StatusWireGaps(ShapeExtend_FAIL2)) cout << "failed to fix some 3D gaps" << endl; } - + sfwf->SetPrecision(tolerance); - + { for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next()) { @@ -350,14 +350,14 @@ namespace netgen if (sfwf->FixSmallEdges()) { - cout << endl << "- fixing wire frames" << endl; + cout << endl << "- fixing wire frames" << endl; if (sfwf->StatusSmallEdges(ShapeExtend_OK)) cout << "no small edges found" << endl; if (sfwf->StatusSmallEdges(ShapeExtend_DONE1)) cout << "some small edges fixed" << endl; if (sfwf->StatusSmallEdges(ShapeExtend_FAIL1)) cout << "failed to fix some small edges" << endl; } - + shape = sfwf->Shape(); //delete sfwf; sfwf = NULL; @@ -368,7 +368,7 @@ namespace netgen - + { for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next()) { @@ -383,18 +383,18 @@ namespace netgen if (fixspotstripfaces) { - + cout << endl << "- fixing spot and strip faces" << endl; Handle(ShapeFix_FixSmallFace) sffsm = new ShapeFix_FixSmallFace(); sffsm -> Init (shape); sffsm -> SetPrecision (tolerance); sffsm -> Perform(); - + shape = sffsm -> FixShape(); //delete sffsm; sffsm = NULL; } - + { for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next()) { @@ -417,9 +417,9 @@ namespace netgen TopoDS_Face face = TopoDS::Face (exp0.Current()); sewedObj.Add (face); } - + sewedObj.Perform(); - + if (!sewedObj.SewedShape().IsNull()) shape = sewedObj.SewedShape(); else @@ -427,7 +427,7 @@ namespace netgen } - + { Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape; rebuild->Apply(shape); @@ -442,7 +442,7 @@ namespace netgen if (makesolids) - { + { cout << endl << "- making solids" << endl; BRepBuilderAPI_MakeSolid ms; @@ -452,7 +452,7 @@ namespace netgen count++; ms.Add (TopoDS::Shell(exp0.Current())); } - + if (!count) { cout << " not possible (no shells)" << endl; @@ -468,7 +468,7 @@ namespace netgen sfs->SetMaxTolerance(tolerance); sfs->Perform(); shape = sfs->Shape(); - + for (exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) { TopoDS_Solid solid = TopoDS::Solid(exp0.Current()); @@ -488,7 +488,7 @@ namespace netgen cout << " not possible" << endl; } } - + if (splitpartitions) @@ -498,30 +498,30 @@ namespace netgen TopExp_Explorer e2; Partition_Spliter ps; int count = 0; - + for (e2.Init (shape, TopAbs_SOLID); e2.More(); e2.Next()) { count++; ps.AddShape (e2.Current()); } - + ps.Compute(); shape = ps.Shape(); - + cout << " before: " << count << " solids" << endl; - + count = 0; for (e2.Init (shape, TopAbs_SOLID); e2.More(); e2.Next()) count++; - + cout << " after : " << count << " solids" << endl; } BuildFMap(); - + { for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next()) { @@ -533,8 +533,8 @@ namespace netgen double newsurfacecont = 0; - - + + for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next()) { TopoDS_Face face = TopoDS::Face(exp0.Current()); @@ -542,7 +542,7 @@ namespace netgen BRepGProp::SurfaceProperties(face, system); newsurfacecont += system.Mass(); } - + int nnrc = 0, nnrcs = 0, nnrso = somap.Extent(), @@ -580,7 +580,7 @@ namespace netgen // } } - + @@ -592,9 +592,9 @@ namespace netgen wmap.Clear(); emap.Clear(); vmap.Clear(); - + TopExp_Explorer exp0, exp1, exp2, exp3, exp4, exp5; - + for (exp0.Init(shape, TopAbs_COMPOUND); exp0.More(); exp0.Next()) { @@ -607,16 +607,16 @@ namespace netgen (*testout) << "shell " << ++i << endl; } } - + for (exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) { TopoDS_Solid solid = TopoDS::Solid (exp0.Current()); - + if (somap.FindIndex(solid) < 1) { somap.Add (solid); - + for (exp1.Init(solid, TopAbs_SHELL); exp1.More(); exp1.Next()) { @@ -625,7 +625,7 @@ namespace netgen if (shmap.FindIndex(shell) < 1) { shmap.Add (shell); - + for (exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()) { @@ -636,7 +636,7 @@ namespace netgen fmap.Add (face); (*testout) << "face " << fmap.FindIndex(face) << " "; (*testout) << ((face.Orientation() == TopAbs_REVERSED) ? "-" : "+") << ", "; - (*testout) << ((exp2.Current().Orientation() == TopAbs_REVERSED) ? "-" : "+") << endl; + (*testout) << ((exp2.Current().Orientation() == TopAbs_REVERSED) ? "-" : "+") << endl; for (exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()) { @@ -645,7 +645,7 @@ namespace netgen if (wmap.FindIndex(wire) < 1) { wmap.Add (wire); - + for (exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()) { @@ -671,7 +671,7 @@ namespace netgen } } } - + // Free Shells for (exp1.Init(shape, TopAbs_SHELL, TopAbs_SOLID); // for (exp1.Init(exp0.Current(), TopAbs_SHELL, TopAbs_SOLID); @@ -682,10 +682,10 @@ namespace netgen if (shmap.FindIndex(shell) < 1) { shmap.Add (shell); - + (*testout) << "shell " << shmap.FindIndex(shell) << " "; (*testout) << ((shell.Orientation() == TopAbs_REVERSED) ? "-" : "+") << ", "; - (*testout) << ((exp1.Current().Orientation() == TopAbs_REVERSED) ? "-" : "+") << endl; + (*testout) << ((exp1.Current().Orientation() == TopAbs_REVERSED) ? "-" : "+") << endl; for (exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()) @@ -695,7 +695,7 @@ namespace netgen if (fmap.FindIndex(face) < 1) { fmap.Add (face); - + for (exp3.Init(face, TopAbs_WIRE); exp3.More(); exp3.Next()) { @@ -704,7 +704,7 @@ namespace netgen if (wmap.FindIndex(wire) < 1) { wmap.Add (wire); - + for (exp4.Init(wire, TopAbs_EDGE); exp4.More(); exp4.Next()) { @@ -729,10 +729,10 @@ namespace netgen } } } - - + + // Free Faces - + for (exp2.Init(shape, TopAbs_FACE, TopAbs_SHELL); exp2.More(); exp2.Next()) { @@ -741,7 +741,7 @@ namespace netgen if (fmap.FindIndex(face) < 1) { fmap.Add (face); - + for (exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()) { @@ -750,7 +750,7 @@ namespace netgen if (wmap.FindIndex(wire) < 1) { wmap.Add (wire); - + for (exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()) { @@ -776,7 +776,7 @@ namespace netgen // Free Wires - + for (exp3.Init(shape, TopAbs_WIRE, TopAbs_FACE); exp3.More(); exp3.Next()) { @@ -784,7 +784,7 @@ namespace netgen if (wmap.FindIndex(wire) < 1) { wmap.Add (wire); - + for (exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()) { @@ -806,7 +806,7 @@ namespace netgen // Free Edges - + for (exp4.Init(shape, TopAbs_EDGE, TopAbs_WIRE); exp4.More(); exp4.Next()) { @@ -826,7 +826,7 @@ namespace netgen // Free Vertices - + for (exp5.Init(shape, TopAbs_VERTEX, TopAbs_EDGE); exp5.More(); exp5.Next()) { @@ -837,10 +837,18 @@ namespace netgen - + facemeshstatus.SetSize (fmap.Extent()); facemeshstatus = 0; + // Philippose - 15/01/2009 + face_maxh.SetSize (fmap.Extent()); + face_maxh = mparam.maxh; + + // Philippose - 17/01/2009 + face_sel_status.SetSize (fmap.Extent()); + face_sel_status = 0; + fvispar.SetSize (fmap.Extent()); evispar.SetSize (emap.Extent()); vvispar.SetSize (vmap.Extent()); @@ -867,9 +875,9 @@ namespace netgen TopoDS_Face face = TopoDS::Face (fmap(i)); sewedObj.Add (face); } - + sewedObj.Perform(); - + if (!sewedObj.SewedShape().IsNull()) { shape = sewedObj.SewedShape(); @@ -877,7 +885,7 @@ namespace netgen } else cout << " not possible"; - + /* ShapeUpgrade_ShellSewing sewing; TopoDS_Shape sh = sewing.ApplySewing (shape); @@ -915,14 +923,14 @@ namespace netgen { Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape; sfs->Init (ms); - + sfs->SetPrecision(1e-5); sfs->SetMaxTolerance(1e-5); - + sfs->Perform(); shape = sfs->Shape(); - + for (exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()) { TopoDS_Solid solid = TopoDS::Solid(exp0.Current()); @@ -932,13 +940,13 @@ namespace netgen // rebuild->Apply(shape); rebuild->Replace(solid, newsolid, Standard_False); // TopoDS_Shape newshape = rebuild->Apply(shape); - + TopoDS_Shape newshape = rebuild->Apply(shape, TopAbs_SHAPE, 1); shape = newshape; } //delete sfs; sfs = NULL; - + cout << " done" << endl; } else @@ -961,17 +969,17 @@ namespace netgen Bnd_Box bb; BRepBndLib::Add (shape, bb); - + double x1,y1,z1,x2,y2,z2; bb.Get (x1,y1,z1,x2,y2,z2); Point<3> p1 = Point<3> (x1,y1,z1); Point<3> p2 = Point<3> (x2,y2,z2); - + (*testout) << "Bounding Box = [" << p1 << " - " << p2 << "]" << endl; boundingbox = Box<3> (p1,p2); SetCenter(); - + } @@ -979,7 +987,7 @@ namespace netgen { static int cnt = 0; if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl; - + gp_Pnt pnt(p(0), p(1), p(2)); //(*testout) << "before " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl; @@ -1000,7 +1008,7 @@ namespace netgen } */ - + double u,v; Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi))); Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf ); @@ -1008,65 +1016,65 @@ namespace netgen suval.Coord( u, v); pnt = thesurf->Value( u, v ); - + p = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); - + } bool OCCGeometry :: FastProject (int surfi, Point<3> & ap, double& u, double& v) const { gp_Pnt p(ap(0), ap(1), ap(2)); - + Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi))); - + gp_Pnt x = surface->Value (u,v); - + if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true; - + gp_Vec du, dv; - + surface->D1(u,v,x,du,dv); - + int count = 0; - + gp_Pnt xold; gp_Vec n; double det, lambda, mu; - + do { count++; - + n = du^dv; - + det = Det3 (n.X(), du.X(), dv.X(), n.Y(), du.Y(), dv.Y(), n.Z(), du.Z(), dv.Z()); - - if (det < 1e-15) return false; - + + if (det < 1e-15) return false; + lambda = Det3 (n.X(), p.X()-x.X(), dv.X(), n.Y(), p.Y()-x.Y(), dv.Y(), n.Z(), p.Z()-x.Z(), dv.Z())/det; - + mu = Det3 (n.X(), du.X(), p.X()-x.X(), n.Y(), du.Y(), p.Y()-x.Y(), n.Z(), du.Z(), p.Z()-x.Z())/det; - + u += lambda; v += mu; - + xold = x; surface->D1(u,v,x,du,dv); - + } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50); // (*testout) << "FastProject count: " << count << endl; - + if (count == 50) return false; - + ap = Point<3> (x.X(), x.Y(), x.Z()); - + return true; } @@ -1082,22 +1090,22 @@ namespace netgen cout << "done" << endl; } - + OCCGeometry * LoadOCC_IGES (const char * filename) { OCCGeometry * occgeo; occgeo = new OCCGeometry; - + IGESControl_Reader reader; - + Standard_Integer stat = reader.ReadFile((char*)filename); // pre OCC52-times: // Standard_Integer stat = reader.LoadFile((char*)filename); // reader.Clear(); - - + + reader.TransferRoots(); // Tranlate IGES -> OCC // pre OCC52-times: @@ -1108,7 +1116,7 @@ namespace netgen occgeo->shape = reader.OneShape(); occgeo->changed = 1; occgeo->BuildFMap(); - + // // Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape; // sfs->Init(occgeo->shape); @@ -1127,69 +1135,145 @@ namespace netgen return occgeo; } - OCCGeometry * LoadOCC_STEP (const char * filename) - { - OCCGeometry * occgeo; - occgeo = new OCCGeometry; - - STEPControl_Reader reader; - Standard_Integer stat = reader.ReadFile((char*)filename); - Standard_Integer nb = reader.NbRootsForTransfer(); - reader.TransferRoots (); // Tranlate STEP -> OCC + + +// Philippose - 29/01/2009 +/* Special STEP File load function including the ability + to extract individual surface colours via the extended + OpenCascade XDE and XCAF Feature set. +*/ +OCCGeometry * LoadOCC_STEP (const char * filename) +{ + OCCGeometry * occgeo; + occgeo = new OCCGeometry; + + // Initiate a dummy XCAF Application to handle the STEP XCAF Document + static Handle_XCAFApp_Application dummy_app = XCAFApp_Application::GetApplication(); + + // Create an XCAF Document to contain the STEP file itself + Handle_TDocStd_Document step_doc; + + // Check if a STEP File is already open under this handle, if so, close it to prevent + // Segmentation Faults when trying to create a new document + if(dummy_app->NbDocuments() > 0) + { + dummy_app->GetDocument(1,step_doc); + dummy_app->Close(step_doc); + } + dummy_app->NewDocument ("STEP-XCAF",step_doc); + + STEPCAFControl_Reader reader; + + Standard_Integer stat = reader.ReadFile((char*)filename); + + // Enable transfer of colours + reader.SetColorMode(Standard_True); + + reader.Transfer(step_doc); + + // Read in the shape(s) and the colours present in the STEP File + Handle_XCAFDoc_ShapeTool step_shape_contents = XCAFDoc_DocumentTool::ShapeTool(step_doc->Main()); + Handle_XCAFDoc_ColorTool step_colour_contents = XCAFDoc_DocumentTool::ColorTool(step_doc->Main()); + + TDF_LabelSequence step_shapes; + step_shape_contents->GetShapes(step_shapes); + +/* + // List out the available colours in the STEP File as Colour Names + TDF_LabelSequence allColours; + stepColourContents->GetColors(allColours); + cout << "Number of colours in STEP = " << allColours.Length() << endl; + for(int i = 1; i <= allColours.Length(); i++) + { + Quantity_Color col; + stepColourContents->GetColor(allColours.Value(i),col); + cout << "Colour [" << i << "] = " << col.StringName(col.Name()) << endl; + } +*/ + + occgeo->shape = step_shape_contents->GetShape(step_shapes.Value(1)); + occgeo->face_colours = step_colour_contents; + occgeo->changed = 1; + occgeo->BuildFMap(); + + occgeo->BuildVisualizationMesh(); + PrintContents (occgeo); + + return occgeo; +} - occgeo->shape = reader.OneShape(); - occgeo->changed = 1; - occgeo->BuildFMap(); - // - //Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape; - //sfs->Init(occgeo->shape); - //sfs->Perform(); - //Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe(occgeo->shape); - //sfwf->FixSmallEdges(); - //sfwf->FixWireGaps(); +// Philippose - 29/01/2009 +// The LOADOCC_STEP Function has been replaced by the one +// above, which also includes support for the OpenCascade +// XDE Features. +// +// OCCGeometry * LoadOCC_STEP (const char * filename) +// { +// OCCGeometry * occgeo; +// occgeo = new OCCGeometry; +// +// STEPControl_Reader reader; +// Standard_Integer stat = reader.ReadFile((char*)filename); +// Standard_Integer nb = reader.NbRootsForTransfer(); +// reader.TransferRoots (); // Tranlate STEP -> OCC +// +// +// +// occgeo->shape = reader.OneShape(); +// occgeo->changed = 1; +// occgeo->BuildFMap(); +// // +// //Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape; +// //sfs->Init(occgeo->shape); +// //sfs->Perform(); +// //Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe(occgeo->shape); +// //sfwf->FixSmallEdges(); +// //sfwf->FixWireGaps(); +// +// +// +// /* +// // JS +// TopoDS_Compound aRes; +// BRep_Builder aBuilder; +// aBuilder.MakeCompound (aRes); +// +// for (TopExp_Explorer exp(occgeo->shape, TopAbs_SOLID); exp.More(); exp.Next()) +// { +// aBuilder.Add (aRes, exp.Current()); +// cout << "solid" << endl; +// } +// +// for (TopExp_Explorer exp(aRes, TopAbs_SOLID); exp.More(); exp.Next()) +// { +// cout << "compound has shapes solid" << endl; +// } +// occgeo->shape = aRes; +// occgeo->changed = 1; +// occgeo->BuildFMap(); +// */ +// +// +// occgeo->BuildVisualizationMesh(); +// PrintContents (occgeo); +// +// return occgeo; +// } - /* - // JS - TopoDS_Compound aRes; - BRep_Builder aBuilder; - aBuilder.MakeCompound (aRes); - - for (TopExp_Explorer exp(occgeo->shape, TopAbs_SOLID); exp.More(); exp.Next()) - { - aBuilder.Add (aRes, exp.Current()); - cout << "solid" << endl; - } - - for (TopExp_Explorer exp(aRes, TopAbs_SOLID); exp.More(); exp.Next()) - { - cout << "compound has shapes solid" << endl; - } - occgeo->shape = aRes; - occgeo->changed = 1; - occgeo->BuildFMap(); - */ - - - // - occgeo->BuildVisualizationMesh(); - PrintContents (occgeo); - - return occgeo; - } OCCGeometry * LoadOCC_BREP (const char * filename) { OCCGeometry * occgeo; occgeo = new OCCGeometry; - + BRep_Builder aBuilder; Standard_Boolean result = BRepTools::Read(occgeo->shape, const_cast (filename),aBuilder); - - + + // cout << "Writing VRML" << endl; // VrmlAPI::Write(occgeo->shape,"test.vrml"); // cout << "Writing STL" << endl; @@ -1202,7 +1286,7 @@ namespace netgen PrintContents (occgeo); return occgeo; - } + } char * shapesname[] = {" ", "CompSolids", "Solids", "Shells", @@ -1256,9 +1340,9 @@ namespace netgen case TopAbs_VERTEX: count2 = vmap.FindIndex(TopoDS::Vertex(e.Current())); break; } - + int nrsubshapes = 0; - + if (l <= TopAbs_WIRE) { TopExp_Explorer e2; @@ -1266,9 +1350,9 @@ namespace netgen e2.More(); e2.Next()) nrsubshapes++; } - + str << "{" << shapename[l] << " " << count2; - + if (l <= TopAbs_EDGE) { str << " (" << orientationstring[e.Current().Orientation()]; @@ -1280,7 +1364,7 @@ namespace netgen RecursiveTopologyTree (e.Current(), str, TopAbs_ShapeEnum (l+1), false, (char*)lname2.str().c_str()); - + } } @@ -1381,7 +1465,7 @@ namespace netgen str << "FaceSplitByVertices/Face" << i << " "; str << "{Face " << i << " (split by " << count << "vertex/vertices)} "; } - + int whatrow, sens; if (int type = csm.CheckPin (face, whatrow, sens)) { @@ -1471,9 +1555,9 @@ namespace netgen for (i = 1; i <= shmap.Extent(); i++) { TopoDS_Shell shell = TopoDS::Shell (shmap(i)); - if (!shell.Closed()) + if (!shell.Closed()) cout << "Shell " << i << " is not closed" << endl; - if (shell.Infinite()) + if (shell.Infinite()) cout << "Shell " << i << " is infinite" << endl; BRepCheck_Analyzer ba(shell); @@ -1484,9 +1568,9 @@ namespace netgen for (i = 1; i <= somap.Extent(); i++) { TopoDS_Solid solid = TopoDS::Solid (somap(i)); - if (!solid.Closed()) + if (!solid.Closed()) cout << "Solid " << i << " is not closed" << endl; - if (solid.Infinite()) + if (solid.Infinite()) cout << "Solid " << i << " is infinite" << endl; BRepCheck_Analyzer ba(solid); diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 085b1459..fda4a0ae 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -11,9 +11,9 @@ #include -#include -#include -#include +#include "BRep_Tool.hxx" +#include "Geom_Curve.hxx" +#include "Geom2d_Curve.hxx" #include "Geom_Surface.hxx" #include "GeomAPI_ProjectPointOnSurf.hxx" #include "GeomAPI_ProjectPointOnCurve.hxx" @@ -37,6 +37,7 @@ #include "TopoDS.hxx" #include "TopoDS_Solid.hxx" #include "TopExp_Explorer.hxx" +#include "TopTools_ListIteratorOfListOfShape.hxx" #include "BRep_Tool.hxx" #include "Geom_Curve.hxx" #include "Geom2d_Curve.hxx" @@ -62,8 +63,6 @@ #include "Poly_Triangle.hxx" #include "GProp_GProps.hxx" #include "BRepGProp.hxx" -#include "IGESControl_Reader.hxx" -#include "STEPControl_Reader.hxx" #include "TopoDS_Shape.hxx" #include "TopoDS_Face.hxx" #include "IGESToBRep_Reader.hxx" @@ -80,8 +79,33 @@ #include "Bnd_Box.hxx" #include "ShapeAnalysis.hxx" #include "ShapeBuild_ReShape.hxx" + +// Philippose - 29/01/2009 +// OpenCascade XDE Support +// Include support for OpenCascade XDE Features +#include "TDocStd_Document.hxx" +#include "Quantity_Color.hxx" +#include "XCAFApp_Application.hxx" +#include "XCAFDoc_ShapeTool.hxx" +#include "XCAFDoc_Color.hxx" +#include "XCAFDoc_ColorTool.hxx" +#include "XCAFDoc_ColorType.hxx" +#include "XCAFDoc_LayerTool.hxx" +#include "XCAFDoc_DimTolTool.hxx" +#include "XCAFDoc_MaterialTool.hxx" +#include "XCAFDoc_DocumentTool.hxx" +#include "TDF_Label.hxx" +#include "TDF_LabelSequence.hxx" +#include "STEPCAFControl_Reader.hxx" +#include "STEPCAFControl_Writer.hxx" +#include "IGESCAFControl_Reader.hxx" +#include "IGESCAFControl_Writer.hxx" + +#include "IGESControl_Reader.hxx" +#include "STEPControl_Reader.hxx" #include "IGESControl_Writer.hxx" #include "STEPControl_Writer.hxx" + #include "StlAPI_Writer.hxx" #include "STEPControl_StepModelType.hxx" @@ -89,175 +113,241 @@ namespace netgen { #include "../visualization/vispar.hpp" - // class VisualizationParameters; - // extern VisualizationParameters vispar; + // class VisualizationParameters; + // extern VisualizationParameters vispar; #include "occmeshsurf.hpp" #define PROJECTION_TOLERANCE 1e-10 - #define ENTITYISVISIBLE 1 #define ENTITYISHIGHLIGHTED 2 #define ENTITYISDRAWABLE 4 -class EntityVisualizationCode -{ - int code; + class EntityVisualizationCode + { + int code; -public: + public: - EntityVisualizationCode() - { code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE; } + EntityVisualizationCode() + { code = ENTITYISVISIBLE + !ENTITYISHIGHLIGHTED + ENTITYISDRAWABLE;} - int IsVisible () - { return code & ENTITYISVISIBLE; } + int IsVisible () + { return code & ENTITYISVISIBLE;} - int IsHighlighted () - { return code & ENTITYISHIGHLIGHTED; } + int IsHighlighted () + { return code & ENTITYISHIGHLIGHTED;} - int IsDrawable () - { return code & ENTITYISDRAWABLE; } + int IsDrawable () + { return code & ENTITYISDRAWABLE;} - void Show () - { code |= ENTITYISVISIBLE; } + void Show () + { code |= ENTITYISVISIBLE;} - void Hide () - { code &= ~ENTITYISVISIBLE; } + void Hide () + { code &= ~ENTITYISVISIBLE;} - void Highlight () - { code |= ENTITYISHIGHLIGHTED; } + void Highlight () + { code |= ENTITYISHIGHLIGHTED;} - void Lowlight () - { code &= ~ENTITYISHIGHLIGHTED; } + void Lowlight () + { code &= ~ENTITYISHIGHLIGHTED;} - void SetDrawable () - { code |= ENTITYISDRAWABLE; } - - void SetNotDrawable () - { code &= ~ENTITYISDRAWABLE; } -}; - - - -inline double Det3 (double a00, double a01, double a02, - double a10, double a11, double a12, - double a20, double a21, double a22) -{ - return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00; -} + void SetDrawable () + { code |= ENTITYISDRAWABLE;} + void SetNotDrawable () + { code &= ~ENTITYISDRAWABLE;} + }; + inline double Det3 (double a00, double a01, double a02, + double a10, double a11, double a12, + double a20, double a21, double a22) + { + return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00; + } #define OCCGEOMETRYVISUALIZATIONNOCHANGE 0 #define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1 - // == compute transformation matrices and redraw + // == compute transformation matrices and redraw #define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2 - // == redraw + // == redraw -class OCCGeometry -{ - Point<3> center; + class OCCGeometry + { + Point<3> center; -public: - TopoDS_Shape shape; - TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; - Array fsingular, esingular, vsingular; - Box<3> boundingbox; + public: + TopoDS_Shape shape; + TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap; + Array fsingular, esingular, vsingular; + Box<3> boundingbox; - int changed; - Array facemeshstatus; + // Philippose - 29/01/2009 + // OpenCascade XDE Support + // XCAF Handle to make the face colours available to the rest of + // the system + Handle_XCAFDoc_ColorTool face_colours; - Array fvispar, evispar, vvispar; + int changed; + Array facemeshstatus; - double tolerance; - bool fixsmalledges; - bool fixspotstripfaces; - bool sewfaces; - bool makesolids; - bool splitpartitions; + // Philippose - 15/01/2009 + // Maximum mesh size for a given face + // (Used to explicitly define mesh size limits on individual faces) + Array face_maxh; + // Philippose - 15/01/2009 + // Indicates which faces have been selected by the user in geometry mode + // (Currently handles only selection of one face at a time, but an array would + // help to extend this to multiple faces) + Array face_sel_status; - OCCGeometry() - { - somap.Clear(); - shmap.Clear(); - fmap.Clear(); - wmap.Clear(); - emap.Clear(); - vmap.Clear(); - } + Array fvispar, evispar, vvispar; + double tolerance; + bool fixsmalledges; + bool fixspotstripfaces; + bool sewfaces; + bool makesolids; + bool splitpartitions; - void BuildFMap(); + OCCGeometry() + { + somap.Clear(); + shmap.Clear(); + fmap.Clear(); + wmap.Clear(); + emap.Clear(); + vmap.Clear(); + } - Box<3> GetBoundingBox() - { return boundingbox; } + void BuildFMap(); - int NrSolids() - { return somap.Extent(); } + Box<3> GetBoundingBox() + { return boundingbox;} - void SetCenter() - { center = boundingbox.Center(); } + int NrSolids() + { return somap.Extent();} - Point<3> Center() - { return center; } + // Philippose - 17/01/2009 + // Total number of faces in the geometry + int NrFaces() + { return fmap.Extent();} - void Project (int surfi, Point<3> & p) const; - bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; + void SetCenter() + { center = boundingbox.Center();} - - OCCSurface GetSurface (int surfi) - { - cout << "OCCGeometry::GetSurface using PLANESPACE" << endl; - return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE); - } - + Point<3> Center() + { return center;} - void BuildVisualizationMesh (); + void Project (int surfi, Point<3> & p) const; + bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; - void RecursiveTopologyTree (const TopoDS_Shape & sh, - stringstream & str, - TopAbs_ShapeEnum l, - bool free, - const char * lname); + OCCSurface GetSurface (int surfi) + { + cout << "OCCGeometry::GetSurface using PLANESPACE" << endl; + return OCCSurface (TopoDS::Face(fmap(surfi)), PLANESPACE); + } - void GetTopologyTree (stringstream & str); + void BuildVisualizationMesh (); - void PrintNrShapes (); + void RecursiveTopologyTree (const TopoDS_Shape & sh, + stringstream & str, + TopAbs_ShapeEnum l, + bool free, + const char * lname); - void CheckIrregularEntities (stringstream & str); + void GetTopologyTree (stringstream & str); - void SewFaces(); + void PrintNrShapes (); - void MakeSolid(); + void CheckIrregularEntities (stringstream & str); - void HealGeometry(); + void SewFaces(); - void LowLightAll() - { - for (int i = 1; i <= fmap.Extent(); i++) - fvispar[i-1].Lowlight(); - for (int i = 1; i <= emap.Extent(); i++) - evispar[i-1].Lowlight(); - for (int i = 1; i <= vmap.Extent(); i++) - vvispar[i-1].Lowlight(); - } + void MakeSolid(); - void GetUnmeshedFaceInfo (stringstream & str); - void GetNotDrawableFaces (stringstream & str); - bool ErrorInSurfaceMeshing (); + void HealGeometry(); - void WriteOCC_STL(char * filename); -}; + // Philippose - 15/01/2009 + // Sets the maximum mesh size for a given face + // (Note: Local mesh size limited by the global max mesh size) + void SetFaceMaxH(int facenr, double faceh) + { + if((facenr> 0) && (facenr <= fmap.Extent())) + { + face_maxh[facenr-1] = min(mparam.maxh,faceh); + } + } + // Philippose - 15/01/2009 + // Returns the local mesh size of a given face + double GetFaceMaxH(int facenr) + { + if((facenr> 0) && (facenr <= fmap.Extent())) + { + return face_maxh[facenr-1]; + } + else + { + return 0.0; + } + } -void PrintContents (OCCGeometry * geom); + // Philippose - 17/01/2009 + // Returns the index of the currently selected face + int SelectedFace() + { + int i; -OCCGeometry * LoadOCC_IGES (const char * filename); -OCCGeometry * LoadOCC_STEP (const char * filename); -OCCGeometry * LoadOCC_BREP (const char * filename); + for(i = 1; i <= fmap.Extent(); i++) + { + if(face_sel_status[i-1]) + { + return i; + } + } + + return 0; + } + + // Philippose - 17/01/2009 + // Sets the currently selected face + void SetSelectedFace(int facenr) + { + face_sel_status = 0; + + if((facenr >= 1) && (facenr <= fmap.Extent())) + { + face_sel_status[facenr-1] = 1; + } + } + + void LowLightAll() + { + for (int i = 1; i <= fmap.Extent(); i++) + fvispar[i-1].Lowlight(); + for (int i = 1; i <= emap.Extent(); i++) + evispar[i-1].Lowlight(); + for (int i = 1; i <= vmap.Extent(); i++) + vvispar[i-1].Lowlight(); + } + + void GetUnmeshedFaceInfo (stringstream & str); + void GetNotDrawableFaces (stringstream & str); + bool ErrorInSurfaceMeshing (); + + void WriteOCC_STL(char * filename); + }; + + void PrintContents (OCCGeometry * geom); + + OCCGeometry * LoadOCC_IGES (const char * filename); + OCCGeometry * LoadOCC_STEP (const char * filename); + OCCGeometry * LoadOCC_BREP (const char * filename); } diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 410dfaf0..84848e4a 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -9,18 +9,25 @@ #include #include -#include - // #include +#ifdef OCCGEOMETRY +// Philippose - 30/01/2009 +// Required for OpenCascade XDE Support +#include +#endif - +#include namespace netgen { - +#ifdef OCCGEOMETRY + // Philippose - 30/01/2009 + // Required for OpenCascade XDE Support + extern OCCGeometry * occgeometry; +#endif extern AutoPtr mesh; @@ -50,12 +57,12 @@ namespace netgen linetimestamp = GetTimeStamp(); edgetimestamp = GetTimeStamp(); pointnumbertimestamp = GetTimeStamp(); - + tettimestamp = GetTimeStamp(); prismtimestamp = GetTimeStamp(); hextimestamp = GetTimeStamp(); pyramidtimestamp = GetTimeStamp(); - + badeltimestamp = GetTimeStamp(); identifiedtimestamp = GetTimeStamp(); domainsurftimestamp = GetTimeStamp(); @@ -74,25 +81,25 @@ namespace netgen ; } - + void VisualSceneMesh :: DrawScene () { - if (!mesh) + if (!mesh) { VisualScene::DrawScene(); return; } - + lock = NULL; - + clock_t starttime, endtime; starttime = clock(); - + BuildScene(); - + glClearColor(backcolor, backcolor, backcolor, 1.0); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - + glEnable (GL_COLOR_MATERIAL); glColor3f (1.0f, 1.0f, 1.0f); glLineWidth (1.0f); @@ -104,8 +111,8 @@ namespace netgen GLdouble projmat[16]; glGetDoublev (GL_PROJECTION_MATRIX, projmat); - - + + #ifdef PARALLEL glEnable (GL_BLEND); glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); @@ -121,23 +128,23 @@ namespace netgen // glDisable (GL_DEPTH_TEST); // glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // glHint (GL_LINE_SMOOTH_HINT, GL_DONT_CARE); - + glDisable (GL_COLOR_MATERIAL); - + GLfloat matcol0[] = { 0, 0, 0, 1 }; GLfloat matcol1[] = { 1, 1, 1, 1 }; GLfloat matcolf[] = { 0, 1, 0, 1 }; GLfloat matcolb[] = { 0.5, 0, 0, 1 }; GLfloat matcolblue[] = { 0, 0, 1, 1 }; - - glMatrixMode (GL_MODELVIEW); - + + glMatrixMode (GL_MODELVIEW); + glMaterialfv(GL_FRONT, GL_EMISSION, matcol0); glMaterialfv(GL_BACK, GL_EMISSION, matcol0); glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, matcol1); glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolf); glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, matcolb); - + glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); @@ -198,15 +205,15 @@ namespace netgen } glDisable (GL_POLYGON_OFFSET_FILL); - + // draw lines - glMatrixMode (GL_MODELVIEW); - + glMatrixMode (GL_MODELVIEW); + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcol0); glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, matcol0); glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, matcol0); - + glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); glLineWidth (1.0f); glColor3f (0.0f, 0.0f, 0.0f); @@ -215,7 +222,7 @@ namespace netgen - + if (vispar.drawoutline) { glPolygonOffset (1, 1); @@ -246,14 +253,14 @@ namespace netgen glCallList (identifiedlist); glDisable (GL_POLYGON_OFFSET_LINE); } - + if (vispar.drawpointnumbers || vispar.drawedgenumbers || vispar.drawfacenumbers || vispar.drawelementnumbers) glCallList (pointnumberlist); - - + + glPopName(); if (vispar.drawedges) @@ -270,20 +277,20 @@ namespace netgen glColor3d (0, 0, 1); glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolblue); glBegin (GL_POINTS); - + const Point3d p = mesh->Point(selpoint); glVertex3f (p.X(), p.Y(), p.Z()); glEnd(); */ glColor3d (0, 0, 1); - - static GLubyte cross[] = + + static GLubyte cross[] = { 0xc6, 0xee, 0x7c, 0x38, 0x7c, 0xee, 0xc6 }; glPixelStorei(GL_UNPACK_ALIGNMENT, 1); - + glDisable (GL_COLOR_MATERIAL); glDisable (GL_LIGHTING); glDisable (GL_CLIP_PLANE0); @@ -311,7 +318,7 @@ namespace netgen lock = NULL; } - glFinish(); + glFinish(); endtime = clock(); @@ -334,11 +341,11 @@ namespace netgen } int i, j; - + Point3d pmin, pmax; static double oldrad = 0; - + Array faces; int meshtimestamp = mesh->GetTimeStamp(); @@ -352,9 +359,9 @@ namespace netgen else { // otherwise strange zooms douring mesh generation - mesh->GetBox (pmin, pmax, SURFACEPOINT); + mesh->GetBox (pmin, pmax, SURFACEPOINT); } - + if (vispar.use_center_coords && zoomall == 2) { center.X() = vispar.centerx; center.Y() = vispar.centery; center.Z() = vispar.centerz; @@ -367,9 +374,9 @@ namespace netgen center = Center (pmin, pmax); rad = 0.5 * Dist (pmin, pmax); if(rad == 0) rad = 1e-6; - - if (rad > 1.2 * oldrad || - mesh->GetMajorTimeStamp() > vstimestamp || + + if (rad > 1.2 * oldrad || + mesh->GetMajorTimeStamp() > vstimestamp || zoomall) { CalcTransformationMatrices(); @@ -381,7 +388,7 @@ namespace netgen if (pointnumberlist) { - glDeleteLists (pointnumberlist, 1); + glDeleteLists (pointnumberlist, 1); pointnumberlist = 0; } @@ -431,7 +438,7 @@ namespace netgen glNormal3d (0, 0, 1); glPushAttrib (GL_LIST_BIT); glListBase (fontbase); - + char buf[30]; if (vispar.drawpointnumbers) @@ -456,12 +463,12 @@ namespace netgen const Point3d & p2 = mesh->Point(seg.p2); const Point3d p = Center (p1, p2); glRasterPos3d (p.X(), p.Y(), p.Z()); - + sprintf (buf, "%d", seg.edgenr); glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); } */ - + const MeshTopology & top = mesh->GetTopology(); for (i = 1; i <= top.GetNEdges(); i++) { @@ -475,8 +482,8 @@ namespace netgen sprintf (buf, "%d", i); glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); } - - } + + } if (vispar.drawfacenumbers) @@ -506,9 +513,9 @@ namespace netgen sprintf (buf, "%d", i); glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); } - } + } + - if (vispar.drawelementnumbers) { @@ -523,12 +530,12 @@ namespace netgen if ( ! el.PNum(5)) // eltype == TET ) { - + pnums.SetSize(4); for( int j = 0; j < pnums.Size(); j++) pnums[j] = mesh->VolumeElement(i).PNum(j+1); - - + + const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p2 = mesh->Point(pnums[1]); const Point3d & p3 = mesh->Point(pnums[2]); @@ -540,14 +547,14 @@ namespace netgen pnums.SetSize(5); for( int j = 0; j < pnums.Size(); j++) pnums[j] = mesh->VolumeElement(i).PNum(j+1); - - + + const Point3d & p1 = mesh->Point(pnums[0]); const Point3d & p2 = mesh->Point(pnums[1]); const Point3d & p3 = mesh->Point(pnums[2]); const Point3d & p4 = mesh->Point(pnums[3]); const Point3d & p5 = mesh->Point(pnums[4]); - + p.X() = 0.3 * p5.X() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . X(); p.Y() = 0.3 * p5.Y() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . Y(); p.Z() = 0.3 * p5.Z() + 0.7 * Center ( Center(p1, p3) , Center(p2, p4) ) . Z(); @@ -583,17 +590,17 @@ namespace netgen const Point3d & p7 = mesh->Point(pnums[6]); const Point3d & p8 = mesh->Point(pnums[7]); - p = Center ( Center ( Center(p1, p3), Center(p2, p4) ) , Center( Center(p5, p7) , Center(p6, p8 ) ) ); + p = Center ( Center ( Center(p1, p3), Center(p2, p4) ) , Center( Center(p5, p7) , Center(p6, p8 ) ) ); } glRasterPos3d (p.X(), p.Y(), p.Z()); sprintf (buf, "%d", i); glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf); - - } - } - + } + } + + glPopAttrib (); // glDisable (GL_COLOR_MATERIAL); } @@ -616,8 +623,8 @@ namespace netgen for (i = 1; i <= mesh->GetNE(); i++) { - if (mesh->VolumeElement(i).flags.badel || - mesh->VolumeElement(i).flags.illegal || + if (mesh->VolumeElement(i).flags.badel || + mesh->VolumeElement(i).flags.illegal || (i == vispar.drawelement)) { // copy to be thread-safe @@ -631,7 +638,7 @@ namespace netgen if (el.PNum(1)) { glBegin (GL_TRIANGLES); - + for (j = 1; j <= faces.Size(); j++) { Element2d & face = faces.Elem(j); @@ -645,7 +652,7 @@ namespace netgen glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); glVertex3d (lp3.X(), lp3.Y(), lp3.Z()); } - + glEnd(); } } @@ -669,7 +676,7 @@ namespace netgen } } } - + for (i = 1; i <= mesh->GetNE(); i++) { @@ -688,7 +695,7 @@ namespace netgen if (el.GetNP() == 4) { - int et[6][2] = + int et[6][2] = { { 1, 2 }, { 1, 3 }, { 1, 4 }, @@ -710,7 +717,7 @@ namespace netgen if (el.GetNP() == 10) { - int et[12][2] = + int et[12][2] = { { 1, 5 }, { 2, 5 }, { 1, 6 }, @@ -758,7 +765,7 @@ namespace netgen case TRIG: { glBegin (GL_TRIANGLES); - + Point3d lp1 = mesh->Point (el.PNum(1)); Point3d lp2 = mesh->Point (el.PNum(2)); Point3d lp3 = mesh->Point (el.PNum(3)); @@ -774,12 +781,12 @@ namespace netgen case QUAD: { glBegin (GL_QUADS); - + const Point3d & lp1 = mesh->Point (el.PNum(1)); const Point3d & lp2 = mesh->Point (el.PNum(2)); const Point3d & lp3 = mesh->Point (el.PNum(4)); const Point3d & lp4 = mesh->Point (el.PNum(3)); - Vec3d n = Cross (Vec3d (lp1, lp2), + Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, Center (lp3, lp4))); n /= (n.Length() + 1e-12); glNormal3d (n.X(), n.Y(), n.Z()); @@ -796,7 +803,7 @@ namespace netgen { 1, 6 }, { 2, 6 }, { 1, 5 }, { 3, 5 }, { 2, 4 }, { 3, 4 } }; - + glBegin (GL_LINES); for (j = 0; j < 6; j++) { @@ -813,14 +820,14 @@ namespace netgen { 1, 5 }, { 2, 5 }, { 3, 6 }, { 4, 6 }, { 1, 4 }, { 2, 3 } }; - + glBegin (GL_LINES); - + for (j = 0; j < 6; j++) { const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0])); const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1])); - + glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); } @@ -828,31 +835,31 @@ namespace netgen break; } default: - PrintSysError ("Cannot draw surface element of type ", + PrintSysError ("Cannot draw surface element of type ", int(el.GetType())); } } - glLoadName (0); + glLoadName (0); } glEndList (); - + if (1) { - + identifiedlist = glGenLists (1); glNewList (identifiedlist, GL_COMPILE); - + GLfloat identifiedcol[] = { 1, 0, 1, 1 }; - + glLineWidth (3); - + // for (i = 1; i <= mesh->GetNSeg(); i++) if (& mesh -> GetIdentifications() ) { - INDEX_2_HASHTABLE & idpts = + INDEX_2_HASHTABLE & idpts = mesh->GetIdentifications().GetIdentifiedPoints(); if (&idpts) { @@ -861,18 +868,18 @@ namespace netgen { INDEX_2 pts; int val; - + idpts.GetData (i, j, pts, val); const Point3d & p1 = mesh->Point(pts.I1()); const Point3d & p2 = mesh->Point(pts.I2()); - - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, + + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, identifiedcol); - + glBegin (GL_LINES); glVertex3f (p1.X(), p1.Y(), p1.Z()); glVertex3f (p2.X(), p2.Y(), p2.Z()); - glEnd(); + glEnd(); } } } @@ -897,7 +904,7 @@ namespace netgen { static int timer = NgProfiler::CreateTimer ("Mesh::BuildFilledList"); NgProfiler::RegionTimer reg (timer); - + // cout << "buildilled" << endl; #ifdef PARALLELGL @@ -925,10 +932,10 @@ namespace netgen filledlist = glGenLists (1); glNewList (filledlist, GL_COMPILE); - + for ( int dest = 1; dest < ntasks; dest++ ) glCallList (par_filledlists[dest]); - + glEndList(); @@ -956,15 +963,15 @@ namespace netgen filledlist = glGenLists (1); glNewList (filledlist, GL_COMPILE); - + // cout << "I am p " << id << " and got filledlist " << filledlist << endl; bool checkvicinity = (stlgeometry != NULL) && stldoctor.showvicinity; glEnable (GL_NORMALIZE); - + glLineWidth (1.0f); - + Vector locms; if (vispar.colormeshsize) @@ -995,7 +1002,7 @@ namespace netgen #endif CurvedElements & curv = mesh->GetCurvedElements(); - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; for (int col = 1; col <= 2; col++) { @@ -1016,7 +1023,7 @@ namespace netgen else glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_coll_transp); } - else + else { if (col == 2) glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcolsel); @@ -1025,6 +1032,23 @@ namespace netgen } #endif +#ifdef OCCGEOMETRY + // Philippose - 30/01/2009 + // OpenCascade XDE Support + // Update the colour of each face based on the STEP File Data + // if the advanced OpenCascade XDE Support has been enabled + if((col == 1) && (occgeometry)) + { + TopoDS_Face face = TopoDS::Face(occgeometry->fmap(el.GetIndex())); + Quantity_Color face_colour; + occgeometry->face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour); + matcol[0] = face_colour.Red(); + matcol[1] = face_colour.Green(); + matcol[2] = face_colour.Blue(); + matcol[3] = 1.0; + glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol); + } +#endif bool drawel = !el.IsDeleted(); @@ -1057,7 +1081,7 @@ namespace netgen for (int i = 0; i < hoplotn; i++) { glBegin (GL_TRIANGLE_STRIP); - + for (int j = 0; j <= hoplotn-i; j++) for (int k = 0; k < 2; k++) { @@ -1093,15 +1117,15 @@ namespace netgen } glEnd(); } - } + } else // not high order { glBegin (GL_TRIANGLES); - + const Point<3> & lp0 = (*mesh) [el[0]]; const Point<3> & lp1 = (*mesh) [el[1]]; const Point<3> & lp2 = (*mesh) [el[2]]; - + Vec<3> n = Cross (lp1-lp0, lp2-lp0); glNormal3dv (n); @@ -1123,7 +1147,7 @@ namespace netgen glEnd(); } - + break; } case QUAD: @@ -1164,15 +1188,15 @@ namespace netgen } } - + glEnd(); - } + } else // not high order { glBegin (GL_QUADS); - + const Point<3> & lp1 = mesh->Point (el.PNum(1)); const Point<3> & lp2 = mesh->Point (el.PNum(2)); const Point<3> & lp3 = mesh->Point (el.PNum(4)); @@ -1186,7 +1210,7 @@ namespace netgen glVertex3dv (lp2); glVertex3dv (lp4); glVertex3dv (lp3); - + glEnd (); } break; @@ -1201,7 +1225,7 @@ namespace netgen { 2, 4, 6 }, { 3, 5, 4 }, { 4, 5, 6 } }; - + for (int j = 0; j < 4; j++) { const Point<3> & lp1 = mesh->Point (el.PNum(trigs[j][0])); @@ -1225,7 +1249,7 @@ namespace netgen static int quads[2][4] = { { 1, 5, 6, 4 }, { 5, 2, 3, 6 } }; - + for (int j = 0; j < 2; j++) { Point3d lp1 = mesh->Point (el.PNum(quads[j][0])); @@ -1247,9 +1271,9 @@ namespace netgen case QUAD8: { glBegin (GL_TRIANGLES); - static int boundary[] = + static int boundary[] = { 1, 5, 2, 8, 3, 6, 4, 7, 1 }; - + Point3d c(0,0,0); for (int j = 0; j < 4; j++) { @@ -1284,7 +1308,7 @@ namespace netgen default: - PrintSysError ("Cannot draw (2) surface element of type ", + PrintSysError ("Cannot draw (2) surface element of type ", int(el.GetType())); } @@ -1333,13 +1357,13 @@ namespace netgen linelist = glGenLists (1); glNewList (linelist, GL_COMPILE); - + for ( int dest = 1; dest < ntasks; dest++ ) glCallList (par_linelists[dest]); - + glEndList(); - + linetimestamp = NextTimeStamp(); return; } @@ -1358,7 +1382,7 @@ namespace netgen if (linelist) glDeleteLists (linelist, 1); - + linelist = glGenLists (1); glNewList (linelist, GL_COMPILE); @@ -1367,13 +1391,13 @@ namespace netgen glLineWidth (1.0f); - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; // PrintMessage (3, "nse = ", mesh->GetNSE()); for (SurfaceElementIndex sei = 0; sei < mesh->GetNSE(); sei++) { const Element2d & el = (*mesh)[sei]; - + bool drawel = !el.IsDeleted(); if (checkvicinity) for (int j = 0; j < el.GetNP(); j++) @@ -1412,8 +1436,8 @@ namespace netgen } glEnd(); - } - else + } + else { glBegin (GL_TRIANGLES); @@ -1430,8 +1454,8 @@ namespace netgen break; - } - + } + case QUAD: { CurvedElements & curv = mesh->GetCurvedElements(); @@ -1492,46 +1516,46 @@ namespace netgen glEnd(); } - + break; - + } - + case TRIG6: { int lines[6][2] = { { 1, 6 }, { 2, 6 }, { 1, 5 }, { 3, 5 }, { 2, 4 }, { 3, 4 } }; - + glBegin (GL_LINES); for (int j = 0; j < 6; j++) { const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0])); const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1])); - + glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); } - + glEnd(); break; } - + case QUAD6: { int lines[6][2] = { { 1, 5 }, { 2, 5 }, { 3, 6 }, { 4, 6 }, { 1, 4 }, { 2, 3 } }; - + glBegin (GL_LINES); - + for (int j = 0; j < 6; j++) { const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0])); const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1])); - + glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); } @@ -1545,14 +1569,14 @@ namespace netgen { 1, 5 }, { 2, 5 }, { 3, 6 }, { 4, 6 }, { 1, 7 }, { 4, 7 }, { 2, 8 }, { 3, 8 } }; - + glBegin (GL_LINES); - + for (int j = 0; j < 8; j++) { const Point3d & lp1 = mesh->Point (el.PNum(lines[j][0])); const Point3d & lp2 = mesh->Point (el.PNum(lines[j][1])); - + glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); } @@ -1563,14 +1587,14 @@ namespace netgen default: - PrintSysError ("Cannot draw (4) surface element of type ", + PrintSysError ("Cannot draw (4) surface element of type ", int(el.GetType())); } } - + glEndList (); - + #ifdef PARALLELGL glFinish(); if (id > 0) @@ -1595,7 +1619,7 @@ namespace netgen if (edgelist) glDeleteLists (edgelist, 1); - + edgelist = glGenLists (1); glNewList (edgelist, GL_COMPILE); @@ -1605,21 +1629,21 @@ namespace netgen glEnable (GL_POLYGON_OFFSET_LINE); glPolygonOffset (1, -1); - + glEnable (GL_COLOR_MATERIAL); glDisable (GL_LIGHTING); - + for (int i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); const Point3d & p1 = (*mesh)[seg.p1]; const Point3d & p2 = (*mesh)[seg.p2]; - + if (seg.singedge_left || seg.singedge_right) - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolsingedge); else - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcoledge); if (seg.singedge_left || seg.singedge_right) @@ -1631,31 +1655,31 @@ namespace netgen glLineWidth(5); else glLineWidth(2); - - if (mesh->GetCurvedElements().IsHighOrder()) + + if (mesh->GetCurvedElements().IsHighOrder()) { int j; - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; // mesh->GetCurvedElements().GetNVisualSubsecs(); - + Point<3> x; glBegin (GL_LINE_STRIP); - - for (int j = 0; j <= hoplotn; j++) + + for (int j = 0; j <= hoplotn; j++) { mesh->GetCurvedElements().CalcSegmentTransformation ((double) j/hoplotn, i-1, x); glVertex3d (x(0), x(1), x(2)); /* - cout << "x = " << x(0) << ", " << x(1) << ", " << x(2) + cout << "x = " << x(0) << ", " << x(1) << ", " << x(2) << ", norm = 1+" << sqrt(x(0)*x(0)+x(1)*x(1))-1 << ", phi = " << atan2(x(1), x(0))/M_PI << endl; */ } - + glEnd(); - - } + + } else { glBegin (GL_LINES); @@ -1672,10 +1696,10 @@ namespace netgen glEnd(); } } - + glLineWidth (2); glDisable (GL_POLYGON_OFFSET_LINE); - + glDisable (GL_COLOR_MATERIAL); glEnable (GL_LIGHTING); @@ -1691,14 +1715,14 @@ namespace netgen } - + // Bernstein Pol B_{n,i}(x) = n! / i! / (n-i)! (1-x)^{n-i} x^i static inline double Bernstein (int n, int i, double x) { double val = 1; - for (int j = 1; j <= i; j++) + for (int j = 1; j <= i; j++) val *= x; - for (int j = 1; j <= n-i; j++) + for (int j = 1; j <= n-i; j++) val *= (1-x) * (j+i) / j; return val; } @@ -1774,7 +1798,7 @@ namespace netgen Array faces; - + BitArray shownode(mesh->GetNP()); if (vispar.clipenable) { @@ -1782,13 +1806,13 @@ namespace netgen for (int i = 1; i <= shownode.Size(); i++) { Point<3> p = mesh->Point(i); - + double val = p[0] * clipplane[0] + p[1] * clipplane[1] + p[2] * clipplane[2] + clipplane[3]; - + if (val > 0) shownode.Set (i); } } @@ -1798,7 +1822,7 @@ namespace netgen #ifdef PARALLEL - static float tetcols[][8] = + static float tetcols[][8] = { { 1.0f, 1.0f, 0.0f, 1.0f }, { 1.0f, 0.0f, 0.0f, 1.0f }, @@ -1811,7 +1835,7 @@ namespace netgen }; #else - static float tetcols[][4] = + static float tetcols[][4] = { { 1.0f, 1.0f, 0.0f, 1.0f }, { 1.0f, 0.0f, 0.0f, 1.0f }, @@ -1847,8 +1871,8 @@ namespace netgen if (!shownode.Test(el[j])) drawtet = 0; if (!drawtet) continue; - - + + int ind = el.GetIndex() % 4; @@ -1868,11 +1892,11 @@ namespace netgen } #else glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, tetcols[ind]); -#endif +#endif + - if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) { const ELEMENT_FACE * faces = MeshTopology :: GetFaces (TET); @@ -1882,51 +1906,51 @@ namespace netgen Point<3> grid[11][11]; Point<3> fpts[3]; int order = vispar.subdivisions+1; - + for (int trig = 0; trig < 4; trig++) - { + { for (int j = 0; j < 3; j++) fpts[j] = vertices[faces[trig][j]-1]; - + static Point<3> c(0.25, 0.25, 0.25); if (vispar.shrink < 1) for (int j = 0; j < 3; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[3] = + double lami[3] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), double(iy)/order }; - + Point<3> xl; for (int l = 0; l < 3; l++) xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l); - + curv.CalcElementTransformation (xl, i-1, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(8, 0.0, 0.999, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } */ @@ -1936,30 +1960,30 @@ namespace netgen Array > ploc ( (order+1)*(order+1) ); Array > pglob ( (order+1)*(order+1) ); Point<3> fpts[3]; - + for (int trig = 0; trig < 4; trig++) - { + { for (int j = 0; j < 3; j++) fpts[j] = vertices[faces[trig][j]-1]; - + static Point<3> c(0.25, 0.25, 0.25); if (vispar.shrink < 1) for (int j = 0; j < 3; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0, ii = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++, ii++) { - double lami[3] = + double lami[3] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), double(iy)/order }; - + Point<3> xl; for (int l = 0; l < 3; l++) xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l); - + ploc[ii] = xl; } @@ -1969,25 +1993,25 @@ namespace netgen for (int ix = 0, ii = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++, ii++) grid[ix][iy] = pglob[ii]; - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(hoplotn, 0.0, 0.9999, hoplotn, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, hoplotn, 0, hoplotn); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } } else @@ -2031,7 +2055,7 @@ namespace netgen } } } - + glEndList (); } @@ -2061,11 +2085,11 @@ namespace netgen glNewList (prismlist, GL_COMPILE); static float prismcol[] = { 0.0f, 1.0f, 1.0f, 1.0f }; - glLineWidth (1.0f); + glLineWidth (1.0f); Array faces; - + glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, prismcol); for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++) @@ -2081,109 +2105,109 @@ namespace netgen { const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PRISM); const Point3d * vertices = MeshTopology :: GetVertices (PRISM); - + Point<3> grid[11][11]; Point<3> fpts[4]; int order = vispar.subdivisions+1; - + for (int trig = 0; trig < 2; trig++) - { + { for (int j = 0; j < 3; j++) fpts[j] = vertices[faces[trig][j]-1]; - + static Point<3> c(1.0/3.0, 1.0/3.0, 0.5); if (vispar.shrink < 1) for (int j = 0; j < 3; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[3] = + double lami[3] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), double(iy)/order }; - + Point<3> xl; for (int l = 0; l < 3; l++) xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l); - + curv.CalcElementTransformation (xl, i-1, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(8, 0.0, 0.999, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } for (int quad = 2; quad < 5; quad++) - { + { for (int j = 0; j < 4; j++) fpts[j] = vertices[faces[quad][j]-1]; - + static Point<3> c(1.0/3.0, 1.0/3.0, 0.5); if (vispar.shrink < 1) for (int j = 0; j < 4; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[4] = + double lami[4] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * ( double(iy)/order), (1-double(ix)/order) * ( double(iy)/order) }; - + Point<3> xl; for (int l = 0; l < 3; l++) - xl(l) = + xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l) + lami[3] * fpts[3](l); - + curv.CalcElementTransformation (xl, ei, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } /* - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; // int hoplotn = curv.GetNVisualSubsecs(); const Point3d * facepoint = MeshTopology :: GetVertices (TRIG); @@ -2193,7 +2217,7 @@ namespace netgen for (int trig = 0; trig<2; trig++) { - + Vec<3> x0,x1,d0,d1; x0 = facepoint[1] - facepoint[2]; x1 = facepoint[0] - facepoint[2]; @@ -2221,7 +2245,7 @@ namespace netgen xr[1](0) = (double)(i1+1)/hoplotn; xr[1](1) = (double)(j1+1)/hoplotn; xr[2](0) = (double) i1/hoplotn; xr[2](1) = (double)(j1+1)/hoplotn; }; - + for (int l=0; l<3; l++) { Mat<3,3> dxdxi; @@ -2241,7 +2265,7 @@ namespace netgen glVertex3d (xg(0), xg(1), xg(2)); } } - + } glEnd (); @@ -2249,7 +2273,7 @@ namespace netgen glBegin (GL_QUADS); for (int quad = 0; quad<3; quad++) - { + { const Point3d * facepoint = MeshTopology :: GetVertices (PRISM); Vec<3> x0,x1; @@ -2290,7 +2314,7 @@ namespace netgen xr[1](xyz) = (double)(i1+1)/hoplotn; xr[1](2) = (double) j1/hoplotn; xr[2](xyz) = (double)(i1+1)/hoplotn; xr[2](2) = (double)(j1+1)/hoplotn; xr[3](xyz) = (double) i1/hoplotn; xr[3](2) = (double)(j1+1)/hoplotn; - + for (int l=0; l<4; l++) { switch (quad) @@ -2319,9 +2343,9 @@ namespace netgen } glEnd (); */ - } + } else - { + { Point3d c(0,0,0); if (vispar.shrink < 1) { @@ -2355,14 +2379,14 @@ namespace netgen glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); glVertex3d (lp3.X(), lp3.Y(), lp3.Z()); } - + glEnd(); } } } glEndList (); } - + @@ -2387,11 +2411,11 @@ namespace netgen static float hexcol[] = { 1.0f, 1.0f, 0.0f, 1.0f }; - glLineWidth (1.0f); + glLineWidth (1.0f); glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, hexcol); Array faces; - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++) { @@ -2401,9 +2425,9 @@ namespace netgen CurvedElements & curv = mesh->GetCurvedElements(); if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) { - /* // classical + /* // classical glBegin (GL_QUADS); - + const ELEMENT_FACE * faces = MeshTopology :: GetFaces (HEX); const Point3d * vertices = MeshTopology :: GetVertices (HEX); @@ -2411,7 +2435,7 @@ namespace netgen Vec<3> gridn[33][33]; Point<3> fpts[4]; for (int quad = 0; quad<6; quad++) - { + { for (int j = 0; j < 4; j++) fpts[j] = vertices[faces[quad][j]-1]; @@ -2428,7 +2452,7 @@ namespace netgen { Point<3> xl; Mat<3,3> dxdxi; - double lami[4] = + double lami[4] = { (1-double(ix)/hoplotn) * (1-double(iy)/hoplotn), ( double(ix)/hoplotn) * (1-double(iy)/hoplotn), ( double(ix)/hoplotn) * ( double(iy)/hoplotn), @@ -2438,12 +2462,12 @@ namespace netgen lami[2] * fpts[2](l) + lami[3] * fpts[3](l); curv.CalcElementTransformation (xl, ei, grid[ix][iy], dxdxi); - + Vec<3> gtaux = dxdxi * taux; Vec<3> gtauy = dxdxi * tauy; gridn[ix][iy] = Cross (gtauy, gtaux).Normalize(); } - + for (int ix = 0; ix < hoplotn; ix++) for (int iy = 0; iy < hoplotn; iy++) { @@ -2460,7 +2484,7 @@ namespace netgen glVertex3dv (grid[ix][iy+1]); } } - + glEnd (); */ @@ -2472,7 +2496,7 @@ namespace netgen int order = vispar.subdivisions+1; for (int quad = 0; quad<6; quad++) - { + { for (int j = 0; j < 4; j++) fpts[j] = vertices[faces[quad][j]-1]; @@ -2484,7 +2508,7 @@ namespace netgen for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[4] = + double lami[4] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * ( double(iy)/order), @@ -2497,18 +2521,18 @@ namespace netgen curv.CalcElementTransformation (xl, ei, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - glMap2d(GL_MAP2_VERTEX_3, + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); + glEnable(GL_AUTO_NORMAL); glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); @@ -2516,9 +2540,9 @@ namespace netgen glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); } - } + } else - { + { Point3d c(0,0,0); if (vispar.shrink < 1) { @@ -2543,10 +2567,10 @@ namespace netgen Point<3> lp1 = mesh->Point (el.PNum(face.PNum(1))); Point<3> lp2 = mesh->Point (el.PNum(face.PNum(2))); Point<3> lp3 = mesh->Point (el.PNum(face.PNum(3))); - Vec<3> n = Cross (lp3-lp1, lp2-lp1); + Vec<3> n = Cross (lp3-lp1, lp2-lp1); n.Normalize(); glNormal3dv (n); - + if (vispar.shrink < 1) { lp1 = c + vispar.shrink * (lp1 - c); @@ -2573,7 +2597,7 @@ namespace netgen - + void VisualSceneMesh :: BuildPyramidList() { if (pyramidtimestamp > mesh->GetTimeStamp () && @@ -2595,7 +2619,7 @@ namespace netgen pyramidlist = glGenLists (1); glNewList (pyramidlist, GL_COMPILE); - + static float pyramidcol[] = { 1.0f, 0.0f, 1.0f, 1.0f }; glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, pyramidcol); @@ -2616,102 +2640,102 @@ namespace netgen const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PYRAMID); const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID); - + Point<3> grid[11][11]; Point<3> fpts[4]; int order = vispar.subdivisions+1; - + for (int trig = 0; trig < 4; trig++) - { + { for (int j = 0; j < 3; j++) fpts[j] = vertices[faces[trig][j]-1]; - + static Point<3> c(0.375, 0.375, 0.25); if (vispar.shrink < 1) for (int j = 0; j < 3; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[3] = + double lami[3] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), double(iy)/order }; - + Point<3> xl; for (int l = 0; l < 3; l++) xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l); - + curv.CalcElementTransformation (xl, i-1, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(8, 0.0, 0.999, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } for (int quad = 4; quad < 5; quad++) - { + { for (int j = 0; j < 4; j++) fpts[j] = vertices[faces[quad][j]-1]; - + static Point<3> c(0.375, 0.375, 0.25); if (vispar.shrink < 1) for (int j = 0; j < 4; j++) fpts[j] += (1-vispar.shrink) * (c-fpts[j]); - + for (int ix = 0; ix <= order; ix++) for (int iy = 0; iy <= order; iy++) { - double lami[4] = + double lami[4] = { (1-double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * (1-double(iy)/order), ( double(ix)/order) * ( double(iy)/order), (1-double(ix)/order) * ( double(iy)/order) }; - + Point<3> xl; for (int l = 0; l < 3; l++) - xl(l) = + xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) + lami[2] * fpts[2](l) + lami[3] * fpts[3](l); - + curv.CalcElementTransformation (xl, ei, grid[ix][iy]); } - + for (int j = 0; j <= order; j++) ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]); for (int j = 0; j <= order; j++) ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]); - - glMap2d(GL_MAP2_VERTEX_3, + + glMap2d(GL_MAP2_VERTEX_3, 0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1, 0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1, &grid[0][0](0)); glEnable(GL_MAP2_VERTEX_3); - glEnable(GL_AUTO_NORMAL); - + glEnable(GL_AUTO_NORMAL); + glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0); glEvalMesh2(GL_FILL, 0, 8, 0, 8); - + glDisable (GL_AUTO_NORMAL); glDisable (GL_MAP2_VERTEX_3); - } + } @@ -2719,7 +2743,7 @@ namespace netgen /* - int hoplotn = 1 << vispar.subdivisions; + int hoplotn = 1 << vispar.subdivisions; const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PYRAMID); const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID); @@ -2729,9 +2753,9 @@ namespace netgen glBegin (GL_TRIANGLES); - + for (int trig = 0; trig < 4; trig++) - { + { Point<3> p0 = vertices[faces[trig][0]-1]; Point<3> p1 = vertices[faces[trig][1]-1]; Point<3> p2 = vertices[faces[trig][2]-1]; @@ -2743,7 +2767,7 @@ namespace netgen p1 = c + vispar.shrink * (p1 - c); p2 = c + vispar.shrink * (p2 - c); } - + Vec<3> taux = p0-p2; Vec<3> tauy = p1-p2; @@ -2756,13 +2780,13 @@ namespace netgen for (int iy = 0; iy <= hoplotn-ix; iy++) { for (int l = 0; l < 3; l++) - xl(l) = + xl(l) = (1-double(ix+iy)/hoplotn) * p2(l) + (double(ix)/hoplotn) * p0(l) + (double(iy)/hoplotn) * p1(l); - + curv.CalcElementTransformation (xl, i-1, grid[ix][iy], dxdxi); - + gtaux = dxdxi * taux; gtauy = dxdxi * tauy; gridn[ix][iy] = Cross (gtauy, gtaux).Normalize(); @@ -2784,7 +2808,7 @@ namespace netgen { glNormal3dv (gridn[ix][iy+1]); glVertex3dv (grid[ix][iy+1]); - + glNormal3dv (gridn[ix+1][iy]); glVertex3dv (grid[ix+1][iy]); @@ -2800,9 +2824,9 @@ namespace netgen glBegin (GL_QUADS); - + for (int quad = 4; quad < 5; quad++) - { + { Point<3> p0 = vertices[faces[quad][0]-1]; Point<3> p1 = vertices[faces[quad][1]-1]; Point<3> p2 = vertices[faces[quad][2]-1]; @@ -2829,14 +2853,14 @@ namespace netgen { Point<3> xl; for (int l = 0; l < 3; l++) - xl(l) = + xl(l) = (1-double(ix)/hoplotn)*(1-double(iy)/hoplotn) * p0(l) + ( double(ix)/hoplotn)*(1-double(iy)/hoplotn) * p1(l) + ( double(ix)/hoplotn)*( double(iy)/hoplotn) * p2(l) + (1-double(ix)/hoplotn)*( double(iy)/hoplotn) * p3(l); - + curv.CalcElementTransformation (xl, i-1, grid[ix][iy], dxdxi); - + gtaux = dxdxi * taux; gtauy = dxdxi * tauy; gridn[ix][iy] = Cross (gtauy, gtaux).Normalize(); @@ -2863,9 +2887,9 @@ namespace netgen */ - } + } else - { + { @@ -2883,11 +2907,11 @@ namespace netgen el.GetSurfaceTriangles (faces); - + if (el.PNum(1)) { glBegin (GL_TRIANGLES); - + for (int j = 1; j <= faces.Size(); j++) { Element2d & face = faces.Elem(j); @@ -2910,7 +2934,7 @@ namespace netgen glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); glVertex3d (lp3.X(), lp3.Y(), lp3.Z()); } - + glEnd(); } } @@ -2928,7 +2952,7 @@ namespace netgen { ; } - + void VisualSceneMesh :: BuildDomainSurfList() { if (domainsurflist) @@ -2939,23 +2963,23 @@ namespace netgen int i, j; glLineWidth (1.0f); - + glDisable (GL_COLOR_MATERIAL); - + for (i = 1; i <= mesh->GetNSE(); i++) { Element2d el = mesh->SurfaceElement (i); - + int drawel = 1; for (j = 1; j <= el.GetNP(); j++) { if (!el.PNum(j)) drawel = 0; } - + if (!drawel) continue; - + if (el.GetIndex() < 1 || el.GetIndex() > mesh->GetNFD()) continue; int domin = mesh->GetFaceDescriptor(el.GetIndex()).DomainIn(); @@ -2968,23 +2992,23 @@ namespace netgen fac = -1; else continue; - - + + GLfloat matcol[] = { 1, 0, 0, 1 }; glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matcol); - - + + if (el.GetNP() == 3) { glBegin (GL_TRIANGLES); - + const Point3d & lp1 = mesh->Point (el.PNum(1)); const Point3d & lp2 = mesh->Point (el.PNum(2)); const Point3d & lp3 = mesh->Point (el.PNum(3)); Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, lp3)); n /= ( fac * (n.Length()+1e-12)); glNormal3d (n.X(), n.Y(), n.Z()); - + if (!vispar.colormeshsize) { glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); @@ -2996,15 +3020,15 @@ namespace netgen else if (el.GetNP() == 4) { glBegin (GL_QUADS); - + const Point3d & lp1 = mesh->Point (el.PNum(1)); const Point3d & lp2 = mesh->Point (el.PNum(2)); const Point3d & lp3 = mesh->Point (el.PNum(4)); const Point3d & lp4 = mesh->Point (el.PNum(3)); - Vec3d n = Cross (Vec3d (lp1, lp2), + Vec3d n = Cross (Vec3d (lp1, lp2), Vec3d (lp1, Center (lp3, lp4))); n /= (fac * (n.Length()+1e-12)); - glNormal3d (n.X(), n.Y(), n.Z()); + glNormal3d (n.X(), n.Y(), n.Z()); glVertex3d (lp1.X(), lp1.Y(), lp1.Z()); glVertex3d (lp2.X(), lp2.Y(), lp2.Z()); glVertex3d (lp4.X(), lp4.Y(), lp4.Z()); @@ -3019,7 +3043,7 @@ namespace netgen { 2, 4, 6 }, { 3, 5, 4 }, { 4, 5, 6 } }; - + for (j = 0; j < 4; j++) { const Point3d & lp1 = mesh->Point (el.PNum(trigs[j][0])); @@ -3053,14 +3077,14 @@ namespace netgen delete lock; lock = NULL; } - + MouseDblClickSelect(px,py,clipplane,backcolor,transformationmat,center,rad, filledlist,selelement,selface,seledge,selpoint,selpoint2,locpi); selecttimestamp = NextTimeStamp(); - + /* int i, hits; @@ -3076,22 +3100,22 @@ namespace netgen glGetIntegerv (GL_VIEWPORT, viewport); - glMatrixMode (GL_PROJECTION); + glMatrixMode (GL_PROJECTION); glPushMatrix(); GLdouble projmat[16]; glGetDoublev (GL_PROJECTION_MATRIX, projmat); - glLoadIdentity(); - gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); + glLoadIdentity(); + gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); glMultMatrixd (projmat); - + glClearColor(backcolor, backcolor, backcolor, 1.0); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glMatrixMode (GL_MODELVIEW); + glMatrixMode (GL_MODELVIEW); glPushMatrix(); glMultMatrixf (transformationmat); @@ -3106,7 +3130,7 @@ namespace netgen glEnable (GL_POLYGON_OFFSET_FILL); glDisable(GL_CLIP_PLANE0); - + if (vispar.clipenable) { Vec<3> n(clipplane[0], clipplane[1], clipplane[2]); @@ -3116,11 +3140,11 @@ namespace netgen n /= len; Vec<3> t1 = n.GetNormal (); Vec<3> t2 = Cross (n, t1); - + double xi1mid = (center - p) * t1; double xi2mid = (center - p) * t2; - - glLoadName (0); + + glLoadName (0); glBegin (GL_QUADS); glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid-rad) * t2); glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid-rad) * t2); @@ -3134,18 +3158,18 @@ namespace netgen glCallList (filledlist); glDisable (GL_POLYGON_OFFSET_FILL); - + glPopName(); - glMatrixMode (GL_PROJECTION); + glMatrixMode (GL_PROJECTION); glPopMatrix(); - glMatrixMode (GL_MODELVIEW); + glMatrixMode (GL_MODELVIEW); glPopMatrix(); - glFlush(); + glFlush(); + - hits = glRenderMode (GL_RENDER); // cout << "hits = " << hits << endl; @@ -3194,21 +3218,21 @@ namespace netgen locpi = (locpi % sel.GetNP()) + 1; selpoint2 = selpoint; selpoint = sel.PNum(locpi); - cout << "selected point " << selpoint - << ", pos = " << mesh->Point (selpoint) + cout << "selected point " << selpoint + << ", pos = " << mesh->Point (selpoint) << endl; for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - if (seg.p1 == selpoint && seg.p2 == selpoint2 || + if (seg.p1 == selpoint && seg.p2 == selpoint2 || seg.p2 == selpoint && seg.p1 == selpoint2) { seledge = seg.edgenr; cout << "seledge = " << seledge << endl; } } - + } else { @@ -3252,22 +3276,22 @@ namespace netgen glGetIntegerv (GL_VIEWPORT, viewport); - glMatrixMode (GL_PROJECTION); + glMatrixMode (GL_PROJECTION); glPushMatrix(); GLdouble projmat[16]; glGetDoublev (GL_PROJECTION_MATRIX, projmat); - glLoadIdentity(); - gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); + glLoadIdentity(); + gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); glMultMatrixd (projmat); - + glClearColor(backcolor, backcolor, backcolor, 1.0); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glMatrixMode (GL_MODELVIEW); + glMatrixMode (GL_MODELVIEW); glPushMatrix(); glMultMatrixf (transformationmat); @@ -3282,7 +3306,7 @@ namespace netgen glEnable (GL_POLYGON_OFFSET_FILL); glDisable(GL_CLIP_PLANE0); - + if (vispar.clipenable) { Vec<3> n(clipplane[0], clipplane[1], clipplane[2]); @@ -3292,11 +3316,11 @@ namespace netgen n /= len; Vec<3> t1 = n.GetNormal (); Vec<3> t2 = Cross (n, t1); - + double xi1mid = (center - p) * t1; double xi2mid = (center - p) * t2; - - glLoadName (0); + + glLoadName (0); glBegin (GL_QUADS); glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid-rad) * t2); glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid-rad) * t2); @@ -3310,18 +3334,18 @@ namespace netgen glCallList (displaylist); glDisable (GL_POLYGON_OFFSET_FILL); - + glPopName(); - glMatrixMode (GL_PROJECTION); + glMatrixMode (GL_PROJECTION); glPopMatrix(); - glMatrixMode (GL_MODELVIEW); + glMatrixMode (GL_MODELVIEW); glPopMatrix(); - glFlush(); + glFlush(); + - hits = glRenderMode (GL_RENDER); //cout << "hits = " << hits << endl; @@ -3343,7 +3367,7 @@ namespace netgen int curname = selbuf[4*i+3]; GLuint curdepth = selbuf[4*i+1]; /* - cout << selbuf[4*i] << " " << selbuf[4*i+1] << " " + cout << selbuf[4*i] << " " << selbuf[4*i+1] << " " << selbuf[4*i+2] << " " << selbuf[4*i+3] << endl; */ if (curname && (curdepth > clipdepth) && @@ -3373,21 +3397,21 @@ namespace netgen locpi = (locpi % sel.GetNP()) + 1; selpoint2 = selpoint; selpoint = sel.PNum(locpi); - cout << "selected point " << selpoint - << ", pos = " << mesh->Point (selpoint) + cout << "selected point " << selpoint + << ", pos = " << mesh->Point (selpoint) << endl; - + for (i = 1; i <= mesh->GetNSeg(); i++) { const Segment & seg = mesh->LineSegment(i); - if (seg.p1 == selpoint && seg.p2 == selpoint2 || + if (seg.p1 == selpoint && seg.p2 == selpoint2 || seg.p2 == selpoint && seg.p1 == selpoint2) { seledge = seg.edgenr; cout << "seledge = " << seledge << endl; } } - + } else { @@ -3403,14 +3427,14 @@ namespace netgen #ifdef PARALLELGL vsmesh.Broadcast (); -#endif +#endif } void VisualSceneMesh :: SetSelectedFace (int asf) - { - selface = asf; - selecttimestamp = GetTimeStamp(); + { + selface = asf; + selecttimestamp = GetTimeStamp(); } @@ -3420,4 +3444,4 @@ namespace netgen -#endif // NOTCL +#endif // NOTCL diff --git a/libsrc/visualization/vsocc.cpp b/libsrc/visualization/vsocc.cpp index f2905830..80877191 100644 --- a/libsrc/visualization/vsocc.cpp +++ b/libsrc/visualization/vsocc.cpp @@ -2,7 +2,6 @@ #ifdef OCCGEOMETRY - #include #include #include @@ -30,723 +29,725 @@ #include "incvis.hpp" - namespace netgen { #include "mvdraw.hpp" + extern OCCGeometry * occgeometry; -extern OCCGeometry * occgeometry; + /* *********************** Draw OCC Geometry **************** */ + VisualSceneOCCGeometry :: VisualSceneOCCGeometry () + : VisualScene() + { + trilists.SetSize(0); + linelists.SetSize(1); + } + VisualSceneOCCGeometry :: ~VisualSceneOCCGeometry () + { + ; + } -/* *********************** Draw OCC Geometry **************** */ - - -VisualSceneOCCGeometry :: VisualSceneOCCGeometry () - : VisualScene() -{ - trilists.SetSize(0); - linelists.SetSize(1); - -} - -VisualSceneOCCGeometry :: ~VisualSceneOCCGeometry () -{ - ; -} - -void VisualSceneOCCGeometry :: DrawScene () -{ - if ( occgeometry->changed ) - { - BuildScene(); - occgeometry -> changed = 0; - } - - glClearColor(backcolor, backcolor, backcolor, 1.0); - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - SetLight(); - - glPushMatrix(); - glMultMatrixf (transformationmat); - - glShadeModel (GL_SMOOTH); - glDisable (GL_COLOR_MATERIAL); - glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); - - glEnable (GL_BLEND); - glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); - - // glEnable (GL_LIGHTING); - - double shine = vispar.shininess; - // double transp = vispar.transp; - - glMaterialf (GL_FRONT_AND_BACK, GL_SHININESS, shine); - glLogicOp (GL_COPY); - - - float mat_col[] = { 0.2f, 0.2f, 0.8f, 1.0f }; - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); - - glPolygonOffset (1, 1); - glEnable (GL_POLYGON_OFFSET_FILL); - - GLfloat matcoledge[] = { 0, 0, 1, 1 }; - GLfloat matcolhiedge[] = { 1, 0, 0, 1 }; - - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, - matcoledge); - glLineWidth (1.0f); - - if (vispar.occshowedges) glCallList (linelists.Get(1)); - if (vispar.occshowsurfaces) glCallList (trilists.Get(1)); - - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, - matcolhiedge); - glLineWidth (5.0f); - - if (vispar.occshowedges) glCallList (linelists.Get(2)); - - for (int i = 1; i <= occgeometry->vmap.Extent(); i++) - if (occgeometry->vvispar[i-1].IsHighlighted()) + void VisualSceneOCCGeometry :: DrawScene () + { + if ( occgeometry->changed ) { - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, - matcolhiedge); - glLineWidth (5.0f); - - glBegin (GL_LINES); - - gp_Pnt p = BRep_Tool::Pnt(TopoDS::Vertex(occgeometry->vmap(i))); - double d = rad/100; - glVertex3f (p.X()-d, p.Y(), p.Z()); - glVertex3f (p.X()+d, p.Y(), p.Z()); - glVertex3f (p.X(), p.Y()-d, p.Z()); - glVertex3f (p.X(), p.Y()+d, p.Z()); - glVertex3f (p.X(), p.Y(), p.Z()-d); - glVertex3f (p.X(), p.Y(), p.Z()+d); - glEnd(); + BuildScene(); + occgeometry -> changed = 0; } + glClearColor(backcolor, backcolor, backcolor, 1.0); + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + SetLight(); - glDisable (GL_POLYGON_OFFSET_FILL); - - glPopMatrix(); - // DrawCoordinateCross (); - // DrawNetgenLogo (); - glFinish(); + glPushMatrix(); + glMultMatrixf (transformationmat); - glDisable (GL_POLYGON_OFFSET_FILL); -} + glShadeModel (GL_SMOOTH); + glDisable (GL_COLOR_MATERIAL); + glPolygonMode (GL_FRONT_AND_BACK, GL_FILL); + glEnable (GL_BLEND); + glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); -/* -void VisualSceneOCCGeometry :: BuildScene (int zoomall) -{ - int i = 0, j, k; - - TopExp_Explorer ex, ex_edge; + // glEnable (GL_LIGHTING); - if (vispar.occvisproblemfaces || (occgeometry -> changed != 2)) - { - Box<3> bb = occgeometry -> GetBoundingBox(); + double shine = vispar.shininess; + // double transp = vispar.transp; - center = bb.Center(); - rad = bb.Diam() / 2; - - - - if (vispar.occvisproblemfaces) - { - for (i = 1; i <= occgeometry->fmap.Extent(); i++) - if (occgeometry->facemeshstatus[i-1] == -1) - { - GProp_GProps system; - BRepGProp::LinearProperties(occgeometry->fmap(i), system); - gp_Pnt pnt = system.CentreOfMass(); - center = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); - cout << "Setting center to mid of face " << i << " = " << center << endl; - } - } - - - CalcTransformationMatrices(); - } - - - for (i = 1; i <= linelists.Size(); i++) - glDeleteLists (linelists.Elem(i), 1); - linelists.SetSize(0); - - linelists.Append (glGenLists (1)); - glNewList (linelists.Last(), GL_COMPILE); - - i = 0; - for (ex_edge.Init(occgeometry -> shape, TopAbs_EDGE); - ex_edge.More(); ex_edge.Next()) - { - if (BRep_Tool::Degenerated(TopoDS::Edge(ex_edge.Current()))) continue; - i++; - - - TopoDS_Edge edge = TopoDS::Edge(ex_edge.Current()); - - Handle(Poly_PolygonOnTriangulation) aEdgePoly; - Handle(Poly_Triangulation) T; - TopLoc_Location aEdgeLoc; - BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); - - if(aEdgePoly.IsNull()) - { - cout << "cannot visualize edge " << i << endl; - continue; - } - - glBegin (GL_LINE_STRIP); - - int nbnodes = aEdgePoly -> NbNodes(); - for (j = 1; j <= nbnodes; j++) - { - gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); - glVertex3f (p.X(), p.Y(), p.Z()); - } - - glEnd (); - - - } - - glEndList (); - - for (i = 1; i <= trilists.Size(); i++) - glDeleteLists (trilists.Elem(i), 1); - trilists.SetSize(0); - - - trilists.Append (glGenLists (1)); - glNewList (trilists.Last(), GL_COMPILE); - - i = 0; - - TopExp_Explorer exp0, exp1, exp2, exp3; - int shapenr = 0; - for (exp0.Init(occgeometry -> shape, TopAbs_SOLID); exp0.More(); exp0.Next()) - { - shapenr++; - - if (vispar.occshowvolumenr != 0 && - vispar.occshowvolumenr != shapenr) continue; - - float mat_col[4]; - mat_col[3] = 1; - switch (shapenr) - { - case 1: - mat_col[0] = 0.2; - mat_col[1] = 0.2; - mat_col[2] = 0.8; - break; - case 2: - mat_col[0] = 0.8; - mat_col[1] = 0.2; - mat_col[2] = 0.8; - break; - case 3: - mat_col[0] = 0.2; - mat_col[1] = 0.8; - mat_col[2] = 0.8; - break; - case 4: - mat_col[0] = 0.8; - mat_col[1] = 0.2; - mat_col[2] = 0.2; - break; - case 5: - mat_col[0] = 0.8; - mat_col[1] = 0.8; - mat_col[2] = 0.8; - break; - case 6: - mat_col[0] = 0.6; - mat_col[1] = 0.6; - mat_col[2] = 0.6; - break; - case 7: - mat_col[0] = 0.2; - mat_col[1] = 0.8; - mat_col[2] = 0.2; - break; - case 8: - mat_col[0] = 0.8; - mat_col[1] = 0.8; - mat_col[2] = 0.2; - break; - default: - // mat_col[0] = 1-(1.0/double(shapenr)); - // mat_col[1] = 0.5; - mat_col[0] = 0.5+double((shapenr*shapenr*shapenr*shapenr) % 10)/20.0; - mat_col[1] = 0.5+double(int(shapenr*shapenr*shapenr*shapenr*sin(double(shapenr))) % 10)/20.0; - mat_col[2] = 0.5+double((shapenr*shapenr*shapenr) % 10)/20.0; - } + glMaterialf (GL_FRONT_AND_BACK, GL_SHININESS, shine); + glLogicOp (GL_COPY); + float mat_col[] = { 0.2f, 0.2f, 0.8f, 1.0f}; glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); - for (exp1.Init(exp0.Current(), TopAbs_SHELL); exp1.More(); exp1.Next()) - for (exp2.Init(exp1.Current().Composed(exp0.Current().Orientation()), TopAbs_FACE); exp2.More(); exp2.Next()) - { - TopoDS_Face face = TopoDS::Face (exp2.Current().Composed(exp1.Current().Orientation())); - - i = occgeometry->fmap.FindIndex(face); + glPolygonOffset (1, 1); + glEnable (GL_POLYGON_OFFSET_FILL); - TopLoc_Location loc; - Handle(Geom_Surface) surf = BRep_Tool::Surface (face); - BRepAdaptor_Surface sf(face, Standard_False); - BRepLProp_SLProps prop(sf, 1, 1e-5); - Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); - - if (triangulation.IsNull()) - { - cout << "cannot visualize face " << i << endl; - continue; - } + // Philippose - 30/01/2009 + // Added clipping planes to Geometry view + SetClippingPlane(); - if (vispar.occvisproblemfaces) - { - switch (occgeometry->facemeshstatus[i-1]) - { - case 0: - mat_col[0] = 0.2; - mat_col[1] = 0.2; - mat_col[2] = 0.8; - break; - case 1: - mat_col[0] = 0.2; - mat_col[1] = 0.8; - mat_col[2] = 0.2; - break; - case -1: - mat_col[0] = 0.8; - mat_col[1] = 0.2; - mat_col[2] = 0.2; - break; - } - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); + GLfloat matcoledge[] = { 0, 0, 1, 1}; + GLfloat matcolhiedge[] = { 1, 0, 0, 1}; - } - glBegin (GL_TRIANGLES); - - int ntriangles = triangulation -> NbTriangles(); - for (j = 1; j <= ntriangles; j++) - { - Poly_Triangle triangle = (triangulation -> Triangles())(j); - for (k = 1; k <= 3; k++) - { - gp_Pnt2d uv = (triangulation -> UVNodes())(triangle(k)); - gp_Pnt pnt; - gp_Vec du, dv; - prop.SetParameters (uv.X(), uv.Y()); - surf->D0 (uv.X(), uv.Y(), pnt); - gp_Vec n; - - if (prop.IsNormalDefined()) - n = prop.Normal(); - else - n = gp_Vec (0,0,0); - - if (face.Orientation() == TopAbs_REVERSED) n *= -1; - glNormal3f (n.X(), n.Y(), n.Z()); - glVertex3f (pnt.X(), pnt.Y(), pnt.Z()); - } - } - glEnd (); + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcoledge); + glLineWidth (1.0f); - } - } + if (vispar.occshowedges) glCallList (linelists.Get(1)); + if (vispar.occshowsurfaces) glCallList (trilists.Get(1)); + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolhiedge); + glLineWidth (5.0f); - glEndList (); + if (vispar.occshowedges) glCallList (linelists.Get(2)); -} -*/ - - -void VisualSceneOCCGeometry :: BuildScene (int zoomall) -{ - if (occgeometry -> changed == OCCGEOMETRYVISUALIZATIONFULLCHANGE) - { - center = occgeometry -> Center(); - rad = occgeometry -> GetBoundingBox().Diam() / 2; - - if (vispar.occzoomtohighlightedentity) - { - bool hilite = false; - bool hiliteonepoint = false; - Bnd_Box bb; - - for (int i = 1; i <= occgeometry->fmap.Extent(); i++) - if (occgeometry->fvispar[i-1].IsHighlighted()) - { - hilite = true; - BRepBndLib::Add (occgeometry->fmap(i), bb); - } - - for (int i = 1; i <= occgeometry->emap.Extent(); i++) - if (occgeometry->evispar[i-1].IsHighlighted()) - { - hilite = true; - BRepBndLib::Add (occgeometry->emap(i), bb); - } - - for (int i = 1; i <= occgeometry->vmap.Extent(); i++) - if (occgeometry->vvispar[i-1].IsHighlighted()) - { - hiliteonepoint = true; - BRepBndLib::Add (occgeometry->vmap(i), bb); - } - - if (hilite || hiliteonepoint) - { - double x1,y1,z1,x2,y2,z2; - bb.Get (x1,y1,z1,x2,y2,z2); - Point<3> p1 = Point<3> (x1,y1,z1); - Point<3> p2 = Point<3> (x2,y2,z2); - Box<3> boundingbox(p1,p2); - - center = boundingbox.Center(); - if (hiliteonepoint) - rad = occgeometry -> GetBoundingBox().Diam() / 100; - else - rad = boundingbox.Diam() / 2; - } - } - - CalcTransformationMatrices(); - } - - - // Clear lists - - for (int i = 1; i <= linelists.Size(); i++) - glDeleteLists (linelists.Elem(i), 1); - linelists.SetSize(0); - - for (int i = 1; i <= trilists.Size(); i++) - glDeleteLists (trilists.Elem(i), 1); - trilists.SetSize(0); - - - // Total wireframe - - linelists.Append (glGenLists (1)); - glNewList (linelists.Last(), GL_COMPILE); - - for (int i = 1; i <= occgeometry->emap.Extent(); i++) - { - TopoDS_Edge edge = TopoDS::Edge(occgeometry->emap(i)); - if (BRep_Tool::Degenerated(edge)) continue; - if (occgeometry->evispar[i-1].IsHighlighted()) continue; - - Handle(Poly_PolygonOnTriangulation) aEdgePoly; - Handle(Poly_Triangulation) T; - TopLoc_Location aEdgeLoc; - BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); - - if(aEdgePoly.IsNull()) - { - (*testout) << "visualizing edge " << occgeometry->emap.FindIndex (edge) - << " without using the occ visualization triangulation" << endl; - - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - - glBegin (GL_LINE_STRIP); - for (int i = 0; i<=50; i++) - { - gp_Pnt p = c->Value (s0 + i*(s1-s0)/50.0); - glVertex3f (p.X(),p.Y(),p.Z()); - } - glEnd (); - - continue; - } - - int nbnodes = aEdgePoly -> NbNodes(); - glBegin (GL_LINE_STRIP); - for (int j = 1; j <= nbnodes; j++) - { - gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); - glVertex3f (p.X(), p.Y(), p.Z()); - } - glEnd (); - } - - glEndList (); - - - // Highlighted edge list - - linelists.Append (glGenLists (1)); - glNewList (linelists.Last(), GL_COMPILE); - - for (int i = 1; i <= occgeometry->emap.Extent(); i++) - if (occgeometry->evispar[i-1].IsHighlighted()) + for (int i = 1; i <= occgeometry->vmap.Extent(); i++) + if (occgeometry->vvispar[i-1].IsHighlighted()) { - TopoDS_Edge edge = TopoDS::Edge(occgeometry->emap(i)); - if (BRep_Tool::Degenerated(edge)) continue; + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcolhiedge); + glLineWidth (5.0f); - Handle(Poly_PolygonOnTriangulation) aEdgePoly; - Handle(Poly_Triangulation) T; - TopLoc_Location aEdgeLoc; - BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); - - if(aEdgePoly.IsNull()) - { - (*testout) << "visualizing edge " << occgeometry->emap.FindIndex (edge) - << " without using the occ visualization triangulation" << endl; + glBegin (GL_LINES); - double s0, s1; - Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); - - glBegin (GL_LINE_STRIP); - for (int i = 0; i<=50; i++) - { - gp_Pnt p = c->Value (s0 + i*(s1-s0)/50.0); - glVertex3f (p.X(),p.Y(),p.Z()); - } - glEnd (); - - continue; - } - - int nbnodes = aEdgePoly -> NbNodes(); - glBegin (GL_LINE_STRIP); - for (int j = 1; j <= nbnodes; j++) - { - gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); - glVertex3f (p.X(), p.Y(), p.Z()); - } - glEnd (); + gp_Pnt p = BRep_Tool::Pnt(TopoDS::Vertex(occgeometry->vmap(i))); + double d = rad/100; + glVertex3f (p.X()-d, p.Y(), p.Z()); + glVertex3f (p.X()+d, p.Y(), p.Z()); + glVertex3f (p.X(), p.Y()-d, p.Z()); + glVertex3f (p.X(), p.Y()+d, p.Z()); + glVertex3f (p.X(), p.Y(), p.Z()-d); + glVertex3f (p.X(), p.Y(), p.Z()+d); + glEnd(); } - - glEndList (); + glDisable (GL_POLYGON_OFFSET_FILL); + glPopMatrix(); + // DrawCoordinateCross (); + // DrawNetgenLogo (); + glFinish(); + glDisable (GL_POLYGON_OFFSET_FILL); + } - - // display faces - - trilists.Append (glGenLists (1)); - glNewList (trilists.Last(), GL_COMPILE); - - for (int i = 1; i <= occgeometry->fmap.Extent(); i++) + /* + void VisualSceneOCCGeometry :: BuildScene (int zoomall) { - if (!occgeometry->fvispar[i-1].IsVisible()) continue; + int i = 0, j, k; - glLoadName (i); - float mat_col[4]; - mat_col[3] = 1; + TopExp_Explorer ex, ex_edge; - if (!occgeometry->fvispar[i-1].IsHighlighted()) - { - mat_col[0] = 0.2; - mat_col[1] = 0.2; - mat_col[2] = 0.8; - } + if (vispar.occvisproblemfaces || (occgeometry -> changed != 2)) + { + Box<3> bb = occgeometry -> GetBoundingBox(); + + center = bb.Center(); + rad = bb.Diam() / 2; + + + + if (vispar.occvisproblemfaces) + { + for (i = 1; i <= occgeometry->fmap.Extent(); i++) + if (occgeometry->facemeshstatus[i-1] == -1) + { + GProp_GProps system; + BRepGProp::LinearProperties(occgeometry->fmap(i), system); + gp_Pnt pnt = system.CentreOfMass(); + center = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); + cout << "Setting center to mid of face " << i << " = " << center << endl; + } + } + + + CalcTransformationMatrices(); + } + + + for (i = 1; i <= linelists.Size(); i++) + glDeleteLists (linelists.Elem(i), 1); + linelists.SetSize(0); + + linelists.Append (glGenLists (1)); + glNewList (linelists.Last(), GL_COMPILE); + + i = 0; + for (ex_edge.Init(occgeometry -> shape, TopAbs_EDGE); + ex_edge.More(); ex_edge.Next()) + { + if (BRep_Tool::Degenerated(TopoDS::Edge(ex_edge.Current()))) continue; + i++; + + + TopoDS_Edge edge = TopoDS::Edge(ex_edge.Current()); + + Handle(Poly_PolygonOnTriangulation) aEdgePoly; + Handle(Poly_Triangulation) T; + TopLoc_Location aEdgeLoc; + BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); + + if(aEdgePoly.IsNull()) + { + cout << "cannot visualize edge " << i << endl; + continue; + } + + glBegin (GL_LINE_STRIP); + + int nbnodes = aEdgePoly -> NbNodes(); + for (j = 1; j <= nbnodes; j++) + { + gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); + glVertex3f (p.X(), p.Y(), p.Z()); + } + + glEnd (); + + + } + + glEndList (); + + for (i = 1; i <= trilists.Size(); i++) + glDeleteLists (trilists.Elem(i), 1); + trilists.SetSize(0); + + + trilists.Append (glGenLists (1)); + glNewList (trilists.Last(), GL_COMPILE); + + i = 0; + + TopExp_Explorer exp0, exp1, exp2, exp3; + int shapenr = 0; + for (exp0.Init(occgeometry -> shape, TopAbs_SOLID); exp0.More(); exp0.Next()) + { + shapenr++; + + if (vispar.occshowvolumenr != 0 && + vispar.occshowvolumenr != shapenr) continue; + + float mat_col[4]; + mat_col[3] = 1; + switch (shapenr) + { + case 1: + mat_col[0] = 0.2; + mat_col[1] = 0.2; + mat_col[2] = 0.8; + break; + case 2: + mat_col[0] = 0.8; + mat_col[1] = 0.2; + mat_col[2] = 0.8; + break; + case 3: + mat_col[0] = 0.2; + mat_col[1] = 0.8; + mat_col[2] = 0.8; + break; + case 4: + mat_col[0] = 0.8; + mat_col[1] = 0.2; + mat_col[2] = 0.2; + break; + case 5: + mat_col[0] = 0.8; + mat_col[1] = 0.8; + mat_col[2] = 0.8; + break; + case 6: + mat_col[0] = 0.6; + mat_col[1] = 0.6; + mat_col[2] = 0.6; + break; + case 7: + mat_col[0] = 0.2; + mat_col[1] = 0.8; + mat_col[2] = 0.2; + break; + case 8: + mat_col[0] = 0.8; + mat_col[1] = 0.8; + mat_col[2] = 0.2; + break; + default: + // mat_col[0] = 1-(1.0/double(shapenr)); + // mat_col[1] = 0.5; + mat_col[0] = 0.5+double((shapenr*shapenr*shapenr*shapenr) % 10)/20.0; + mat_col[1] = 0.5+double(int(shapenr*shapenr*shapenr*shapenr*sin(double(shapenr))) % 10)/20.0; + mat_col[2] = 0.5+double((shapenr*shapenr*shapenr) % 10)/20.0; + } + + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); + + for (exp1.Init(exp0.Current(), TopAbs_SHELL); exp1.More(); exp1.Next()) + for (exp2.Init(exp1.Current().Composed(exp0.Current().Orientation()), TopAbs_FACE); exp2.More(); exp2.Next()) + { + TopoDS_Face face = TopoDS::Face (exp2.Current().Composed(exp1.Current().Orientation())); + + i = occgeometry->fmap.FindIndex(face); + + TopLoc_Location loc; + Handle(Geom_Surface) surf = BRep_Tool::Surface (face); + BRepAdaptor_Surface sf(face, Standard_False); + BRepLProp_SLProps prop(sf, 1, 1e-5); + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); + + if (triangulation.IsNull()) + { + cout << "cannot visualize face " << i << endl; + continue; + } + + if (vispar.occvisproblemfaces) + { + switch (occgeometry->facemeshstatus[i-1]) + { + case 0: + mat_col[0] = 0.2; + mat_col[1] = 0.2; + mat_col[2] = 0.8; + break; + case 1: + mat_col[0] = 0.2; + mat_col[1] = 0.8; + mat_col[2] = 0.2; + break; + case -1: + mat_col[0] = 0.8; + mat_col[1] = 0.2; + mat_col[2] = 0.2; + break; + } + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); + + } + glBegin (GL_TRIANGLES); + + int ntriangles = triangulation -> NbTriangles(); + for (j = 1; j <= ntriangles; j++) + { + Poly_Triangle triangle = (triangulation -> Triangles())(j); + for (k = 1; k <= 3; k++) + { + gp_Pnt2d uv = (triangulation -> UVNodes())(triangle(k)); + gp_Pnt pnt; + gp_Vec du, dv; + prop.SetParameters (uv.X(), uv.Y()); + surf->D0 (uv.X(), uv.Y(), pnt); + gp_Vec n; + + if (prop.IsNormalDefined()) + n = prop.Normal(); + else + n = gp_Vec (0,0,0); + + if (face.Orientation() == TopAbs_REVERSED) n *= -1; + glNormal3f (n.X(), n.Y(), n.Z()); + glVertex3f (pnt.X(), pnt.Y(), pnt.Z()); + } + } + glEnd (); + + } + } + + + glEndList (); + + } + */ + + void VisualSceneOCCGeometry :: BuildScene (int zoomall) + { + if (occgeometry -> changed == OCCGEOMETRYVISUALIZATIONFULLCHANGE) + { + center = occgeometry -> Center(); + rad = occgeometry -> GetBoundingBox().Diam() / 2; + + if (vispar.occzoomtohighlightedentity) + { + bool hilite = false; + bool hiliteonepoint = false; + Bnd_Box bb; + + for (int i = 1; i <= occgeometry->fmap.Extent(); i++) + if (occgeometry->fvispar[i-1].IsHighlighted()) + { + hilite = true; + BRepBndLib::Add (occgeometry->fmap(i), bb); + } + + for (int i = 1; i <= occgeometry->emap.Extent(); i++) + if (occgeometry->evispar[i-1].IsHighlighted()) + { + hilite = true; + BRepBndLib::Add (occgeometry->emap(i), bb); + } + + for (int i = 1; i <= occgeometry->vmap.Extent(); i++) + if (occgeometry->vvispar[i-1].IsHighlighted()) + { + hiliteonepoint = true; + BRepBndLib::Add (occgeometry->vmap(i), bb); + } + + if (hilite || hiliteonepoint) + { + double x1,y1,z1,x2,y2,z2; + bb.Get (x1,y1,z1,x2,y2,z2); + Point<3> p1 = Point<3> (x1,y1,z1); + Point<3> p2 = Point<3> (x2,y2,z2); + Box<3> boundingbox(p1,p2); + + center = boundingbox.Center(); + if (hiliteonepoint) + rad = occgeometry -> GetBoundingBox().Diam() / 100; + else + rad = boundingbox.Diam() / 2; + } + } + + CalcTransformationMatrices(); + } + + // Clear lists + + for (int i = 1; i <= linelists.Size(); i++) + glDeleteLists (linelists.Elem(i), 1); + linelists.SetSize(0); + + for (int i = 1; i <= trilists.Size(); i++) + glDeleteLists (trilists.Elem(i), 1); + trilists.SetSize(0); + + // Total wireframe + + linelists.Append (glGenLists (1)); + glNewList (linelists.Last(), GL_COMPILE); + + for (int i = 1; i <= occgeometry->emap.Extent(); i++) + { + TopoDS_Edge edge = TopoDS::Edge(occgeometry->emap(i)); + if (BRep_Tool::Degenerated(edge)) continue; + if (occgeometry->evispar[i-1].IsHighlighted()) continue; + + Handle(Poly_PolygonOnTriangulation) aEdgePoly; + Handle(Poly_Triangulation) T; + TopLoc_Location aEdgeLoc; + BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); + + if(aEdgePoly.IsNull()) + { + (*testout) << "visualizing edge " << occgeometry->emap.FindIndex (edge) + << " without using the occ visualization triangulation" << endl; + + double s0, s1; + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + + glBegin (GL_LINE_STRIP); + for (int i = 0; i<=50; i++) + { + gp_Pnt p = c->Value (s0 + i*(s1-s0)/50.0); + glVertex3f (p.X(),p.Y(),p.Z()); + } + glEnd (); + + continue; + } + + int nbnodes = aEdgePoly -> NbNodes(); + glBegin (GL_LINE_STRIP); + for (int j = 1; j <= nbnodes; j++) + { + gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); + glVertex3f (p.X(), p.Y(), p.Z()); + } + glEnd (); + } + + glEndList (); + + // Highlighted edge list + + linelists.Append (glGenLists (1)); + glNewList (linelists.Last(), GL_COMPILE); + + for (int i = 1; i <= occgeometry->emap.Extent(); i++) + if (occgeometry->evispar[i-1].IsHighlighted()) + { + TopoDS_Edge edge = TopoDS::Edge(occgeometry->emap(i)); + if (BRep_Tool::Degenerated(edge)) continue; + + Handle(Poly_PolygonOnTriangulation) aEdgePoly; + Handle(Poly_Triangulation) T; + TopLoc_Location aEdgeLoc; + BRep_Tool::PolygonOnTriangulation(edge, aEdgePoly, T, aEdgeLoc); + + if(aEdgePoly.IsNull()) + { + (*testout) << "visualizing edge " << occgeometry->emap.FindIndex (edge) + << " without using the occ visualization triangulation" << endl; + + double s0, s1; + Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); + + glBegin (GL_LINE_STRIP); + for (int i = 0; i<=50; i++) + { + gp_Pnt p = c->Value (s0 + i*(s1-s0)/50.0); + glVertex3f (p.X(),p.Y(),p.Z()); + } + glEnd (); + + continue; + } + + int nbnodes = aEdgePoly -> NbNodes(); + glBegin (GL_LINE_STRIP); + for (int j = 1; j <= nbnodes; j++) + { + gp_Pnt p = (T -> Nodes())(aEdgePoly->Nodes()(j)).Transformed(aEdgeLoc); + glVertex3f (p.X(), p.Y(), p.Z()); + } + glEnd (); + } + + glEndList (); + + // display faces + + trilists.Append (glGenLists (1)); + glNewList (trilists.Last(), GL_COMPILE); + + for (int i = 1; i <= occgeometry->fmap.Extent(); i++) + { + if (!occgeometry->fvispar[i-1].IsVisible()) continue; + + glLoadName (i); + float mat_col[4]; + mat_col[3] = 1; + + TopoDS_Face face = TopoDS::Face(occgeometry->fmap(i)); + + if (!occgeometry->fvispar[i-1].IsHighlighted()) + { + // Philippose - 30/01/2009 + // OpenCascade XDE Support + Quantity_Color face_colour; + occgeometry->face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour); + mat_col[0] = face_colour.Red(); + mat_col[1] = face_colour.Green(); + mat_col[2] = face_colour.Blue(); + } + else + { + mat_col[0] = 0.8; + mat_col[1] = 0.2; + mat_col[2] = 0.2; + } + + glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); + + TopLoc_Location loc; + Handle(Geom_Surface) surf = BRep_Tool::Surface (face); + BRepAdaptor_Surface sf(face, Standard_False); + BRepLProp_SLProps prop(sf, 1, 1e-5); + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); + + if (triangulation.IsNull()) + { + cout << "cannot visualize face " << i << endl; + occgeometry->fvispar[i-1].SetNotDrawable(); + continue; + } + + gp_Pnt2d uv; + gp_Pnt pnt; + gp_Vec n; + + glBegin (GL_TRIANGLES); + + int ntriangles = triangulation -> NbTriangles(); + for (int j = 1; j <= ntriangles; j++) + { + Poly_Triangle triangle = (triangulation -> Triangles())(j); + gp_Pnt p[3]; + for (int k = 1; k <= 3; k++) + p[k-1] = (triangulation -> Nodes())(triangle(k)).Transformed(loc); + + for (int k = 1; k <= 3; k++) + { + uv = (triangulation -> UVNodes())(triangle(k)); + prop.SetParameters (uv.X(), uv.Y()); + + // surf->D0 (uv.X(), uv.Y(), pnt); + + if (prop.IsNormalDefined()) + n = prop.Normal(); + else + { + (*testout) << "Visualization of face " << i + << ": Normal vector not defined" << endl; + // n = gp_Vec (0,0,0); + gp_Vec a(p[0],p[1]); + gp_Vec b(p[0],p[2]); + n = b^a; + } + + if (face.Orientation() == TopAbs_REVERSED) n *= -1; + glNormal3f (n.X(), n.Y(), n.Z()); + glVertex3f (p[k-1].X(), p[k-1].Y(), p[k-1].Z()); + } + } + glEnd (); + + } + glEndList (); + + } + + void SelectFaceInOCCDialogTree (int facenr); + + void VisualSceneOCCGeometry :: MouseDblClick (int px, int py) + { + int hits; + + // select surface triangle by mouse click + + GLuint selbuf[10000]; + glSelectBuffer (10000, selbuf); + + glRenderMode (GL_SELECT); + + GLint viewport[4]; + glGetIntegerv (GL_VIEWPORT, viewport); + + glMatrixMode (GL_PROJECTION); + glPushMatrix(); + + GLdouble projmat[16]; + glGetDoublev (GL_PROJECTION_MATRIX, projmat); + + glLoadIdentity(); + gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); + glMultMatrixd (projmat); + + glClearColor(backcolor, backcolor, backcolor, 1.0); + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + glMatrixMode (GL_MODELVIEW); + + glPushMatrix(); + glMultMatrixf (transformationmat); + + glInitNames(); + glPushName (1); + + glPolygonOffset (1, 1); + glEnable (GL_POLYGON_OFFSET_FILL); + + glDisable(GL_CLIP_PLANE0); + + // Philippose - 30/01/2009 + // Enable clipping planes for Selection mode in OCC Geometry + if (vispar.clipenable) + { + Vec<3> n(clipplane[0], clipplane[1], clipplane[2]); + double len = Abs(n); + double mu = -clipplane[3] / (len*len); + Point<3> p (mu * n); + n /= len; + Vec<3> t1 = n.GetNormal (); + Vec<3> t2 = Cross (n, t1); + + double xi1mid = (center - p) * t1; + double xi2mid = (center - p) * t2; + + glLoadName (0); + glBegin (GL_QUADS); + glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid-rad) * t2); + glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid-rad) * t2); + glVertex3dv (p + (xi1mid+rad) * t1 + (xi2mid+rad) * t2); + glVertex3dv (p + (xi1mid-rad) * t1 + (xi2mid+rad) * t2); + glEnd (); + } + + glCallList (trilists.Get(1)); + + glDisable (GL_POLYGON_OFFSET_FILL); + + glPopName(); + + glMatrixMode (GL_PROJECTION); + glPopMatrix(); + + glMatrixMode (GL_MODELVIEW); + glPopMatrix(); + + glFlush(); + + hits = glRenderMode (GL_RENDER); + + int minname = 0; + GLuint mindepth = 0; + + // find clippingplane + GLuint clipdepth = 0; // GLuint(-1); + + for (int i = 0; i < hits; i++) + { + int curname = selbuf[4*i+3]; + if (!curname) clipdepth = selbuf[4*i+1]; + } + + for (int i = 0; i < hits; i++) + { + int curname = selbuf[4*i+3]; + GLuint curdepth = selbuf[4*i+1]; + if (curname && (curdepth> clipdepth) && + (curdepth < mindepth || !minname)) + { + mindepth = curdepth; + minname = curname; + } + } + + occgeometry->LowLightAll(); + + if (minname) + { + occgeometry->fvispar[minname-1].Highlight(); + + if (vispar.occzoomtohighlightedentity) + occgeometry->changed = OCCGEOMETRYVISUALIZATIONFULLCHANGE; + else + occgeometry->changed = OCCGEOMETRYVISUALIZATIONHALFCHANGE; + cout << "Selected face: " << minname << endl; + } else - { - mat_col[0] = 0.8; - mat_col[1] = 0.2; - mat_col[2] = 0.2; - } + { + occgeometry->changed = OCCGEOMETRYVISUALIZATIONHALFCHANGE; + } - glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_col); + glDisable(GL_CLIP_PLANE0); - TopoDS_Face face = TopoDS::Face(occgeometry->fmap(i)); - TopLoc_Location loc; - Handle(Geom_Surface) surf = BRep_Tool::Surface (face); - BRepAdaptor_Surface sf(face, Standard_False); - BRepLProp_SLProps prop(sf, 1, 1e-5); - Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); - - if (triangulation.IsNull()) - { - cout << "cannot visualize face " << i << endl; - occgeometry->fvispar[i-1].SetNotDrawable(); - continue; - } - - gp_Pnt2d uv; - gp_Pnt pnt; - gp_Vec n; + SelectFaceInOCCDialogTree (minname); - glBegin (GL_TRIANGLES); - - int ntriangles = triangulation -> NbTriangles(); - for (int j = 1; j <= ntriangles; j++) - { - Poly_Triangle triangle = (triangulation -> Triangles())(j); - gp_Pnt p[3]; - for (int k = 1; k <= 3; k++) - p[k-1] = (triangulation -> Nodes())(triangle(k)).Transformed(loc); + // Philippose - 30/01/2009 + // Set the currently selected face in the array + // for local face mesh size definition + occgeometry->SetSelectedFace(minname); - for (int k = 1; k <= 3; k++) - { - uv = (triangulation -> UVNodes())(triangle(k)); - prop.SetParameters (uv.X(), uv.Y()); - - // surf->D0 (uv.X(), uv.Y(), pnt); - - if (prop.IsNormalDefined()) - n = prop.Normal(); - else - { - (*testout) << "Visualization of face " << i - << ": Normal vector not defined" << endl; - // n = gp_Vec (0,0,0); - gp_Vec a(p[0],p[1]); - gp_Vec b(p[0],p[2]); - n = b^a; - } - - if (face.Orientation() == TopAbs_REVERSED) n *= -1; - glNormal3f (n.X(), n.Y(), n.Z()); - glVertex3f (p[k-1].X(), p[k-1].Y(), p[k-1].Z()); - } - } - glEnd (); - - } - glEndList (); + // selecttimestamp = NextTimeStamp(); + } } -void SelectFaceInOCCDialogTree (int facenr); - -void VisualSceneOCCGeometry :: MouseDblClick (int px, int py) -{ - int hits; - - // select surface triangle by mouse click - - GLuint selbuf[10000]; - glSelectBuffer (10000, selbuf); - - - glRenderMode (GL_SELECT); - - GLint viewport[4]; - glGetIntegerv (GL_VIEWPORT, viewport); - - - glMatrixMode (GL_PROJECTION); - glPushMatrix(); - - GLdouble projmat[16]; - glGetDoublev (GL_PROJECTION_MATRIX, projmat); - - glLoadIdentity(); - gluPickMatrix (px, viewport[3] - py, 1, 1, viewport); - glMultMatrixd (projmat); - - - - glClearColor(backcolor, backcolor, backcolor, 1.0); - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - - glMatrixMode (GL_MODELVIEW); - - glPushMatrix(); - glMultMatrixf (transformationmat); - - - - glInitNames(); - glPushName (1); - - glPolygonOffset (1, 1); - glEnable (GL_POLYGON_OFFSET_FILL); - - glDisable(GL_CLIP_PLANE0); - - glCallList (trilists.Get(1)); - - glDisable (GL_POLYGON_OFFSET_FILL); - - glPopName(); - - glMatrixMode (GL_PROJECTION); - glPopMatrix(); - - glMatrixMode (GL_MODELVIEW); - glPopMatrix(); - - glFlush(); - - - hits = glRenderMode (GL_RENDER); - - int minname = 0; - GLuint mindepth = 0; - - // find clippingplane - GLuint clipdepth = 0; // GLuint(-1); - - for (int i = 0; i < hits; i++) - { - int curname = selbuf[4*i+3]; - if (!curname) clipdepth = selbuf[4*i+1]; - } - - for (int i = 0; i < hits; i++) - { - int curname = selbuf[4*i+3]; - GLuint curdepth = selbuf[4*i+1]; - if (curname && (curdepth > clipdepth) && - (curdepth < mindepth || !minname)) - { - mindepth = curdepth; - minname = curname; - } - } - - occgeometry->LowLightAll(); - - if (minname) - { - occgeometry->fvispar[minname-1].Highlight(); - - if (vispar.occzoomtohighlightedentity) - occgeometry->changed = OCCGEOMETRYVISUALIZATIONFULLCHANGE; - else - occgeometry->changed = OCCGEOMETRYVISUALIZATIONHALFCHANGE; - cout << "Selected face: " << minname << endl; - } - else - { - occgeometry->changed = OCCGEOMETRYVISUALIZATIONHALFCHANGE; - } - - glDisable(GL_CLIP_PLANE0); - - SelectFaceInOCCDialogTree (minname); - - - // selecttimestamp = NextTimeStamp(); -} - - - - - - -} - - - #endif #endif // NOTCL -