ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-15 18:38:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/master' into boundarylayer_fixes

Browse Source
This commit is contained in:
Matthias Hochsteger 2024-02-29 09:45:11 +01:00
commit 4a95414ec8
13 changed files with 1377 additions and 1204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
.gitignore vendored
View File

@ -1,3 +1,11 @@
*.whl *.whl
dist dist
build build
*.vol.gz
*.vol
*.ini
__pycache__
*.json
*.zip
.cache
*.patch

3
libsrc/meshing/basegeom.cpp
View File

@ -1271,9 +1271,6 @@ namespace netgen
if (res != MESHING3_OK) return 1; if (res != MESHING3_OK) return 1;
if (multithread.terminate) return 0; if (multithread.terminate) return 0;
RemoveIllegalElements (*mesh);
if (multithread.terminate) return 0;
MeshQuality3d (*mesh); MeshQuality3d (*mesh);
} }

5
libsrc/meshing/delaunay.cpp
View File

@ -1625,7 +1625,8 @@ namespace netgen
// tempmesh.Save ("tempmesh.vol"); // tempmesh.Save ("tempmesh.vol");
{ {
MeshOptimize3d meshopt(mp); MeshOptimize3d meshopt(tempmesh, mp);
meshopt.SetGoal(OPT_CONFORM);
tempmesh.Compress(); tempmesh.Compress();
tempmesh.FindOpenElements (); tempmesh.FindOpenElements ();
#ifndef EMSCRIPTEN #ifndef EMSCRIPTEN
@ -1638,7 +1639,7 @@ namespace netgen
if(i%5==0) if(i%5==0)
tempmesh.FreeOpenElementsEnvironment (1); tempmesh.FreeOpenElementsEnvironment (1);
meshopt.SwapImprove(tempmesh, OPT_CONFORM); meshopt.SwapImprove();
} }
tempmesh.Compress(); tempmesh.Compress();
} }

323
libsrc/meshing/improve3.cpp
View File

