ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'small_improvements' into 'master'

Browse Source 
Small improvements

See merge request jschoeberl/netgen!390
This commit is contained in:
Joachim Schöberl 2021-06-10 10:43:13 +00:00
commit d922d3abdd
12 changed files with 300 additions and 279 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
libsrc/core/CMakeLists.txt
View File

@ -69,7 +69,7 @@ install(TARGETS ngcore DESTINATION ${NG_INSTALL_DIR} COMPONENT netgen)
target_link_libraries(ngcore PUBLIC netgen_mpi PRIVATE "$<BUILD_INTERFACE:netgen_python>" ${CMAKE_THREAD_LIBS_INIT})
install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp logging.hpp
install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp logging.hpp memtracer.hpp
exception.hpp symboltable.hpp paje_trace.hpp utils.hpp profiler.hpp mpi_wrapper.hpp
array.hpp taskmanager.hpp concurrentqueue.h localheap.hpp python_ngcore.hpp flags.hpp
xbool.hpp signal.hpp bitarray.hpp table.hpp hashtable.hpp ranges.hpp

2
libsrc/core/array.hpp
View File

@ -11,7 +11,7 @@
#include "archive.hpp"
#include "exception.hpp"
#include "localheap.hpp"
#include "profiler.hpp"
#include "memtracer.hpp"
#include "utils.hpp"
namespace ngcore

174
libsrc/core/memtracer.hpp Normal file
View File

@ -0,0 +1,174 @@
#ifndef NETGEN_CORE_MEMTRACER_HPP
#define NETGEN_CORE_MEMTRACER_HPP
#include <array>
#include <chrono>
#include <functional>
#include <string>
#include "array.hpp"
#include "logging.hpp"
#include "paje_trace.hpp"
#include "utils.hpp"
namespace ngcore
{
class MemoryTracer;
namespace detail
{
//Type trait to check if a class implements a 'void SetMemoryTacing(int)' function
template<typename T>
struct has_StartMemoryTracing
{
private:
template<typename T2>
static constexpr auto check(T2*) ->
typename std::is_same<decltype(std::declval<T2>().StartMemoryTracing()),void>::type;
template<typename>
static constexpr std::false_type check(...);
using type = decltype(check<T>(nullptr)); // NOLINT
public:
static constexpr bool value = type::value;
};
} // namespace detail
class MemoryTracer
{
#ifdef NETGEN_TRACE_MEMORY
NGCORE_API static std::vector<std::string> names;
NGCORE_API static std::vector<int> parents;
static int CreateId(const std::string& name)
{
int id = names.size();
names.push_back(name);
parents.push_back(0);
if(id==10*8*1024)
std::cerr << "Allocated " << id << " MemoryTracer objects" << std::endl;
return id;
}
int id;
public:
MemoryTracer( std::string name )
{
id = CreateId(name);
}
// not tracing
MemoryTracer() : id(0) {}
template <typename... TRest>
MemoryTracer( std::string name, TRest & ... rest )
{
id = CreateId(name);
Track(rest...);
}
NETGEN_INLINE void Alloc(size_t size) const
{
if(id && trace)
trace->AllocMemory(id, size);
}
void Free(size_t size) const
{
if(id && trace)
trace->FreeMemory(id, size);
}
void Swap(size_t mysize, MemoryTracer& other, size_t other_size) const
{
if(!trace || (id == 0 && other.id == 0))
return;
if(id == 0)
return trace->ChangeMemory(other.id, mysize - other_size);
if(other.id == 0)
return trace->ChangeMemory(id, other_size - mysize);
// first decrease memory, otherwise have artificial/wrong high peak memory usage
if(mysize<other_size)
{
trace->ChangeMemory(other.id, mysize-other_size);
trace->ChangeMemory(id, other_size-mysize);
}
else
{
trace->ChangeMemory(id, other_size-mysize);
trace->ChangeMemory(other.id, mysize-other_size);
}
}
int GetId() const { return id; }
template <typename T1, typename... TRest>
void Track( T1 & obj, const std::string& name, TRest & ... rest ) const
{
Track(obj, name);
Track(rest...);
}
template<typename T>
void Track( T & obj, const std::string& name ) const
{
obj.GetMemoryTracer().Activate(obj, name);
parents[obj.GetMemoryTracer().GetId()] = id;
}
static std::string GetName(int id)
{
return names[id];
}
std::string GetName() const
{
return names[id];
}
template<typename T>
void Activate(T& me, const std::string& name) const
{
if(!id)
{
const_cast<MemoryTracer*>(this)->id = CreateId(name);
if constexpr(detail::has_StartMemoryTracing<T>::value)
me.StartMemoryTracing();
}
else
SetName(name);
}
void SetName(const std::string& name) const
{
names[id] = name;
}
static const std::vector<std::string> & GetNames() { return names; }
static const std::vector<int> & GetParents() { return parents; }
#else // NETGEN_TRACE_MEMORY
public:
MemoryTracer() {}
MemoryTracer( std::string /* name */ ) {}
template <typename... TRest>
MemoryTracer( std::string /* name */, TRest & ... ) {}
void Alloc(size_t /* size */) const {}
void Free(size_t /* size */) const {}
void Swap(...) const {}
int GetId() const { return 0; }
template <typename... TRest>
void Track(TRest&...) const {}
static std::string GetName(int /* id */) { return ""; }
std::string GetName() const { return ""; }
void SetName(std::string /* name */) const {}
#endif // NETGEN_TRACE_MEMORY
};
} // namespace ngcore
#endif // NETGEN_CORE_MEMTRACER_HPP

