From 83e8b1ec53e934675fef236e677574c7fa1c58db Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Mon, 16 Nov 2009 08:18:00 +0000 Subject: [PATCH] mesh optimization improvements --- libsrc/csg/algprim.cpp | 8 +- libsrc/csg/algprim.hpp | 4 +- libsrc/csg/genmesh.cpp | 21 +- libsrc/csg/solid.cpp | 7 +- libsrc/csg/specpoin.cpp | 297 +++++++- libsrc/csg/specpoin.hpp | 10 + libsrc/meshing/adfront2.cpp | 27 +- libsrc/meshing/improve2.cpp | 1317 +++++++++++++++++---------------- libsrc/meshing/meshclass.cpp | 1 - libsrc/meshing/meshing2.cpp | 74 +- libsrc/meshing/meshtype.hpp | 3 +- libsrc/meshing/smoothing2.cpp | 23 +- 12 files changed, 1093 insertions(+), 699 deletions(-) diff --git a/libsrc/csg/algprim.cpp b/libsrc/csg/algprim.cpp index 1dde9179..a5b07441 100644 --- a/libsrc/csg/algprim.cpp +++ b/libsrc/csg/algprim.cpp @@ -351,7 +351,8 @@ namespace netgen { c = ac; r = ar; - + invr = 1.0/r; + cxx = cyy = czz = 0.5 / r; cxy = cxz = cyz = 0; cx = - c(0) / r; @@ -377,6 +378,7 @@ namespace netgen c(2) = coeffs.Elem(3); r = coeffs.Elem(4); + invr = 1.0/r; cxx = cyy = czz = 0.5 / r; cxy = cxz = cyz = 0; cx = - c(0) / r; @@ -412,6 +414,10 @@ namespace netgen } + double Sphere :: CalcFunctionValue (const Point<3> & point) const + { + return 0.5* (invr * Abs2 (point-c) - r); + } int Sphere :: IsIdentic (const Surface & s2, int & inv, double eps) const diff --git a/libsrc/csg/algprim.hpp b/libsrc/csg/algprim.hpp index e178a881..c16d52ba 100644 --- a/libsrc/csg/algprim.hpp +++ b/libsrc/csg/algprim.hpp @@ -122,7 +122,7 @@ namespace netgen /// Point<3> c; /// - double r; + double r, invr; public: /// Sphere (const Point<3> & ac, double ar); @@ -135,6 +135,8 @@ namespace netgen virtual Primitive * Copy () const; virtual void Transform (Transformation<3> & trans); + virtual double CalcFunctionValue (const Point<3> & point) const; + virtual int IsIdentic (const Surface & s2, int & inv, double eps) const; diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index a04bd05e..cec90311 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -221,16 +221,6 @@ namespace netgen static void MeshSurface (CSGeometry & geom, Mesh & mesh) { - /* - Point3d pmin, pmax; - mesh.GetBox(pmin, pmax); - cout << "box = " << pmin << " - " << pmax << endl; - cout << "localhbox = " - << mesh.LocalHFunction().GetBoundingBox().PMin() << " - " - << mesh.LocalHFunction().GetBoundingBox().PMax() - << endl; - */ - const char * savetask = multithread.task; multithread.task = "Surface meshing"; @@ -429,7 +419,7 @@ namespace netgen double eps = 1e-8 * geom.MaxSize(); for (PointIndex pi = PointIndex::BASE; pi < noldp+PointIndex::BASE; pi++) { - //if(surf->PointOnSurface(mesh[pi])) + // if(surf->PointOnSurface(mesh[pi])) meshing.AddPoint (mesh[pi], pi, NULL, (surf->PointOnSurface(mesh[pi], eps) != 0)); } @@ -511,18 +501,23 @@ namespace netgen if (segments.Size()) { // surface was meshed, not copied + + static int timer = NgProfiler::CreateTimer ("total surface mesh optimization"); + NgProfiler::RegionTimer reg (timer); + + PrintMessage (2, "Optimize Surface"); for (int i = 1; i <= mparam.optsteps2d; i++) { if (multithread.terminate) return; - + { MeshOptimize2dSurfaces meshopt(geom); meshopt.SetFaceIndex (k); meshopt.SetImproveEdges (0); meshopt.SetMetricWeight (mparam.elsizeweight); meshopt.SetWriteStatus (0); - + meshopt.EdgeSwapping (mesh, (i > mparam.optsteps2d/2)); } diff --git a/libsrc/csg/solid.cpp b/libsrc/csg/solid.cpp index 6eea88f2..70f7049f 100644 --- a/libsrc/csg/solid.cpp +++ b/libsrc/csg/solid.cpp @@ -689,7 +689,7 @@ namespace netgen case UNION: { int in1, in2, strin1, strin2; - Solid * tansol1, * tansol2; + Solid * tansol1 = 0, * tansol2 = 0; s1 -> RecTangentialSolid (p, tansol1, surfids, in1, strin1, eps); s2 -> RecTangentialSolid (p, tansol2, surfids, in2, strin2, eps); @@ -703,6 +703,11 @@ namespace netgen else if (tansol2) tansol = tansol2; } + else + { + delete tansol1; + delete tansol2; + } in = (in1 || in2); strin = (strin1 || strin2); break; diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp index b86666e0..f2462b7f 100644 --- a/libsrc/csg/specpoin.cpp +++ b/libsrc/csg/specpoin.cpp @@ -76,6 +76,10 @@ namespace netgen CalcSpecialPoints (const CSGeometry & ageometry, Array & apoints) { + static int timer = NgProfiler::CreateTimer ("CSG: find special points"); + NgProfiler::RegionTimer reg (timer); + + geometry = &ageometry; points = &apoints; @@ -231,7 +235,7 @@ namespace netgen // explicit solution for planes only and at most one quadratic if (numprim <= 5) { - int nplane = 0, nquad = 0, quadi = -1; + int nplane = 0, nquad = 0, quadi = -1, nsphere = 0; const QuadraticSurface *qsurf = 0, *qsurfi; for (int i = 0; i < numprim; i++) @@ -247,6 +251,9 @@ namespace netgen quadi = i; qsurf = qsurfi; } + + if (dynamic_cast (qsurfi)) + nsphere++; } /* @@ -374,6 +381,82 @@ namespace netgen return; } + + + + if (nsphere == numprim) // && calccp == false) + { + Array > pts; + Array surfids; + + + + for (int k1 = 0; k1 < numprim; k1++) + for (int k2 = 0; k2 < k1; k2++) + for (int k3 = 0; k3 < k2; k3++) + { + ComputeCrossPoints (dynamic_cast (geometry->GetSurface(locsurf[k1])), + dynamic_cast (geometry->GetSurface(locsurf[k2])), + dynamic_cast (geometry->GetSurface(locsurf[k3])), + pts); + + for (int j = 0; j < pts.Size(); j++) + if (Dist (pts[j], box.Center()) < box.Diam()/2) + { + Solid * tansol; + sol -> TangentialSolid (pts[j], tansol, surfids, 1e-9*size); + + if(!tansol) + continue; + + bool ok1 = false, ok2 = false, ok3 = false; + int rep1 = geometry->GetSurfaceClassRepresentant(locsurf[k1]); + int rep2 = geometry->GetSurfaceClassRepresentant(locsurf[k2]); + int rep3 = geometry->GetSurfaceClassRepresentant(locsurf[k3]); + for(int jj=0; jjGetSurfaceClassRepresentant(surfids[jj]); + if(actrep == rep1) ok1 = true; + if(actrep == rep2) ok2 = true; + if(actrep == rep3) ok3 = true; + } + + + if (tansol && ok1 && ok2 && ok3) + // if (sol -> IsIn (pts[j], 1e-6*size) && !sol->IsStrictIn (pts[j], 1e-6*size)) + { + if (AddPoint (pts[j], layer)) + (*testout) << "cross point found, 1: " << pts[j] << endl; + } + delete tansol; + } + } + + + for (int k1 = 0; k1 < numprim; k1++) + for (int k2 = 0; k2 < k1; k2++) + { + ComputeExtremalPoints (dynamic_cast (geometry->GetSurface(locsurf[k1])), + dynamic_cast (geometry->GetSurface(locsurf[k2])), + pts); + + for (int j = 0; j < pts.Size(); j++) + if (Dist (pts[j], box.Center()) < box.Diam()/2) + { + Solid * tansol; + sol -> TangentialSolid (pts[j], tansol, surfids, 1e-9*size); + if (tansol) + // sol -> IsIn (pts[j], 1e-6*size) && !sol->IsStrictIn (pts[j], 1e-6*size) ) + { + if (AddPoint (pts[j], layer)) + (*testout) << "extremal point found, spheres: " << pts[j] << endl; + } + delete tansol; + } + } + + return; + } } @@ -1143,6 +1226,83 @@ namespace netgen + void SpecialPointCalculation :: + ComputeCrossPoints (const Sphere * sphere1, + const Sphere * sphere2, + const Sphere * sphere3, + Array > & pts) + { + Mat<2,3> mat; + Mat<3,2> inv; + Vec<2> rhs; + Vec<3> sol, t; + Point<3> p0(0,0,0); + + pts.SetSize (0); + + + Point<3> c1 = sphere1 -> Center(); + Point<3> c2 = sphere2 -> Center(); + Point<3> c3 = sphere3 -> Center(); + double r1 = sphere1 -> Radius(); + double r2 = sphere2 -> Radius(); + double r3 = sphere3 -> Radius(); + + + Vec<3> a1 = c2-c1; + double b1 = 0.5 * (sqr(r1) - sqr(r2) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c2)) ); + + Vec<3> a2 = c3-c1; + double b2 = 0.5 * (sqr(r1) - sqr(r3) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c3)) ); + + + for (int j = 0; j < 3; j++) + { + mat(0,j) = a1(j); + mat(1,j) = a2(j); + } + + rhs(0) = b1; + rhs(1) = b2; + + + CalcInverse (mat, inv); + sol = inv * rhs; + t = Cross (mat.Row(0), mat.Row(1)); + + if (t.Length() > 1e-8) + { + Point<3> p (sol); + // quadratic on p + s t = 0 + double quad_a; + Vec<3> quad_b; + Mat<3> quad_c; + + quad_a = sphere1 -> CalcFunctionValue(p); + sphere1 -> CalcGradient (p, quad_b); + sphere1 -> CalcHesse (p, quad_c); + + double a, b, c; + a = quad_a; + b = quad_b * t; + c = 0.5 * t * (quad_c * t); + + // a + s b + s^2 c = 0; + double disc = b*b-4*a*c; + if (disc > 1e-10 * fabs (b)) + { + disc = sqrt (disc); + double s1 = (-b-disc) / (2*c); + double s2 = (-b+disc) / (2*c); + + pts.Append (p + s1 * t); + pts.Append (p + s2 * t); + } + } + } + + + @@ -1238,6 +1398,138 @@ namespace netgen + + void SpecialPointCalculation :: + ComputeExtremalPoints (const Sphere * sphere1, + const Sphere * sphere2, + Array > & pts) + { + // 3 equations: + // surf1 = 0 <===> |x-c1|^2 - r1^2 = 0; + // surf2 = 0 <===> |x-c2|^2 - r2^2 = 0; + // (grad 1 x grad 2)(i) = 0 <====> (x-p1) x (p1-p2) . e_i = 0; + + pts.SetSize (0); + + Point<3> c1 = sphere1 -> Center(); + Point<3> c2 = sphere2 -> Center(); + double r1 = sphere1 -> Radius(); + double r2 = sphere2 -> Radius(); + + /* + *testout << "\n\ncompute extremalpoint, sphere-sphere" << endl; + *testout << "c1 = " << c1 << ", r1 = " << r1 << endl; + *testout << "c2 = " << c2 << ", r2 = " << r2 << endl; + *testout << "dist = " << Abs (c2-c1) << ", r1+r2 = " << r1+r2 << endl; + */ + + Vec<3> v12 = c2 - c1; + + Vec<3> a1, a2; + double b1, b2; + + // eqn: ai . x = bi + + a1 = v12; + b1 = 0.5 * (sqr(r1) - sqr(r2) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c2)) ); + + int dir = 0; + for (int j = 1; j < 3; j++) + if (fabs (v12(j)) > v12(dir)) + dir = j; + + // *testout << "dir = " << dir << endl; + + Vec<3> ei = 0.0; + ei(dir) = 1; + a2 = Cross (v12, ei); + b2 = Vec<3>(c1) * a2; + + + Point<3> p0 (0,0,0); + double quad_a; + Vec<3> quad_b; + Mat<3> quad_c; + + quad_a = sphere1 -> CalcFunctionValue(p0); + sphere1 -> CalcGradient (p0, quad_b); + sphere1 -> CalcHesse (p0, quad_c); + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + quad_c(i,j) *= 0.5; + + + // find line of two linear equations: + + Vec<2> rhs; + Vec<3> sol; + Mat<2,3> mat; + + for (int j = 0; j < 3; j++) + { + mat(0,j) = a1(j); + mat(1,j) = a2(j); + } + rhs(0) = b1; + rhs(1) = b2; + + + // *testout << "mat = " << endl << mat << endl; + // *testout << "rhs = " << endl << rhs << endl; + + Vec<3> t = Cross (a1, a2); + if (Abs2(t) > 0) + { + mat.Solve (rhs, sol); + + /* + *testout << "sol = " << endl << sol << endl; + + *testout << "a * sol = " << mat * sol << endl; + + *testout << "c1-sol = " << Abs (Vec<3>(c1)-sol) << endl; + *testout << "c2-sol = " << Abs (Vec<3>(c2)-sol) << endl; + */ + + // solve quadratic equation along line sol + alpha t .... + double a = quad_a + quad_b * sol + sol * (quad_c * sol); + double b = quad_b * t + 2 * (sol * (quad_c * t)); + double c = t * (quad_c * t); + + // solve a + b alpha + c alpha^2: + + if (fabs (c) > 1e-32) + { + double disc = sqr (0.5*b/c) - a/c; + if (disc > 0) + { + disc = sqrt (disc); + double alpha1 = -0.5*b/c + disc; + double alpha2 = -0.5*b/c - disc; + + pts.Append (Point<3> (sol+alpha1*t)); + pts.Append (Point<3> (sol+alpha2*t)); + + // *testout << "pts = " << endl << pts << endl; + + /* + cout << "sol1 = " << sol + alpha1 * t + << ", sol2 = " << sol + alpha2 * t << endl; + */ + } + } + } + } + + + + + + + + + + /* bool SpecialPointCalculation :: ExtremalPointPossible (const Surface * f1, const Surface * f2, @@ -1445,11 +1737,11 @@ namespace netgen if (tlo->GetLayer() != apoints[i].GetLayer()) continue; + Solid * locsol; sol -> TangentialSolid (p, locsol, surfind, ideps*geomsize); - rep_surfind.SetSize (surfind.Size()); int num_indep_surfs = 0; @@ -1470,6 +1762,7 @@ namespace netgen #endif if (!locsol) continue; + // get all surface indices, if (surf) diff --git a/libsrc/csg/specpoin.hpp b/libsrc/csg/specpoin.hpp index bc36903c..8284ab3b 100644 --- a/libsrc/csg/specpoin.hpp +++ b/libsrc/csg/specpoin.hpp @@ -161,6 +161,11 @@ namespace netgen const QuadraticSurface * quadric, Array > & pts); + void ComputeExtremalPoints (const Sphere * sphere1, + const Sphere * sphere2, + Array > & pts); + + void ComputeCrossPoints (const Plane * plane1, const Plane * plane2, const Plane * plane3, @@ -170,6 +175,11 @@ namespace netgen const Plane * plane2, const QuadraticSurface * quadratic, Array > & pts); + + void ComputeCrossPoints (const Sphere * sphere1, + const Sphere * sphere2, + const Sphere * sphere3, + Array > & pts); }; } diff --git a/libsrc/meshing/adfront2.cpp b/libsrc/meshing/adfront2.cpp index fbe048fc..3b273bc7 100644 --- a/libsrc/meshing/adfront2.cpp +++ b/libsrc/meshing/adfront2.cpp @@ -94,8 +94,9 @@ namespace netgen if (mgi) cpointsearchtree.Insert (p, pi); - pointsearchtree.Insert (p, pi); - + if (pointonsurface) + pointsearchtree.Insert (p, pi); + return pi; } @@ -287,7 +288,9 @@ namespace netgen Array & lindex, double xh) { - // baselineindex += 1-lines.Begin(); + static int timer = NgProfiler::CreateTimer ("adfront2::GetLocals"); + NgProfiler::RegionTimer reg (timer); + int pstind; Point<3> midp, p0; @@ -296,13 +299,12 @@ namespace netgen p0 = points[pstind].P(); loclines.Append(lines[baselineindex].L()); - lindex.Append(baselineindex); // +1-lines.Begin()); + lindex.Append(baselineindex); - static Array nearlines; - nearlines.SetSize(0); - static Array nearpoints; - nearpoints.SetSize(0); + ArrayMem nearlines(0); + ArrayMem nearpoints(0); + // dominating costs !! linesearchtree.GetIntersecting (p0 - Vec3d(xh, xh, xh), p0 + Vec3d(xh, xh, xh), nearlines); @@ -310,12 +312,11 @@ namespace netgen pointsearchtree.GetIntersecting (p0 - Vec3d(xh, xh, xh), p0 + Vec3d(xh, xh, xh), nearpoints); - - - for (int ii = 1; ii <= nearlines.Size(); ii++) + + for (int ii = 0; ii < nearlines.Size(); ii++) { - int i = nearlines.Get(ii); - if (lines[i].Valid() && i != baselineindex) // + 1-lines.Begin()) + int i = nearlines[ii]; + if (lines[i].Valid() && i != baselineindex) { loclines.Append(lines[i].L()); lindex.Append(i); diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index 9626c05d..eea6c511 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -12,406 +12,413 @@ namespace netgen { -class Neighbour -{ - int nr[3]; - int orient[3]; + class Neighbour + { + int nr[3]; + int orient[3]; -public: - Neighbour () { nr[0] = nr[1] = nr[2] = -1; orient[0] = orient[1] = orient[2] = 0; } + public: + Neighbour () { ; } - void SetNr (int side, int anr) { nr[side-1] = anr; } - int GetNr (int side) { return nr[side-1]; } + void SetNr (int side, int anr) { nr[side] = anr; } + int GetNr (int side) { return nr[side]; } - void SetOrientation (int side, int aorient) { orient[side-1] = aorient; } - int GetOrientation (int side) { return orient[side-1]; } -}; + void SetOrientation (int side, int aorient) { orient[side] = aorient; } + int GetOrientation (int side) { return orient[side]; } + + + + /* + void SetNr1 (int side, int anr) { nr[side-1] = anr; } + int GetNr1 (int side) { return nr[side-1]; } + + void SetOrientation1 (int side, int aorient) { orient[side-1] = aorient; } + int GetOrientation1 (int side) { return orient[side-1]; } + */ + }; -class trionedge -{ -public: - int tnr; - int sidenr; + class trionedge + { + public: + int tnr; + int sidenr; - trionedge () { tnr = 0; sidenr = 0; } - trionedge (int atnr, int asidenr) + trionedge () { tnr = 0; sidenr = 0; } + trionedge (int atnr, int asidenr) { tnr = atnr; sidenr = asidenr; } -}; + }; -void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric) -{ - // return; - - if (!faceindex) - { - if (usemetric) - PrintMessage (3, "Edgeswapping, metric"); - else - PrintMessage (3, "Edgeswapping, topological"); - - for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) - { - EdgeSwapping (mesh, usemetric); - - if (multithread.terminate) - throw NgException ("Meshing stopped"); - } - - faceindex = 0; - mesh.CalcSurfacesOfNode(); - return; - } - - - static int timer = NgProfiler::CreateTimer ("EdgeSwapping 2D"); - NgProfiler::RegionTimer reg1 (timer); - - - int i, i2, j, j2; - bool should; - PointIndex pi; - - Array seia; - mesh.GetSurfaceElementsOfFace (faceindex, seia); - - for (i = 0; i < seia.Size(); i++) - if (mesh[seia[i]].GetNP() != 3) + void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric) + { + if (!faceindex) { - GenericImprove (mesh); + if (usemetric) + PrintMessage (3, "Edgeswapping, metric"); + else + PrintMessage (3, "Edgeswapping, topological"); + + for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) + { + EdgeSwapping (mesh, usemetric); + + if (multithread.terminate) + throw NgException ("Meshing stopped"); + } + + faceindex = 0; + mesh.CalcSurfacesOfNode(); return; } - int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); - Array neighbors(mesh.GetNSE()); - INDEX_2_HASHTABLE other(seia.Size() + 2); + static int timer = NgProfiler::CreateTimer ("EdgeSwapping 2D"); + NgProfiler::RegionTimer reg1 (timer); + + static int timerstart = NgProfiler::CreateTimer ("EdgeSwapping 2D start"); + NgProfiler::StartTimer (timerstart); - Array swapped(mesh.GetNSE()); - Array pdef(mesh.GetNP()); - Array pangle(mesh.GetNP()); + Array seia; + mesh.GetSurfaceElementsOfFace (faceindex, seia); - SurfaceElementIndex t1, t2; - int o1, o2; - - PointIndex pi1, pi2, pi3, pi4; - PointGeomInfo gi1, gi2, gi3, gi4; - - - int nswaps = 0; - int e, done; - double d; - Vec3d nv1, nv2; - double horder; - double loch(-1); - static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 }; - - pangle = 0; - - for (i = 0; i < seia.Size(); i++) - { - const Element2d & sel = mesh[seia[i]]; - for (j = 0; j < 3; j++) + for (int i = 0; i < seia.Size(); i++) + if (mesh[seia[i]].GetNP() != 3) { - POINTTYPE typ = mesh[sel[j]].Type(); - if (typ == FIXEDPOINT || typ == EDGEPOINT) - { - pangle[sel[j]] += - Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]], - mesh[sel[(j+2)%3]] - mesh[sel[j]]); - } + GenericImprove (mesh); + return; } - } - for (pi = PointIndex::BASE; - pi < mesh.GetNP()+PointIndex::BASE; pi++) - { - if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT) - pdef[pi] = -6; - else - for (j = 0; j < 8; j++) - if (pangle[pi] >= minangle[j]) - pdef[pi] = -1-j; - } + int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); - for (i = 0; i < seia.Size(); i++) - { - const Element2d & sel = mesh[seia[i]]; - for (j = 0; j < 3; j++) - pdef[sel[j]]++; - } + Array neighbors(mesh.GetNSE()); + INDEX_2_HASHTABLE other(seia.Size() + 2); - for (i = 0; i < seia.Size(); i++) - { - //const Element2d & sel = mesh[seia[i]]; - for (j = 0; j < 3; j++) - { - neighbors[seia[i]].SetNr (j+1, -1); - neighbors[seia[i]].SetOrientation (j+1, 0); - } - } - /* - Array normals(mesh.GetNP()); - for (i = 1; i <= mesh.GetNSE(); i++) - { + Array swapped(mesh.GetNSE()); + Array pdef(mesh.GetNP()); + Array pangle(mesh.GetNP()); + + + // int e; + // double d; + // Vec3d nv1, nv2; + + // double loch(-1); + static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 }; + + + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + pangle[sel[j]] = 0.0; + } + // pangle = 0; + + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + { + POINTTYPE typ = mesh[sel[j]].Type(); + if (typ == FIXEDPOINT || typ == EDGEPOINT) + { + pangle[sel[j]] += + Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]], + mesh[sel[(j+2)%3]] - mesh[sel[j]]); + } + } + } + + // for (PointIndex pi = PointIndex::BASE; + // pi < mesh.GetNP()+PointIndex::BASE; pi++) + + // pdef = 0; + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + { + PointIndex pi = sel[j]; + if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT) + pdef[pi] = -6; + else + for (int j = 0; j < 8; j++) + if (pangle[pi] >= minangle[j]) + pdef[pi] = -1-j; + } + } + + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + pdef[sel[j]]++; + } + + for (int i = 0; i < seia.Size(); i++) + { + for (int j = 0; j < 3; j++) + { + neighbors[seia[i]].SetNr (j, -1); + neighbors[seia[i]].SetOrientation (j, 0); + } + } + + /* + Array normals(mesh.GetNP()); + for (i = 1; i <= mesh.GetNSE(); i++) + { Element2d & hel = mesh.SurfaceElement(i); if (hel.GetIndex() == faceindex) - for (k = 1; k <= 3; k++) + for (k = 1; k <= 3; k++) + { + int pi = hel.PNum(k); + SelectSurfaceOfPoint (mesh.Point(pi), hel.GeomInfoPi(k)); + int surfi = mesh.GetFaceDescriptor(faceindex).SurfNr(); + GetNormalVector (surfi, mesh.Point(pi), normals.Elem(pi)); + normals.Elem(pi) /= normals.Elem(pi).Length(); + } + } + */ + + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + + for (int j = 0; j < 3; j++) { - int pi = hel.PNum(k); - SelectSurfaceOfPoint (mesh.Point(pi), hel.GeomInfoPi(k)); - int surfi = mesh.GetFaceDescriptor(faceindex).SurfNr(); - GetNormalVector (surfi, mesh.Point(pi), normals.Elem(pi)); - normals.Elem(pi) /= normals.Elem(pi).Length(); - } - } - */ - - - for (i = 0; i < seia.Size(); i++) - { - const Element2d & sel = mesh[seia[i]]; - - for (j = 1; j <= 3; j++) - { - pi1 = sel.PNumMod(j+1); - pi2 = sel.PNumMod(j+2); + PointIndex pi1 = sel.PNumMod(j+2); + PointIndex pi2 = sel.PNumMod(j+3); - loch = mesh.GetH(mesh[pi1]); + // double loch = mesh.GetH(mesh[pi1]); - INDEX_2 edge(pi1, pi2); - edge.Sort(); + INDEX_2 edge(pi1, pi2); + edge.Sort(); - if (mesh.IsSegment (pi1, pi2)) - continue; + if (mesh.IsSegment (pi1, pi2)) + continue; - /* - if (segments.Used (edge)) - continue; - */ - INDEX_2 ii2 (pi1, pi2); - if (other.Used (ii2)) - { - // INDEX_2 i2s(ii2); - // i2s.Sort(); + /* + if (segments.Used (edge)) + continue; + */ + INDEX_2 ii2 (pi1, pi2); + if (other.Used (ii2)) + { + // INDEX_2 i2s(ii2); + // i2s.Sort(); - i2 = other.Get(ii2).tnr; - j2 = other.Get(ii2).sidenr; + int i2 = other.Get(ii2).tnr; + int j2 = other.Get(ii2).sidenr; - neighbors[seia[i]].SetNr (j, i2); - neighbors[seia[i]].SetOrientation (j, j2); - neighbors[i2].SetNr (j2, seia[i]); - neighbors[i2].SetOrientation (j2, j); - } - else - { - other.Set (INDEX_2 (pi2, pi1), trionedge (seia[i], j)); - } - } - } + neighbors[seia[i]].SetNr (j, i2); + neighbors[seia[i]].SetOrientation (j, j2); + neighbors[i2].SetNr (j2, seia[i]); + neighbors[i2].SetOrientation (j2, j); + } + else + { + other.Set (INDEX_2 (pi2, pi1), trionedge (seia[i], j)); + } + } + } - for (i = 0; i < seia.Size(); i++) - swapped[seia[i]] = 0; + for (int i = 0; i < seia.Size(); i++) + swapped[seia[i]] = 0; + + NgProfiler::StopTimer (timerstart); + - int t = 4; - done = 0; - while (!done && t >= 2) - { - for (i = 0; i < seia.Size(); i++) - { - t1 = seia[i]; + int t = 4; + int done = 0; + while (!done && t >= 2) + { + for (int i = 0; i < seia.Size(); i++) + { + SurfaceElementIndex t1 = seia[i]; - if (mesh[t1].IsDeleted()) - continue; + if (mesh[t1].IsDeleted()) + continue; - if (mesh[t1].GetIndex() != faceindex) - continue; + if (mesh[t1].GetIndex() != faceindex) + continue; - if (multithread.terminate) - throw NgException ("Meshing stopped"); + if (multithread.terminate) + throw NgException ("Meshing stopped"); - for (o1 = 1; o1 <= 3; o1++) - { - t2 = neighbors[t1].GetNr (o1); - o2 = neighbors[t1].GetOrientation (o1); + for (int o1 = 0; o1 < 3; o1++) + { + bool should; - if (t2 == -1) continue; - if (swapped[t1] || swapped[t2]) continue; + + SurfaceElementIndex t2 = neighbors[t1].GetNr (o1); + int o2 = neighbors[t1].GetOrientation (o1); + + if (t2 == -1) continue; + if (swapped[t1] || swapped[t2]) continue; - pi1 = mesh[t1].PNumMod(o1+1); - pi2 = mesh[t1].PNumMod(o1+2); - pi3 = mesh[t1].PNumMod(o1); - pi4 = mesh[t2].PNumMod(o2); + PointIndex pi1 = mesh[t1].PNumMod(o1+1+1); + PointIndex pi2 = mesh[t1].PNumMod(o1+1+2); + PointIndex pi3 = mesh[t1].PNumMod(o1+1); + PointIndex pi4 = mesh[t2].PNumMod(o2+1); - gi1 = mesh[t1].GeomInfoPiMod(o1+1); - gi2 = mesh[t1].GeomInfoPiMod(o1+2); - gi3 = mesh[t1].GeomInfoPiMod(o1); - gi4 = mesh[t2].GeomInfoPiMod(o2); + PointGeomInfo gi1 = mesh[t1].GeomInfoPiMod(o1+1+1); + PointGeomInfo gi2 = mesh[t1].GeomInfoPiMod(o1+1+2); + PointGeomInfo gi3 = mesh[t1].GeomInfoPiMod(o1+1); + PointGeomInfo gi4 = mesh[t2].GeomInfoPiMod(o2+1); - bool allowswap = true; + bool allowswap = true; - Vec3d auxvec1,auxvec2; + Vec<3> auxvec1 = mesh[pi3]-mesh[pi4]; + Vec<3> auxvec2 = mesh[pi1]-mesh[pi4]; - auxvec1 = mesh.Point(pi3)-mesh.Point(pi4); - auxvec2 = mesh.Point(pi1)-mesh.Point(pi4); - allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4; + allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4; - if(!allowswap) - continue; + if(!allowswap) + continue; - // normal of new - nv1 = Cross (auxvec1, - auxvec2); + // normal of new + Vec<3> nv1 = Cross (auxvec1, auxvec2); - auxvec1 = mesh.Point(pi4)-mesh.Point(pi3); - auxvec2 = mesh.Point(pi2)-mesh.Point(pi3); - allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4; + auxvec1 = mesh.Point(pi4)-mesh.Point(pi3); + auxvec2 = mesh.Point(pi2)-mesh.Point(pi3); + allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4; - if(!allowswap) - continue; + if(!allowswap) + continue; - nv2 = Cross (auxvec1, - auxvec2); + Vec<3> nv2 = Cross (auxvec1, auxvec2); - // normals of original - Vec3d nv3, nv4; - nv3 = Cross (mesh.Point(pi1)-mesh.Point(pi4), - mesh.Point(pi2)-mesh.Point(pi4)); - nv4 = Cross (mesh.Point(pi2)-mesh.Point(pi3), - mesh.Point(pi1)-mesh.Point(pi3)); + // normals of original + Vec<3> nv3 = Cross (mesh[pi1]-mesh[pi4], mesh[pi2]-mesh[pi4]); + Vec<3> nv4 = Cross (mesh[pi2]-mesh[pi3], mesh[pi1]-mesh[pi3]); - nv3 *= -1; - nv4 *= -1; - nv3.Normalize(); - nv4.Normalize(); + nv3 *= -1; + nv4 *= -1; + nv3.Normalize(); + nv4.Normalize(); - nv1.Normalize(); - nv2.Normalize(); + nv1.Normalize(); + nv2.Normalize(); - Vec<3> nvp3, nvp4; - SelectSurfaceOfPoint (mesh.Point(pi3), gi3); - GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3); + Vec<3> nvp3, nvp4; + SelectSurfaceOfPoint (mesh.Point(pi3), gi3); + GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3); - nvp3.Normalize(); + nvp3.Normalize(); - SelectSurfaceOfPoint (mesh.Point(pi4), gi4); - GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4); + SelectSurfaceOfPoint (mesh.Point(pi4), gi4); + GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4); - nvp4.Normalize(); + nvp4.Normalize(); - double critval = cos (M_PI / 6); // 30 degree - allowswap = allowswap && - (nv1 * nvp3 > critval) && - (nv1 * nvp4 > critval) && - (nv2 * nvp3 > critval) && - (nv2 * nvp4 > critval) && - (nvp3 * nv3 > critval) && - (nvp4 * nv4 > critval); + double critval = cos (M_PI / 6); // 30 degree + allowswap = allowswap && + (nv1 * nvp3 > critval) && + (nv1 * nvp4 > critval) && + (nv2 * nvp3 > critval) && + (nv2 * nvp4 > critval) && + (nvp3 * nv3 > critval) && + (nvp4 * nv4 > critval); - horder = Dist (mesh.Point(pi1), mesh.Point(pi2)); + double horder = Dist (mesh.Point(pi1), mesh.Point(pi2)); - if ( // nv1 * nv2 >= 0 && - nv1.Length() > 1e-3 * horder * horder && - nv2.Length() > 1e-3 * horder * horder && - allowswap ) - { - if (!usemetric) - { - e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4]; - d = - Dist2 (mesh.Point(pi1), mesh.Point(pi2)) - - Dist2 (mesh.Point(pi3), mesh.Point(pi4)); + if ( // nv1 * nv2 >= 0 && + nv1.Length() > 1e-3 * horder * horder && + nv2.Length() > 1e-3 * horder * horder && + allowswap ) + { + if (!usemetric) + { + int e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4]; + double d = + Dist2 (mesh.Point(pi1), mesh.Point(pi2)) - + Dist2 (mesh.Point(pi3), mesh.Point(pi4)); - should = e >= t && (e > 2 || d > 0); - } - else - { - should = - CalcTriangleBadness (mesh.Point(pi4), mesh.Point(pi3), mesh.Point(pi1), - metricweight, loch) + - CalcTriangleBadness (mesh.Point(pi3), mesh.Point(pi4), mesh.Point(pi2), - metricweight, loch) < - CalcTriangleBadness (mesh.Point(pi1), mesh.Point(pi2), mesh.Point(pi3), - metricweight, loch) + - CalcTriangleBadness (mesh.Point(pi2), mesh.Point(pi1), mesh.Point(pi4), - metricweight, loch); + should = e >= t && (e > 2 || d > 0); + } + else + { + double loch = mesh.GetH(mesh[pi1]); + should = + CalcTriangleBadness (mesh.Point(pi4), mesh.Point(pi3), mesh.Point(pi1), + metricweight, loch) + + CalcTriangleBadness (mesh.Point(pi3), mesh.Point(pi4), mesh.Point(pi2), + metricweight, loch) < + CalcTriangleBadness (mesh.Point(pi1), mesh.Point(pi2), mesh.Point(pi3), + metricweight, loch) + + CalcTriangleBadness (mesh.Point(pi2), mesh.Point(pi1), mesh.Point(pi4), + metricweight, loch); - } + } - if (allowswap) - { - Element2d sw1 (pi4, pi3, pi1); - Element2d sw2 (pi3, pi4, pi2); + if (allowswap) + { + Element2d sw1 (pi4, pi3, pi1); + Element2d sw2 (pi3, pi4, pi2); - int legal1 = - mesh.LegalTrig (mesh.SurfaceElement (t1)) + - mesh.LegalTrig (mesh.SurfaceElement (t2)); - int legal2 = - mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2); + int legal1 = + mesh.LegalTrig (mesh.SurfaceElement (t1)) + + mesh.LegalTrig (mesh.SurfaceElement (t2)); + int legal2 = + mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2); - if (legal1 < legal2) should = 1; - if (legal2 < legal1) should = 0; - } + if (legal1 < legal2) should = 1; + if (legal2 < legal1) should = 0; + } - if (should) - { - // do swapping ! + if (should) + { + // do swapping ! - // cout << "swap " << endl; + done = 1; + + mesh[t1].PNum(1) = pi1; + mesh[t1].PNum(2) = pi4; + mesh[t1].PNum(3) = pi3; + + mesh[t2].PNum(1) = pi2; + mesh[t2].PNum(2) = pi3; + mesh[t2].PNum(3) = pi4; + + mesh[t1].GeomInfoPi(1) = gi1; + mesh[t1].GeomInfoPi(2) = gi4; + mesh[t1].GeomInfoPi(3) = gi3; + + mesh[t2].GeomInfoPi(1) = gi2; + mesh[t2].GeomInfoPi(2) = gi3; + mesh[t2].GeomInfoPi(3) = gi4; + + pdef[pi1]--; + pdef[pi2]--; + pdef[pi3]++; + pdef[pi4]++; + + swapped[t1] = 1; + swapped[t2] = 1; + } + } + } + } + t--; + } - nswaps ++; - - // testout << "nv1 = " << nv1 << " nv2 = " << nv2 << endl; - - - done = 1; - - mesh[t1].PNum(1) = pi1; - mesh[t1].PNum(2) = pi4; - mesh[t1].PNum(3) = pi3; - - mesh[t2].PNum(1) = pi2; - mesh[t2].PNum(2) = pi3; - mesh[t2].PNum(3) = pi4; - - mesh[t1].GeomInfoPi(1) = gi1; - mesh[t1].GeomInfoPi(2) = gi4; - mesh[t1].GeomInfoPi(3) = gi3; - - mesh[t2].GeomInfoPi(1) = gi2; - mesh[t2].GeomInfoPi(2) = gi3; - mesh[t2].GeomInfoPi(3) = gi4; - - pdef[pi1]--; - pdef[pi2]--; - pdef[pi3]++; - pdef[pi4]++; - - swapped[t1] = 1; - swapped[t2] = 1; - } - } - } - } - t--; - } - - mesh.SetNextTimeStamp(); -} + mesh.SetNextTimeStamp(); + } @@ -420,412 +427,428 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric) -void MeshOptimize2d :: CombineImprove (Mesh & mesh) -{ - if (!faceindex) - { - PrintMessage (3, "Combine improve"); + void MeshOptimize2d :: CombineImprove (Mesh & mesh) + { + if (!faceindex) + { + PrintMessage (3, "Combine improve"); - for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) - { - CombineImprove (mesh); + for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++) + { + CombineImprove (mesh); - if (multithread.terminate) - throw NgException ("Meshing stopped"); - } - faceindex = 0; - return; - } + if (multithread.terminate) + throw NgException ("Meshing stopped"); + } + faceindex = 0; + return; + } - static int timer = NgProfiler::CreateTimer ("Combineimprove 2D"); - NgProfiler::RegionTimer reg (timer); + static int timer = NgProfiler::CreateTimer ("Combineimprove 2D"); + NgProfiler::RegionTimer reg (timer); + + static int timerstart = NgProfiler::CreateTimer ("Combineimprove 2D start"); + NgProfiler::StartTimer (timerstart); + + + static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1"); + NgProfiler::StartTimer (timerstart1); - int i, j, k, l; - PointIndex pi; - SurfaceElementIndex sei; + // int i, j, k, l; + // PointIndex pi; + // SurfaceElementIndex sei; - Array seia; - mesh.GetSurfaceElementsOfFace (faceindex, seia); + Array seia; + mesh.GetSurfaceElementsOfFace (faceindex, seia); - for (i = 0; i < seia.Size(); i++) - if (mesh[seia[i]].GetNP() != 3) - return; + for (int i = 0; i < seia.Size(); i++) + if (mesh[seia[i]].GetNP() != 3) + return; - int surfnr = 0; - if (faceindex) - surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); + int surfnr = 0; + if (faceindex) + surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr(); - PointIndex pi1, pi2; - MeshPoint p1, p2, pnew; - double bad1, bad2; - Vec<3> nv; + // PointIndex pi1, pi2; + // MeshPoint p1, p2, pnew; + double bad1, bad2; + Vec<3> nv; - int np = mesh.GetNP(); - //int nse = mesh.GetNSE(); + int np = mesh.GetNP(); + //int nse = mesh.GetNSE(); - TABLE elementsonnode(np); - Array hasonepi, hasbothpi; + TABLE elementsonnode(np); + Array hasonepi, hasbothpi; - for (i = 0; i < seia.Size(); i++) - { - Element2d & el = mesh[seia[i]]; - for (j = 0; j < el.GetNP(); j++) - { + for (int i = 0; i < seia.Size(); i++) + { + Element2d & el = mesh[seia[i]]; + for (int j = 0; j < el.GetNP(); j++) elementsonnode.Add (el[j], seia[i]); - } - } + } + Array fixed(np); + fixed = false; - Array fixed(np); - fixed = false; + NgProfiler::StopTimer (timerstart1); - SegmentIndex si; - for (si = 0; si < mesh.GetNSeg(); si++) - { - INDEX_2 i2(mesh[si][0], mesh[si][1]); - fixed[i2.I1()] = true; - fixed[i2.I2()] = true; - } + /* + for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++) + { + INDEX_2 i2(mesh[si][0], mesh[si][1]); + fixed[i2.I1()] = true; + fixed[i2.I2()] = true; + } + */ - for(i = 0; i,PointIndex::BASE> normals(np); - - for (pi = PointIndex::BASE; - pi < np + PointIndex::BASE; pi++) - { - if (elementsonnode[pi].Size()) - { - Element2d & hel = mesh[elementsonnode[pi][0]]; - for (k = 0; k < 3; k++) - if (hel[k] == pi) - { - SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1)); - GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]); - break; + for (int i = 0; i < seia.Size(); i++) + { + Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < sel.GetNP(); j++) + { + PointIndex pi1 = sel.PNumMod(j+2); + PointIndex pi2 = sel.PNumMod(j+3); + if (mesh.IsSegment (pi1, pi2)) + { + fixed[pi1] = true; + fixed[pi2] = true; } - if (k == 3) - { - cerr << "Neuer Fehler von Joachim, code 17121" << endl; - } - } - } + } + } - for (i = 0; i < seia.Size(); i++) - { - sei = seia[i]; - Element2d & elem = mesh[sei]; - if (elem.IsDeleted()) continue; + for(int i = 0; i < mesh.LockedPoints().Size(); i++) + fixed[mesh.LockedPoints()[i]] = true; - for (j = 0; j < 3; j++) - { - pi1 = elem[j]; - pi2 = elem[(j+1) % 3]; - if (pi1 < PointIndex::BASE || - pi2 < PointIndex::BASE) - continue; - /* - INDEX_2 i2(pi1, pi2); - i2.Sort(); - if (segmentht.Used(i2)) - continue; - */ + Array,PointIndex::BASE> normals(np); - bool debugflag = 0; - - if (debugflag) - { - (*testout) << "Combineimprove, face = " << faceindex - << "pi1 = " << pi1 << " pi2 = " << pi2 << endl; - } - - /* - // save version: - if (fixed.Get(pi1) || fixed.Get(pi2)) - continue; - if (pi2 < pi1) swap (pi1, pi2); - */ - - // more general - if (fixed[pi2]) - Swap (pi1, pi2); - - if (fixed[pi2]) - continue; - - double loch = mesh.GetH (mesh[pi1]); - - INDEX_2 si2 (pi1, pi2); - si2.Sort(); - - /* - if (edgetested.Used (si2)) - continue; - edgetested.Set (si2, 1); - */ - - hasonepi.SetSize(0); - hasbothpi.SetSize(0); - - for (k = 0; k < elementsonnode[pi1].Size(); k++) - { - const Element2d & el2 = mesh[elementsonnode[pi1][k]]; - - if (el2.IsDeleted()) continue; - - if (el2[0] == pi2 || el2[1] == pi2 || el2[2] == pi2) + for (PointIndex pi = PointIndex::BASE; + pi < np + PointIndex::BASE; pi++) + { + if (elementsonnode[pi].Size()) + { + Element2d & hel = mesh[elementsonnode[pi][0]]; + for (int k = 0; k < 3; k++) + if (hel[k] == pi) { - hasbothpi.Append (elementsonnode[pi1][k]); - nv = Cross (Vec3d (mesh[el2[0]], mesh[el2[1]]), - Vec3d (mesh[el2[0]], mesh[el2[2]])); + SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1)); + GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]); + break; } - else - { - hasonepi.Append (elementsonnode[pi1][k]); - } - } + } + } + NgProfiler::StopTimer (timerstart); - Element2d & hel = mesh[hasbothpi[0]]; - for (k = 0; k < 3; k++) - if (hel[k] == pi1) + for (int i = 0; i < seia.Size(); i++) + { + SurfaceElementIndex sei = seia[i]; + Element2d & elem = mesh[sei]; + if (elem.IsDeleted()) continue; + + for (int j = 0; j < 3; j++) + { + PointIndex pi1 = elem[j]; + PointIndex pi2 = elem[(j+1) % 3]; + + if (pi1 < PointIndex::BASE || + pi2 < PointIndex::BASE) + continue; + + /* + INDEX_2 i2(pi1, pi2); + i2.Sort(); + if (segmentht.Used(i2)) + continue; + */ + + bool debugflag = 0; + + if (debugflag) { - SelectSurfaceOfPoint (mesh[pi1], - hel.GeomInfoPi(k+1)); - GetNormalVector (surfnr, mesh[pi1], hel.GeomInfoPi(k+1), nv); - break; + (*testout) << "Combineimprove, face = " << faceindex + << "pi1 = " << pi1 << " pi2 = " << pi2 << endl; } - if (k == 3) - { - cerr << "Neuer Fehler von Joachim, code 32434" << endl; - } + + /* + // save version: + if (fixed.Get(pi1) || fixed.Get(pi2)) + continue; + if (pi2 < pi1) swap (pi1, pi2); + */ + + // more general + if (fixed[pi2]) + Swap (pi1, pi2); + + if (fixed[pi2]) + continue; + + double loch = mesh.GetH (mesh[pi1]); + + INDEX_2 si2 (pi1, pi2); + si2.Sort(); + + /* + if (edgetested.Used (si2)) + continue; + edgetested.Set (si2, 1); + */ + + hasonepi.SetSize(0); + hasbothpi.SetSize(0); + + for (int k = 0; k < elementsonnode[pi1].Size(); k++) + { + const Element2d & el2 = mesh[elementsonnode[pi1][k]]; + + if (el2.IsDeleted()) continue; + + if (el2[0] == pi2 || el2[1] == pi2 || el2[2] == pi2) + { + hasbothpi.Append (elementsonnode[pi1][k]); + nv = Cross (Vec3d (mesh[el2[0]], mesh[el2[1]]), + Vec3d (mesh[el2[0]], mesh[el2[2]])); + } + else + { + hasonepi.Append (elementsonnode[pi1][k]); + } + } - // nv = normals.Get(pi1); + Element2d & hel = mesh[hasbothpi[0]]; + for (int k = 0; k < 3; k++) + if (hel[k] == pi1) + { + SelectSurfaceOfPoint (mesh[pi1], + hel.GeomInfoPi(k+1)); + GetNormalVector (surfnr, mesh[pi1], hel.GeomInfoPi(k+1), nv); + break; + } + + // nv = normals.Get(pi1); - for (k = 0; k < elementsonnode[pi2].Size(); k++) - { - const Element2d & el2 = mesh[elementsonnode[pi2][k]]; - if (el2.IsDeleted()) continue; + for (int k = 0; k < elementsonnode[pi2].Size(); k++) + { + const Element2d & el2 = mesh[elementsonnode[pi2][k]]; + if (el2.IsDeleted()) continue; - if (el2[0] == pi1 || el2[1] == pi1 || el2[2] == pi1) - ; - else - hasonepi.Append (elementsonnode[pi2][k]); - } + if (el2[0] == pi1 || el2[1] == pi1 || el2[2] == pi1) + ; + else + hasonepi.Append (elementsonnode[pi2][k]); + } - bad1 = 0; - int illegal1 = 0, illegal2 = 0; - for (k = 0; k < hasonepi.Size(); k++) - { - const Element2d & el = mesh[hasonepi[k]]; - bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], - nv, -1, loch); - illegal1 += 1-mesh.LegalTrig(el); - } + bad1 = 0; + int illegal1 = 0, illegal2 = 0; + for (int k = 0; k < hasonepi.Size(); k++) + { + const Element2d & el = mesh[hasonepi[k]]; + bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], + nv, -1, loch); + illegal1 += 1-mesh.LegalTrig(el); + } - for (k = 0; k < hasbothpi.Size(); k++) - { - const Element2d & el = mesh[hasbothpi[k]]; - bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], - nv, -1, loch); - illegal1 += 1-mesh.LegalTrig(el); - } - bad1 /= (hasonepi.Size()+hasbothpi.Size()); + for (int k = 0; k < hasbothpi.Size(); k++) + { + const Element2d & el = mesh[hasbothpi[k]]; + bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], + nv, -1, loch); + illegal1 += 1-mesh.LegalTrig(el); + } + bad1 /= (hasonepi.Size()+hasbothpi.Size()); - p1 = mesh[pi1]; - p2 = mesh[pi2]; + MeshPoint p1 = mesh[pi1]; + MeshPoint p2 = mesh[pi2]; - pnew = p1; - mesh[pi1] = pnew; - mesh[pi2] = pnew; + MeshPoint pnew = p1; + mesh[pi1] = pnew; + mesh[pi2] = pnew; - bad2 = 0; - for (k = 0; k < hasonepi.Size(); k++) - { - Element2d & el = mesh[hasonepi[k]]; - double err = - CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], - nv, -1, loch); - bad2 += err; + bad2 = 0; + for (int k = 0; k < hasonepi.Size(); k++) + { + Element2d & el = mesh[hasonepi[k]]; + double err = + CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], + nv, -1, loch); + bad2 += err; - Vec<3> hnv = Cross (Vec3d (mesh[el[0]], - mesh[el[1]]), - Vec3d (mesh[el[0]], - mesh[el[2]])); - if (hnv * nv < 0) - bad2 += 1e10; - - for (l = 0; l < 3; l++) - if ( (normals[el[l]] * nv) < 0.5) + Vec<3> hnv = Cross (Vec3d (mesh[el[0]], + mesh[el[1]]), + Vec3d (mesh[el[0]], + mesh[el[2]])); + if (hnv * nv < 0) bad2 += 1e10; + + for (int l = 0; l < 3; l++) + if ( (normals[el[l]] * nv) < 0.5) + bad2 += 1e10; - illegal2 += 1-mesh.LegalTrig(el); - } - bad2 /= hasonepi.Size(); + illegal2 += 1-mesh.LegalTrig(el); + } + bad2 /= hasonepi.Size(); - mesh[pi1] = p1; - mesh[pi2] = p2; + mesh[pi1] = p1; + mesh[pi2] = p2; - if (debugflag) - { - (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl; - } + if (debugflag) + { + (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl; + } - bool should = (bad2 < bad1 && bad2 < 1e4); - if (bad2 < 1e4) - { - if (illegal1 > illegal2) should = 1; - if (illegal2 > illegal1) should = 0; - } + bool should = (bad2 < bad1 && bad2 < 1e4); + if (bad2 < 1e4) + { + if (illegal1 > illegal2) should = 1; + if (illegal2 > illegal1) should = 0; + } - if (should) - { - // (*testout) << "combine !" << endl; - // (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl; + if (should) + { + // (*testout) << "combine !" << endl; + // (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl; - mesh[pi1] = pnew; - PointGeomInfo gi; - bool gi_set(false); + mesh[pi1] = pnew; + PointGeomInfo gi; + bool gi_set(false); - Element2d *el1p(NULL); - l=0; - while(mesh[elementsonnode[pi1][l]].IsDeleted() && lGetNP(); l++) - if ((*el1p)[l] == pi1) + for (l = 0; l < el1p->GetNP(); l++) + if ((*el1p)[l] == pi1) + { + gi = el1p->GeomInfoPi (l+1); + gi_set = true; + } + + // (*testout) << "Connect point " << pi2 << " to " << pi1 << "\n"; + for (int k = 0; k < elementsonnode[pi2].Size(); k++) { - gi = el1p->GeomInfoPi (l+1); - gi_set = true; + Element2d & el = mesh[elementsonnode[pi2][k]]; + if(el.IsDeleted()) continue; + elementsonnode.Add (pi1, elementsonnode[pi2][k]); + + bool haspi1 = 0; + for (l = 0; l < el.GetNP(); l++) + if (el[l] == pi1) + haspi1 = 1; + if (haspi1) continue; + + for (int l = 0; l < el.GetNP(); l++) + { + if (el[l] == pi2) + { + el[l] = pi1; + el.GeomInfoPi (l+1) = gi; + } + + fixed[el[l]] = true; + } } - // (*testout) << "Connect point " << pi2 << " to " << pi1 << "\n"; - for (k = 0; k < elementsonnode[pi2].Size(); k++) - { - Element2d & el = mesh[elementsonnode[pi2][k]]; - if(el.IsDeleted()) continue; - elementsonnode.Add (pi1, elementsonnode[pi2][k]); - - bool haspi1 = 0; - for (l = 0; l < el.GetNP(); l++) - if (el[l] == pi1) - haspi1 = 1; - if (haspi1) continue; - - for (l = 0; l < el.GetNP(); l++) - { - if (el[l] == pi2) - { - el[l] = pi1; - el.GeomInfoPi (l+1) = gi; - } - - fixed[el[l]] = true; - } - } - - /* - for (k = 0; k < hasbothpi.Size(); k++) - { + /* + for (k = 0; k < hasbothpi.Size(); k++) + { cout << mesh[hasbothpi[k]] << endl; for (l = 0; l < 3; l++) - cout << mesh[mesh[hasbothpi[k]][l]] << " "; + cout << mesh[mesh[hasbothpi[k]][l]] << " "; cout << endl; - } - */ + } + */ - for (k = 0; k < hasbothpi.Size(); k++) - { - mesh[hasbothpi[k]].Delete(); - /* - for (l = 0; l < 4; l++) - mesh[hasbothpi[k]][l] = PointIndex::BASE-1; - */ - } + for (int k = 0; k < hasbothpi.Size(); k++) + { + mesh[hasbothpi[k]].Delete(); + /* + for (l = 0; l < 4; l++) + mesh[hasbothpi[k]][l] = PointIndex::BASE-1; + */ + } - } - } - } + } + } + } - // mesh.Compress(); - mesh.SetNextTimeStamp(); -} + // mesh.Compress(); + mesh.SetNextTimeStamp(); + } -void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh) -{ - // Check angles between elements and normals at corners - /* + void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh) + { + // Check angles between elements and normals at corners + /* - int i, j; - int ne = mesh.GetNSE(); - int surfnr; + int i, j; + int ne = mesh.GetNSE(); + int surfnr; - Vec3d n, ng; - Array ngs(3); + Vec3d n, ng; + Array ngs(3); - (*mycout) << "Check Surface Approxiamtion" << endl; - (*testout) << "Check Surface Approxiamtion" << endl; + (*mycout) << "Check Surface Approxiamtion" << endl; + (*testout) << "Check Surface Approxiamtion" << endl; - for (i = 1; i <= ne; i++) + for (i = 1; i <= ne; i++) { - const Element2d & el = mesh.SurfaceElement(i); - surfnr = mesh.GetFaceDescriptor (el.GetIndex()).SurfNr(); - Vec3d n = Cross (mesh.Point (el.PNum(1)) - mesh.Point (el.PNum(2)), - mesh.Point (el.PNum(1)) - mesh.Point (el.PNum(3))); - n /= n.Length(); + const Element2d & el = mesh.SurfaceElement(i); + surfnr = mesh.GetFaceDescriptor (el.GetIndex()).SurfNr(); + Vec3d n = Cross (mesh.Point (el.PNum(1)) - mesh.Point (el.PNum(2)), + mesh.Point (el.PNum(1)) - mesh.Point (el.PNum(3))); + n /= n.Length(); - for (j = 1; j <= el.GetNP(); j++) - { - SelectSurfaceOfPoint (mesh.Point(el.PNum(j)), el.GeomInfoPi(j)); - GetNormalVector (surfnr, mesh.Point(el.PNum(j)), ng); - ng /= ng.Length(); - ngs.Elem(j) = ng; + for (j = 1; j <= el.GetNP(); j++) + { + SelectSurfaceOfPoint (mesh.Point(el.PNum(j)), el.GeomInfoPi(j)); + GetNormalVector (surfnr, mesh.Point(el.PNum(j)), ng); + ng /= ng.Length(); + ngs.Elem(j) = ng; - double angle = (180.0 / M_PI) * Angle (n, ng); - if (angle > 60) - { - (*testout) << "el " << i << " node " << el.PNum(j) - << "has angle = " << angle << endl; - } - } - - for (j = 1; j <= 3; j++) - { - double angle = (180.0 / M_PI) * Angle (ngs.Get(j), ngs.Get(j%3+1)); - if (angle > 60) - { - (*testout) << "el " << i << " node-node " - << ngs.Get(j) << " - " << ngs.Get(j%3+1) - << " has angle = " << angle << endl; - } - } + double angle = (180.0 / M_PI) * Angle (n, ng); + if (angle > 60) + { + (*testout) << "el " << i << " node " << el.PNum(j) + << "has angle = " << angle << endl; } - */ -} + } + + for (j = 1; j <= 3; j++) + { + double angle = (180.0 / M_PI) * Angle (ngs.Get(j), ngs.Get(j%3+1)); + if (angle > 60) + { + (*testout) << "el " << i << " node-node " + << ngs.Get(j) << " - " << ngs.Get(j%3+1) + << " has angle = " << angle << endl; + } + } + } + */ + } } diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 35495ae7..e3f9c689 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -896,7 +896,6 @@ namespace netgen } } - if (strcmp (str, "volumeelements") == 0) { infile >> n; diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp index 5e4112c6..84c216c6 100644 --- a/libsrc/meshing/meshing2.cpp +++ b/libsrc/meshing/meshing2.cpp @@ -196,6 +196,14 @@ namespace netgen MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, double gh, int facenr) { + static int timer = NgProfiler::CreateTimer ("surface meshing"); + + static int timer1 = NgProfiler::CreateTimer ("surface meshing1"); + static int timer2 = NgProfiler::CreateTimer ("surface meshing2"); + static int timer3 = NgProfiler::CreateTimer ("surface meshing3"); + NgProfiler::RegionTimer reg (timer); + + Array pindex, lindex; Array delpoints, dellines; @@ -205,7 +213,6 @@ namespace netgen Array locelements; int z1, z2, oldnp(-1); - SurfaceElementIndex sei; bool found; int rulenr(-1); int globind; @@ -254,7 +261,9 @@ namespace netgen double totalarea = Area (); double meshedarea = 0; + // search tree for surface elements: + /* for (sei = 0; sei < mesh.GetNSE(); sei++) { const Element2d & sel = mesh[sei]; @@ -269,18 +278,45 @@ namespace netgen box.Add ( mesh[sel[2]] ); surfeltree.Insert (box, sei); } - - double trigarea = Cross ( mesh[sel[1]]-mesh[sel[0]], - mesh[sel[2]]-mesh[sel[0]] ).Length() / 2; - - - if (sel.GetNP() == 4) - trigarea += Cross (Vec3d (mesh.Point (sel.PNum(1)), - mesh.Point (sel.PNum(3))), - Vec3d (mesh.Point (sel.PNum(1)), - mesh.Point (sel.PNum(4)))).Length() / 2;; - meshedarea += trigarea; } + */ + Array seia; + mesh.GetSurfaceElementsOfFace (facenr, seia); + for (int i = 0; i < seia.Size(); i++) + { + const Element2d & sel = mesh[seia[i]]; + + if (sel.IsDeleted()) continue; + + Box<3> box; + box.Set ( mesh[sel[0]] ); + box.Add ( mesh[sel[1]] ); + box.Add ( mesh[sel[2]] ); + surfeltree.Insert (box, i); + } + + + + + if (totalarea > 0 || maxarea > 0) + for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) + { + const Element2d & sel = mesh[sei]; + if (sel.IsDeleted()) continue; + + double trigarea = Cross ( mesh[sel[1]]-mesh[sel[0]], + mesh[sel[2]]-mesh[sel[0]] ).Length() / 2; + + + if (sel.GetNP() == 4) + trigarea += Cross (Vec3d (mesh.Point (sel.PNum(1)), + mesh.Point (sel.PNum(3))), + Vec3d (mesh.Point (sel.PNum(1)), + mesh.Point (sel.PNum(4)))).Length() / 2;; + meshedarea += trigarea; + } + + const char * savetask = multithread.task; @@ -293,8 +329,12 @@ namespace netgen double meshedarea_before = meshedarea; + while (!adfront ->Empty() && !multithread.terminate) { + NgProfiler::RegionTimer reg1 (timer1); + + if (multithread.terminate) throw NgException ("Meshing stopped"); @@ -370,6 +410,10 @@ namespace netgen adfront ->GetLocals (baselineindex, locpoints, mpgeominfo, loclines, pindex, lindex, 2*hinner); + + + NgProfiler::RegionTimer reg2 (timer2); + //(*testout) << "h for locals: " << 2*hinner << endl; @@ -707,6 +751,9 @@ namespace netgen } } + NgProfiler::RegionTimer reg3 (timer3); + + for (int i = 1; i <= locelements.Size() && found; i++) { const Element2d & el = locelements.Get(i); @@ -1451,14 +1498,15 @@ namespace netgen PrintMessage (3, "Surface meshing done"); + adfront->PrintOpenSegments (*testout); multithread.task = savetask; - // cout << "surfeltree.depth = " << surfeltree.Tree().Depth() << endl; EndMesh (); + if (!adfront->Empty()) return MESHING2_GIVEUP; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index eb83dbfc..05ac8cac 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -232,7 +232,8 @@ class MeshPoint : public Point<3> #endif public: - MeshPoint () : layer(1), singular(0.), type(INNERPOINT) + MeshPoint () + // : layer(1), singular(0.), type(INNERPOINT) // would unnecessarily initialize large arrays ! { #ifdef PARALLEL isghost = 0; diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp index 1794b840..0d25efbd 100644 --- a/libsrc/meshing/smoothing2.cpp +++ b/libsrc/meshing/smoothing2.cpp @@ -667,14 +667,15 @@ namespace netgen } - static int timer = NgProfiler::CreateTimer ("MeshSmoothing 2D"); - NgProfiler::RegionTimer reg (timer); + static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start"); + NgProfiler::RegionTimer reg (timer); + NgProfiler::StartTimer (timer1); CheckMeshApproximation (mesh); - SurfaceElementIndex sei; + // SurfaceElementIndex sei; Array seia; mesh.GetSurfaceElementsOfFace (faceindex, seia); @@ -694,14 +695,17 @@ namespace netgen PointGeomInfo ngi; Point3d origp; + + Vec3d n1, n2; Vector x(2), xedge(1); Array savepoints(mesh.GetNP()); + uselocalh = mparam.uselocalh; - Array nelementsonpoint(mesh.GetNP()); + Array nelementsonpoint(mesh.GetNP()); nelementsonpoint = 0; for (int i = 0; i < seia.Size(); i++) { @@ -711,6 +715,7 @@ namespace netgen } TABLE elementsonpoint(nelementsonpoint); + for (int i = 0; i < seia.Size(); i++) { const Element2d & el = mesh[seia[i]]; @@ -718,6 +723,7 @@ namespace netgen elementsonpoint.Add (el[j], seia[i]); } + loch = mparam.maxh; locmetricweight = metricweight; meshthis = this; @@ -813,6 +819,10 @@ namespace netgen } int cnt = 0; + + NgProfiler::StopTimer (timer1); + + for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++) if (mesh[pi].Type() == SURFACEPOINT) { @@ -825,9 +835,10 @@ namespace netgen printeddot = 1; PrintDot (plotchar); } - + if (elementsonpoint[pi].Size() == 0) continue; + sp1 = mesh[pi]; @@ -850,7 +861,7 @@ namespace netgen for (int j = 0; j < elementsonpoint[pi].Size(); j++) { - sei = elementsonpoint[pi][j]; + SurfaceElementIndex sei = elementsonpoint[pi][j]; const Element2d & bel = mesh[sei]; surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();