@ -41,7 +41,7 @@ static ArrayMem<Element, 3> SplitElement (Element old, PointIndex pi0, PointInde
ArrayMem<Element, 3> new_elements; ArrayMem<Element, 3> new_elements;
// split element by cutting edge pi0,pi1 at pinew // split element by cutting edge pi0,pi1 at pinew
auto np = old.GetNP(); auto np = old.GetNP();
old.Flags().illegal_valid = 0; old.Touch();
if(np == 4) if(np == 4)
{ {
// Split tet into two tets // Split tet into two tets
@ -89,7 +89,7 @@ static ArrayMem<Element, 3> SplitElement (Element old, PointIndex pi0, PointInde
} }
return new_elements; return new_elements;
}; }
static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingParameters & mp, Element old, PointIndex pi0, PointIndex pi1, MeshPoint & pnew) static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingParameters & mp, Element old, PointIndex pi0, PointIndex pi1, MeshPoint & pnew)
{ {
@ -133,7 +133,63 @@ static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingP
} }
return badness; return badness;
}; }
tuple<double, double, int> MeshOptimize3d :: UpdateBadness()
{
static Timer tbad("UpdateBadness");
RegionTimer reg(tbad);
double totalbad = 0.0;
double maxbad = 0.0;
atomic<int> bad_elements = 0;
ParallelForRange(Range(mesh.GetNE()), [&] (auto myrange) {
double totalbad_local = 0.0;
double maxbad_local = 0.0;
int bad_elements_local = 0;
for (ElementIndex ei : myrange)
{
auto & el = mesh[ei];
if(mp.only3D_domain_nr && mp.only3D_domain_nr != el.GetIndex()) continue;
if(!el.BadnessValid())
el.SetBadness(CalcBad(mesh.Points(), el, 0));
double bad = el.GetBadness();
totalbad_local += bad;
maxbad_local = max(maxbad_local, bad);
if(bad > min_badness)
bad_elements_local++;
}
AtomicAdd(totalbad, totalbad_local);
AtomicMax(maxbad, maxbad_local);
bad_elements += bad_elements_local;
});
return {totalbad, maxbad, bad_elements};
}
bool MeshOptimize3d :: HasBadElement(FlatArray<ElementIndex> els)
{
for(auto ei : els)
if(mesh[ei].GetBadness()>min_badness)
return true;
return false;
}
bool MeshOptimize3d :: HasIllegalElement(FlatArray<ElementIndex> els)
{
for(auto ei : els)
if(!mesh.LegalTet(mesh[ei]))
return true;
return false;
}
bool MeshOptimize3d :: NeedsOptimization(FlatArray<ElementIndex> els)
{
if(goal == OPT_LEGAL) return HasIllegalElement(els);
if(goal == OPT_QUALITY) return HasBadElement(els);
return true;
}
/* /*
@ -143,10 +199,8 @@ static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingP
Connect inner point to boundary point, if one Connect inner point to boundary point, if one
point is inner point. point is inner point.
*/ */
double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh, double MeshOptimize3d :: CombineImproveEdge (
const MeshingParameters & mp,
Table<ElementIndex, PointIndex> & elements_of_point, Table<ElementIndex, PointIndex> & elements_of_point,
Array<double> & elerrs,
PointIndex pi0, PointIndex pi1, PointIndex pi0, PointIndex pi1,
FlatArray<bool, PointIndex> is_point_removed, FlatArray<bool, PointIndex> is_point_removed,
bool check_only) bool check_only)
@ -199,9 +253,9 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
double badness_old = 0.0; double badness_old = 0.0;
for (auto ei : has_one_point) for (auto ei : has_one_point)
badness_old += elerrs[ei]; badness_old += mesh[ei].GetBadness();
for (auto ei : has_both_points) for (auto ei : has_both_points)
badness_old += elerrs[ei]; badness_old += mesh[ei].GetBadness();
MeshPoint pnew = p0; MeshPoint pnew = p0;
if (p0.Type() == INNERPOINT) if (p0.Type() == INNERPOINT)
@ -232,7 +286,7 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
break; break;
} }
elem.Flags().illegal_valid = 0; elem.Touch();
if (!mesh.LegalTet(elem)) if (!mesh.LegalTet(elem))
badness_new += 1e4; badness_new += 1e4;
} }
@ -255,17 +309,17 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
if (elem[l] == pi1) if (elem[l] == pi1)
elem[l] = pi0; elem[l] = pi0;
elem.Flags().illegal_valid = 0; elem.Touch();
if (!mesh.LegalTet (elem)) if (!mesh.LegalTet (elem))
(*testout) << "illegal tet " << ei << endl; (*testout) << "illegal tet " << ei << endl;
} }
for (auto i : Range(has_one_point)) for (auto i : Range(has_one_point))
elerrs[has_one_point[i]] = one_point_badness[i]; mesh[has_one_point[i]].SetBadness(one_point_badness[i]);
for (auto ei : has_both_points) for (auto ei : has_both_points)
{ {
mesh[ei].Flags().illegal_valid = 0; mesh[ei].Touch();
mesh[ei].Delete(); mesh[ei].Delete();
} }
} }
@ -273,14 +327,12 @@ double MeshOptimize3d :: CombineImproveEdge (Mesh & mesh,
return d_badness; return d_badness;
} }
void MeshOptimize3d :: CombineImprove (Mesh & mesh, void MeshOptimize3d :: CombineImprove ()
OPTIMIZEGOAL goal)
{ {
static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t); static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t);
static Timer topt("Optimize"); static Timer topt("Optimize");
static Timer tsearch("Search"); static Timer tsearch("Search");
static Timer tbuild_elements_table("Build elements table"); static Timer tbuild_elements_table("Build elements table");
static Timer tbad("CalcBad");
mesh.BuildBoundaryEdges(false); mesh.BuildBoundaryEdges(false);
@ -288,7 +340,6 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
int ne = mesh.GetNE(); int ne = mesh.GetNE();
int ntasks = 4*ngcore::TaskManager::GetNumThreads(); int ntasks = 4*ngcore::TaskManager::GetNumThreads();
Array<double> elerrs (ne);
Array<bool, PointIndex> is_point_removed (np); Array<bool, PointIndex> is_point_removed (np);
is_point_removed = false; is_point_removed = false;
@ -300,26 +351,11 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
multithread.task = "Optimize Volume: Combine Improve"; multithread.task = "Optimize Volume: Combine Improve";
tbad.Start(); UpdateBadness();
double totalbad = 0.0;
ParallelForRange(Range(ne), [&] (auto myrange)
{
double totalbad_local = 0.0;
for (ElementIndex ei : myrange)
{
if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
continue;
double elerr = CalcBad (mesh.Points(), mesh[ei], 0);
totalbad_local += elerr;
elerrs[ei] = elerr;
}
AtomicAdd(totalbad, totalbad_local);
}, ntasks);
tbad.Stop();
if (goal == OPT_QUALITY) if (goal == OPT_QUALITY && testout->good())
{ {
totalbad = mesh.CalcTotalBad (mp); double totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl; (*testout) << "Total badness = " << totalbad << endl;
} }
@ -338,7 +374,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
for(auto i : myrange) for(auto i : myrange)
{ {
auto [p0,p1] = edges[i]; auto [p0,p1] = edges[i];
double d_badness = CombineImproveEdge (mesh, mp, elementsonnode, elerrs, p0, p1, is_point_removed, true); double d_badness = CombineImproveEdge (elementsonnode, p0, p1, is_point_removed, true);
if(d_badness<0.0) if(d_badness<0.0)
{ {
int index = improvement_counter++; int index = improvement_counter++;
@ -360,7 +396,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
for(auto [d_badness, ei] : edges_with_improvement) for(auto [d_badness, ei] : edges_with_improvement)
{ {
auto [p0,p1] = edges[ei]; auto [p0,p1] = edges[ei];
if (CombineImproveEdge (mesh, mp, elementsonnode, elerrs, p0, p1, is_point_removed, false) < 0.0) if (CombineImproveEdge (elementsonnode, p0, p1, is_point_removed, false) < 0.0)
cnt++; cnt++;
} }
topt.Stop(); topt.Stop();
@ -371,9 +407,9 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
PrintMessage (5, cnt, " elements combined"); PrintMessage (5, cnt, " elements combined");
(*testout) << "CombineImprove done" << "\n"; (*testout) << "CombineImprove done" << "\n";
if (goal == OPT_QUALITY) if (goal == OPT_QUALITY && testout->good())
{ {
totalbad = mesh.CalcTotalBad (mp); double totalbad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << totalbad << endl; (*testout) << "Total badness = " << totalbad << endl;
int cntill = 0; int cntill = 0;
@ -389,7 +425,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only) double MeshOptimize3d :: SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only)
{ {
double d_badness = 0.0; double d_badness = 0.0;
// int cnt = 0; // int cnt = 0;
@ -418,21 +454,14 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
if(mp.only3D_domain_nr != mesh[ei].GetIndex()) if(mp.only3D_domain_nr != mesh[ei].GetIndex())
return 0.0; return 0.0;
if (goal == OPT_LEGAL) if (!NeedsOptimization(hasbothpoints))
{
bool all_tets_legal = true;
for(auto ei : hasbothpoints)
if( !mesh.LegalTet (mesh[ei]) || elerrs[ei] > 1e3)
all_tets_legal = false;
if(all_tets_legal)
return 0.0; return 0.0;
}
double bad1 = 0.0; double bad1 = 0.0;
double bad1_max = 0.0; double bad1_max = 0.0;
for (ElementIndex ei : hasbothpoints) for (ElementIndex ei : hasbothpoints)
{ {
double bad = elerrs[ei]; double bad = mesh[ei].GetBadness();
bad1 += bad; bad1 += bad;
bad1_max = max(bad1_max, bad); bad1_max = max(bad1_max, bad);
} }
@ -505,9 +534,8 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
Element newel1 = oldel; Element newel1 = oldel;
Element newel2 = oldel; Element newel2 = oldel;
oldel.Flags().illegal_valid = 0; newel1.Touch();
newel1.Flags().illegal_valid = 0; newel2.Touch();
newel2.Flags().illegal_valid = 0;
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
{ {
@ -536,11 +564,11 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
Element newel1 = oldel; Element newel1 = oldel;
Element newel2 = oldel; Element newel2 = oldel;
oldel.Flags().illegal_valid = 0; oldel.Touch();
oldel.Delete(); oldel.Delete();
newel1.Flags().illegal_valid = 0; newel1.Touch();
newel2.Flags().illegal_valid = 0; newel2.Touch();
for (int l = 0; l < 4; l++) for (int l = 0; l < 4; l++)
{ {
@ -555,8 +583,7 @@ double MeshOptimize3d :: SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table
return d_badness; return d_badness;
} }
void MeshOptimize3d :: SplitImprove (Mesh & mesh, void MeshOptimize3d :: SplitImprove ()
OPTIMIZEGOAL goal)
{ {
static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t); static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t);
static Timer topt("Optimize"); static Timer topt("Optimize");
@ -569,8 +596,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
auto elementsonnode = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr); auto elementsonnode = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr);
Array<double> elerrs(ne);
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Optimize Volume: Split Improve"; multithread.task = "Optimize Volume: Split Improve";
@ -578,17 +603,9 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
(*testout) << "start SplitImprove" << "\n"; (*testout) << "start SplitImprove" << "\n";
mesh.BuildBoundaryEdges(false); mesh.BuildBoundaryEdges(false);
ParallelFor( mesh.VolumeElements().Range(), [&] (ElementIndex ei) NETGEN_LAMBDA_INLINE UpdateBadness();
{
if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
return;
elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0); if (goal == OPT_QUALITY && testout->good())
bad += elerrs[ei];
AtomicMax(badmax, elerrs[ei]);
});
if (goal == OPT_QUALITY)
{ {
bad = mesh.CalcTotalBad (mp); bad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad << endl; (*testout) << "Total badness = " << bad << endl;
@ -610,7 +627,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
for(auto i : myrange) for(auto i : myrange)
{ {
auto [p0,p1] = edges[i]; auto [p0,p1] = edges[i];
double d_badness = SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, true); double d_badness = SplitImproveEdge (elementsonnode, locfaces, badmax, p0, p1, ptmp, true);
if(d_badness<0.0) if(d_badness<0.0)
{ {
int index = improvement_counter++; int index = improvement_counter++;
@ -633,7 +650,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
for(auto [d_badness, ei] : edges_with_improvement) for(auto [d_badness, ei] : edges_with_improvement)
{ {
auto [p0,p1] = edges[ei]; auto [p0,p1] = edges[ei];
if (SplitImproveEdge (mesh, goal, elementsonnode, elerrs, locfaces, badmax, p0, p1, ptmp, false) < 0.0) if (SplitImproveEdge (elementsonnode, locfaces, badmax, p0, p1, ptmp, false) < 0.0)
cnt++; cnt++;
} }
topt.Stop(); topt.Stop();
@ -642,9 +659,12 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
(*testout) << "Splitt - Improve done" << "\n"; (*testout) << "Splitt - Improve done" << "\n";
if (goal == OPT_QUALITY) if (goal == OPT_QUALITY)
{
if(testout->good())
{ {
bad = mesh.CalcTotalBad (mp); bad = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad << endl; (*testout) << "Total badness = " << bad << endl;
}
[[maybe_unused]] int cntill = 0; [[maybe_unused]] int cntill = 0;
ne = mesh.GetNE(); ne = mesh.GetNE();
@ -658,7 +678,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
} }
double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, double MeshOptimize3d :: SwapImproveEdge (
const NgBitArray * working_elements, const NgBitArray * working_elements,
Table<ElementIndex, PointIndex> & elementsonnode, Table<ElementIndex, PointIndex> & elementsonnode,
INDEX_3_HASHTABLE<int> & faces, INDEX_3_HASHTABLE<int> & faces,
@ -706,8 +726,6 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
} }
} }
bool have_bad_element = false;
for (ElementIndex ei : hasbothpoints) for (ElementIndex ei : hasbothpoints)
{ {
if (mesh[ei].GetType () != TET) if (mesh[ei].GetType () != TET)
@ -727,14 +745,9 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
if (mesh[ei].IsDeleted()) if (mesh[ei].IsDeleted())
return 0.0; return 0.0;
if ((goal == OPT_LEGAL) &&
mesh.LegalTet (mesh[ei]) &&
CalcBad (mesh.Points(), mesh[ei], 0) >= 1e3)
have_bad_element = true;
} }
if ((goal == OPT_LEGAL) && !have_bad_element) if(!NeedsOptimization(hasbothpoints))
return 0.0; return 0.0;
int nsuround = hasbothpoints.Size(); int nsuround = hasbothpoints.Size();
@ -799,9 +812,9 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el32, 0) + CalcBad (mesh.Points(), el32, 0) +
CalcBad (mesh.Points(), el33, 0); CalcBad (mesh.Points(), el33, 0);
el31.Flags().illegal_valid = 0; el31.Touch();
el32.Flags().illegal_valid = 0; el32.Touch();
el33.Flags().illegal_valid = 0; el33.Touch();
if (!mesh.LegalTet(el31) || if (!mesh.LegalTet(el31) ||
!mesh.LegalTet(el32) || !mesh.LegalTet(el32) ||
@ -823,15 +836,15 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
bad2 = CalcBad (mesh.Points(), el21, 0) + bad2 = CalcBad (mesh.Points(), el21, 0) +
CalcBad (mesh.Points(), el22, 0); CalcBad (mesh.Points(), el22, 0);
el21.Flags().illegal_valid = 0; el21.Touch();
el22.Flags().illegal_valid = 0; el22.Touch();
if (!mesh.LegalTet(el21) || if (!mesh.LegalTet(el21) ||
!mesh.LegalTet(el22)) !mesh.LegalTet(el22))
bad2 += 1e4; bad2 += 1e4;
if (goal == OPT_CONFORM && NotTooBad(bad1, bad2)) if ((goal == OPT_CONFORM) && NotTooBad(bad1, bad2))
{ {
INDEX_3 face(pi3, pi4, pi5); INDEX_3 face(pi3, pi4, pi5);
face.Sort(); face.Sort();
@ -866,8 +879,8 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
mesh[hasbothpoints[1]].Delete(); mesh[hasbothpoints[1]].Delete();
mesh[hasbothpoints[2]].Delete(); mesh[hasbothpoints[2]].Delete();
el21.Flags().illegal_valid = 0; el21.Touch();
el22.Flags().illegal_valid = 0; el22.Touch();
mesh.AddVolumeElement(el21); mesh.AddVolumeElement(el21);
mesh.AddVolumeElement(el22); mesh.AddVolumeElement(el22);
} }
@ -942,10 +955,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el4, 0); CalcBad (mesh.Points(), el4, 0);
el1.Flags().illegal_valid = 0; el1.Touch();
el2.Flags().illegal_valid = 0; el2.Touch();
el3.Flags().illegal_valid = 0; el3.Touch();
el4.Flags().illegal_valid = 0; el4.Touch();
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
@ -978,10 +991,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el3, 0) + CalcBad (mesh.Points(), el3, 0) +
CalcBad (mesh.Points(), el4, 0); CalcBad (mesh.Points(), el4, 0);
el1.Flags().illegal_valid = 0; el1.Touch();
el2.Flags().illegal_valid = 0; el2.Touch();
el3.Flags().illegal_valid = 0; el3.Touch();
el4.Flags().illegal_valid = 0; el4.Touch();
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
{ {
@ -1014,10 +1027,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
CalcBad (mesh.Points(), el3b, 0) + CalcBad (mesh.Points(), el3b, 0) +
CalcBad (mesh.Points(), el4b, 0); CalcBad (mesh.Points(), el4b, 0);
el1b.Flags().illegal_valid = 0; el1b.Touch();
el2b.Flags().illegal_valid = 0; el2b.Touch();
el3b.Flags().illegal_valid = 0; el3b.Touch();
el4b.Flags().illegal_valid = 0; el4b.Touch();
if (goal != OPT_CONFORM) if (goal != OPT_CONFORM)
{ {
@ -1054,10 +1067,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
for (auto i : IntRange(4)) for (auto i : IntRange(4))
mesh[hasbothpoints[i]].Delete(); mesh[hasbothpoints[i]].Delete();
el1.Flags().illegal_valid = 0; el1.Touch();
el2.Flags().illegal_valid = 0; el2.Touch();
el3.Flags().illegal_valid = 0; el3.Touch();
el4.Flags().illegal_valid = 0; el4.Touch();
mesh.AddVolumeElement (el1); mesh.AddVolumeElement (el1);
mesh.AddVolumeElement (el2); mesh.AddVolumeElement (el2);
mesh.AddVolumeElement (el3); mesh.AddVolumeElement (el3);
@ -1068,10 +1081,10 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
for (auto i : IntRange(4)) for (auto i : IntRange(4))
mesh[hasbothpoints[i]].Delete(); mesh[hasbothpoints[i]].Delete();
el1b.Flags().illegal_valid = 0; el1b.Touch();
el2b.Flags().illegal_valid = 0; el2b.Touch();
el3b.Flags().illegal_valid = 0; el3b.Touch();
el4b.Flags().illegal_valid = 0; el4b.Touch();
mesh.AddVolumeElement (el1b); mesh.AddVolumeElement (el1b);
mesh.AddVolumeElement (el2b); mesh.AddVolumeElement (el2b);
mesh.AddVolumeElement (el3b); mesh.AddVolumeElement (el3b);
@ -1174,7 +1187,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
hel[3] = pi2; hel[3] = pi2;
bad2 += CalcBad (mesh.Points(), hel, 0); bad2 += CalcBad (mesh.Points(), hel, 0);
hel.Flags().illegal_valid = 0; hel.Touch();
if (!mesh.LegalTet(hel)) bad2 += 1e4; if (!mesh.LegalTet(hel)) bad2 += 1e4;
hel[2] = suroundpts[k % nsuround]; hel[2] = suroundpts[k % nsuround];
@ -1183,7 +1196,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
bad2 += CalcBad (mesh.Points(), hel, 0); bad2 += CalcBad (mesh.Points(), hel, 0);
hel.Flags().illegal_valid = 0; hel.Touch();
if (!mesh.LegalTet(hel)) bad2 += 1e4; if (!mesh.LegalTet(hel)) bad2 += 1e4;
} }
// (*testout) << "bad2," << l << " = " << bad2 << endl; // (*testout) << "bad2," << l << " = " << bad2 << endl;
@ -1253,7 +1266,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
hel[1] = suroundpts[k % nsuround]; hel[1] = suroundpts[k % nsuround];
hel[2] = suroundpts[(k+1) % nsuround]; hel[2] = suroundpts[(k+1) % nsuround];
hel[3] = pi2; hel[3] = pi2;
hel.Flags().illegal_valid = 0; hel.Touch();
/* /*
(*testout) << nsuround << "-swap, new el,top = " (*testout) << nsuround << "-swap, new el,top = "
@ -1289,8 +1302,7 @@ double MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
return d_badness; return d_badness;
} }
void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, void MeshOptimize3d :: SwapImprove (const NgBitArray * working_elements)
const NgBitArray * working_elements)
{ {
static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t); static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
static Timer tloop("MeshOptimize3d::SwapImprove loop"); static Timer tloop("MeshOptimize3d::SwapImprove loop");
@ -1344,7 +1356,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
} }
// Calculate total badness // Calculate total badness
if (goal == OPT_QUALITY) if (goal == OPT_QUALITY && testout->good())
{ {
double bad1 = mesh.CalcTotalBad (mp); double bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl; (*testout) << "Total badness = " << bad1 << endl;
@ -1356,6 +1368,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
Array<std::tuple<double, int>> candidate_edges(edges.Size()); Array<std::tuple<double, int>> candidate_edges(edges.Size());
std::atomic<int> improvement_counter(0); std::atomic<int> improvement_counter(0);
UpdateBadness();
tloop.Start(); tloop.Start();
auto num_elements_before = mesh.VolumeElements().Range().Next(); auto num_elements_before = mesh.VolumeElements().Range().Next();
@ -1368,7 +1382,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
break; break;
auto [pi0, pi1] = edges[i]; auto [pi0, pi1] = edges[i];
double d_badness = SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true); double d_badness = SwapImproveEdge (working_elements, elementsonnode, faces, pi0, pi1, true);
if(d_badness<0.0) if(d_badness<0.0)
{ {
int index = improvement_counter++; int index = improvement_counter++;
@ -1383,7 +1397,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
for(auto [d_badness, ei] : edges_with_improvement) for(auto [d_badness, ei] : edges_with_improvement)
{ {
auto [pi0,pi1] = edges[ei]; auto [pi0,pi1] = edges[ei];
if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false) < 0.0) if(SwapImproveEdge (working_elements, elementsonnode, faces, pi0, pi1, false) < 0.0)
cnt++; cnt++;
} }
@ -1425,7 +1439,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal, void MeshOptimize3d :: SwapImproveSurface (
const NgBitArray * working_elements, const NgBitArray * working_elements,
const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps) const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps)
{ {
@ -1517,9 +1531,7 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
if (mesh[ei].IsDeleted()) if (mesh[ei].IsDeleted())
continue; continue;
if ((goal == OPT_LEGAL) && if (goal == OPT_LEGAL && mesh.LegalTet (mesh[ei]))
mesh.LegalTet (mesh[ei]) &&
CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
continue; continue;
const Element & elemi = mesh[ei]; const Element & elemi = mesh[ei];
@ -2287,7 +2299,7 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
2 -> 3 conversion 2 -> 3 conversion
*/ */
double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face,
Table<ElementIndex, PointIndex> & elementsonnode, Table<ElementIndex, PointIndex> & elementsonnode,
TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only ) TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only )
{ {
@ -2370,6 +2382,10 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
if (elem2.GetType() != TET) if (elem2.GetType() != TET)
continue; continue;
ArrayMem<ElementIndex, 2> elis = {eli1, eli2};
if(!NeedsOptimization(elis))
continue;
int comnodes=0; int comnodes=0;
for (int l = 1; l <= 4; l++) for (int l = 1; l <= 4; l++)
if (elem2.PNum(l) == pi1 || elem2.PNum(l) == pi2 || if (elem2.PNum(l) == pi1 || elem2.PNum(l) == pi2 ||
@ -2384,8 +2400,7 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
if (comnodes == 3) if (comnodes == 3)
{ {
bad1 = CalcBad (mesh.Points(), elem, 0) + bad1 = elem.GetBadness() + elem2.GetBadness();
CalcBad (mesh.Points(), elem2, 0);
if (!mesh.LegalTet(elem) || if (!mesh.LegalTet(elem) ||
!mesh.LegalTet(elem2)) !mesh.LegalTet(elem2))
@ -2415,9 +2430,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
CalcBad (mesh.Points(), el33, 0); CalcBad (mesh.Points(), el33, 0);
el31.Flags().illegal_valid = 0; el31.Touch();
el32.Flags().illegal_valid = 0; el32.Touch();
el33.Flags().illegal_valid = 0; el33.Touch();
if (!mesh.LegalTet(el31) || if (!mesh.LegalTet(el31) ||
!mesh.LegalTet(el32) || !mesh.LegalTet(el32) ||
@ -2436,9 +2451,9 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
if (d_badness<0.0) if (d_badness<0.0)
{ {
el31.Flags().illegal_valid = 0; el31.Touch();
el32.Flags().illegal_valid = 0; el32.Touch();
el33.Flags().illegal_valid = 0; el33.Touch();
mesh[eli1].Delete(); mesh[eli1].Delete();
mesh[eli2].Delete(); mesh[eli2].Delete();
@ -2457,7 +2472,7 @@ double MeshOptimize3d :: SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementI
2 -> 3 conversion 2 -> 3 conversion
*/ */
void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) void MeshOptimize3d :: SwapImprove2 ()
{ {
static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t); static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t);
@ -2478,8 +2493,11 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
PrintMessage (3, "SwapImprove2 "); PrintMessage (3, "SwapImprove2 ");
(*testout) << "\n" << "Start SwapImprove2" << "\n"; (*testout) << "\n" << "Start SwapImprove2" << "\n";
if(testout->good())
{
double bad1 = mesh.CalcTotalBad (mp); double bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl; (*testout) << "Total badness = " << bad1 << endl;
}
// find elements on node // find elements on node
@ -2495,6 +2513,8 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
Array<std::tuple<double, ElementIndex, int>> faces_with_improvement; Array<std::tuple<double, ElementIndex, int>> faces_with_improvement;
Array<Array<std::tuple<double, ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads); Array<Array<std::tuple<double, ElementIndex, int>>> faces_with_improvement_threadlocal(num_threads);
UpdateBadness();
ParallelForRange( Range(ne), [&]( auto myrange ) ParallelForRange( Range(ne), [&]( auto myrange )
{ {
int tid = ngcore::TaskManager::GetThreadId(); int tid = ngcore::TaskManager::GetThreadId();
@ -2510,9 +2530,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
if (mesh[eli1].GetType() != TET) if (mesh[eli1].GetType() != TET)
continue; continue;
if ((goal == OPT_LEGAL) && if (goal == OPT_LEGAL && mesh.LegalTet (mesh[eli1]))
mesh.LegalTet (mesh[eli1]) &&
CalcBad (mesh.Points(), mesh[eli1], 0) < 1e3)
continue; continue;
if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(eli1).GetIndex()) if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(eli1).GetIndex())
@ -2520,7 +2538,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
{ {
double d_badness = SwapImprove2( mesh, goal, eli1, j, elementsonnode, belementsonnode, true); double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, true);
if(d_badness<0.0) if(d_badness<0.0)
my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) ); my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) );
} }
@ -2536,22 +2554,24 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
{ {
if(mesh[eli].IsDeleted()) if(mesh[eli].IsDeleted())
continue; continue;
if(SwapImprove2( mesh, goal, eli, j, elementsonnode, belementsonnode, false) < 0.0) if(SwapImprove2( eli, j, elementsonnode, belementsonnode, false) < 0.0)
cnt++; cnt++;
} }
PrintMessage (5, cnt, " swaps performed"); PrintMessage (5, cnt, " swaps performed");
mesh.Compress(); mesh.Compress();
bad1 = mesh.CalcTotalBad (mp); if(testout->good())
{
double bad1 = mesh.CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl; (*testout) << "Total badness = " << bad1 << endl;
(*testout) << "swapimprove2 done" << "\n"; (*testout) << "swapimprove2 done" << "\n";
} }
}
double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh, double MeshOptimize3d :: SplitImprove2Element (
ElementIndex ei, ElementIndex ei,
const Table<ElementIndex, PointIndex> & elements_of_point, const Table<ElementIndex, PointIndex> & elements_of_point,
const Array<double> & el_badness,
bool check_only) bool check_only)
{ {
auto & el = mesh[ei]; auto & el = mesh[ei];
@ -2559,7 +2579,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
return false; return false;
// Optimize only bad elements // Optimize only bad elements
if(el_badness[ei] < 100) if(el.GetBadness() < 100)
return false; return false;
// search for very flat tets, with two disjoint edges nearly crossing, like a rectangle with diagonals // search for very flat tets, with two disjoint edges nearly crossing, like a rectangle with diagonals
@ -2635,21 +2655,21 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
has_both_points1.Append (ei1); has_both_points1.Append (ei1);
} }
double badness_before = el_badness[ei]; double badness_before = mesh[ei].GetBadness();
double badness_after = 0.0; double badness_after = 0.0;
for (auto ei0 : has_both_points0) for (auto ei0 : has_both_points0)
{ {
if(mesh[ei0].GetType()!=TET) if(mesh[ei0].GetType()!=TET)
return false; return false;
badness_before += el_badness[ei0]; badness_before += mesh[ei0].GetBadness();
badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei0], pi0, pi1, pnew); badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei0], pi0, pi1, pnew);
} }
for (auto ei1 : has_both_points1) for (auto ei1 : has_both_points1)
{ {
if(mesh[ei1].GetType()!=TET) if(mesh[ei1].GetType()!=TET)
return false; return false;
badness_before += el_badness[ei1]; badness_before += mesh[ei1].GetBadness();
badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei1], pi2, pi3, pnew); badness_after += SplitElementBadness (mesh.Points(), mp, mesh[ei1], pi2, pi3, pnew);
} }
@ -2659,7 +2679,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
if(badness_after<badness_before) if(badness_after<badness_before)
{ {
PointIndex pinew = mesh.AddPoint (center); PointIndex pinew = mesh.AddPoint (center);
el.Flags().illegal_valid = 0; el.Touch();
el.Delete(); el.Delete();
for (auto ei1 : has_both_points0) for (auto ei1 : has_both_points0)
@ -2682,7 +2702,7 @@ double MeshOptimize3d :: SplitImprove2Element (Mesh & mesh,
// Split two opposite edges of very flat tet and let all 4 new segments have one common vertex // Split two opposite edges of very flat tet and let all 4 new segments have one common vertex
// Imagine a square with 2 diagonals -> new point where diagonals cross, remove the flat tet // Imagine a square with 2 diagonals -> new point where diagonals cross, remove the flat tet
void MeshOptimize3d :: SplitImprove2 (Mesh & mesh) void MeshOptimize3d :: SplitImprove2 ()
{ {
static Timer t("MeshOptimize3d::SplitImprove2"); RegionTimer reg(t); static Timer t("MeshOptimize3d::SplitImprove2"); RegionTimer reg(t);
static Timer tsearch("Search"); static Timer tsearch("Search");
@ -2695,18 +2715,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh)
const char * savetask = multithread.task; const char * savetask = multithread.task;
multithread.task = "Optimize Volume: Split Improve 2"; multithread.task = "Optimize Volume: Split Improve 2";
Array<double> el_badness (ne); UpdateBadness();
ParallelForRange(Range(ne), [&] (auto myrange)
{
for (ElementIndex ei : myrange)
{
if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
continue;
el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
}
});
mesh.BuildBoundaryEdges(false); mesh.BuildBoundaryEdges(false);
Array<std::tuple<double, ElementIndex>> split_candidates(ne); Array<std::tuple<double, ElementIndex>> split_candidates(ne);
@ -2719,7 +2728,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh)
{ {
if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex()) if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
continue; continue;
double d_badness = SplitImprove2Element(mesh, ei, elements_of_point, el_badness, true); double d_badness = SplitImprove2Element(ei, elements_of_point, true);
if(d_badness<0.0) if(d_badness<0.0)
{ {
int index = improvement_counter++; int index = improvement_counter++;
@ -2736,7 +2745,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh)
topt.Start(); topt.Start();
for(auto [d_badness, ei] : elements_with_improvement) for(auto [d_badness, ei] : elements_with_improvement)
{ {
if( SplitImprove2Element(mesh, ei, elements_of_point, el_badness, false) < 0.0) if( SplitImprove2Element(ei, elements_of_point, false) < 0.0)
cnt++; cnt++;
} }
topt.Stop(); topt.Stop();

44
libsrc/meshing/improve3.hpp
View File

@ -13,32 +13,46 @@ extern double CalcTotalBad (const Mesh::T_POINTS & points,
class MeshOptimize3d class MeshOptimize3d
{ {
const MeshingParameters & mp; const MeshingParameters & mp;
Mesh & mesh;
OPTIMIZEGOAL goal = OPT_QUALITY;
double min_badness = 0;
bool HasBadElement(FlatArray<ElementIndex> els);
bool HasIllegalElement(FlatArray<ElementIndex> els);
bool NeedsOptimization(FlatArray<ElementIndex> els);
public: public:
MeshOptimize3d (const MeshingParameters & amp) : mp(amp) { ; }
double CombineImproveEdge (Mesh & mesh, const MeshingParameters & mp, MeshOptimize3d (Mesh & m, const MeshingParameters & amp, OPTIMIZEGOAL agoal = OPT_QUALITY) :
mesh(m), mp(amp), goal(agoal) { ; }
void SetGoal(OPTIMIZEGOAL agoal) { goal = agoal; }
void SetMinBadness(double badness) { min_badness = badness; }
tuple<double, double, int> UpdateBadness();
double CombineImproveEdge (
Table<ElementIndex, PointIndex> & elements_of_point, Table<ElementIndex, PointIndex> & elements_of_point,
Array<double> & elerrs, PointIndex pi0, PointIndex pi1, PointIndex pi0, PointIndex pi1,
FlatArray<bool, PointIndex> is_point_removed, bool check_only=false); FlatArray<bool, PointIndex> is_point_removed, bool check_only=false);
void CombineImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY); void CombineImprove ();
void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY); void SplitImprove ();
double SplitImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, Table<ElementIndex,PointIndex> & elementsonnode, Array<double> &elerrs, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false); double SplitImproveEdge (Table<ElementIndex,PointIndex> & elementsonnode, NgArray<INDEX_3> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false);
void SplitImprove2 (Mesh & mesh); void SplitImprove2 ();
double SplitImprove2Element (Mesh & mesh, ElementIndex ei, const Table<ElementIndex, PointIndex> & elements_of_point, const Array<double> & elerrs, bool check_only); double SplitImprove2Element (ElementIndex ei, const Table<ElementIndex, PointIndex> & elements_of_point, bool check_only);
double SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false); double SwapImproveEdge (const NgBitArray * working_elements, Table<ElementIndex,PointIndex> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY, void SwapImprove (const NgBitArray * working_elements = NULL);
const NgBitArray * working_elements = NULL); void SwapImproveSurface (const NgBitArray * working_elements = NULL,
void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
const NgBitArray * working_elements = NULL,
const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL); const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY); void SwapImprove2 ();
double SwapImprove2 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false ); double SwapImprove2 (ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );
void ImproveMesh() { mesh.ImproveMesh(mp, goal); }
double double
CalcBad (const Mesh::T_POINTS & points, const Element & elem, double h) CalcBad (const Mesh::T_POINTS & points, const Element & elem, double h)

4
libsrc/meshing/meshclass.cpp
View File

@ -603,7 +603,7 @@ namespace netgen
{ {
volelements.Append (el); volelements.Append (el);
} }
volelements.Last().Flags().illegal_valid = 0; volelements.Last().Touch();
volelements.Last().Flags().fixed = 0; volelements.Last().Flags().fixed = 0;
volelements.Last().Flags().deleted = 0; volelements.Last().Flags().deleted = 0;
@ -626,7 +626,7 @@ namespace netgen
*/ */
volelements[ei] = el; volelements[ei] = el;
volelements[ei].Flags().illegal_valid = 0; volelements[ei].Touch();
volelements[ei].Flags().fixed = 0; volelements[ei].Flags().fixed = 0;
volelements[ei].Flags().deleted = 0; volelements[ei].Flags().deleted = 0;
} }

54
libsrc/meshing/meshfunc.cpp
View File

@ -480,7 +480,7 @@ namespace netgen
meshed = 0; meshed = 0;
PrintMessage (5, mesh.GetNOpenElements(), " open faces found"); PrintMessage (5, mesh.GetNOpenElements(), " open faces found");
MeshOptimize3d optmesh(mp); MeshOptimize3d optmesh(mesh, mp, OPT_REST);
const char * optstr = "mcmstmcmstmcmstmcm"; const char * optstr = "mcmstmcmstmcmstmcm";
for (size_t j = 1; j <= strlen(optstr); j++) for (size_t j = 1; j <= strlen(optstr); j++)
@ -492,11 +492,11 @@ namespace netgen
switch (optstr[j-1]) switch (optstr[j-1])
{ {
case 'c': optmesh.CombineImprove(mesh, OPT_REST); break; case 'c': optmesh.CombineImprove(); break;
case 'd': optmesh.SplitImprove(mesh, OPT_REST); break; case 'd': optmesh.SplitImprove(); break;
case 's': optmesh.SwapImprove(mesh, OPT_REST); break; case 's': optmesh.SwapImprove(); break;
case 't': optmesh.SwapImprove2(mesh, OPT_REST); break; case 't': optmesh.SwapImprove2(); break;
case 'm': mesh.ImproveMesh(mp, OPT_REST); break; case 'm': optmesh.ImproveMesh(); break;
} }
} }
@ -528,6 +528,7 @@ namespace netgen
throw NgException ("Stop meshing since surface mesh not consistent"); throw NgException ("Stop meshing since surface mesh not consistent");
} }
} }
RemoveIllegalElements (mesh);
} }
void MergeMeshes( Mesh & mesh, Array<MeshingData> & md ) void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
@ -682,13 +683,28 @@ namespace netgen
*/ */
mesh3d.CalcSurfacesOfNode(); mesh3d.CalcSurfacesOfNode();
MeshOptimize3d optmesh(mesh3d, mp);
// optimize only bad elements first
optmesh.SetMinBadness(1000.);
for(auto i : Range(mp.optsteps3d))
{
auto [total_badness, max_badness, bad_els] = optmesh.UpdateBadness();
if(bad_els==0) break;
optmesh.SplitImprove();
optmesh.SwapImprove();
optmesh.SwapImprove2();
}
// Now optimize all elements
optmesh.SetMinBadness(0);
for (auto i : Range(mp.optsteps3d)) for (auto i : Range(mp.optsteps3d))
{ {
if (multithread.terminate) if (multithread.terminate)
break; break;
MeshOptimize3d optmesh(mp);
// teterrpow = mp.opterrpow; // teterrpow = mp.opterrpow;
// for (size_t j = 1; j <= strlen(mp.optimize3d); j++) // for (size_t j = 1; j <= strlen(mp.optimize3d); j++)
for (auto j : Range(mp.optimize3d.size())) for (auto j : Range(mp.optimize3d.size()))
@ -699,12 +715,16 @@ namespace netgen
switch (mp.optimize3d[j]) switch (mp.optimize3d[j])
{ {
case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break; case 'c':
case 'd': optmesh.SplitImprove(mesh3d); break; optmesh.SetGoal(OPT_REST);
case 'D': optmesh.SplitImprove2(mesh3d); break; optmesh.CombineImprove();
case 's': optmesh.SwapImprove(mesh3d); break; optmesh.SetGoal(OPT_QUALITY);
break;
case 'd': optmesh.SplitImprove(); break;
case 'D': optmesh.SplitImprove2(); break;
case 's': optmesh.SwapImprove(); break;
// case 'u': optmesh.SwapImproveSurface(mesh3d); break; // case 'u': optmesh.SwapImproveSurface(mesh3d); break;
case 't': optmesh.SwapImprove2(mesh3d); break; case 't': optmesh.SwapImprove2(); break;
#ifdef SOLIDGEOM #ifdef SOLIDGEOM
case 'm': mesh3d.ImproveMesh(*geometry); break; case 'm': mesh3d.ImproveMesh(*geometry); break;
case 'M': mesh3d.ImproveMesh(*geometry); break; case 'M': mesh3d.ImproveMesh(*geometry); break;
@ -744,19 +764,19 @@ namespace netgen
nillegal = mesh3d.MarkIllegalElements(); nillegal = mesh3d.MarkIllegalElements();
MeshingParameters dummymp; MeshingParameters dummymp;
MeshOptimize3d optmesh(dummymp); MeshOptimize3d optmesh(mesh3d, dummymp, OPT_LEGAL);
while (nillegal && (it--) > 0) while (nillegal && (it--) > 0)
{ {
if (multithread.terminate) if (multithread.terminate)
break; break;
PrintMessage (5, nillegal, " illegal tets"); PrintMessage (5, nillegal, " illegal tets");
optmesh.SplitImprove (mesh3d, OPT_LEGAL); optmesh.SplitImprove ();
mesh3d.MarkIllegalElements(); // test mesh3d.MarkIllegalElements(); // test
optmesh.SwapImprove (mesh3d, OPT_LEGAL); optmesh.SwapImprove ();
mesh3d.MarkIllegalElements(); // test mesh3d.MarkIllegalElements(); // test
optmesh.SwapImprove2 (mesh3d, OPT_LEGAL); optmesh.SwapImprove2 ();
oldn = nillegal; oldn = nillegal;
nillegal = mesh3d.MarkIllegalElements(); nillegal = mesh3d.MarkIllegalElements();

26
libsrc/meshing/meshtype.hpp
View File

@ -13,6 +13,7 @@
#include <gprim/geom3d.hpp> #include <gprim/geom3d.hpp>
#include <linalg.hpp> #include <linalg.hpp>
#include "core/exception.hpp"
#include "msghandler.hpp" #include "msghandler.hpp"
namespace netgen namespace netgen
@ -994,7 +995,10 @@ namespace netgen
{ return flags.strongrefflag; } { return flags.strongrefflag; }
int Illegal () const int Illegal () const
{ return flags.illegal; } {
NETGEN_CHECK_SAME(flags.illegal_valid, true);
return flags.illegal;
}
int IllegalValid () const int IllegalValid () const
{ return flags.illegal_valid; } { return flags.illegal_valid; }
void SetIllegal (int aillegal) void SetIllegal (int aillegal)
@ -1008,6 +1012,26 @@ namespace netgen
flags.illegal_valid = 1; flags.illegal_valid = 1;
} }
bool BadnessValid()
{ return flags.badness_valid; }
float GetBadness()
{
NETGEN_CHECK_SAME(flags.badness_valid, true);
return badness;
}
void SetBadness(float value)
{
badness = value;
flags.badness_valid = 1;
}
void Touch() {
flags.illegal_valid = 0;
flags.badness_valid = 0;
}
void Delete () { flags.deleted = 1; } void Delete () { flags.deleted = 1; }
bool IsDeleted () const bool IsDeleted () const
{ {

12
libsrc/meshing/validate.cpp
View File

@ -470,11 +470,11 @@ namespace netgen
multithread.terminate != 1) multithread.terminate != 1)
{ {
MeshingParameters dummymp; MeshingParameters dummymp;
MeshOptimize3d optmesh(dummymp); MeshOptimize3d optmesh(mesh, dummymp, OPT_QUALITY);
for(int i=0; i<numtopimprove; i++) for(int i=0; i<numtopimprove; i++)
{ {
optmesh.SwapImproveSurface(mesh,OPT_QUALITY,&working_elements,&idmaps); optmesh.SwapImproveSurface(&working_elements,&idmaps);
optmesh.SwapImprove(mesh,OPT_QUALITY,&working_elements); optmesh.SwapImprove(&working_elements);
} }
@ -518,11 +518,11 @@ namespace netgen
} }
MeshingParameters dummymp; MeshingParameters dummymp;
MeshOptimize3d optmesh(dummymp); MeshOptimize3d optmesh(mesh, dummymp, OPT_QUALITY);
for(int i=0; i<numtopimprove && multithread.terminate != 1; i++) for(int i=0; i<numtopimprove && multithread.terminate != 1; i++)
{ {
optmesh.SwapImproveSurface(mesh,OPT_QUALITY,NULL,&idmaps); optmesh.SwapImproveSurface(NULL,&idmaps);
optmesh.SwapImprove(mesh,OPT_QUALITY); optmesh.SwapImprove();
//mesh.UpdateTopology(); //mesh.UpdateTopology();
} }
mesh.UpdateTopology(); mesh.UpdateTopology();

21
tests/pytest/compare_results.py
View File

@ -89,11 +89,15 @@ for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data
print(f"{GREEN}ntrigs {f1} got better: {bad1} -> {bad2}".ljust(w) + diff + RESET) print(f"{GREEN}ntrigs {f1} got better: {bad1} -> {bad2}".ljust(w) + diff + RESET)
n = len(data) + 1 n = len(data) + 1
fig,ax = plt.subplots(figsize=(10,7)) fig, ax = plt.subplots(figsize=(15, 7)) # Adjust figsize as needed
plt.xticks([])
plt.yticks([])
ax.yaxis.grid(False)
ax.xaxis.grid(False)
for i, d in enumerate(['min trig angle', 'min tet angle', 'max trig angle', 'max tet angle']): for i, d in enumerate(['min trig angle', 'min tet angle', 'max trig angle', 'max tet angle']):
ax = plt.subplot(2,5,i+1) plt.subplot(2, 4, i + 1) # Remove ax =
plt.title(d) plt.title(d)
ax.set_xticks([1,2]) # plt.xticks([1, 2])
if len(data[d]) == 0 or len(data2[d]) == 0: if len(data[d]) == 0 or len(data2[d]) == 0:
continue continue
@ -101,17 +105,18 @@ for i,d in enumerate(['min trig angle','min tet angle','max trig angle','max tet
med = statistics.median(data[d]) med = statistics.median(data[d])
plt.hlines(med, 1, 2, linestyle='dotted') plt.hlines(med, 1, 2, linestyle='dotted')
if d == 'badness': if d == 'badness':
ax.set_yscale('log') plt.yscale('log')
ax.set_xticklabels([ref, ref2]) plt.xticks([1, 2], [ref, ref2])
for i, d in enumerate(['badness', '#edges', '#trigs', '#tets']): for i, d in enumerate(['badness', '#edges', '#trigs', '#tets']):
ax = plt.subplot(2,5,6+i) plt.xticks([])
plt.subplot(2, 4, 5 + i)
plt.title('difference ' + d + ' (in %)') plt.title('difference ' + d + ' (in %)')
# plt.violinplot([(y-x)/x for x,y in zip(data[d],data2[d])], showmedians=True)
plt.boxplot([100.0 * (y - x) / x for x, y in zip(data[d], data2[d])]) plt.boxplot([100.0 * (y - x) / x for x, y in zip(data[d], data2[d])])
plt.hlines(0.0, 0.5, 1.5, linestyle='dotted') plt.hlines(0.0, 0.5, 1.5, linestyle='dotted')
plt.tight_layout() # Adjust layout
# plt.savefig('comparison.png', dpi=100) # plt.savefig('comparison.png', dpi=100)
plt.show() plt.show()

1954
tests/pytest/results.json
View File

File diff suppressed because it is too large Load Diff

93
tests/pytest/test_nonnative_master Normal file
View File

@ -0,0 +1,93 @@
tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214tstart (must be vertically): (0, 1.82574)
tend (must be vertically): (0, -1.82574)
generate boundarycondition.geo
needed 0.12342000007629395 seconds
generate boxcyl.geo
needed 0.801245927810669 seconds
generate circle_on_cube.geo
needed 0.48620080947875977 seconds
generate cone.geo
needed 1.186929702758789 seconds
generate cube.geo
needed 0.09043073654174805 seconds
generate cubeandring.geo
needed 2.1420021057128906 seconds
generate cubeandspheres.geo
needed 0.26427721977233887 seconds
generate cubemcyl.geo
needed 18.373886108398438 seconds
generate cubemsphere.geo
needed 3.6954052448272705 seconds
generate cylinder.geo
needed 0.44164204597473145 seconds
generate cylsphere.geo
needed 0.8774328231811523 seconds
generate ellipsoid.geo
needed 1.4510962963104248 seconds
generate ellipticcone.geo
needed 3.0906074047088623 seconds
generate ellipticcyl.geo
needed 2.0780415534973145 seconds
generate extrusion.geo
needed 0.40680599212646484 seconds
generate fichera.geo
needed 0.1265270709991455 seconds
generate hinge.stl
needed 5.519087553024292 seconds
generate lense.in2d
needed 0.05641365051269531 seconds
generate lshape3d.geo
needed 0.09937620162963867 seconds
generate manyholes.geo
needed 9.43863034248352 seconds
generate manyholes2.geo
needed 7.40019965171814 seconds
generate matrix.geo
needed 3.734792709350586 seconds
generate ortho.geo
needed 0.09516167640686035 seconds
generate part1.stl
needed 2.6940107345581055 seconds
generate period.geo
needed 2.709449291229248 seconds
generate revolution.geo
needed 7.7064368724823 seconds
generate sculpture.geo
needed 0.6283819675445557 seconds
generate shaft.geo
needed 2.921243190765381 seconds
generate shell.geo
generate sphere.geo
needed 0.18424725532531738 seconds
generate sphereincube.geo
needed 0.6060984134674072 seconds
generate square.in2d
needed 0.021883487701416016 seconds
generate squarecircle.in2d
needed 0.04081606864929199 seconds
generate squarehole.in2d
needed 0.03681302070617676 seconds
generate torus.geo
needed 6.590093612670898 seconds
generate trafo.geo
needed 3.2712368965148926 seconds
generate twobricks.geo
needed 0.13849091529846191 seconds
generate twocubes.geo
needed 0.13361692428588867 seconds
generate twocyl.geo
needed 0.8036918640136719 seconds
generate plane.stl
needed 15.712460041046143 seconds
done
sin (must not be 0) = 0.447214sin (must not be 0) = 0.447214

2
tests/pytest/test_tutorials.py
View File

@ -72,6 +72,8 @@ def getMeshingparameters(filename):
return standard[:-1] return standard[:-1]
if filename == "screw.step": if filename == "screw.step":
return standard[3:] # coarser meshes don't work here return standard[3:] # coarser meshes don't work here
if filename == "cylinder.geo":
return standard[0:-1] # very fine gives inconsistent reults
if filename == "cylsphere.geo": if filename == "cylsphere.geo":
return standard[0:2] + standard[3:] # coarse gives inconsistent reults (other mesh on MacOS) return standard[0:2] + standard[3:] # coarse gives inconsistent reults (other mesh on MacOS)
if filename == "part1.stl": if filename == "part1.stl":