157
libsrc/core/profiler.hpp
View File

@ -6,8 +6,10 @@
#include <functional>
#include <string>
#include "array.hpp"
#include "logging.hpp"
#include "paje_trace.hpp"
#include "taskmanager.hpp"
#include "utils.hpp"
namespace ngcore
@ -300,161 +302,6 @@ namespace ngcore
return tres;
}
class MemoryTracer;
namespace detail
{
//Type trait to check if a class implements a 'void SetMemoryTacing(int)' function
template<typename T>
struct has_StartMemoryTracing
{
private:
template<typename T2>
static constexpr auto check(T2*) ->
typename std::is_same<decltype(std::declval<T2>().StartMemoryTracing()),void>::type;
template<typename>
static constexpr std::false_type check(...);
using type = decltype(check<T>(nullptr)); // NOLINT
public:
static constexpr bool value = type::value;
};
} // namespace detail
class MemoryTracer
{
#ifdef NETGEN_TRACE_MEMORY
NGCORE_API static std::vector<std::string> names;
NGCORE_API static std::vector<int> parents;
static int CreateId(const std::string& name)
{
int id = names.size();
names.push_back(name);
parents.push_back(0);
if(id==10*NgProfiler::SIZE)
std::cerr << "Allocated " << id << " MemoryTracer objects" << std::endl;
return id;
}
int id;
public:
MemoryTracer( std::string name )
{
id = CreateId(name);
}
// not tracing
MemoryTracer() : id(0) {}
template <typename... TRest>
MemoryTracer( std::string name, TRest & ... rest )
{
id = CreateId(name);
Track(rest...);
}
NETGEN_INLINE void Alloc(size_t size) const
{
if(id && trace)
trace->AllocMemory(id, size);
}
void Free(size_t size) const
{
if(id && trace)
trace->FreeMemory(id, size);
}
void Swap(size_t mysize, MemoryTracer& other, size_t other_size) const
{
if(!trace || (id == 0 && other.id == 0))
return;
if(id == 0)
return trace->ChangeMemory(other.id, mysize - other_size);
if(other.id == 0)
return trace->ChangeMemory(id, other_size - mysize);
// first decrease memory, otherwise have artificial/wrong high peak memory usage
if(mysize<other_size)
{
trace->ChangeMemory(other.id, mysize-other_size);
trace->ChangeMemory(id, other_size-mysize);
}
else
{
trace->ChangeMemory(id, other_size-mysize);
trace->ChangeMemory(other.id, mysize-other_size);
}
}
int GetId() const { return id; }
template <typename T1, typename... TRest>
void Track( T1 & obj, const std::string& name, TRest & ... rest ) const
{
Track(obj, name);
Track(rest...);
}
template<typename T>
void Track( T & obj, const std::string& name ) const
{
obj.GetMemoryTracer().Activate(obj, name);
parents[obj.GetMemoryTracer().GetId()] = id;
}
static std::string GetName(int id)
{
return names[id];
}
std::string GetName() const
{
return names[id];
}
template<typename T>
void Activate(T& me, const std::string& name) const
{
if(!id)
{
const_cast<MemoryTracer*>(this)->id = CreateId(name);
if constexpr(detail::has_StartMemoryTracing<T>::value)
me.StartMemoryTracing();
}
else
SetName(name);
}
void SetName(const std::string& name) const
{
names[id] = name;
}
static const std::vector<std::string> & GetNames() { return names; }
static const std::vector<int> & GetParents() { return parents; }
#else // NETGEN_TRACE_MEMORY
public:
MemoryTracer() {}
MemoryTracer( std::string /* name */ ) {}
template <typename... TRest>
MemoryTracer( std::string /* name */, TRest & ... ) {}
void Alloc(size_t /* size */) const {}
void Free(size_t /* size */) const {}
void Swap(...) const {}
int GetId() const { return 0; }
template <typename... TRest>
void Track(TRest&...) const {}
static std::string GetName(int /* id */) { return ""; }
std::string GetName() const { return ""; }
void SetName(std::string /* name */) const {}
#endif // NETGEN_TRACE_MEMORY
};
} // namespace ngcore
// Helper macro to easily add multiple timers in a function for profiling

