diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp index b0cfe87a..9c537df6 100644 --- a/libsrc/stlgeom/stlgeom.cpp +++ b/libsrc/stlgeom/stlgeom.cpp @@ -2104,22 +2104,21 @@ void STLGeometry :: TopologyChanged() int STLGeometry :: CheckGeometryOverlapping() { - int i, j, k; + PrintMessageCR(3,"Check overlapping geometry ..."); Box<3> geombox = GetBoundingBox(); Point<3> pmin = geombox.PMin(); Point<3> pmax = geombox.PMax(); Box3dTree setree(pmin, pmax); - Array inters; int oltrigs = 0; markedtrigs.SetSize(GetNT()); - for (i = 1; i <= GetNT(); i++) + for (int i = 1; i <= GetNT(); i++) SetMarkedTrig(i, 0); - for (i = 1; i <= GetNT(); i++) + for (int i = 1; i <= GetNT(); i++) { const STLTriangle & tri = GetTriangle(i); @@ -2133,48 +2132,57 @@ int STLGeometry :: CheckGeometryOverlapping() setree.Insert (tpmin, tpmax, i); } - for (i = 1; i <= GetNT(); i++) - { - const STLTriangle & tri = GetTriangle(i); - - Point<3> tpmin = tri.box.PMin(); - Point<3> tpmax = tri.box.PMax(); - setree.GetIntersecting (tpmin, tpmax, inters); - - for (j = 1; j <= inters.Size(); j++) - { - const STLTriangle & tri2 = GetTriangle(inters.Get(j)); - - const Point<3> *trip1[3], *trip2[3]; - Point<3> hptri1[3], hptri2[3]; - /* - for (k = 1; k <= 3; k++) - { +#pragma omp parallel + { + Array inters; + +#pragma omp for + for (int i = 1; i <= GetNT(); i++) + { + const STLTriangle & tri = GetTriangle(i); + + Point<3> tpmin = tri.box.PMin(); + Point<3> tpmax = tri.box.PMax(); + + setree.GetIntersecting (tpmin, tpmax, inters); + + for (int j = 1; j <= inters.Size(); j++) + { + const STLTriangle & tri2 = GetTriangle(inters.Get(j)); + + const Point<3> *trip1[3], *trip2[3]; + Point<3> hptri1[3], hptri2[3]; + /* + for (k = 1; k <= 3; k++) + { trip1[k-1] = &GetPoint (tri.PNum(k)); trip2[k-1] = &GetPoint (tri2.PNum(k)); - } - */ + } + */ + + for (int k = 0; k < 3; k++) + { + hptri1[k] = GetPoint (tri[k]); + hptri2[k] = GetPoint (tri2[k]); + trip1[k] = &hptri1[k]; + trip2[k] = &hptri2[k]; + } - for (k = 0; k < 3; k++) - { - hptri1[k] = GetPoint (tri[k]); - hptri2[k] = GetPoint (tri2[k]); - trip1[k] = &hptri1[k]; - trip2[k] = &hptri2[k]; - } - - if (IntersectTriangleTriangle (&trip1[0], &trip2[0])) - { - oltrigs++; - PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!"); - SetMarkedTrig(i, 1); - SetMarkedTrig(inters.Get(j), 1); - } - } - } - - PrintMessage(3,"Check Geometry Overlapping: overlapping triangles = ",oltrigs); + if (IntersectTriangleTriangle (&trip1[0], &trip2[0])) + { +#pragma omp critical + { + oltrigs++; + PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!"); + SetMarkedTrig(i, 1); + SetMarkedTrig(inters.Get(j), 1); + } + } + } + } + } + PrintMessage(3,"Check overlapping geometry ... ", oltrigs, " triangles overlap"); return oltrigs; } @@ -2408,12 +2416,11 @@ void STLGeometry :: CalcEdgeData() PushStatus("Calc Edge Data"); int np1, np2; - int i; int ecnt = 0; edgedata->SetSize(GetNT()/2*3); - for (i = 1; i <= GetNT(); i++) + for (int i = 1; i <= GetNT(); i++) { SetThreadPercent((double)i/(double)GetNT()*100.); @@ -2452,11 +2459,9 @@ void STLGeometry :: CalcEdgeData() void STLGeometry :: CalcEdgeDataAngles() { - PrintMessage(5,"calc edge data angles"); + PrintMessageCR (5,"calc edge data angles ... "); - int i; - - for (i = 1; i <= GetNTE(); i++) + for (int i = 1; i <= GetNTE(); i++) { STLTopEdge & edge = GetTopEdge (i); double cosang = @@ -2465,7 +2470,7 @@ void STLGeometry :: CalcEdgeDataAngles() edge.SetCosAngle (cosang); } - for (i = 1; i <= edgedata->Size(); i++) + for (int i = 1; i <= edgedata->Size(); i++) { /* const STLEdgeData& e = edgedata->Get(i); @@ -2474,7 +2479,7 @@ void STLGeometry :: CalcEdgeDataAngles() edgedata->Elem(i).angle = fabs(ang); */ } - + PrintMessage (5,"calc edge data angles ... done"); } void STLGeometry :: FindEdgesFromAngles() @@ -3044,11 +3049,10 @@ void STLGeometry :: BuildSmoothEdges () PushStatusF("Build Smooth Edges"); - int i, j;//, k, l; int nt = GetNT(); Vec3d ng1, ng2; - for (i = 1; i <= nt; i++) + for (int i = 1; i <= nt; i++) { if (multithread.terminate) {PopStatus();return;} @@ -3060,7 +3064,7 @@ void STLGeometry :: BuildSmoothEdges () ng1 = trig.GeomNormal(points); ng1 /= (ng1.Length() + 1e-24); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { int nbt = NeighbourTrig (i, j); @@ -3084,7 +3088,6 @@ void STLGeometry :: BuildSmoothEdges () } } } - PopStatus(); } diff --git a/libsrc/stlgeom/stlgeomchart.cpp b/libsrc/stlgeom/stlgeomchart.cpp index 105ba8c6..867c0cc2 100644 --- a/libsrc/stlgeom/stlgeomchart.cpp +++ b/libsrc/stlgeom/stlgeomchart.cpp @@ -19,30 +19,27 @@ int chartdebug = 0; void STLGeometry :: MakeAtlas(Mesh & mesh) { - - double h, h2; - - h = mparam.maxh; - + int timer1 = NgProfiler::CreateTimer ("makeatlas"); + int timer2 = NgProfiler::CreateTimer ("makeatlas - part 2"); + int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3"); + int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4"); PushStatusF("Make Atlas"); - int i,j,k,l; - + double h = mparam.maxh; + double atlasminh = 5e-3 * Dist (boundingbox.PMin(), boundingbox.PMax()); PrintMessage(5, "atlasminh = ", atlasminh); //speedup for make atlas if (GetNT() > 50000) - { - mesh.SetGlobalH(0.05*Dist (boundingbox.PMin(), boundingbox.PMax())); - } - + mesh.SetGlobalH(0.05*Dist (boundingbox.PMin(), boundingbox.PMax())); atlas.SetSize(0); ClearSpiralPoints(); BuildSmoothEdges(); + NgProfiler::StartTimer (timer1); double chartangle = stlparam.chartangle; double outerchartangle = stlparam.outerchartangle; @@ -63,20 +60,16 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) Array chartpoints; //point in chart becomes chartnum Array innerchartpoints; Array dirtycharttrigs; - Array chartpointchecked; - Array chartdistacttrigs; //outercharttrigs - chartdistacttrigs.SetSize(GetNT()); - for (i = 1; i <= GetNT(); i++) - { - chartdistacttrigs.Elem(i) = 0; - } - + Array chartdistacttrigs (GetNT()); //outercharttrigs + chartdistacttrigs = 0; + + STLBoundary chartbound(this); //knows the actual chart boundary //int chartboundarydivisions = 10; markedsegs.SetSize(0); //for testing!!! - chartpointchecked.SetSize(GetNP()); //for dirty-chart-trigs + Array chartpointchecked(GetNP()); //for dirty-chart-trigs outermark.SetSize(GetNT()); outertested.SetSize(GetNT()); @@ -84,7 +77,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) innerpointstochart.SetSize(GetNP()); chartmark.SetSize(GetNT()); - for (i = 1; i <= GetNP(); i++) + for (int i = 1; i <= GetNP(); i++) { innerpointstochart.Elem(i) = 0; pointstochart.Elem(i) = 0; @@ -96,70 +89,70 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) int spiralcheckon = stldoctor.spiralcheck; if (!spiralcheckon) {PrintWarning("++++++++++++\nspiral deactivated by user!!!!\n+++++++++++++++"); } - for (i = 1; i <= GetNT(); i++) - { - chartmark.Elem(i) = 0; - } + chartmark = 0; + outermark = 0; + outertested = 0; - for (i = 1; i <= GetNT(); i++) - { - outermark.Elem(i) = 0; - outertested.Elem(i) = 0; - } - - int markedtrigcnt = 0; - int found = 1; double atlasarea = Area(); double workedarea = 0; double showinc = 100.*5000./(double)GetNT(); double nextshow = 0; - Point<3> startp; + // Point<3> startp; int lastunmarked = 1; - int prelastunmarked; PrintMessage(5,"one dot per 5000 triangles: "); - while(markedtrigcnt < GetNT() && found) + int markedtrigcnt = 0; + while (markedtrigcnt < GetNT()) { - if (multithread.terminate) - {PopStatus();return;} + if (multithread.terminate) { PopStatus(); return; } if (workedarea / atlasarea*100. >= nextshow) {PrintDot(); nextshow+=showinc;} SetThreadPercent(100.0 * workedarea / atlasarea); - /* - for (j = 1; j <= GetNT(); j++) - { - outermark.Elem(j) = 0; - } - */ STLChart * chart = new STLChart(this); atlas.Append(chart); //find unmarked trig - prelastunmarked = lastunmarked; - j = lastunmarked; - found = 0; - while (!found && j <= GetNT()) + int prelastunmarked = lastunmarked; + /* + int j = lastunmarked; + bool found = 0; + while (j <= GetNT()) { - if (!GetMarker(j)) {found = 1; lastunmarked = j;} - else {j++;} + if (!GetMarker(j)) + { + found = 1; + lastunmarked = j; + break; + } + else + j++; } + */ + bool found = false; + for (int j = lastunmarked; j <= GetNT(); j++) + if (!GetMarker(j)) + { + found = true; + lastunmarked = j; + break; + } chartpoints.SetSize(0); innerchartpoints.SetSize(0); chartbound.Clear(); chartbound.SetChart(chart); - if (!found) {PrintSysError("Make Atlas, no starttrig found"); return;} + if (!found) { PrintSysError("Make Atlas, no starttrig found"); return; } //find surrounding trigs - int starttrig = j; + // int starttrig = j; + int starttrig = lastunmarked; - double tdist; - startp = GetPoint(GetTriangle(starttrig).PNum(1)); + Point<3> startp = GetPoint(GetTriangle(starttrig).PNum(1)); int accepted; int chartnum = GetNOCharts(); @@ -175,7 +168,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) workedarea += GetTriangle(starttrig).Area(points); - for (i = 1; i <= 3; i++) + for (int i = 1; i <= 3; i++) { innerpointstochart.Elem(GetTriangle(starttrig).PNum(i)) = chartnum; pointstochart.Elem(GetTriangle(starttrig).PNum(i)) = chartnum; @@ -183,32 +176,32 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) innerchartpoints.Append(GetTriangle(starttrig).PNum(i)); } - Vec<3> n2, n3; - int changed = 1; - int nt; - int ic; + bool changed = true; int oldstartic = 1; int oldstartic2; - int np1, np2; + + NgProfiler::StartTimer (timer2); + while (changed) { - changed = 0; + changed = false; oldstartic2 = oldstartic; oldstartic = chart->GetNT(); // for (ic = oldstartic2; ic <= chart->GetNT(); ic++) - for (ic = oldstartic2; ic <= oldstartic; ic++) + for (int ic = oldstartic2; ic <= oldstartic; ic++) { - i = chart->GetTrig(ic); + int i = chart->GetTrig(ic); if (GetMarker(i) == chartnum) { - for (j = 1; j <= NONeighbourTrigs(i); j++) + for (int j = 1; j <= NONeighbourTrigs(i); j++) { - nt = NeighbourTrig(i,j); + int nt = NeighbourTrig(i,j); + int np1, np2; GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2); if (GetMarker(nt) == 0 && !IsEdge(np1,np2)) { - n2 = GetTriangle(nt).Normal(); + Vec<3> n2 = GetTriangle(nt).Normal(); if ( (n2 * sn) >= coschartangle ) { @@ -240,14 +233,13 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) } */ - int nnp1, nnp2; - int nnt; - //find overlapping charts exacter: - for (k = 1; k <= 3; k++) + //find overlapping charts exacter (fast, too): + for (int k = 1; k <= 3; k++) { - nnt = NeighbourTrig(nt,k); + int nnt = NeighbourTrig(nt,k); if (GetMarker(nnt) != chartnum) { + int nnp1, nnp2; GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); accepted = chartbound.TestSeg(GetPoint(nnp1), @@ -255,7 +247,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) sn,sinchartangle,1 /*chartboundarydivisions*/ ,points, eps); - n3 = GetTriangle(nnt).Normal(); + Vec<3> n3 = GetTriangle(nnt).Normal(); if ( (n3 * sn) >= coschartangle && IsSmoothEdge (nnp1, nnp2) ) accepted = 1; @@ -276,14 +268,14 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) if (accepted) { SetMarker(nt, chartnum); - changed = 1; + changed = true; markedtrigcnt++; workedarea += GetTriangle(nt).Area(points); chart->AddChartTrig(nt); chartbound.AddTriangle(GetTriangle(nt)); - for (k = 1; k <= 3; k++) + for (int k = 1; k <= 3; k++) { if (innerpointstochart.Get(GetTriangle(nt).PNum(k)) != chartnum) @@ -302,6 +294,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) } } + NgProfiler::StopTimer (timer2); + NgProfiler::StartTimer (timer3); //find outertrigs @@ -311,155 +305,153 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) outermark.Elem(starttrig) = chartnum; //chart->AddOuterTrig(starttrig); - changed = 1; + changed = true; oldstartic = 1; while (changed) { - changed = 0; + changed = false; oldstartic2 = oldstartic; oldstartic = chart->GetNT(); - //for (ic = oldstartic2; ic <= chart->GetNT(); ic++) - for (ic = oldstartic2; ic <= oldstartic; ic++) + + for (int ic = oldstartic2; ic <= oldstartic; ic++) { - i = chart->GetTrig(ic); - - if (outermark.Get(i) == chartnum) + int i = chart->GetTrig(ic); + if (outermark.Get(i) != chartnum) continue; + + for (int j = 1; j <= NONeighbourTrigs(i); j++) { - for (j = 1; j <= NONeighbourTrigs(i); j++) + int nt = NeighbourTrig(i,j); + if (outermark.Get(nt) == chartnum) continue; + + const STLTriangle & ntrig = GetTriangle(nt); + int np1, np2; + GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2); + + if (IsEdge (np1, np2)) continue; + + + /* + if (outertested.Get(nt) == chartnum) + continue; + */ + outertested.Elem(nt) = chartnum; + + Vec<3> n2 = GetTriangle(nt).Normal(); + /* + double ang; + ang = Angle(n2,sn); + if (ang < -M_PI*0.5) {ang += 2*M_PI;} + + (*testout) << "ang < ocharang = " << (fabs(ang) <= outerchartangle); + (*testout) << " = " << ( (n2 * sn) >= cosouterchartangle) << endl; + + // if (fabs(ang) <= outerchartangle) + */ + //abfragen, ob noch im tolerierten Winkel + if ( (n2 * sn) >= cosouterchartangle ) { - nt = NeighbourTrig(i,j); - if (outermark.Get(nt) == chartnum) - continue; + accepted = 1; + + bool isdirtytrig = false; + Vec<3> gn = GetTriangle(nt).GeomNormal(points); + double gnlen = gn.Length(); + + if (n2 * gn <= cosouterchartanglehalf * gnlen) + isdirtytrig = true; + + //zurueckweisen, falls eine Spiralartige outerchart entsteht + + // NgProfiler::StartTimer (timer4); - const STLTriangle & ntrig = GetTriangle(nt); - GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2); - - if (IsEdge (np1, np2)) - continue; - - - /* - if (outertested.Get(nt) == chartnum) - continue; - */ - outertested.Elem(nt) = chartnum; - - - n2 = GetTriangle(nt).Normal(); - /* - double ang; - ang = Angle(n2,sn); - if (ang < -M_PI*0.5) {ang += 2*M_PI;} + //find overlapping charts exacter: + //do not check dirty trigs! + if (spiralcheckon && !isdirtytrig) + for (int k = 1; k <= 3; k++) + { + int nnt = NeighbourTrig(nt,k); - (*testout) << "ang < ocharang = " << (fabs(ang) <= outerchartangle); - (*testout) << " = " << ( (n2 * sn) >= cosouterchartangle) << endl; - - // if (fabs(ang) <= outerchartangle) - */ - //abfragen, ob noch im tolerierten Winkel - if ( (n2 * sn) >= cosouterchartangle ) - { - accepted = 1; - - int isdirtytrig = 0; - Vec<3> gn = GetTriangle(nt).GeomNormal(points); - double gnlen = gn.Length(); - - if (n2 * gn <= cosouterchartanglehalf * gnlen) - {isdirtytrig = 1;} - - //zurueckweisen, falls eine Spiralartige outerchart entsteht - int nnp1, nnp2; - int nnt; - //find overlapping charts exacter: - //do not check dirty trigs! - - - if (spiralcheckon && !isdirtytrig) - for (k = 1; k <= 3; k++) - { - nnt = NeighbourTrig(nt,k); + if (outermark.Elem(nnt) != chartnum) + { + int nnp1, nnp2; + GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); - if (outermark.Elem(nnt) != chartnum) - { - GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); - - accepted = - chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2), - sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps); - - - n3 = GetTriangle(nnt).Normal(); - if ( (n3 * sn) >= cosouterchartangle && - IsSmoothEdge (nnp1, nnp2) ) - accepted = 1; - } - if (!accepted) {break;} + accepted = + chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2), + sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps); + + + Vec<3> n3 = GetTriangle(nnt).Normal(); + if ( (n3 * sn) >= cosouterchartangle && + IsSmoothEdge (nnp1, nnp2) ) + accepted = 1; + } + if (!accepted) break; + } + + // NgProfiler::StopTimer (timer4); + + // outer chart is only small environment of + // inner chart: + + if (accepted) + { + accepted = 0; + + for (int k = 1; k <= 3; k++) + if (innerpointstochart.Get(ntrig.PNum(k)) == chartnum) + { + accepted = 1; + break; } - //} - - - // outer chart is only small environment of - // inner chart: - if (accepted) - { - accepted = 0; - - for (k = 1; k <= 3; k++) - { - if (innerpointstochart.Get(ntrig.PNum(k)) == chartnum) - { - accepted = 1; - break; - } - } - - if (!accepted) - for (k = 1; k <= 3; k++) + if (!accepted) + for (int k = 1; k <= 3; k++) + { + Point<3> pt = GetPoint(ntrig.PNum(k)); + double h2 = sqr(mesh.GetH(pt)); + + for (int l = 1; l <= innerchartpoints.Size(); l++) { - Point<3> pt = GetPoint(ntrig.PNum(k)); - h2 = sqr(mesh.GetH(pt)); - - for (l = 1; l <= innerchartpoints.Size(); l++) + double tdist = + Dist2(pt, GetPoint (innerchartpoints.Get(l))); + if (tdist < 4 * h2) { - tdist = Dist2(pt, GetPoint (innerchartpoints.Get(l))); - if (tdist < 4 * h2) - { - accepted = 1; - break; - } + accepted = 1; + break; } - if (accepted) {break;} } - } - - - if (accepted) + if (accepted) break; + } + } + + + if (accepted) + { + changed = true; + outermark.Elem(nt) = chartnum; + + if (GetMarker(nt) != chartnum) { - changed = 1; - outermark.Elem(nt) = chartnum; - - if (GetMarker(nt) != chartnum) + chartbound.AddTriangle(GetTriangle(nt)); + chart->AddOuterTrig(nt); + for (int k = 1; k <= 3; k++) { - chartbound.AddTriangle(GetTriangle(nt)); - chart->AddOuterTrig(nt); - for (k = 1; k <= 3; k++) + if (pointstochart.Get(GetTriangle(nt).PNum(k)) + != chartnum) { - if (pointstochart.Get(GetTriangle(nt).PNum(k)) - != chartnum) - { - pointstochart.Elem(GetTriangle(nt).PNum(k)) = chartnum; - chartpoints.Append(GetTriangle(nt).PNum(k)); - } + pointstochart.Elem(GetTriangle(nt).PNum(k)) = chartnum; + chartpoints.Append(GetTriangle(nt).PNum(k)); } } } - } - } + } + } } - } - } + } + } + + NgProfiler::StopTimer (timer3); + //end of while loop for outer chart GetDirtyChartTrigs(chartnum, *chart, outermark, chartpointchecked, dirtycharttrigs); //dirtycharttrigs are local (chart) point numbers!!!!!!!!!!!!!!!! @@ -472,7 +464,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) //if all trigs would be eliminated -> leave 1 trig! dirtycharttrigs.SetSize(dirtycharttrigs.Size() - 1); } - for (k = 1; k <= dirtycharttrigs.Size(); k++) + for (int k = 1; k <= dirtycharttrigs.Size(); k++) { int tn = chart->GetChartTrig(dirtycharttrigs.Get(k)); outermark.Elem(tn) = 0; //not necessary, for later use @@ -489,21 +481,21 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) RestrictHChartDistOneChart(chartnum, chartdistacttrigs, mesh, h, 0.5, atlasminh); } + NgProfiler::StopTimer (timer1); + NgProfiler::Print(stdout); + + PrintMessage(5,""); PrintMessage(5,"NO charts=", atlas.Size()); int cnttrias = 0; - //int found2; outerchartspertrig.SetSize(GetNT()); - - for (i = 1; i <= atlas.Size(); i++) + for (int i = 1; i <= atlas.Size(); i++) { - //found2 = 1; - for (j = 1; j <= GetChart(i).GetNT(); j++) + for (int j = 1; j <= GetChart(i).GetNT(); j++) { int tn = GetChart(i).GetTrig(j); AddOCPT(tn,i); - } cnttrias += GetChart(i).GetNT(); @@ -511,28 +503,22 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) PrintMessage(5, "NO outer chart trias=", cnttrias); //sort outerchartspertrig - for (i = 1; i <= GetNT(); i++) + for (int i = 1; i <= GetNT(); i++) { - int swap; - for (k = 1; k < GetNOCPT(i); k++) - { - - for (j = 1; j < GetNOCPT(i); j++) - { - swap = GetOCPT(i,j); - if (GetOCPT(i,j+1) < swap) - { - SetOCPT(i,j,GetOCPT(i,j+1)); - SetOCPT(i,j+1,swap); - } - } - } + for (int k = 1; k < GetNOCPT(i); k++) + for (int j = 1; j < GetNOCPT(i); j++) + { + int swap = GetOCPT(i,j); + if (GetOCPT(i,j+1) < swap) + { + SetOCPT(i,j,GetOCPT(i,j+1)); + SetOCPT(i,j+1,swap); + } + } // check make atlas if (GetChartNr(i) <= 0 || GetChartNr(i) > GetNOCharts()) - { - PrintSysError("Make Atlas: chartnr(", i, ")=0!!"); - }; + PrintSysError("Make Atlas: chartnr(", i, ")=0!!"); } mesh.SetGlobalH(mparam.maxh); diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index ee06157b..80ed764b 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -15,11 +15,9 @@ namespace netgen //add a point into a pointlist, return pointnumber int AddPointIfNotExists(Array& ap, const Point3d& p, double eps) { - int i; - for (i = 1; i <= ap.Size(); i++) - { - if (Dist(ap.Get(i),p) <= eps ) {return i;} - } + for (int i = 1; i <= ap.Size(); i++) + if (Dist(ap.Get(i),p) <= eps ) + return i; return ap.Append(p); } @@ -74,11 +72,8 @@ void FIOReadInt(istream& ios, int& i) const int ilen = sizeof(int); char buf[ilen]; - int j; - for (j = 0; j < ilen; j++) - { - ios.get(buf[j]); - } + for (int j = 0; j < ilen; j++) + ios.get(buf[j]); memcpy(&i, &buf, ilen); } @@ -89,11 +84,8 @@ void FIOWriteInt(ostream& ios, const int& i) char buf[ilen]; memcpy(&buf, &i, ilen); - int j; - for (j = 0; j < ilen; j++) - { - ios << buf[j]; - } + for (int j = 0; j < ilen; j++) + ios << buf[j]; } void FIOReadDouble(istream& ios, double& i) @@ -101,11 +93,9 @@ void FIOReadDouble(istream& ios, double& i) const int ilen = sizeof(double); char buf[ilen]; - int j; - for (j = 0; j < ilen; j++) - { - ios.get(buf[j]); - } + for (int j = 0; j < ilen; j++) + ios.get(buf[j]); + memcpy(&i, &buf, ilen); } @@ -116,11 +106,8 @@ void FIOWriteDouble(ostream& ios, const double& i) char buf[ilen]; memcpy(&buf, &i, ilen); - int j; - for (j = 0; j < ilen; j++) - { - ios << buf[j]; - } + for (int j = 0; j < ilen; j++) + ios << buf[j]; } void FIOReadFloat(istream& ios, float& i) @@ -143,40 +130,28 @@ void FIOWriteFloat(ostream& ios, const float& i) char buf[ilen]; memcpy(&buf, &i, ilen); - int j; - for (j = 0; j < ilen; j++) - { - ios << buf[j]; - } + for (int j = 0; j < ilen; j++) + ios << buf[j]; } void FIOReadString(istream& ios, char* str, int len) { - int j; - for (j = 0; j < len; j++) - { - ios.get(str[j]); - } + for (int j = 0; j < len; j++) + ios.get(str[j]); } //read string and add terminating 0 void FIOReadStringE(istream& ios, char* str, int len) { - int j; - for (j = 0; j < len; j++) - { - ios.get(str[j]); - } + for (int j = 0; j < len; j++) + ios.get(str[j]); str[len] = 0; } void FIOWriteString(ostream& ios, char* str, int len) { - int j; - for (j = 0; j < len; j++) - { - ios << str[j]; - } + for (int j = 0; j < len; j++) + ios << str[j]; } @@ -317,62 +292,58 @@ STLTriangle :: STLTriangle(const int * apts) int STLTriangle :: IsNeighbourFrom(const STLTriangle& t) const { //triangles must have same orientation!!! - int i, j; - for(i = 0; i <= 2; i++) - { - for(j = 0; j <= 2; j++) - { - if (t.pts[(i+1)%3] == pts[j] && - t.pts[i] == pts[(j+1)%3]) - {return 1;} - } - } + + for(int i = 0; i <= 2; i++) + for(int j = 0; j <= 2; j++) + if (t.pts[(i+1)%3] == pts[j] && + t.pts[i] == pts[(j+1)%3]) + + return 1; + return 0; } int STLTriangle :: IsWrongNeighbourFrom(const STLTriangle& t) const { //triangles have not same orientation!!! - int i, j; - for(i = 0; i <= 2; i++) - { - for(j = 0; j <= 2; j++) - { - if (t.pts[(i+1)%3] == pts[(j+1)%3] && - t.pts[i] == pts[j]) - {return 1;} - } - } + for(int i = 0; i <= 2; i++) + for(int j = 0; j <= 2; j++) + if (t.pts[(i+1)%3] == pts[(j+1)%3] && + t.pts[i] == pts[j]) + + return 1; + return 0; } void STLTriangle :: GetNeighbourPoints(const STLTriangle& t, int& p1, int& p2) const { - int i, j; - for(i = 1; i <= 3; i++) - { - for(j = 1; j <= 3; j++) + for(int i = 1; i <= 3; i++) + for(int j = 1; j <= 3; j++) + if (t.PNumMod(i+1) == PNumMod(j) && + t.PNumMod(i) == PNumMod(j+1)) { - if (t.PNumMod(i+1) == PNumMod(j) && - t.PNumMod(i) == PNumMod(j+1)) - {p1 = PNumMod(j); p2 = PNumMod(j+1); return;} + p1 = PNumMod(j); + p2 = PNumMod(j+1); + return; } - } + PrintSysError("Get neighbourpoints failed!"); } int STLTriangle :: GetNeighbourPointsAndOpposite(const STLTriangle& t, int& p1, int& p2, int& po) const { - int i, j; - for(i = 1; i <= 3; i++) - { - for(j = 1; j <= 3; j++) + for(int i = 1; i <= 3; i++) + for(int j = 1; j <= 3; j++) + if (t.PNumMod(i+1) == PNumMod(j) && + t.PNumMod(i) == PNumMod(j+1)) { - if (t.PNumMod(i+1) == PNumMod(j) && - t.PNumMod(i) == PNumMod(j+1)) - {p1 = PNumMod(j); p2 = PNumMod(j+1); po = PNumMod(j+2); return 1;} + p1 = PNumMod(j); + p2 = PNumMod(j+1); + po = PNumMod(j+2); + return 1; } - } + return 0; } @@ -690,15 +661,12 @@ void STLChart :: AddOuterTrig(int i) int STLChart :: IsInWholeChart(int nr) const { - int i; - for (i = 1; i <= charttrigs->Size(); i++) - { - if (charttrigs->Get(i) == nr) {return 1;} - } - for (i = 1; i <= outertrigs->Size(); i++) - { - if (outertrigs->Get(i) == nr) {return 1;} - } + for (int i = 1; i <= charttrigs->Size(); i++) + if (charttrigs->Get(i) == nr) return 1; + + for (int i = 1; i <= outertrigs->Size(); i++) + if (outertrigs->Get(i) == nr) return 1; + return 0; } @@ -712,7 +680,6 @@ void STLChart :: GetTrianglesInBox (const Point3d & pmin, searchtree -> GetIntersecting (pmin, pmax, trias); else { - int i; Box3d box1(pmin, pmax); box1.Increase (1e-4); Box3d box2; @@ -720,9 +687,8 @@ void STLChart :: GetTrianglesInBox (const Point3d & pmin, trias.SetSize(0); int nt = GetNT(); - for (i = 1; i <= nt; i++) + for (int i = 1; i <= nt; i++) { - int trignum = GetTrig(i); const STLTriangle & trig = geometry->GetTriangle(trignum); box2.SetPoint (geometry->GetPoint (trig.PNum(1))); @@ -730,9 +696,7 @@ void STLChart :: GetTrianglesInBox (const Point3d & pmin, box2.AddPoint (geometry->GetPoint (trig.PNum(3))); if (box1.Intersect (box2)) - { - trias.Append (trignum); - } + trias.Append (trignum); } } } @@ -741,8 +705,7 @@ void STLChart :: GetTrianglesInBox (const Point3d & pmin, void STLChart :: MoveToOuterChart(const Array& trigs) { if (!trigs.Size()) {return;} - int i; - for (i = 1; i <= trigs.Size(); i++) + for (int i = 1; i <= trigs.Size(); i++) { if (charttrigs->Get(trigs.Get(i)) != -1) {AddOuterTrig(charttrigs->Get(trigs.Get(i)));} @@ -756,14 +719,13 @@ void STLChart :: DelChartTrigs(const Array& trigs) { if (!trigs.Size()) {return;} - int i; - for (i = 1; i <= trigs.Size(); i++) + for (int i = 1; i <= trigs.Size(); i++) { charttrigs->Elem(trigs.Get(i)) = -1; } int cnt = 0; - for (i = 1; i <= charttrigs->Size(); i++) + for (int i = 1; i <= charttrigs->Size(); i++) { if (charttrigs->Elem(i) == -1) { @@ -774,7 +736,7 @@ void STLChart :: DelChartTrigs(const Array& trigs) charttrigs->Elem(i-cnt+1) = charttrigs->Get(i+1); } } - i = charttrigs->Size() - trigs.Size(); + int i = charttrigs->Size() - trigs.Size(); charttrigs->SetSize(i); if (!geomsearchtreeon && stlparam.usesearchtree == 1) @@ -784,7 +746,7 @@ void STLChart :: DelChartTrigs(const Array& trigs) searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1), geometry->GetBoundingBox().PMax() + Vec3d(1,1,1)); - for (i = 1; i <= charttrigs->Size(); i++) + for (int i = 1; i <= charttrigs->Size(); i++) { const STLTriangle & trig = geometry->GetTriangle(i); const Point3d & p1 = geometry->GetPoint (trig.PNum(1)); @@ -1124,7 +1086,6 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, const Vec3d& sn) { - int j; int nseg = NOSegments(); Point<2> p2d1 = chart->Project2d (p1); @@ -1133,57 +1094,35 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, Box<2> box2d; box2d.Set (p2d1); box2d.Add (p2d2); - /* - Point2d pmin(p2d1); - pmin.SetToMin (p2d2); - Point2d pmax(p2d1); - pmax.SetToMax (p2d2); - */ Line2d l1 (p2d1, p2d2); - double lam1, lam2; double eps = 1e-3; - - for (j = 1; j <= nseg; j++) + bool ok = true; + + for (int j = 1; j <= nseg; j++) { + if (!ok) continue; const STLBoundarySeg & seg = GetSegment(j); - if (!box2d.Intersect (seg.BoundingBox())) - continue; - /* - if (seg.P2DMin()(0) > pmax(0)) continue; - if (seg.P2DMin()(1) > pmax(1)) continue; - if (seg.P2DMax()(0) < pmin(0)) continue; - if (seg.P2DMax()(1) < pmin(1)) continue; - */ - + if (!box2d.Intersect (seg.BoundingBox())) continue; if (seg.IsSmoothEdge()) continue; const Point<2> & sp1 = seg.P2D1(); const Point<2> & sp2 = seg.P2D2(); - Line2d l2 (sp1, sp2); + double lam1, lam2; int err = CrossPointBarycentric (l1, l2, lam1, lam2); - /* - if (chartdebug) - { - - (*testout) << "lam1 = " << lam1 << ", lam2 = " << lam2 << endl; - (*testout) << "p2d = " << p2d1 << ", " << p2d2 << endl; - (*testout) << "sp2d = " << sp1 << ", " << sp2 << endl; - (*testout) << "i1,2 = " << seg.I1() << ", " << seg.I2() << endl; - - } - */ + if (!err && lam1 > eps && lam1 < 1-eps && lam2 > eps && lam2 < 1-eps) - return 0; + ok = false; } - return 1; + + return ok; } diff --git a/libsrc/stlgeom/stltopology.cpp b/libsrc/stlgeom/stltopology.cpp index 8547034f..441b35a2 100644 --- a/libsrc/stlgeom/stltopology.cpp +++ b/libsrc/stlgeom/stltopology.cpp @@ -56,20 +56,20 @@ STLGeometry * STLTopology :: LoadBinary (istream & ist) Point<3> pts[3]; Vec<3> normal; - int cntface, j; - //int vertex = 0; - float f; char spaces[nospaces+1]; - for (cntface = 0; cntface < nofacets; cntface++) + for (int cntface = 0; cntface < nofacets; cntface++) { - if (cntface % 10000 == 9999) { PrintDot(); } + if (cntface % 10000 == 0) + // { PrintDot(); } + PrintMessageCR (3, cntface, " triangles loaded\r"); + float f; FIOReadFloat(ist,f); normal(0) = f; FIOReadFloat(ist,f); normal(1) = f; FIOReadFloat(ist,f); normal(2) = f; - for (j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { FIOReadFloat(ist,f); pts[j](0) = f; FIOReadFloat(ist,f); pts[j](1) = f; @@ -79,7 +79,7 @@ STLGeometry * STLTopology :: LoadBinary (istream & ist) readtrigs.Append (STLReadTriangle (pts, normal)); FIOReadString(ist,spaces,nospaces); } - + PrintMessage (3, nofacets, " triangles loaded\r"); geom->InitSTLGeometry(readtrigs); @@ -337,7 +337,6 @@ void STLTopology :: Save (const char* filename) const STLGeometry * STLTopology ::Load (istream & ist) { - size_t i; STLGeometry * geom = new STLGeometry(); Array readtrigs; @@ -348,14 +347,14 @@ STLGeometry * STLTopology ::Load (istream & ist) int cntface = 0; int vertex = 0; - bool badnormals = 0; + bool badnormals = false; while (ist.good()) { ist >> buf; - size_t n = strlen (buf); - for (i = 0; i < n; i++) + int n = strlen (buf); + for (int i = 0; i < n; i++) buf[i] = tolower (buf[i]); if (strcmp (buf, "facet") == 0) @@ -391,14 +390,11 @@ STLGeometry * STLTopology ::Load (istream & ist) else { - Vec<3> hnormal; - hnormal = Cross (pts[1]-pts[0], pts[2]-pts[0]); + Vec<3> hnormal = Cross (pts[1]-pts[0], pts[2]-pts[0]); hnormal.Normalize(); if (normal * hnormal < 0.5) - { - badnormals = 1; - } + badnormals = true; } vertex = 0; @@ -407,11 +403,18 @@ STLGeometry * STLTopology ::Load (istream & ist) (Dist2 (pts[0], pts[2]) > 1e-16) && (Dist2 (pts[1], pts[2]) > 1e-16) ) - readtrigs.Append (STLReadTriangle (pts, normal)); + { + readtrigs.Append (STLReadTriangle (pts, normal)); + + if (readtrigs.Size() % 100000 == 0) + PrintMessageCR (3, readtrigs.Size(), " triangles loaded\r"); + } + } } } - + PrintMessage (3, readtrigs.Size(), " triangles loaded"); + if (badnormals) { PrintWarning("File has normal vectors which differ extremly from geometry->correct with stldoctor!!!"); @@ -435,8 +438,6 @@ STLGeometry * STLTopology ::Load (istream & ist) void STLTopology :: InitSTLGeometry(const Array & readtrigs) { - int i, k; - // const double geometry_tol_fact = 1E6; // distances lower than max_box_size/tol are ignored @@ -445,13 +446,12 @@ void STLTopology :: InitSTLGeometry(const Array & readtrigs) PrintMessage(3,"number of triangles = ", readtrigs.Size()); - if (!readtrigs.Size()) - return; + if (!readtrigs.Size()) return; boundingbox.Set (readtrigs[0][0]); - for (i = 0; i < readtrigs.Size(); i++) - for (k = 0; k < 3; k++) + for (int i = 0; i < readtrigs.Size(); i++) + for (int k = 0; k < 3; k++) boundingbox.Add (readtrigs[i][k]); PrintMessage(5,"boundingbox: ", Point3d(boundingbox.PMin()), " - ", @@ -462,21 +462,20 @@ void STLTopology :: InitSTLGeometry(const Array & readtrigs) pointtree = new Point3dTree (bb.PMin(), bb.PMax()); - - Array pintersect; pointtol = boundingbox.Diam() * stldoctor.geom_tol_fact; PrintMessage(5,"point tolerance = ", pointtol); + PrintMessageCR(5,"identify points ..."); - for(i = 0; i < readtrigs.Size(); i++) + for(int i = 0; i < readtrigs.Size(); i++) { const STLReadTriangle & t = readtrigs[i]; + STLTriangle st; - // Vec<3> n = t.Normal(); st.SetNormal (t.Normal()); - for (k = 0; k < 3; k++) + for (int k = 0; k < 3; k++) { Point<3> p = t[k]; @@ -511,7 +510,7 @@ void STLTopology :: InitSTLGeometry(const Array & readtrigs) } } - + PrintMessage(5,"identify points ... done"); FindNeighbourTrigs(); } @@ -540,23 +539,22 @@ void STLTopology :: FindNeighbourTrigs() PushStatusF("Find Neighbour Triangles"); - int i, j, k, l; + PrintMessage(5,"build topology ..."); // build up topology tables - //int np = GetNP(); int nt = GetNT(); INDEX_2_HASHTABLE * oldedges = ht_topedges; ht_topedges = new INDEX_2_HASHTABLE (GetNP()+1); topedges.SetSize(0); - for (i = 1; i <= nt; i++) + for (int i = 1; i <= nt; i++) { STLTriangle & trig = GetTriangle(i); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { int pi1 = trig.PNumMod (j+1); int pi2 = trig.PNumMod (j+2); @@ -577,7 +575,7 @@ void STLTopology :: FindNeighbourTrigs() trig.NBTrigNum(j) = othertn; trig.EdgeNum(j) = enr; - for (k = 1; k <= 3; k++) + for (int k = 1; k <= 3; k++) if (othertrig.EdgeNum(k) == enr) othertrig.NBTrigNum(k) = i; } @@ -596,11 +594,11 @@ void STLTopology :: FindNeighbourTrigs() topology_ok = 1; int ne = GetNTE(); - for (i = 1; i <= nt; i++) + for (int i = 1; i <= nt; i++) GetTriangle(i).flags.toperror = 0; - for (i = 1; i <= nt; i++) - for (j = 1; j <= 3; j++) + for (int i = 1; i <= nt; i++) + for (int j = 1; j <= 3; j++) { const STLTopEdge & edge = GetTopEdge (GetTriangle(i).EdgeNum(j)); if (edge.TrigNum(1) != i && edge.TrigNum(2) != i) @@ -610,7 +608,7 @@ void STLTopology :: FindNeighbourTrigs() } } - for (i = 1; i <= ne; i++) + for (int i = 1; i <= ne; i++) { const STLTopEdge & edge = GetTopEdge (i); if (!edge.TrigNum(2)) @@ -623,10 +621,10 @@ void STLTopology :: FindNeighbourTrigs() if (topology_ok) { orientation_ok = 1; - for (i = 1; i <= nt; i++) + for (int i = 1; i <= nt; i++) { const STLTriangle & t = GetTriangle (i); - for (j = 1; j <= 3; j++) + for (int j = 1; j <= 3; j++) { const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j)); if (!t.IsNeighbourFrom (nbt)) @@ -658,8 +656,8 @@ void STLTopology :: FindNeighbourTrigs() // generate point -> trig table trigsperpoint.SetSize(GetNP()); - for (i = 1; i <= GetNT(); i++) - for (j = 1; j <= 3; j++) + for (int i = 1; i <= GetNT(); i++) + for (int j = 1; j <= 3; j++) trigsperpoint.Add1(GetTriangle(i).PNum(j),i); @@ -674,8 +672,8 @@ void STLTopology :: FindNeighbourTrigs() } */ topedgesperpoint.SetSize (GetNP()); - for (i = 1; i <= ne; i++) - for (j = 1; j <= 2; j++) + for (int i = 1; i <= ne; i++) + for (int j = 1; j <= 2; j++) topedgesperpoint.Add1 (GetTopEdge (i).PNum(j), i); PrintMessage(5,"point -> trig table generated"); @@ -691,7 +689,7 @@ void STLTopology :: FindNeighbourTrigs() for (STLTrigIndex ti = 0; ti < GetNT(); ti++) { STLTriangle & trig = trias[ti]; - for (k = 0; k < 3; k++) + for (int k = 0; k < 3; k++) { STLPointIndex pi = trig[k] - STLBASE; STLPointIndex pi2 = trig[(k+1)%3] - STLBASE; @@ -711,7 +709,7 @@ void STLTopology :: FindNeighbourTrigs() double phimin = 10, phimax = -1; // out of (0, 2 pi) - for (j = 0; j < trigsperpoint[pi].Size(); j++) + for (int j = 0; j < trigsperpoint[pi].Size(); j++) { STLTrigIndex ti2 = trigsperpoint[pi][j] - STLBASE; const STLTriangle & trig2 = trias[ti2]; @@ -719,7 +717,7 @@ void STLTopology :: FindNeighbourTrigs() if (ti == ti2) continue; bool hasboth = 0; - for (l = 0; l < 3; l++) + for (int l = 0; l < 3; l++) if (trig2[l] - STLBASE == pi2) { hasboth = 1; @@ -728,7 +726,7 @@ void STLTopology :: FindNeighbourTrigs() if (!hasboth) continue; STLPointIndex pi4(0); - for (l = 0; l < 3; l++) + for (int l = 0; l < 3; l++) if (trig2[l] - STLBASE != pi && trig2[l] - STLBASE != pi2) pi4 = trig2[l] - STLBASE; @@ -758,8 +756,8 @@ void STLTopology :: FindNeighbourTrigs() { // for compatibility: neighbourtrigs.SetSize(GetNT()); - for (i = 1; i <= GetNT(); i++) - for (k = 1; k <= 3; k++) + for (int i = 1; i <= GetNT(); i++) + for (int k = 1; k <= 3; k++) AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k)); } else @@ -770,7 +768,7 @@ void STLTopology :: FindNeighbourTrigs() int tr, found; int wrongneighbourfound = 0; - for (i = 1; i <= GetNT(); i++) + for (int i = 1; i <= GetNT(); i++) { SetThreadPercent((double)i/(double)GetNT()*100.); if (multithread.terminate) @@ -779,9 +777,9 @@ void STLTopology :: FindNeighbourTrigs() return; } - for (k = 1; k <= 3; k++) + for (int k = 1; k <= 3; k++) { - for (j = 1; j <= trigsperpoint.EntrySize(GetTriangle(i).PNum(k)); j++) + for (int j = 1; j <= trigsperpoint.EntrySize(GetTriangle(i).PNum(k)); j++) { tr = trigsperpoint.Get(GetTriangle(i).PNum(k),j); if (i != tr && (GetTriangle(i).IsNeighbourFrom(GetTriangle(tr))