diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 40a4a7d0..b9ee7585 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -323,7 +323,7 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh, (*testout) << "Total badness = " << totalbad << endl; } - auto elementsonnode = mesh.CreatePoint2ElementTable(); + auto elementsonnode = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr); Array> edges; BuildEdgeList(mesh, elementsonnode, edges); @@ -567,7 +567,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, double bad = 0.0; double badmax = 0.0; - auto elementsonnode = mesh.CreatePoint2ElementTable(); + auto elementsonnode = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr); Array elerrs(ne); @@ -1312,13 +1312,16 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, if(el.Flags().fixed) continue; + if(mp.only3D_domain_nr && mp.only3D_domain_nr != el.GetIndex()) + continue; + for (auto pi : el.PNums()) if(!free_points[pi]) free_points.SetBitAtomic(pi); } }); - auto elementsonnode = mesh.CreatePoint2ElementTable(free_points); + auto elementsonnode = mesh.CreatePoint2ElementTable(free_points, mp.only3D_domain_nr ); NgArray hasbothpoints; @@ -2480,7 +2483,8 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) // find elements on node - auto elementsonnode = mesh.CreatePoint2ElementTable(); + auto elementsonnode = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr); + // todo: respect mp.only3D_domain_nr for (SurfaceElementIndex sei = 0; sei < nse; sei++) for (int j = 0; j < 3; j++) @@ -2685,7 +2689,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh) static Timer topt("Optimize"); int ne = mesh.GetNE(); - auto elements_of_point = mesh.CreatePoint2ElementTable(); + auto elements_of_point = mesh.CreatePoint2ElementTable(nullopt, mp.only3D_domain_nr); int ntasks = 4*ngcore::TaskManager::GetNumThreads(); const char * savetask = multithread.task; diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index dac4e90a..48f10fb8 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -6814,7 +6814,7 @@ namespace netgen } - Table Mesh :: CreatePoint2ElementTable(std::optional points) const + Table Mesh :: CreatePoint2ElementTable(std::optional points, int domain) const { if(points) { @@ -6826,6 +6826,9 @@ namespace netgen if(el.IsDeleted()) return; + if(domain && el.GetIndex() != domain) + return; + for (PointIndex pi : el.PNums()) if(free_points[pi]) table.Add (pi, ei); @@ -6839,6 +6842,9 @@ namespace netgen if(el.IsDeleted()) return; + if(domain && el.GetIndex() != domain) + return; + for (PointIndex pi : el.PNums()) table.Add (pi, ei); }, GetNP()); diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index ff0e6e26..22ac8960 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -797,7 +797,7 @@ namespace netgen - DLL_HEADER Table CreatePoint2ElementTable(std::optional points = std::nullopt) const; + DLL_HEADER Table CreatePoint2ElementTable(std::optional points = std::nullopt, int domain = 0) const; DLL_HEADER Table CreatePoint2SurfaceElementTable( int faceindex=0 ) const; DLL_HEADER bool PureTrigMesh (int faceindex = 0) const; diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 617d5de7..a1e3bf67 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -304,22 +304,20 @@ namespace netgen public: Mesh::T_POINTS & points; const Array & elements; - Table &elementsonpoint; + Table &elementsonpoint; bool own_elementsonpoint; const MeshingParameters & mp; PointIndex actpind; double h; public: - PointFunction (Mesh::T_POINTS & apoints, - const Array & aelements, - const MeshingParameters & amp); + PointFunction (Mesh & mesh, const MeshingParameters & amp); PointFunction (const PointFunction & pf); virtual ~PointFunction () { if(own_elementsonpoint) delete &elementsonpoint; } virtual void SetPointIndex (PointIndex aactpind); void SetLocalH (double ah) { h = ah; } double GetLocalH () const { return h; } - const Table & GetPointToElementTable() { return elementsonpoint; }; + const Table & GetPointToElementTable() { return elementsonpoint; }; virtual double PointFunctionValue (const Point<3> & pp) const; virtual double PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const; virtual double PointFunctionValueDeriv (const Point<3> & pp, const Vec<3> & dir, double & deriv) const; @@ -332,23 +330,26 @@ namespace netgen : points(pf.points), elements(pf.elements), elementsonpoint(pf.elementsonpoint), own_elementsonpoint(false), mp(pf.mp) { } - PointFunction :: PointFunction (Mesh::T_POINTS & apoints, - const Array & aelements, - const MeshingParameters & amp) - : points(apoints), elements(aelements), elementsonpoint(* new Table()), own_elementsonpoint(true), mp(amp) + PointFunction :: PointFunction (Mesh & mesh, const MeshingParameters & amp) + : points(mesh.Points()), elements(mesh.VolumeElements()), elementsonpoint(* new Table()), own_elementsonpoint(true), mp(amp) { static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim); - elementsonpoint = ngcore::CreateSortedTable( elements.Range(), - [&](auto & table, ElementIndex ei) - { - const auto & el = elements[ei]; + BitArray free_points(points.Size()+PointIndex::BASE); + free_points.Clear(); - if(el.NP()!=4) - return; - - for (PointIndex pi : el.PNums()) - table.Add (pi, ei); - }, points.Size()); + ParallelForRange(elements.Range(), [&] (auto myrange) + { + for (ElementIndex eli : myrange) + { + const auto & el = elements[eli]; + if(el.Flags().fixed || el.NP()!=4 || (mp.only3D_domain_nr && mp.only3D_domain_nr != el.GetIndex()) ) + for (auto pi : el.PNums()) + if(free_points[pi]) + free_points.SetBitAtomic(pi); + } + }); + free_points.Invert(); + elementsonpoint = mesh.CreatePoint2ElementTable(free_points, mp.only3D_domain_nr); } void PointFunction :: SetPointIndex (PointIndex aactpind) @@ -493,19 +494,15 @@ namespace netgen { DenseMatrix m; public: - CheapPointFunction (Mesh::T_POINTS & apoints, - const Array & aelements, - const MeshingParameters & amp); + CheapPointFunction (Mesh & mesh, const MeshingParameters & amp); virtual void SetPointIndex (PointIndex aactpind); virtual double PointFunctionValue (const Point<3> & pp) const; virtual double PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const; }; - CheapPointFunction :: CheapPointFunction (Mesh::T_POINTS & apoints, - const Array & aelements, - const MeshingParameters & amp) - : PointFunction (apoints, aelements, amp) + CheapPointFunction :: CheapPointFunction (Mesh & mesh, const MeshingParameters & amp) + : PointFunction (mesh, amp) { ; } @@ -1348,14 +1345,14 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) int np = GetNP(); int ne = GetNE(); - PointFunction pf_glob(points, volelements, mp); + PointFunction pf_glob(*this, mp); auto & elementsonpoint = pf_glob.GetPointToElementTable(); const auto & getDofs = [&] (int i) { i += PointIndex::BASE; - return FlatArray(elementsonpoint[i].Size(), elementsonpoint[i].Data()); + return FlatArray(elementsonpoint[i].Size(), elementsonpoint[i].Data()); }; Array colors(points.Size());