67
libsrc/core/table.hpp
View File

@ -13,8 +13,10 @@
#include "array.hpp"
#include "bitarray.hpp"
#include "taskmanager.hpp"
#include "memtracer.hpp"
#include "ngcore_api.hpp"
#include "profiler.hpp"
#include "taskmanager.hpp"
namespace ngcore
{
@ -672,6 +674,69 @@ namespace ngcore
return s;
}
// Helper function to calculate coloring of a set of indices for parallel processing of independent elements/points/etc.
// Assigns a color to each of colors.Size() elements, such that two elements with the same color don't share a common 'dof',
// the mapping from element to dofs is provided by the function getDofs(int) -> iterable<int>
//
// Returns the number of used colors
template <typename Tmask>
int ComputeColoring( FlatArray<int> colors, size_t ndofs, Tmask const & getDofs)
{
static Timer timer("ComputeColoring - "+Demangle(typeid(Tmask).name())); RegionTimer rt(timer);
static_assert(sizeof(unsigned int)==4, "Adapt type of mask array");
size_t n = colors.Size();
Array<unsigned int> mask(ndofs);
size_t colored_blocks = 0;
// We are coloring with 32 colors at once and use each bit to mask conflicts
unsigned int check = 0;
unsigned int checkbit = 0;
int current_color = 0;
colors = -1;
int maxcolor = 0;
while(colored_blocks<n)
{
mask = 0;
for (auto i : Range(n) )
{
if(colors[i]>-1) continue;
check = 0;
const auto & dofs = getDofs(i);
// Check if adjacent dofs are already marked by current color
for (auto dof : dofs)
check|=mask[dof];
// Did we find a free color?
if(check != 0xFFFFFFFF)
{
checkbit = 1;
int color = current_color;
// find the actual color, which is free (out of 32)
while (check & checkbit)
{
color++;
checkbit *= 2;
}
colors[i] = color;
maxcolor = color > maxcolor ? color : maxcolor;
colored_blocks++;
// mask all adjacent dofs with the found color
for (auto dof : dofs)
mask[dof] |= checkbit;
}
}
current_color+=32;
}
return maxcolor+1;
}
typedef DynamicTable<int> IntTable;
} // namespace ngcore

64
libsrc/core/taskmanager.hpp
View File

