From df97e45bd14421f88cdb7bcc61bf59d41bd5ffbd Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 20 Jul 2020 18:56:36 +0200 Subject: [PATCH] SplitImprove2 - cleanup, new point at min dist of edges --- libsrc/gprim/geomtest3d.cpp | 64 +++++++++++++++++++++++++++++++----- libsrc/gprim/geomtest3d.hpp | 3 ++ libsrc/meshing/improve3.cpp | 27 ++++----------- libsrc/meshing/improve3.hpp | 8 +++++ libsrc/meshing/meshclass.cpp | 2 -- 5 files changed, 73 insertions(+), 31 deletions(-) diff --git a/libsrc/gprim/geomtest3d.cpp b/libsrc/gprim/geomtest3d.cpp index 0c117a10..5c3a827a 100644 --- a/libsrc/gprim/geomtest3d.cpp +++ b/libsrc/gprim/geomtest3d.cpp @@ -1035,6 +1035,39 @@ double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p) } +double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p, double & lam) +{ + Vec3d v(lp1, lp2); + Vec3d vlp(lp1, p); + + // dist(lam) = \| vlp \|^2 - 2 lam (v1p, v) + lam^2 \| v \|^2 + + // lam = (v * vlp) / (v * v); + // if (lam < 0) lam = 0; + // if (lam > 1) lam = 1; + + double num = v*vlp; + double den = v*v; + + if (num <= 0) + { + lam = 0.0; + return Dist2 (lp1, p); + } + + if (num >= den) + { + lam = 1.0; + return Dist2 (lp2, p); + } + + lam = num/den; + if (den > 0) + return vlp.Length2() - num * num /den; + else + return vlp.Length2(); +} + double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, const Point3d & tp3, const Point3d & p) @@ -1102,7 +1135,7 @@ double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, // 0 checks !!! double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, - const Point3d & l2p1, const Point3d & l2p2) + const Point3d & l2p1, const Point3d & l2p2, double & lam1, double & lam2 ) { // dist(lam1,lam2) = \| l2p1+lam2v2 - (l1p1+lam1 v1) \| // min ! @@ -1112,7 +1145,7 @@ double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, Vec3d v2 (l2p1, l2p2); double a11, a12, a22, rs1, rs2; - double lam1, lam2, det; + double det; a11 = v1*v1; a12 = -(v1*v2); @@ -1138,14 +1171,27 @@ double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, } double minv, hv; - minv = MinDistLP2 (l1p1, l1p2, l2p1); - hv = MinDistLP2 (l1p1, l1p2, l2p2); - if (hv < minv) minv = hv; + minv = MinDistLP2 (l1p1, l1p2, l2p1, lam1); + lam2 = 0.; + hv = MinDistLP2 (l1p1, l1p2, l2p2, lam1); + if (hv < minv) + { + lam2 = 1.; + minv = hv; + } - hv = MinDistLP2 (l2p1, l2p2, l1p1); - if (hv < minv) minv = hv; - hv = MinDistLP2 (l2p1, l2p2, l1p2); - if (hv < minv) minv = hv; + hv = MinDistLP2 (l2p1, l2p2, l1p1, lam2); + if (hv < minv) + { + lam1 = 0.; + minv = hv; + } + hv = MinDistLP2 (l2p1, l2p2, l1p2, lam2); + if (hv < minv) + { + lam1 = 1.; + minv = hv; + } return minv; } diff --git a/libsrc/gprim/geomtest3d.hpp b/libsrc/gprim/geomtest3d.hpp index a8bb0266..c1fb1186 100644 --- a/libsrc/gprim/geomtest3d.hpp +++ b/libsrc/gprim/geomtest3d.hpp @@ -91,6 +91,9 @@ extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, const Point3d & l2p1, const Point3d & l2p2); +extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, + const Point3d & l2p1, const Point3d & l2p2, double & lam1, double & lam2 ); + } #endif diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 33af9de3..e08e9b7a 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -3920,12 +3920,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) { if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex()) continue; -// el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0); - const auto & el = mesh[ei]; - if(el.GetNP()==4) - el_badness[ei] = CalcTetBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]], 0, mp); - else - el_badness[ei] = 0.0; + el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0); } }); @@ -3953,6 +3948,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) int minedge = -1; double mindist = 1e99; + double minlam0, minlam1; for (int i : Range(3)) { @@ -3967,11 +3963,14 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) // Point3d p2 = mesh[pi2]; // Point3d p3 = mesh[pi3]; - double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3]); + double lam0, lam1; + double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3], lam0, lam1 ); if(dist has_both_points1; Point3d p[4] = { mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]] }; - auto center = Center(p[0], p[1], p[2], p[3]); + auto center = Center(p[0]+minlam0*(p[1]-p[0]), p[2]+minlam1*(p[3]-p[2])); MeshPoint pnew; pnew(0) = center.X(); @@ -4074,13 +4073,8 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) if(newel1[i] == pi0) newel1[i] = pinew; if(newel2[i] == pi1) newel2[i] = pinew; } - if (!mesh.LegalTet (newel1)) cout << "newel1 illegal tet" << endl; mesh.AddVolumeElement (newel1); - if (!mesh.LegalTet (newel2)) cout << "newel2 illegal tet" << endl; mesh.AddVolumeElement (newel2); -// cout << "del element " << ei1 << ':' << oldel[0] << ',' << oldel[1] << ',' << oldel[2] << ',' << oldel[3] << endl; -// cout << "add element " << newel1[0] << ',' << newel1[1] << ',' << newel1[2] << ',' << newel1[3] << endl; -// cout << "add element " << newel2[0] << ',' << newel2[1] << ',' << newel2[2] << ',' << newel2[3] << endl; } for (auto ei1 : has_both_points1) { @@ -4096,18 +4090,11 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) if(newel1[i] == pi2) newel1[i] = pinew; if(newel2[i] == pi3) newel2[i] = pinew; } -// cout << "del element " << ei1 << ':' << oldel[0] << ',' << oldel[1] << ',' << oldel[2] << ',' << oldel[3] << endl; -// cout << "add element " << newel1[0] << ',' << newel1[1] << ',' << newel1[2] << ',' << newel1[3] << endl; -// cout << "add element " << newel2[0] << ',' << newel2[1] << ',' << newel2[2] << ',' << newel2[3] << endl; - if (!mesh.LegalTet (newel1)) cout << "newel1 illegal tet" << endl; mesh.AddVolumeElement (newel1); - if (!mesh.LegalTet (newel2)) cout << "newel2 illegal tet" << endl; mesh.AddVolumeElement (newel2); } } } - MeshTopology mt(mesh); - mt.Update(); mesh.Compress(); bad1 = mesh.CalcTotalBad (mp); cout << "Total badness after splitimprove2= " << bad1 << endl; diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp index 9da6f548..25de63f3 100644 --- a/libsrc/meshing/improve3.hpp +++ b/libsrc/meshing/improve3.hpp @@ -50,6 +50,14 @@ public: return 0; } + double + CalcBadPow (const Mesh::T_POINTS & points, const Element & elem, double h) + { + if (elem.GetType() == TET) + return pow (max2(CalcTetBadness (points[elem[0]], points[elem[1]], + points[elem[2]], points[elem[3]], h, mp), 1e-10) , 1.0/mp.opterrpow); + return 0; + } double CalcTotalBad (const Mesh::T_POINTS & points, const Array & elements) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 0efbc756..76a22ffc 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -3620,8 +3620,6 @@ namespace netgen // FindOpenElements(); timestamp = NextTimeStamp(); lock.UnLock(); - FindOpenElements(); - CheckVolumeMesh (); } void Mesh :: OrderElements()