From 22222d4cbced467c2bd6d93a1efb3d4bcf92c8e1 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Mon, 28 Apr 2014 07:07:36 +0000 Subject: [PATCH] second order mesh merge --- libsrc/meshing/meshclass.cpp | 32 ++++++++++------ libsrc/meshing/meshtype.hpp | 14 ++++++- libsrc/meshing/secondorder.cpp | 67 ++++++++++++++++++++++++++++----- libsrc/stlgeom/stlgeomchart.cpp | 35 +++++++++++++---- libsrc/stlgeom/stltool.cpp | 20 ++++++---- 5 files changed, 133 insertions(+), 35 deletions(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 70e85eee..94d0e3bb 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -1470,15 +1470,28 @@ namespace netgen // int si = sel.GetIndex(); - for (int j = 0; j < sel.GetNP(); j++) - { - INDEX_2 i2; - i2.I1() = sel.PNumMod(j+1); - i2.I2() = sel.PNumMod(j+2); - i2.Sort(); - if (sel.GetNP() <= 4) + if (sel.GetNP() <= 4) + for (int j = 0; j < sel.GetNP(); j++) + { + INDEX_2 i2; + i2.I1() = sel.PNumMod(j+1); + i2.I2() = sel.PNumMod(j+2); + i2.Sort(); boundaryedges->Set (i2, 1); + } + else if (sel.GetType()==TRIG6) + { + for (int j = 0; j < 3; j++) + { + INDEX_2 i2; + i2.I1() = sel[j]; + i2.I2() = sel[(j+1)%3]; + i2.Sort(); + boundaryedges->Set (i2, 1); + } } + else + cerr << "illegal elemenet for buildboundaryedges" << endl; } @@ -1825,6 +1838,7 @@ namespace netgen for (ii = 0; ii < row.Size(); ii++) { hel = SurfaceElement(row[ii]); + if (hel.GetType() == TRIG6) hel.SetType(TRIG); int ind = hel.GetIndex(); if (GetFaceDescriptor(ind).DomainIn() && @@ -4347,11 +4361,9 @@ namespace netgen int i = 0; const int maxits = 30; - while(delta > 1e-16 && iCalcElementTransformation(lam,element-1,x,Jac); - rhs = p-x; Jac.Solve(rhs,deltalam); @@ -4367,7 +4379,6 @@ namespace netgen if(i==maxits) return false; - for(i=0; i<3; i++) lami[i] = lam(i); @@ -4568,7 +4579,6 @@ namespace netgen ii = locels.Get(i); else ii = i; - if(ii == ps_startelement) continue; if(indices != NULL && indices->Size() > 0) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index ceaed59d..2e1ee4fd 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -354,11 +354,22 @@ namespace netgen /// int GetNV() const { + if (typ == TRIG || typ == TRIG6) + return 3; + else + { +#ifdef DEBUG + if (typ != QUAD && typ != QUAD6 && typ != QUAD8) + PrintSysError ("element2d::GetNV not implemented for typ", typ) +#endif + return 4; + } + /* switch (typ) { case TRIG: case TRIG6: return 3; - + case QUAD: case QUAD8: case QUAD6: return 4; @@ -369,6 +380,7 @@ namespace netgen ; } return np; + */ } /// diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 7cd2205a..8f419b01 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -1,4 +1,4 @@ -#include + #include #include "meshing.hpp" @@ -13,13 +13,64 @@ namespace netgen void Refinement :: MakeSecondOrder (Mesh & mesh) { - int nseg, nse, ne; + /* + Berlin, 2014: if we have curved surface elements, keep them ! + */ mesh.ComputeNVertices(); - mesh.SetNP(mesh.GetNV()); - + // mesh.SetNP(mesh.GetNV()); + mesh.SetNP(mesh.GetNP()); // setup multilevel-table + INDEX_2_HASHTABLE between(mesh.GetNP() + 5); + for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) + { + const Element2d & el = mesh[sei]; + + static int betw_trig[3][3] = + { { 1, 2, 3 }, { 0, 2, 4 }, { 0, 1, 5 } }; + static int betw_quad6[2][3] = + { { 0, 1, 4 }, { 3, 2, 5 } }; + static int betw_quad8[4][3] = + { { 0, 1, 4 }, { 3, 2, 5 }, + { 0, 3, 6 }, { 1, 2, 7 } }; + + int onp = 0; + int (*betw)[3] = NULL; + switch (el.GetType()) + { + case TRIG6: + { + betw = betw_trig; + onp = 3; + break; + } + case QUAD6: + { + betw = betw_quad6; + onp = 4; + break; + } + case QUAD8: + { + betw = betw_quad8; + onp = 4; + break; + } + default: + ; + } + + if (betw) + for (int j = 0; j < el.GetNP()-onp; j++) + { + int pi1 = el[betw[j][0]]; + int pi2 = el[betw[j][1]]; + INDEX_2 i2 = INDEX_2::Sort (pi1, pi2); + between.Set (i2, el[onp+j]); + } + } + bool thinlayers = 0; for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++) @@ -28,7 +79,7 @@ namespace netgen thinlayers = 1; - nseg = mesh.GetNSeg(); + int nseg = mesh.GetNSeg(); for (SegmentIndex si = 0; si < nseg; si++) { Segment & el = mesh.LineSegment(si); @@ -54,8 +105,7 @@ namespace netgen } // refine surface elements - nse = mesh.GetNSE(); - for (SurfaceElementIndex sei = 0; sei < nse; sei++) + for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) { int j; const Element2d & el = mesh.SurfaceElement(sei); @@ -149,8 +199,7 @@ namespace netgen // refine volume elements - ne = mesh.GetNE(); - for (int i = 1; i <= ne; i++) + for (int i = 1; i <= mesh.GetNE(); i++) { const Element & el = mesh.VolumeElement(i); int onp = 0; diff --git a/libsrc/stlgeom/stlgeomchart.cpp b/libsrc/stlgeom/stlgeomchart.cpp index 55f800c7..ca061ba4 100644 --- a/libsrc/stlgeom/stlgeomchart.cpp +++ b/libsrc/stlgeom/stlgeomchart.cpp @@ -22,6 +22,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) 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"); + int timer4a = NgProfiler::CreateTimer ("makeatlas - part 4a"); + int timer5 = NgProfiler::CreateTimer ("makeatlas - part 5"); PushStatusF("Make Atlas"); @@ -113,6 +116,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) STLChart * chart = new STLChart(this); atlas.Append(chart); + // *testout << "Chart " << atlas.Size() << endl; + //find unmarked trig int prelastunmarked = lastunmarked; @@ -143,7 +148,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) Vec<3> sn = GetTriangle(starttrig).Normal(); chart->SetNormal (startp, sn); - + + // *testout << "first trig " << starttrig << ", n = " << sn << endl; SetMarker(starttrig, chartnum); markedtrigcnt++; @@ -181,14 +187,16 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) for (int j = 1; j <= NONeighbourTrigs(i); j++) { int nt = NeighbourTrig(i,j); + // *testout << "check trig " << nt << endl; int np1, np2; GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2); if (GetMarker(nt) == 0 && !IsEdge(np1,np2)) { Vec<3> n2 = GetTriangle(nt).Normal(); + // *testout << "acos = " << 180/M_PI*acos (n2*sn) << endl; if ( (n2 * sn) >= coschartangle ) { - + // *testout << "good angle " << endl; accepted = 1; /* //alter spiralentest, schnell, aber ungenau @@ -229,19 +237,22 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) accepted = chartbound.TestSeg(GetPoint(nnp1), GetPoint(nnp2), sn,sinchartangle,1 /*chartboundarydivisions*/ ,points, eps); - + + // if (!accepted) *testout << "not acc due to testseg" << endl; Vec<3> n3 = GetTriangle(nnt).Normal(); if ( (n3 * sn) >= coschartangle && IsSmoothEdge (nnp1, nnp2) ) accepted = 1; } - if (!accepted) {break;} + if (!accepted) + break; } - + if (accepted) { + // *testout << "trig accepted" << endl; SetMarker(nt, chartnum); changed = true; markedtrigcnt++; @@ -318,6 +329,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) { accepted = 1; + NgProfiler::StartTimer (timer4); + bool isdirtytrig = false; Vec<3> gn = GetTriangle(nt).GeomNormal(points); double gnlen = gn.Length(); @@ -339,11 +352,15 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) int nnp1, nnp2; GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); + // NgProfiler::StartTimer (timer4a); + accepted = chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2), sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps); - + // NgProfiler::StopTimer (timer4a); + + Vec<3> n3 = GetTriangle(nnt).Normal(); if ( (n3 * sn) >= cosouterchartangle && IsSmoothEdge (nnp1, nnp2) ) @@ -352,7 +369,11 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) if (!accepted) break; } - + NgProfiler::StopTimer (timer4); + + NgProfiler::RegionTimer reg5(timer5); + + // outer chart is only small environment of // inner chart: diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index 20cb2c61..96cc384a 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -894,10 +894,10 @@ void STLBoundary ::AddTriangle(const STLTriangle & t) int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> & sn, double sinchartangle, int divisions, Array >& points, double eps) { - if (usechartnormal) return TestSegChartNV (p1, p2, sn); +#ifdef NONE // for statistics { int i; @@ -936,7 +936,7 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> */ } } - +#endif int i,j,k; Point<3> seg1p/*, seg2p*/; @@ -1087,6 +1087,9 @@ 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 timer = NgProfiler::CreateTimer ("TestSegChartNV"); + NgProfiler::StartTimer (timer); + int nseg = NOSegments(); Point<2> p2d1 = chart->Project2d (p1); @@ -1103,11 +1106,10 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, for (int j = 1; j <= nseg; j++) { - if (!ok) continue; const STLBoundarySeg & seg = GetSegment(j); - if (!box2d.Intersect (seg.BoundingBox())) continue; if (seg.IsSmoothEdge()) continue; + if (!box2d.Intersect (seg.BoundingBox())) continue; const Point<2> & sp1 = seg.P2D1(); const Point<2> & sp2 = seg.P2D2(); @@ -1115,14 +1117,18 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, Line2d l2 (sp1, sp2); double lam1, lam2; - int err = - CrossPointBarycentric (l1, l2, lam1, lam2); + int err = CrossPointBarycentric (l1, l2, lam1, lam2); if (!err && lam1 > eps && lam1 < 1-eps && lam2 > eps && lam2 < 1-eps) - ok = false; + { + ok = false; + break; + } } + NgProfiler::StopTimer (timer); + return ok; }