@ -15,7 +15,7 @@
#include "array.hpp"
#include "paje_trace.hpp"
#include "profiler.hpp"
#include "taskmanager.hpp"
#ifdef USE_NUMA
#include <numa.h>
@ -1058,68 +1058,6 @@ public:
#endif // USE_NUMA
// Helper function to calculate coloring of a set of indices for parallel processing of independent elements/points/etc.
// Assigns a color to each of colors.Size() elements, such that two elements with the same color don't share a common 'dof',
// the mapping from element to dofs is provided by the function getDofs(int) -> iterable<int>
//
// Returns the number of used colors
template <typename Tmask>
int ComputeColoring( FlatArray<int> colors, size_t ndofs, Tmask const & getDofs)
{
static Timer timer("ComputeColoring - "+Demangle(typeid(Tmask).name())); RegionTimer rt(timer);
static_assert(sizeof(unsigned int)==4, "Adapt type of mask array");
size_t n = colors.Size();
Array<unsigned int> mask(ndofs);
size_t colored_blocks = 0;
// We are coloring with 32 colors at once and use each bit to mask conflicts
unsigned int check = 0;
unsigned int checkbit = 0;
int current_color = 0;
colors = -1;
int maxcolor = 0;
while(colored_blocks<n)
{
mask = 0;
for (auto i : Range(n) )
{
if(colors[i]>-1) continue;
check = 0;
const auto & dofs = getDofs(i);
// Check if adjacent dofs are already marked by current color
for (auto dof : dofs)
check|=mask[dof];
// Did we find a free color?
if(check != 0xFFFFFFFF)
{
checkbit = 1;
int color = current_color;
// find the actual color, which is free (out of 32)
while (check & checkbit)
{
color++;
checkbit *= 2;
}
colors[i] = color;
maxcolor = color > maxcolor ? color : maxcolor;
colored_blocks++;
// mask all adjacent dofs with the found color
for (auto dof : dofs)
mask[dof] |= checkbit;
}
}
current_color+=32;
}
return maxcolor+1;
}
}

7
libsrc/meshing/delaunay.cpp
View File

@ -233,7 +233,7 @@ namespace netgen
NgArray<Point<3> > & centers, NgArray<double> & radi2,
NgArray<int> & connected, NgArray<int> & treesearch,
NgArray<int> & freelist, SphereList & list,
IndexSet & insphere, IndexSet & closesphere)
IndexSet & insphere, IndexSet & closesphere, Array<DelaunayTet> & newels)
{
static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
// static Timer tsearch("addpoint, search");
@ -399,8 +399,6 @@ namespace netgen
}
} // while (changed)
// NgArray<Element> newels;
static NgArray<DelaunayTet> newels;
newels.SetSize(0);
Element2d face(TRIG);
@ -684,6 +682,7 @@ namespace netgen
for (PointIndex pi : mesh.Points().Range().Modify(0, -4))
mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE );
Array<DelaunayTet> newels;
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
for (PointIndex pi : mesh.Points().Range().Modify(0, -4))
{
@ -710,7 +709,7 @@ namespace netgen
AddDelaunayPoint (newpi, newp, tempels, mesh,
tettree, meshnb, centers, radi2,
connected, treesearch, freelist, list, insphere, closesphere);
connected, treesearch, freelist, list, insphere, closesphere, newels);
}

4
libsrc/meshing/localh.cpp
View File

@ -11,10 +11,6 @@ namespace netgen
for (int i = 0; i < 3; i++)
xmid[i] = 0.5 * (ax1[i] + ax2[i]);
for (int i = 0; i < 8; i++)
childs[i] = NULL;
father = NULL;
flags.cutboundary = 0;
flags.isinner = 0;
flags.oldcell = 0;

4
libsrc/meshing/localh.hpp
View File

@ -20,9 +20,9 @@ namespace netgen
/// half edgelength
float h2;
///
GradingBox * childs[8];
GradingBox * childs[8] = {nullptr};
///
GradingBox * father;
GradingBox * father = nullptr;
///
double hopt;
///

4
libsrc/meshing/meshtype.hpp
View File

