ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-11 21:50:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

stl meshing

Browse Source
This commit is contained in:
Joachim Schoeberl 2012-10-19 08:42:19 +00:00
parent 0166fca8d4
commit 052ef5db85
4 changed files with 394 additions and 468 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
libsrc/stlgeom/stlgeom.cpp
View File

@ -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<int> 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<int> 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();
}

434
libsrc/stlgeom/stlgeomchart.cpp
View File

@ -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<int> chartpoints; //point in chart becomes chartnum
Array<int> innerchartpoints;
Array<int> dirtycharttrigs;
Array<int> chartpointchecked;
Array<int> chartdistacttrigs; //outercharttrigs
chartdistacttrigs.SetSize(GetNT());
for (i = 1; i <= GetNT(); i++)
{
chartdistacttrigs.Elem(i) = 0;
}
Array<int> 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<int> 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);

211
libsrc/stlgeom/stltool.cpp
View File

@ -15,11 +15,9 @@ namespace netgen
//add a point into a pointlist, return pointnumber
int AddPointIfNotExists(Array<Point3d>& 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<int>& 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<int>& 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<int>& 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<int>& 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;
}

106
libsrc/stlgeom/stltopology.cpp
View File

@ -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<STLReadTriangle> 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<STLReadTriangle> & 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<STLReadTriangle> & 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<STLReadTriangle> & readtrigs)
pointtree = new Point3dTree (bb.PMin(), bb.PMax());
Array<int> 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<STLReadTriangle> & 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<int> * oldedges = ht_topedges;
ht_topedges = new INDEX_2_HASHTABLE<int> (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))