@ -49,7 +49,7 @@ namespace netgen
#define ELEMENT2D_MAXPOINTS 8
enum POINTTYPE { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 };
enum POINTTYPE : unsigned char { FIXEDPOINT = 1, EDGEPOINT = 2, SURFACEPOINT = 3, INNERPOINT = 4 };
enum ELEMENTTYPE { FREEELEMENT, FIXEDELEMENT };
enum OPTIMIZEGOAL { OPT_QUALITY, OPT_CONFORM, OPT_REST, OPT_WORSTCASE, OPT_LEGAL };
@ -336,8 +336,8 @@ namespace netgen
*/
class MeshPoint : public Point<3>
{
int layer;
double singular; // singular factor for hp-refinement
int layer;
POINTTYPE type;

2
libsrc/meshing/ruler3.cpp
View File

@ -685,7 +685,7 @@ int Meshing3 :: ApplyRules
for (int i = 1; i <= lfaces.Size() && ok; i++)
{
static NgArray<int> lpi(4);
NgArrayMem<int, 10> lpi(4);
if (!fused.Get(i))
{

90
libsrc/meshing/smoothing3.cpp
View File

@ -304,7 +304,7 @@ namespace netgen
public:
Mesh::T_POINTS & points;
const Array<Element> & elements;
TABLE<int,PointIndex::BASE> &elementsonpoint;
Table<int, PointIndex> &elementsonpoint;
bool own_elementsonpoint;
const MeshingParameters & mp;
PointIndex actpind;
@ -319,7 +319,7 @@ namespace netgen
virtual void SetPointIndex (PointIndex aactpind);
void SetLocalH (double ah) { h = ah; }
double GetLocalH () const { return h; }
const TABLE<int,PointIndex::BASE> & GetPointToElementTable() { return elementsonpoint; };
const Table<int, PointIndex> & 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;
@ -335,13 +335,20 @@ namespace netgen
PointFunction :: PointFunction (Mesh::T_POINTS & apoints,
const Array<Element> & aelements,
const MeshingParameters & amp)
: points(apoints), elements(aelements), elementsonpoint(* new TABLE<int,PointIndex::BASE>(apoints.Size())), own_elementsonpoint(true), mp(amp)
: points(apoints), elements(aelements), elementsonpoint(* new Table<int,PointIndex>()), own_elementsonpoint(true), mp(amp)
{
static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim);
for (int i = 0; i < elements.Size(); i++)
if (elements[i].NP() == 4)
for (int j = 0; j < elements[i].NP(); j++)
elementsonpoint.Add (elements[i][j], i);
elementsonpoint = std::move(ngcore::CreateSortedTable<int, PointIndex>( elements.Range(),
[&](auto & table, ElementIndex ei)
{
const auto & el = elements[ei];
if(el.NP()!=4)
return;
for (PointIndex pi : el.PNums())
table.Add (pi, ei);
}, points.Size()));
}
void PointFunction :: SetPointIndex (PointIndex aactpind)
@ -359,9 +366,9 @@ namespace netgen
hp = points[actpind];
points[actpind] = Point<3> (pp);
for (int j = 0; j < elementsonpoint[actpind].Size(); j++)
for (auto ei : elementsonpoint[actpind])
{
const Element & el = elements[elementsonpoint[actpind][j]];
const Element & el = elements[ei];
badness += CalcTetBadness (points[el[0]], points[el[1]],
points[el[2]], points[el[3]], -1, mp);
}
@ -379,9 +386,9 @@ namespace netgen
Vec<3> vgradi, vgrad(0,0,0);
points[actpind] = Point<3> (pp);
for (int j = 0; j < elementsonpoint[actpind].Size(); j++)
for (auto ei : elementsonpoint[actpind])
{
const Element & el = elements[elementsonpoint[actpind][j]];
const Element & el = elements[ei];
for (int k = 0; k < 4; k++)
if (el[k] == actpind)
{
@ -409,9 +416,9 @@ namespace netgen
points[actpind] = pp;
double f = 0;
for (int j = 0; j < elementsonpoint[actpind].Size(); j++)
for (auto ei : elementsonpoint[actpind])
{
const Element & el = elements[elementsonpoint[actpind][j]];
const Element & el = elements[ei];
for (int k = 1; k <= 4; k++)
if (el.PNum(k) == actpind)
@ -435,10 +442,9 @@ namespace netgen
// try point movement
NgArray<Element2d> faces;
for (int j = 0; j < elementsonpoint[actpind].Size(); j++)
for (auto ei : elementsonpoint[actpind])
{
const Element & el =
elements[elementsonpoint[actpind][j]];
const Element & el = elements[ei];
for (int k = 1; k <= 4; k++)
if (el.PNum(k) == actpind)
@ -1013,11 +1019,8 @@ double JacobianPointFunction :: Func (const Vector & v) const
points[actpind] -= (v(0)*nv(0)+v(1)*nv(1)+v(2)*nv(2)) * nv;
for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
{
int eli = elementsonpoint.Get(actpind, j);
badness += elements[eli-1].CalcJacobianBadness (points);
}
for (auto eli : elementsonpoint[actpind])
badness += elements[eli].CalcJacobianBadness (points);
points[actpind] = hp;
@ -1046,10 +1049,9 @@ FuncGrad (const Vector & x, Vector & g) const
g.SetSize(3);
g = 0;
for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
for (auto ei : elementsonpoint[actpind])
{
int eli = elementsonpoint.Get(actpind, j);
const Element & el = elements[eli-1];
const Element & el = elements[ei];
lpi = 0;
for (k = 1; k <= el.GetNP(); k++)
@ -1057,7 +1059,7 @@ FuncGrad (const Vector & x, Vector & g) const
lpi = k;
if (!lpi) cerr << "loc point not found" << endl;
badness += elements[eli-1].
badness += elements[ei].
CalcJacobianBadnessGradient (points, lpi, hderiv);
for(k=0; k<3; k++)
@ -1119,10 +1121,9 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
vdir -= scal*nv;
}
for (j = 1; j <= elementsonpoint.EntrySize(actpind); j++)
for (auto ei : elementsonpoint[actpind])
{
int eli = elementsonpoint.Get(actpind, j);
const Element & el = elements[eli-1];
const Element & el = elements[ei];
lpi = 0;
for (k = 1; k <= el.GetNP(); k++)
@ -1130,7 +1131,7 @@ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
lpi = k;
if (!lpi) cerr << "loc point not found" << endl;
badness += elements[eli-1].
badness += elements[ei].
CalcJacobianBadnessDirDeriv (points, lpi, vdir, hderiv);
deriv += hderiv;
}
@ -1458,6 +1459,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
static Timer tcalcbadmax("Calc badmax");
static Timer topt("optimize");
static Timer trange("range");
static Timer tloch("loch");
// return ImproveMeshSequential(mp, goal);
BuildBoundaryEdges(false);
@ -1475,24 +1477,20 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
const auto & getDofs = [&] (int i)
{
i += PointIndex::BASE;
return FlatArray<int>(elementsonpoint[i].Size(), &elementsonpoint[i][0]);
return FlatArray<int>(elementsonpoint[i].Size(), elementsonpoint[i].Data());
};
Array<int> colors(points.Size());
tcoloring.Start();
int ncolors = ngcore::ComputeColoring( colors, ne, getDofs );
TableCreator<int> creator(ncolors);
for ( ; !creator.Done(); creator++)
auto color_table = CreateTable<PointIndex, int>( points.Size(),
[&] ( auto & table, int i )
{
ParallelForRange( Range(colors), [&](auto myrange)
{
for(auto i : myrange)
creator.Add(colors[i], i);
});
}
PointIndex pi = i+static_cast<int>(PointIndex::BASE);
table.Add(colors[i], pi);
}, ncolors);
auto color_table = creator.MoveTable();
tcoloring.Stop();
if (goal == OPT_QUALITY)
@ -1505,12 +1503,16 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
(*testout) << setprecision(8);
NgArray<double, PointIndex::BASE> pointh (points.Size());
Array<double, PointIndex> pointh (points.Size());
if(lochfunc)
{
for (PointIndex pi : points.Range())
RegionTimer rt(tloch);
ParallelForRange(points.Range(), [&] (auto myrange)
{
for(auto pi : myrange)
pointh[pi] = GetH(points[pi]);
});
}
else
{
@ -1529,12 +1531,12 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
topt.Start();
int counter = 0;
for (int color : Range(color_table.Size()))
for (auto icolor : Range(ncolors))
{
if (multithread.terminate)
throw NgException ("Meshing stopped");
ParallelForRange( Range(color_table[color].Size()), [&](auto myrange)
ParallelForRange( color_table[icolor].Range(), [&](auto myrange)
{
RegionTracer reg(ngcore::TaskManager::GetThreadId(), trange, myrange.Size());
Vector x(3);
@ -1549,7 +1551,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
for (auto i : myrange)
{
PointIndex pi(color_table[color][i]+PointIndex::BASE);
PointIndex pi = color_table[icolor][i];
if ( (*this)[pi].Type() == INNERPOINT )
{
counter++;