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

Merge branch 'NGSolve:master' into master

Browse Source
This commit is contained in:
Julius Zimmermann 2022-04-20 00:01:15 +02:00 committed by GitHub
commit eb8276c132
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
21 changed files with 1723 additions and 3300 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

319
libsrc/core/python_ngcore.hpp
View File

@ -156,165 +156,7 @@ namespace ngcore
{ return std::string("sp_")+GetPyName<T>(); } { return std::string("sp_")+GetPyName<T>(); }
}; };
template<typename T> // *************** Archiving functionality **************
Array<T> makeCArray(const py::object& obj)
{
Array<T> arr;
if(py::isinstance<py::list>(obj))
for(auto& val : py::cast<py::list>(obj))
arr.Append(py::cast<T>(val));
else if(py::isinstance<py::tuple>(obj))
for(auto& val : py::cast<py::tuple>(obj))
arr.Append(py::cast<T>(val));
else
throw py::type_error("Cannot convert Python object to C Array");
return arr;
}
template <typename T, typename TIND=typename FlatArray<T>::index_type>
void ExportArray (py::module &m)
{
using TFlat = FlatArray<T, TIND>;
using TArray = Array<T, TIND>;
std::string suffix = GetPyName<T>() + "_" +
GetPyName<TIND>();
std::string fname = std::string("FlatArray_") + suffix;
auto flatarray_class = py::class_<TFlat>(m, fname.c_str(),
py::buffer_protocol())
.def ("__len__", [] ( TFlat &self ) { return self.Size(); } )
.def ("__getitem__",
[](TFlat & self, TIND i) -> T&
{
static constexpr int base = IndexBASE<TIND>();
if (i < base || i >= self.Size()+base)
throw py::index_error();
return self[i];
},
py::return_value_policy::reference)
.def ("__setitem__",
[](TFlat & self, TIND i, T val) -> T&
{
static constexpr int base = IndexBASE<TIND>();
if (i < base || i >= self.Size()+base)
throw py::index_error();
self[i] = val;
return self[i];
},
py::return_value_policy::reference)
.def ("__setitem__",
[](TFlat & self, py::slice slice, T val)
{
size_t start, stop, step, slicelength;
if (!slice.compute(self.Size(), &start, &stop, &step, &slicelength))
throw py::error_already_set();
static constexpr int base = IndexBASE<TIND>();
if (start < base || start+(slicelength-1)*step >= self.Size()+base)
throw py::index_error();
for (size_t i = 0; i < slicelength; i++, start+=step)
self[start] = val;
})
.def("__iter__", [] ( TFlat & self) {
return py::make_iterator (self.begin(),self.end());
}, py::keep_alive<0,1>()) // keep array alive while iterator is used
.def("__str__", [](TFlat& self)
{
return ToString(self);
})
;
if constexpr (detail::HasPyFormat<T>::value)
{
if(ngcore_have_numpy && !py::detail::npy_format_descriptor<T>::dtype().is_none())
{
flatarray_class
.def_buffer([](TFlat& self)
{
return py::buffer_info(
self.Addr(0),
sizeof(T),
py::format_descriptor<T>::format(),
1,
{ self.Size() },
{ sizeof(T) * (self.Addr(1) - self.Addr(0)) });
})
.def("NumPy", [](py::object self)
{
return py::module::import("numpy")
.attr("frombuffer")(self, py::detail::npy_format_descriptor<T>::dtype());
})
;
}
}
std::string aname = std::string("Array_") + suffix;
py::class_<TArray, TFlat>(m, aname.c_str())
.def(py::init([] (size_t n) { return new TArray(n); }),py::arg("n"), "Makes array of given length")
.def(py::init([] (std::vector<T> const & x)
{
size_t s = x.size();
TArray tmp(s);
for (size_t i : Range(tmp))
tmp[TIND(i)] = x[i];
return tmp;
}), py::arg("vec"), "Makes array with given list of elements")
;
py::implicitly_convertible<std::vector<T>, TArray>();
}
template <typename T>
void ExportTable (py::module &m)
{
py::class_<ngcore::Table<T>, std::shared_ptr<ngcore::Table<T>>> (m, ("Table_"+GetPyName<T>()).c_str())
.def(py::init([] (py::list blocks)
{
size_t size = py::len(blocks);
Array<int> cnt(size);
size_t i = 0;
for (auto block : blocks)
cnt[i++] = py::len(block);
i = 0;
Table<T> blocktable(cnt);
for (auto block : blocks)
{
auto row = blocktable[i++];
size_t j = 0;
for (auto val : block)
row[j++] = val.cast<T>();
}
// cout << "blocktable = " << *blocktable << endl;
return blocktable;
}), py::arg("blocks"), "a list of lists")
.def ("__len__", [] (Table<T> &self ) { return self.Size(); } )
.def ("__getitem__",
[](Table<T> & self, size_t i) -> FlatArray<T>
{
if (i >= self.Size())
throw py::index_error();
return self[i];
})
.def("__str__", [](Table<T> & self)
{
return ToString(self);
})
;
}
void NGCORE_API SetFlag(Flags &flags, std::string s, py::object value);
// Parse python kwargs to flags
Flags NGCORE_API CreateFlagsFromKwArgs(const py::kwargs& kwargs, py::object pyclass = py::none(),
py::list info = py::list());
// Create python dict from kwargs
py::dict NGCORE_API CreateDictFromFlags(const Flags& flags);
// *************** Archiving functionality **************
template<typename T> template<typename T>
Archive& Archive :: Shallow(T& val) Archive& Archive :: Shallow(T& val)
@ -429,6 +271,165 @@ namespace ngcore
}); });
} }
template<typename T>
Array<T> makeCArray(const py::object& obj)
{
Array<T> arr;
if(py::isinstance<py::list>(obj))
for(auto& val : py::cast<py::list>(obj))
arr.Append(py::cast<T>(val));
else if(py::isinstance<py::tuple>(obj))
for(auto& val : py::cast<py::tuple>(obj))
arr.Append(py::cast<T>(val));
else
throw py::type_error("Cannot convert Python object to C Array");
return arr;
}
template <typename T, typename TIND=typename FlatArray<T>::index_type>
void ExportArray (py::module &m)
{
using TFlat = FlatArray<T, TIND>;
using TArray = Array<T, TIND>;
std::string suffix = GetPyName<T>() + "_" +
GetPyName<TIND>();
std::string fname = std::string("FlatArray_") + suffix;
auto flatarray_class = py::class_<TFlat>(m, fname.c_str(),
py::buffer_protocol())
.def ("__len__", [] ( TFlat &self ) { return self.Size(); } )
.def ("__getitem__",
[](TFlat & self, TIND i) -> T&
{
static constexpr int base = IndexBASE<TIND>();
if (i < base || i >= self.Size()+base)
throw py::index_error();
return self[i];
},
py::return_value_policy::reference)
.def ("__setitem__",
[](TFlat & self, TIND i, T val) -> T&
{
static constexpr int base = IndexBASE<TIND>();
if (i < base || i >= self.Size()+base)
throw py::index_error();
self[i] = val;
return self[i];
},
py::return_value_policy::reference)
.def ("__setitem__",
[](TFlat & self, py::slice slice, T val)
{
size_t start, stop, step, slicelength;
if (!slice.compute(self.Size(), &start, &stop, &step, &slicelength))
throw py::error_already_set();
static constexpr int base = IndexBASE<TIND>();
if (start < base || start+(slicelength-1)*step >= self.Size()+base)
throw py::index_error();
for (size_t i = 0; i < slicelength; i++, start+=step)
self[start] = val;
})
.def("__iter__", [] ( TFlat & self) {
return py::make_iterator (self.begin(),self.end());
}, py::keep_alive<0,1>()) // keep array alive while iterator is used
.def("__str__", [](TFlat& self)
{
return ToString(self);
})
;
if constexpr (detail::HasPyFormat<T>::value)
{
if(ngcore_have_numpy && !py::detail::npy_format_descriptor<T>::dtype().is_none())
{
flatarray_class
.def_buffer([](TFlat& self)
{
return py::buffer_info(
self.Addr(0),
sizeof(T),
py::format_descriptor<T>::format(),
1,
{ self.Size() },
{ sizeof(T) * (self.Addr(1) - self.Addr(0)) });
})
.def("NumPy", [](py::object self)
{
return py::module::import("numpy")
.attr("frombuffer")(self, py::detail::npy_format_descriptor<T>::dtype());
})
;
}
}
std::string aname = std::string("Array_") + suffix;
auto arr = py::class_<TArray, TFlat> (m, aname.c_str())
.def(py::init([] (size_t n) { return new TArray(n); }),py::arg("n"), "Makes array of given length")
.def(py::init([] (std::vector<T> const & x)
{
size_t s = x.size();
TArray tmp(s);
for (size_t i : Range(tmp))
tmp[TIND(i)] = x[i];
return tmp;
}), py::arg("vec"), "Makes array with given list of elements")
;
if constexpr(is_archivable<TArray>)
arr.def(NGSPickle<TArray>());
py::implicitly_convertible<std::vector<T>, TArray>();
}
template <typename T>
void ExportTable (py::module &m)
{
py::class_<ngcore::Table<T>, std::shared_ptr<ngcore::Table<T>>> (m, ("Table_"+GetPyName<T>()).c_str())
.def(py::init([] (py::list blocks)
{
size_t size = py::len(blocks);
Array<int> cnt(size);
size_t i = 0;
for (auto block : blocks)
cnt[i++] = py::len(block);
i = 0;
Table<T> blocktable(cnt);
for (auto block : blocks)
{
auto row = blocktable[i++];
size_t j = 0;
for (auto val : block)
row[j++] = val.cast<T>();
}
// cout << "blocktable = " << *blocktable << endl;
return blocktable;
}), py::arg("blocks"), "a list of lists")
.def ("__len__", [] (Table<T> &self ) { return self.Size(); } )
.def ("__getitem__",
[](Table<T> & self, size_t i) -> FlatArray<T>
{
if (i >= self.Size())
throw py::index_error();
return self[i];
})
.def("__str__", [](Table<T> & self)
{
return ToString(self);
})
;
}
void NGCORE_API SetFlag(Flags &flags, std::string s, py::object value);
// Parse python kwargs to flags
Flags NGCORE_API CreateFlagsFromKwArgs(const py::kwargs& kwargs, py::object pyclass = py::none(),
py::list info = py::list());
// Create python dict from kwargs
py::dict NGCORE_API CreateDictFromFlags(const Flags& flags);
} // namespace ngcore } // namespace ngcore

1
libsrc/core/simd_arm64.hpp
View File

@ -42,6 +42,7 @@ namespace ngcore
SIMD (const SIMD &) = default; SIMD (const SIMD &) = default;
// SIMD (double v0, double v1) : data{v0,v1} { } // SIMD (double v0, double v1) : data{v0,v1} { }
SIMD (double v0, double v1) : data{vcombine_f64(float64x1_t{v0}, float64x1_t{v1})} { } SIMD (double v0, double v1) : data{vcombine_f64(float64x1_t{v0}, float64x1_t{v1})} { }
SIMD (SIMD<double,1> v0, SIMD<double,1> v1) : data{vcombine_f64(float64x1_t{v0.Data()}, float64x1_t{v1.Data()})} { }
SIMD (std::array<double, 2> arr) : data{arr[0], arr[1]} { } SIMD (std::array<double, 2> arr) : data{arr[0], arr[1]} { }
SIMD & operator= (const SIMD &) = default; SIMD & operator= (const SIMD &) = default;

9
libsrc/core/simd_avx.hpp
View File

@ -143,8 +143,6 @@ namespace ngcore
NETGEN_INLINE double & operator[] (int i) { return ((double*)(&data))[i]; } NETGEN_INLINE double & operator[] (int i) { return ((double*)(&data))[i]; }
// [[deprecated("don't write to individual elements of SIMD")]] // [[deprecated("don't write to individual elements of SIMD")]]
// NETGEN_INLINE double & operator[] (int i) { return ((double*)(&data))[i]; } // NETGEN_INLINE double & operator[] (int i) { return ((double*)(&data))[i]; }
template <int I>
double Get() const { return ((double*)(&data))[I]; }
NETGEN_INLINE __m256d Data() const { return data; } NETGEN_INLINE __m256d Data() const { return data; }
NETGEN_INLINE __m256d & Data() { return data; } NETGEN_INLINE __m256d & Data() { return data; }
@ -153,6 +151,13 @@ namespace ngcore
operator std::tuple<double&,double&,double&,double&> () operator std::tuple<double&,double&,double&,double&> ()
{ return std::tuple<double&,double&,double&,double&>((*this)[0], (*this)[1], (*this)[2], (*this)[3]); } { return std::tuple<double&,double&,double&,double&>((*this)[0], (*this)[1], (*this)[2], (*this)[3]); }
template <int I>
double Get() const
{
static_assert(I>=0 && I<4, "Index out of range");
return (*this)[I];
}
}; };
NETGEN_INLINE auto Unpack (SIMD<double,4> a, SIMD<double,4> b) NETGEN_INLINE auto Unpack (SIMD<double,4> a, SIMD<double,4> b)

12
libsrc/core/simd_avx512.hpp
View File

@ -92,6 +92,12 @@ namespace ngcore
SIMD (double const * p, SIMD<mask64,8> mask) SIMD (double const * p, SIMD<mask64,8> mask)
{ data = _mm512_mask_loadu_pd(_mm512_setzero_pd(), mask.Data(), p); } { data = _mm512_mask_loadu_pd(_mm512_setzero_pd(), mask.Data(), p); }
SIMD (__m512d _data) { data = _data; } SIMD (__m512d _data) { data = _data; }
SIMD (SIMD<double,4> v0, SIMD<double,4> v1)
: data(_mm512_set_pd(v1[3], v1[2], v1[1], v1[0], v0[3], v0[2], v0[1], v0[0]))
{}
SIMD (SIMD<double,6> v0, SIMD<double,2> v1)
: data(_mm512_set_pd(v1[1], v1[0], v0[5], v0[4], v0[3], v0[2], v0[1], v0[0]))
{}
template<typename T, typename std::enable_if<std::is_convertible<T, std::function<double(int)>>::value, int>::type = 0> template<typename T, typename std::enable_if<std::is_convertible<T, std::function<double(int)>>::value, int>::type = 0>
SIMD (const T & func) SIMD (const T & func)
@ -129,6 +135,12 @@ namespace ngcore
NETGEN_INLINE __m512d Data() const { return data; } NETGEN_INLINE __m512d Data() const { return data; }
NETGEN_INLINE __m512d & Data() { return data; } NETGEN_INLINE __m512d & Data() { return data; }
template <int I>
double Get() const
{
static_assert(I>=0 && I<8, "Index out of range");
return (*this)[I];
}
}; };
NETGEN_INLINE SIMD<double,8> operator- (SIMD<double,8> a) { return -a.Data(); } NETGEN_INLINE SIMD<double,8> operator- (SIMD<double,8> a) { return -a.Data(); }

36
libsrc/core/simd_generic.hpp
View File

@ -28,6 +28,28 @@ namespace ngcore
#endif #endif
} }
constexpr bool IsNativeSIMDSize(int n) {
if(n==1) return true;
#if defined NETGEN_ARCH_AMD64 || defined __SSE__ || defined __aarch64__
if(n==2) return true;
#endif
#if defined __AVX__
if(n==4) return true;
#endif
#if defined __AVX512F__
if(n==8) return true;
#endif
return false;
}
// split n = k+l such that k is the largest natively supported simd size < n
constexpr int GetLargestNativeSIMDPart(int n) {
int k = n-1;
while(!IsNativeSIMDSize(k))
k--;
return k;
}
template <typename T, int N=GetDefaultSIMDSize()> class SIMD; template <typename T, int N=GetDefaultSIMDSize()> class SIMD;
@ -67,9 +89,9 @@ namespace ngcore
template <int N> template <int N>
class alignas(GetDefaultSIMDSize()*sizeof(int64_t)) SIMD<mask64,N> class alignas(GetLargestNativeSIMDPart(N)*sizeof(int64_t)) SIMD<mask64,N>
{ {
static constexpr int N1 = std::min(GetDefaultSIMDSize(), N/2); static constexpr int N1 = GetLargestNativeSIMDPart(N);
static constexpr int N2 = N-N1; static constexpr int N2 = N-N1;
SIMD<mask64,N1> lo; SIMD<mask64,N1> lo;
@ -123,9 +145,9 @@ namespace ngcore
}; };
template<int N> template<int N>
class alignas(GetDefaultSIMDSize()*sizeof(int64_t)) SIMD<int64_t,N> class alignas(GetLargestNativeSIMDPart(N)*sizeof(int64_t)) SIMD<int64_t,N>
{ {
static constexpr int N1 = std::min(GetDefaultSIMDSize(), N/2); static constexpr int N1 = GetLargestNativeSIMDPart(N);
static constexpr int N2 = N-N1; static constexpr int N2 = N-N1;
SIMD<int64_t,N1> lo; SIMD<int64_t,N1> lo;
@ -240,9 +262,9 @@ namespace ngcore
template<int N> template<int N>
class alignas(GetDefaultSIMDSize()*sizeof(double)) SIMD<double, N> class alignas(GetLargestNativeSIMDPart(N)*sizeof(double)) SIMD<double, N>
{ {
static constexpr int N1 = std::min(GetDefaultSIMDSize(), N/2); static constexpr int N1 = GetLargestNativeSIMDPart(N);
static constexpr int N2 = N-N1; static constexpr int N2 = N-N1;
SIMD<double, N1> lo; SIMD<double, N1> lo;
@ -543,7 +565,7 @@ namespace ngcore
template <int i, typename T, int N> template <int i, typename T, int N>
T get(SIMD<T,N> a) { return a[i]; } T get(SIMD<T,N> a) { return a.template Get<i>(); }
template <int NUM, typename FUNC> template <int NUM, typename FUNC>
NETGEN_INLINE void Iterate2 (FUNC f) NETGEN_INLINE void Iterate2 (FUNC f)

10
libsrc/core/simd_sse.hpp
View File

@ -86,6 +86,9 @@ NETGEN_INLINE SIMD<int64_t,2> operator- (SIMD<int64_t,2> a, SIMD<int64_t,2> b) {
SIMD () {} SIMD () {}
SIMD (const SIMD &) = default; SIMD (const SIMD &) = default;
SIMD (double v0, double v1) { data = _mm_set_pd(v1,v0); } SIMD (double v0, double v1) { data = _mm_set_pd(v1,v0); }
SIMD (SIMD<double,1> v0, SIMD<double,1> v1)
: data{_mm_set_pd(v0.Data(), v1.Data())}
{ }
SIMD (std::array<double, 2> arr) SIMD (std::array<double, 2> arr)
: data{_mm_set_pd(arr[1], arr[0])} : data{_mm_set_pd(arr[1], arr[0])}
{} {}
@ -137,6 +140,13 @@ NETGEN_INLINE SIMD<int64_t,2> operator- (SIMD<int64_t,2> a, SIMD<int64_t,2> b) {
NETGEN_INLINE __m128d Data() const { return data; } NETGEN_INLINE __m128d Data() const { return data; }
NETGEN_INLINE __m128d & Data() { return data; } NETGEN_INLINE __m128d & Data() { return data; }
template <int I>
double Get()
{
static_assert(I>=0 && I<2, "Index out of range");
return (*this)[I];
}
operator std::tuple<double&,double&> () operator std::tuple<double&,double&> ()
{ {
auto pdata = (double*)&data; auto pdata = (double*)&data;

17
libsrc/meshing/boundarylayer.cpp
View File

@ -696,8 +696,6 @@ namespace netgen
else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut())) else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut()))
{ {
type = 3; type = 3;
if(fd.DomainIn() == 0 || fd.DomainOut() == 0)
is_boundary_projected.SetBit(segj.si);
is_boundary_moved.SetBit(segj.si); is_boundary_moved.SetBit(segj.si);
} }
else else
@ -742,6 +740,8 @@ namespace netgen
else else
continue; continue;
if(!params.project_boundaries.Contains(sel.GetIndex()))
continue;
auto& g = growthvectors[pi]; auto& g = growthvectors[pi];
auto ng = n * g; auto ng = n * g;
auto gg = g * g; auto gg = g * g;
@ -818,11 +818,22 @@ namespace netgen
point_found = true; point_found = true;
break; break;
} }
else if(seg[1] == points.Last() &&
points[points.Size()-2] != seg[0])
{
edge_len += (mesh[points.Last()] - mesh[seg[0]]).Length();
points.Append(seg[0]);
point_found = true;
break;
}
} }
if(is_end_point(points.Last())) if(is_end_point(points.Last()))
break; break;
if(!point_found) if(!point_found)
throw Exception(string("Could not find connected list of line segments for edge ") + edgenr); {
cout << "points = " << points << endl;
throw Exception(string("Could not find connected list of line segments for edge ") + edgenr);
}
} }
// tangential part of growth vectors // tangential part of growth vectors

22
libsrc/meshing/improve2.cpp
View File

@ -145,7 +145,7 @@ namespace netgen
} }
else else
{ {
double loch = mesh.GetH(mesh[pi1]); double loch = 0.25*(mesh.GetH(pi1) + mesh.GetH(pi2) + mesh.GetH(pi3) + mesh.GetH(pi4));
should = should =
CalcTriangleBadness (mesh[pi4], mesh[pi3], mesh[pi1], metricweight, loch) + CalcTriangleBadness (mesh[pi4], mesh[pi3], mesh[pi1], metricweight, loch) +
CalcTriangleBadness (mesh[pi3], mesh[pi4], mesh[pi2], metricweight, loch) < CalcTriangleBadness (mesh[pi3], mesh[pi4], mesh[pi2], metricweight, loch) <
@ -383,21 +383,10 @@ namespace netgen
<< "pi1 = " << pi1 << " pi2 = " << pi2 << endl; << "pi1 = " << pi1 << " pi2 = " << pi2 << endl;
} }
/*
// save version:
if (fixed.Get(pi1) || fixed.Get(pi2))
return 0.0;
if (pi2 < pi1) swap (pi1, pi2);
*/
// more general
if (fixed[pi2])
Swap (pi1, pi2);
if (fixed[pi2]) if (fixed[pi2])
return 0.0; return 0.0;
double loch = mesh.GetH (mesh[pi1]); double loch = 0.5*(mesh.GetH(pi1) + mesh.GetH(pi2));
int faceindex = -1; int faceindex = -1;
for (SurfaceElementIndex sei2 : elementsonnode[pi1]) for (SurfaceElementIndex sei2 : elementsonnode[pi1])
@ -655,6 +644,9 @@ namespace netgen
double d_badness = CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, metricweight, true); double d_badness = CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, metricweight, true);
if(d_badness < 0.0) if(d_badness < 0.0)
candidate_edges[improvement_counter++] = make_tuple(d_badness, i); candidate_edges[improvement_counter++] = make_tuple(d_badness, i);
d_badness = CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi2, pi1, metricweight, true);
if(d_badness < 0.0)
candidate_edges[improvement_counter++] = make_tuple(d_badness, -i);
}, TasksPerThread(4)); }, TasksPerThread(4));
auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load()); auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
@ -662,7 +654,9 @@ namespace netgen
for(auto [d_badness, ei] : edges_with_improvement) for(auto [d_badness, ei] : edges_with_improvement)
{ {
auto [pi1, pi2] = edges[ei]; auto [pi1, pi2] = edges[ei < 0 ? -ei : ei];
if(ei<0)
Swap(pi1,pi2);
CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, metricweight, false); CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, metricweight, false);
} }

1521
libsrc/meshing/improve3.cpp
View File

File diff suppressed because it is too large Load Diff

5
libsrc/meshing/improve3.hpp
View File

@ -21,10 +21,8 @@ public:
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 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void CombineImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY); void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SplitImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
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 (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);
void SplitImprove2 (Mesh & mesh); void SplitImprove2 (Mesh & mesh);
@ -34,12 +32,9 @@ public:
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 (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);
void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY, void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
const NgBitArray * working_elements = NULL); const NgBitArray * working_elements = NULL);
void SwapImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
const NgBitArray * working_elements = NULL);
void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY, void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
const NgBitArray * working_elements = NULL, const NgBitArray * working_elements = NULL,
const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL); const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
void SwapImprove2Sequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY); void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
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 ( Mesh & mesh, OPTIMIZEGOAL goal, ElementIndex eli1, int face, Table<ElementIndex, PointIndex> & elementsonnode, TABLE<SurfaceElementIndex, PointIndex::BASE> & belementsonnode, bool check_only=false );

1
libsrc/meshing/meshclass.cpp
View File

@ -6407,6 +6407,7 @@ namespace netgen
} }
Compress(); Compress();
SetNextMajorTimeStamp();
} }
void Mesh :: RebuildSurfaceElementLists () void Mesh :: RebuildSurfaceElementLists ()

1
libsrc/meshing/meshclass.hpp
View File

@ -566,7 +566,6 @@ namespace netgen
DLL_HEADER void DoArchive (Archive & archive); DLL_HEADER void DoArchive (Archive & archive);
/// ///
DLL_HEADER void ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY); DLL_HEADER void ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY);
DLL_HEADER void ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY);
/// ///
void ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY, const NgBitArray * usepoint = NULL); void ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY, const NgBitArray * usepoint = NULL);

9
libsrc/meshing/meshfunc.cpp
View File

@ -451,6 +451,7 @@ namespace netgen
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++)
{ {
mesh.FindOpenElements();
mesh.CalcSurfacesOfNode(); mesh.CalcSurfacesOfNode();
mesh.FreeOpenElementsEnvironment(2); mesh.FreeOpenElementsEnvironment(2);
mesh.CalcSurfacesOfNode(); mesh.CalcSurfacesOfNode();
@ -466,12 +467,10 @@ namespace netgen
} }
mesh.FindOpenElements(); mesh.FindOpenElements(domain);
PrintMessage (3, "Call remove problem"); PrintMessage (3, "Call remove problem");
// mesh.Save("before_remove.vol");
RemoveProblem (mesh, domain); RemoveProblem (mesh, domain);
// mesh.Save("after_remove.vol"); mesh.FindOpenElements(domain);
mesh.FindOpenElements();
} }
else else
{ {
@ -581,7 +580,7 @@ namespace netgen
FillCloseSurface( md[i] ); FillCloseSurface( md[i] );
CloseOpenQuads( md[i] ); CloseOpenQuads( md[i] );
MeshDomain(md[i]); MeshDomain(md[i]);
}); }, md.Size());
MergeMeshes(mesh3d, md); MergeMeshes(mesh3d, md);

11
libsrc/meshing/python_mesh.cpp
View File

@ -854,6 +854,17 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points), static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points),
py::return_value_policy::reference) py::return_value_policy::reference)
.def("Coordinates", [](Mesh & self) {
return py::array
(
py::memoryview::from_buffer
(&self.Points()[PointIndex::BASE](0), sizeof(double),
py::format_descriptor<double>::value,
{ self.Points().Size(), size_t(self.GetDimension()) },
{ sizeof(self.Points()[PointIndex::BASE]), sizeof(double) } )
);
})
.def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor), .def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
py::return_value_policy::reference) py::return_value_policy::reference)
.def("GetNFaceDescriptors", &Mesh::GetNFD) .def("GetNFaceDescriptors", &Mesh::GetNFD)

122
libsrc/meshing/smoothing3.cpp
View File

@ -1331,127 +1331,6 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL goal)
{
static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t);
(*testout) << "Improve Mesh" << "\n";
PrintMessage (3, "ImproveMesh");
int np = GetNP();
int ne = GetNE();
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
Vector x(3);
(*testout) << setprecision(8);
//int uselocalh = mparam.uselocalh;
PointFunction pf(points, volelements, mp);
Opti3FreeMinFunction freeminf(pf);
OptiParameters par;
par.maxit_linsearch = 20;
par.maxit_bfgs = 20;
NgArray<double, PointIndex::BASE> pointh (points.Size());
if(HasLocalHFunction())
{
for (PointIndex pi : points.Range())
pointh[pi] = GetH(pi);
}
else
{
pointh = 0;
for (Element & el : VolumeElements())
{
double h = pow(el.Volume(points),1./3.);
for (PointIndex pi : el.PNums())
if (h > pointh[pi])
pointh[pi] = h;
}
}
int printmod = 1;
char printdot = '.';
if (points.Size() > 1000)
{
printmod = 10;
printdot = '+';
}
if (points.Size() > 10000)
{
printmod = 100;
printdot = '*';
}
const char * savetask = multithread.task;
multithread.task = "Optimize Volume: Smooth Mesh";
for (PointIndex pi : points.Range())
if ( (*this)[pi].Type() == INNERPOINT )
{
if (multithread.terminate)
throw NgException ("Meshing stopped");
multithread.percent = 100.0 * (pi+1-PointIndex::BASE) / points.Size();
if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot);
double lh = pointh[pi];
pf.SetLocalH (lh);
par.typx = lh;
freeminf.SetPoint (points[pi]);
pf.SetPointIndex (pi);
x = 0;
int pok;
pok = freeminf.Func (x) < 1e10;
if (!pok)
{
pok = pf.MovePointToInner ();
freeminf.SetPoint (points[pi]);
pf.SetPointIndex (pi);
}
if (pok)
{
//*testout << "start BFGS, pok" << endl;
BFGS (x, freeminf, par);
//*testout << "BFGS complete, pok" << endl;
points[pi](0) += x(0);
points[pi](1) += x(1);
points[pi](2) += x(2);
}
}
PrintDot ('\n');
multithread.task = savetask;
if (goal == OPT_QUALITY)
{
double bad1 = CalcTotalBad (mp);
(*testout) << "Total badness = " << bad1 << endl;
PrintMessage (5, "Total badness = ", bad1);
}
}
void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
{ {
static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t); static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t);
@ -1461,7 +1340,6 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
static Timer trange("range"); static Timer trange("range");
static Timer tloch("loch"); static Timer tloch("loch");
// return ImproveMeshSequential(mp, goal);
BuildBoundaryEdges(false); BuildBoundaryEdges(false);
(*testout) << "Improve Mesh" << "\n"; (*testout) << "Improve Mesh" << "\n";

2
libsrc/visualization/CMakeLists.txt
View File

@ -14,6 +14,6 @@ install( TARGETS visual ${NG_INSTALL_DIR})
install(FILES install(FILES
meshdoc.hpp mvdraw.hpp meshdoc.hpp mvdraw.hpp
vispar.hpp visual.hpp vssolution.hpp vispar.hpp visual.hpp vssolution.hpp vsfieldlines.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/visualization COMPONENT netgen_devel DESTINATION ${NG_INSTALL_DIR_INCLUDE}/visualization COMPONENT netgen_devel
) )

176
libsrc/visualization/vsfieldlines.cpp
View File

@ -71,7 +71,7 @@ namespace netgen
K.SetSize(steps); K.SetSize(steps);
} }
void RKStepper :: StartNextValCalc(const Point3d & astartval, const double astartt, const double ah, const bool aadaptive) void RKStepper :: StartNextValCalc(const Point<3> & astartval, const double astartt, const double ah, const bool aadaptive)
{ {
//cout << "Starting RK-Step with h=" << ah << endl; //cout << "Starting RK-Step with h=" << ah << endl;
@ -83,12 +83,9 @@ namespace netgen
adrun = 0; adrun = 0;
} }
bool RKStepper :: GetNextData(Point3d & val, double & t, double & ah) bool RKStepper :: GetNextData(Point<3> & val, double & t, double & ah)
{ {
bool finished(false); bool finished = false;
//cout << "stepcount " << stepcount << endl;
if(stepcount <= steps) if(stepcount <= steps)
{ {
@ -125,9 +122,9 @@ namespace netgen
} }
else if (adrun == 2) else if (adrun == 2)
{ {
Point3d valh2 = val; Point<3> valh2 = val;
val = valh2 + 1./(pow(2.,order)-1.) * (valh2 - valh); val = valh2 + 1./(pow(2.,order)-1.) * (valh2 - valh);
Vec3d errvec = val - valh; auto errvec = val - valh;
double err = errvec.Length(); double err = errvec.Length();
@ -172,7 +169,7 @@ namespace netgen
} }
bool RKStepper :: FeedNextF(const Vec3d & f) bool RKStepper :: FeedNextF(const Vec<3> & f)
{ {
K[stepcount] = f; K[stepcount] = f;
stepcount++; stepcount++;
@ -181,19 +178,17 @@ namespace netgen
void FieldLineCalc :: GenerateFieldLines(NgArray<Point3d> & potential_startpoints, const int numlines, const int gllist, void FieldLineCalc :: GenerateFieldLines(Array<Point<3>> & potential_startpoints, const int numlines)
const double minval, const double maxval, const int logscale, double phaser, double phasei)
{ {
NgArray<Point3d> points; Array<Point<3>> line_points;
NgArray<double> values; Array<double> line_values;
NgArray<bool> drawelems; Array<bool> drawelems;
NgArray<int> dirstart; Array<int> dirstart;
pstart.SetSize0();
pend.SetSize0();
if(vsol -> iscomplex) values.SetSize0();
SetPhase(phaser,phasei);
double crit = 1.0; double crit = 1.0;
@ -201,8 +196,7 @@ namespace netgen
{ {
double sum = 0; double sum = 0;
double lami[3]; double lami[3];
double values[6]; Vec<3> v;
Vec3d v;
for(int i=0; i<potential_startpoints.Size(); i++) for(int i=0; i<potential_startpoints.Size(); i++)
{ {
@ -212,15 +206,8 @@ namespace netgen
mesh.SetPointSearchStartElement(elnr); mesh.SetPointSearchStartElement(elnr);
if (mesh.GetDimension()==3) func(elnr, lami, v);
vss.GetValues ( vsol, elnr, lami[0], lami[1], lami[2], values); sum += v.Length();
else
vss.GetSurfValues ( vsol, elnr, -1, lami[0], lami[1], values);
VisualSceneSolution::RealVec3d ( values, v, vsol->iscomplex, phaser, phasei);
sum += v.Length();
} }
crit = sum/double(numlines); crit = sum/double(numlines);
@ -232,8 +219,6 @@ namespace netgen
cout << endl; cout << endl;
for(int i=0; i<potential_startpoints.Size(); i++) for(int i=0; i<potential_startpoints.Size(); i++)
{ {
cout << "\rFieldline Calculation " << int(100.*i/potential_startpoints.Size()) << "%"; cout.flush(); cout << "\rFieldline Calculation " << int(100.*i/potential_startpoints.Size()) << "%"; cout.flush();
@ -243,7 +228,7 @@ namespace netgen
if(calculated >= numlines) break; if(calculated >= numlines) break;
Calc(potential_startpoints[i],points,values,drawelems,dirstart); Calc(potential_startpoints[i],line_points,line_values,drawelems,dirstart);
bool usable = false; bool usable = false;
@ -253,16 +238,9 @@ namespace netgen
if(!drawelems[k] || !drawelems[k+1]) continue; if(!drawelems[k] || !drawelems[k+1]) continue;
usable = true; usable = true;
pstart.Append(line_points[k]);
// vss.SetOpenGlColor (0.5*(values[k]+values[k+1]), minval, maxval, logscale); pend.Append(line_points[k+1]);
values.Append( 0.5*(line_values[k]+line_values[k+1]) );
/*
if (vss.usetexture == 1)
glTexCoord1f ( 0.5*(values[k]+values[k+1]) );
else
*/
vss.SetOpenGlColor (0.5*(values[k]+values[k+1]) );
vss.DrawCylinder (points[k], points[k+1], thickness);
} }
if(usable) calculated++; if(usable) calculated++;
@ -273,10 +251,10 @@ namespace netgen
FieldLineCalc :: FieldLineCalc(const Mesh & amesh, VisualSceneSolution & avss, const VisualSceneSolution::SolData * solution, FieldLineCalc :: FieldLineCalc(const Mesh & amesh, const VectorFunction & afunc,
const double rel_length, const int amaxpoints, const double rel_length, const int amaxpoints,
const double rel_thickness, const double rel_tolerance, const int rk_type, const int adirection) : const double rel_thickness, const double rel_tolerance, const int rk_type, const int adirection) :
mesh(amesh), vss(avss), vsol(solution), stepper(rk_type) mesh(amesh), func(afunc), stepper(rk_type)
{ {
mesh.GetBox (pmin, pmax); mesh.GetBox (pmin, pmax);
rad = 0.5 * Dist (pmin, pmax); rad = 0.5 * Dist (pmin, pmax);
@ -305,9 +283,6 @@ namespace netgen
} }
phaser = 1;
phasei = 0;
critical_value = -1; critical_value = -1;
randomized = false; randomized = false;
@ -317,24 +292,10 @@ namespace netgen
void FieldLineCalc :: Calc(const Point3d & startpoint, NgArray<Point3d> & points, NgArray<double> & vals, NgArray<bool> & drawelems, NgArray<int> & dirstart) void FieldLineCalc :: Calc(const Point<3> & startpoint, Array<Point<3>> & points, Array<double> & vals, Array<bool> & drawelems, Array<int> & dirstart)
{ {
double lami[3], startlami[3]; Vec<3> v = 0.0;
double values[6]; double startlami[3] = {0.0, 0.0, 0.0};
double dummyt(0);
Vec3d v;
Vec3d startv;
Point3d newp;
double h;
double startval;
bool startdraw;
bool drawelem = false;
int elnr;
for (int i=0; i<6; i++) values[i]=0.0;
for (int i=0; i<3; i++) lami[i]=0.0;
for (int i=0; i<3; i++) startlami[i]=0.0;
points.SetSize(0); points.SetSize(0);
vals.SetSize(0); vals.SetSize(0);
@ -351,14 +312,10 @@ namespace netgen
mesh.SetPointSearchStartElement(startelnr); mesh.SetPointSearchStartElement(startelnr);
if (mesh.GetDimension()==3) Vec<3> startv;
startdraw = vss.GetValues ( vsol, startelnr, startlami[0], startlami[1], startlami[2], values); bool startdraw = func(startelnr, startlami, startv);
else
startdraw = vss.GetSurfValues ( vsol, startelnr, -1, startlami[0], startlami[1], values);
VisualSceneSolution::RealVec3d ( values, startv, vsol->iscomplex, phaser, phasei); double startval = startv.Length();
startval = startv.Length();
if(critical_value > 0 && fabs(startval) < critical_value) if(critical_value > 0 && fabs(startval) < critical_value)
return; return;
@ -375,13 +332,13 @@ namespace netgen
vals.Append(startval); vals.Append(startval);
drawelems.Append(startdraw); drawelems.Append(startdraw);
h = 0.001*rad/startval; // otherwise no nice lines; should be made accessible from outside double h = 0.001*rad/startval; // otherwise no nice lines; should be made accessible from outside
v = startv; v = startv;
if(dir == -1) v *= -1.; if(dir == -1) v *= -1.;
elnr = startelnr; int elnr = startelnr;
lami[0] = startlami[0]; lami[1] = startlami[1]; lami[2] = startlami[2]; double lami[3] = { startlami[0], startlami[1], startlami[2]};
for(double length = 0; length < maxlength; length += h*vals.Last()) for(double length = 0; length < maxlength; length += h*vals.Last())
@ -392,21 +349,19 @@ namespace netgen
break; break;
} }
double dummyt;
stepper.StartNextValCalc(points.Last(),dummyt,h,true); stepper.StartNextValCalc(points.Last(),dummyt,h,true);
stepper.FeedNextF(v); stepper.FeedNextF(v);
bool drawelem = false;
Point<3> newp;
while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1) while(!stepper.GetNextData(newp,dummyt,h) && elnr != -1)
{ {
elnr = mesh.GetElementOfPoint(newp,lami,true) - 1; elnr = mesh.GetElementOfPoint(newp,lami,true) - 1;
if(elnr != -1) if(elnr != -1)
{ {
mesh.SetPointSearchStartElement(elnr); mesh.SetPointSearchStartElement(elnr);
if (mesh.GetDimension()==3) drawelem = func(elnr, lami, v);
drawelem = vss.GetValues (vsol, elnr, lami[0], lami[1], lami[2], values);
else
drawelem = vss.GetSurfValues (vsol, elnr, -1, lami[0], lami[1], values);
VisualSceneSolution::RealVec3d (values, v, vsol->iscomplex, phaser, phasei);
if(dir == -1) v *= -1.; if(dir == -1) v *= -1.;
stepper.FeedNextF(v); stepper.FeedNextF(v);
} }
@ -440,7 +395,7 @@ namespace netgen
void VisualSceneSolution :: BuildFieldLinesFromBox(NgArray<Point3d> & startpoints) void VisualSceneSolution :: BuildFieldLinesFromBox(Array<Point<3>> & startpoints)
{ {
shared_ptr<Mesh> mesh = GetMesh(); shared_ptr<Mesh> mesh = GetMesh();
if (!mesh) return; if (!mesh) return;
@ -462,7 +417,7 @@ namespace netgen
for (int i = 1; i <= startpoints.Size(); i++) for (int i = 1; i <= startpoints.Size(); i++)
{ {
Point3d p (fieldlines_startarea_parameter[0] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[3]-fieldlines_startarea_parameter[0]), Point<3> p (fieldlines_startarea_parameter[0] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[3]-fieldlines_startarea_parameter[0]),
fieldlines_startarea_parameter[1] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[4]-fieldlines_startarea_parameter[1]), fieldlines_startarea_parameter[1] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[4]-fieldlines_startarea_parameter[1]),
fieldlines_startarea_parameter[2] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[5]-fieldlines_startarea_parameter[2])); fieldlines_startarea_parameter[2] + double (rand()) / RAND_MAX * (fieldlines_startarea_parameter[5]-fieldlines_startarea_parameter[2]));
@ -470,7 +425,7 @@ namespace netgen
} }
} }
void VisualSceneSolution :: BuildFieldLinesFromLine(NgArray<Point3d> & startpoints) void VisualSceneSolution :: BuildFieldLinesFromLine(Array<Point<3>> & startpoints)
{ {
shared_ptr<Mesh> mesh = GetMesh(); shared_ptr<Mesh> mesh = GetMesh();
if (!mesh) return; if (!mesh) return;
@ -480,7 +435,7 @@ namespace netgen
{ {
double s = double (rand()) / RAND_MAX; double s = double (rand()) / RAND_MAX;
Point3d p (fieldlines_startarea_parameter[0] + s * (fieldlines_startarea_parameter[3]-fieldlines_startarea_parameter[0]), Point<3> p (fieldlines_startarea_parameter[0] + s * (fieldlines_startarea_parameter[3]-fieldlines_startarea_parameter[0]),
fieldlines_startarea_parameter[1] + s * (fieldlines_startarea_parameter[4]-fieldlines_startarea_parameter[1]), fieldlines_startarea_parameter[1] + s * (fieldlines_startarea_parameter[4]-fieldlines_startarea_parameter[1]),
fieldlines_startarea_parameter[2] + s * (fieldlines_startarea_parameter[5]-fieldlines_startarea_parameter[2])); fieldlines_startarea_parameter[2] + s * (fieldlines_startarea_parameter[5]-fieldlines_startarea_parameter[2]));
@ -489,7 +444,7 @@ namespace netgen
} }
void VisualSceneSolution :: BuildFieldLinesFromFile(NgArray<Point3d> & startpoints) void VisualSceneSolution :: BuildFieldLinesFromFile(Array<Point<3>> & startpoints)
{ {
shared_ptr<Mesh> mesh = GetMesh(); shared_ptr<Mesh> mesh = GetMesh();
if (!mesh) return; if (!mesh) return;
@ -538,7 +493,9 @@ namespace netgen
if (keyword == "point") if (keyword == "point")
{ {
(*infile) >> startpoints[numpoints].X(); (*infile) >> startpoints[numpoints].Y(); (*infile) >> startpoints[numpoints].Z(); (*infile) >> startpoints[numpoints][0];
(*infile) >> startpoints[numpoints][1];
(*infile) >> startpoints[numpoints][2];
numpoints++; numpoints++;
} }
else if (keyword == "line" || keyword == "box") else if (keyword == "line" || keyword == "box")
@ -546,7 +503,7 @@ namespace netgen
for(int i=0; i<6; i++) (*infile) >> fieldlines_startarea_parameter[i]; for(int i=0; i<6; i++) (*infile) >> fieldlines_startarea_parameter[i];
(*infile) >> iparam; (*infile) >> iparam;
NgArray<Point3d> auxpoints(iparam); Array<Point<3>> auxpoints(iparam);
if (keyword == "box") if (keyword == "box")
BuildFieldLinesFromBox(auxpoints); BuildFieldLinesFromBox(auxpoints);
@ -571,7 +528,7 @@ namespace netgen
} }
void VisualSceneSolution :: BuildFieldLinesFromFace(NgArray<Point3d> & startpoints) void VisualSceneSolution :: BuildFieldLinesFromFace(Array<Point<3>> & startpoints)
{ {
shared_ptr<Mesh> mesh = GetMesh(); shared_ptr<Mesh> mesh = GetMesh();
if (!mesh) return; if (!mesh) return;
@ -678,8 +635,25 @@ namespace netgen
num_fieldlineslists = (vsol -> iscomplex && !fieldlines_fixedphase) ? 100 : 1; num_fieldlineslists = (vsol -> iscomplex && !fieldlines_fixedphase) ? 100 : 1;
double phaser=1.0;
double phasei=0.0;
std::function eval_func = [&](int elnr, const double * lami, Vec<3> & vec)
{
double values[6] = {0., 0., 0., 0., 0., 0.};
bool drawelem;
auto mesh = GetMesh();
if (mesh->GetDimension()==3)
drawelem = GetValues (vsol, elnr, lami[0], lami[1], lami[2], values);
else
drawelem = GetSurfValues (vsol, elnr, -1, lami[0], lami[1], values);
FieldLineCalc linecalc(*mesh,*this,vsol, Vec3d v;
RealVec3d (values, v, vsol->iscomplex, phaser, phasei);
vec = v;
return drawelem;
};
FieldLineCalc linecalc(*mesh, eval_func,
fieldlines_rellength,fieldlines_maxpoints,fieldlines_relthickness,fieldlines_reltolerance,fieldlines_rktype); fieldlines_rellength,fieldlines_maxpoints,fieldlines_relthickness,fieldlines_reltolerance,fieldlines_rktype);
if(fieldlines_randomstart) if(fieldlines_randomstart)
@ -694,7 +668,7 @@ namespace netgen
num_startpoints *= 10; num_startpoints *= 10;
NgArray<Point3d> startpoints(num_startpoints); Array<Point<3>> startpoints(num_startpoints);
for (int ln = 0; ln < num_fieldlineslists; ln++) for (int ln = 0; ln < num_fieldlineslists; ln++)
@ -722,17 +696,27 @@ namespace netgen
cout << "phi = " << phi << endl; cout << "phi = " << phi << endl;
double phaser = cos(phi), phasei = sin(phi); phaser = cos(phi);
phasei = sin(phi);
linecalc.GenerateFieldLines(startpoints,num_fieldlines / num_fieldlineslists+1);
auto & pstart = linecalc.GetPStart();
auto & pend = linecalc.GetPEnd();
auto & values = linecalc.GetValues();
auto nlines = values.Size();
glNewList(fieldlineslist+ln, GL_COMPILE); glNewList(fieldlineslist+ln, GL_COMPILE);
SetTextureMode (usetexture); SetTextureMode (usetexture);
linecalc.GenerateFieldLines(startpoints,num_fieldlines / num_fieldlineslists+1,
fieldlineslist+ln,minval,maxval,logscale,phaser,phasei); for(auto i : Range(nlines))
{
SetOpenGlColor (values[i]);
DrawCylinder (pstart[i], pend[i], fieldlines_relthickness);
}
glEndList (); glEndList ();
} }
} }

101
libsrc/visualization/vsfieldlines.hpp Normal file
View File

@ -0,0 +1,101 @@
#ifndef VSFIELDLINES_HPP_INCLUDED
#define VSFIELDLINES_HPP_INCLUDED
namespace netgen
{
class RKStepper
{
private:
Array<double> c,b;
TABLE<double> *a;
int steps;
int order;
double tolerance;
Array<Vec<3>> K;
int stepcount;
double h;
double startt;
double startt_bak;
Point<3> startval;
Point<3> startval_bak;
bool adaptive;
int adrun;
Point<3> valh;
int notrestarted;
public:
~RKStepper();
RKStepper(int type = 0);
void SetTolerance(const double tol){tolerance = tol;}
void StartNextValCalc(const Point<3> & astartval, const double astartt, const double ah, const bool aadaptive = false);
bool GetNextData(Point<3> & val, double & t, double & ah);
bool FeedNextF(const Vec<3> & f);
};
class FieldLineCalc
{
private:
const Mesh & mesh;
typedef std::function<bool (int elnr, const double *, Vec<3> &)> VectorFunction;
const VectorFunction & func;
RKStepper stepper;
Array<double> values;
Array<Point<3>> pstart, pend;
double maxlength;
int maxpoints;
int direction;
Point3d pmin, pmax;
double rad;
double critical_value;
bool randomized;
double thickness;
public:
FieldLineCalc(const Mesh & amesh, const VectorFunction & afunc,
const double rel_length, const int amaxpoints = -1,
const double rel_thickness = -1, const double rel_tolerance = -1, const int rk_type = 0, const int adirection = 0);
void SetCriticalValue(const double val) { critical_value = val; }
void Randomized(void) { randomized = true; }
void NotRandomized(void) { randomized = false; }
void Calc(const Point<3> & startpoint, Array<Point<3>> & points, Array<double> & vals, Array<bool> & drawelems, Array<int> & dirstart);
void GenerateFieldLines(Array<Point<3>> & potential_startpoints, const int numlines);
const auto & GetPStart() const { return pstart; }
const auto & GetPEnd() const { return pend; }
const auto & GetValues() const { return values; }
const auto GetThickness() const { return thickness; }
};
} // namespace netgen
#endif // VSFIELDLINES_HPP_INCLUDED

100
libsrc/visualization/vssolution.hpp
View File

@ -1,6 +1,7 @@
#ifndef FILE_VSSOLUTION #ifndef FILE_VSSOLUTION
#define FILE_VSSOLUTION #define FILE_VSSOLUTION
#include "vsfieldlines.hpp"
typedef void * ClientData; typedef void * ClientData;
struct Tcl_Interp; struct Tcl_Interp;
@ -12,8 +13,6 @@ namespace netgen
DLL_HEADER extern void ImportSolution (const char * filename); DLL_HEADER extern void ImportSolution (const char * filename);
class FieldLineCalc;
extern int Ng_Vis_Set (ClientData clientData, extern int Ng_Vis_Set (ClientData clientData,
Tcl_Interp * interp, Tcl_Interp * interp,
int argc, const char *argv[]); int argc, const char *argv[]);
@ -183,10 +182,10 @@ public:
bool imag_part; bool imag_part;
private: private:
void BuildFieldLinesFromFile(NgArray<Point3d> & startpoints); void BuildFieldLinesFromFile(Array<Point<3>> & startpoints);
void BuildFieldLinesFromFace(NgArray<Point3d> & startpoints); void BuildFieldLinesFromFace(Array<Point<3>> & startpoints);
void BuildFieldLinesFromBox(NgArray<Point3d> & startpoints); void BuildFieldLinesFromBox(Array<Point<3>> & startpoints);
void BuildFieldLinesFromLine(NgArray<Point3d> & startpoints); void BuildFieldLinesFromLine(Array<Point<3>> & startpoints);
// weak_ptr<Mesh> wp_mesh; // weak_ptr<Mesh> wp_mesh;
public: public:
VisualSceneSolution (); VisualSceneSolution ();
@ -359,95 +358,6 @@ public:
class RKStepper
{
private:
NgArray<double> c,b;
TABLE<double> *a;
int steps;
int order;
double tolerance;
NgArray<Vec3d> K;
int stepcount;
double h;
double startt;
double startt_bak;
Point3d startval;
Point3d startval_bak;
bool adaptive;
int adrun;
Point3d valh;
int notrestarted;
public:
~RKStepper();
RKStepper(int type = 0);
void SetTolerance(const double tol){tolerance = tol;}
void StartNextValCalc(const Point3d & astartval, const double astartt, const double ah, const bool aadaptive = false);
bool GetNextData(Point3d & val, double & t, double & ah);
bool FeedNextF(const Vec3d & f);
};
class FieldLineCalc
{
private:
const Mesh & mesh;
VisualSceneSolution & vss;
const VisualSceneSolution::SolData * vsol;
RKStepper stepper;
double maxlength;
int maxpoints;
int direction;
Point3d pmin, pmax;
double rad;
double phaser, phasei;
double critical_value;
bool randomized;
double thickness;
public:
FieldLineCalc(const Mesh & amesh, VisualSceneSolution & avss, const VisualSceneSolution::SolData * solution,
const double rel_length, const int amaxpoints = -1,
const double rel_thickness = -1, const double rel_tolerance = -1, const int rk_type = 0, const int adirection = 0);
void SetPhase(const double real, const double imag) { phaser = real; phasei = imag; }
void SetCriticalValue(const double val) { critical_value = val; }
void Randomized(void) { randomized = true; }
void NotRandomized(void) { randomized = false; }
void Calc(const Point3d & startpoint, NgArray<Point3d> & points, NgArray<double> & vals, NgArray<bool> & drawelems, NgArray<int> & dirstart);
void GenerateFieldLines(NgArray<Point3d> & potential_startpoints, const int numlines, const int gllist,
const double minval, const double maxval, const int logscale, double phaser, double phasei);
};

24
tests/pytest/compare_results.py
View File

@ -66,19 +66,27 @@ data2 = readData(s2, filenames)
assert(len(data) == len(data2)) assert(len(data) == len(data2))
w = 90
GREEN = '\033[92m'
RED = '\033[91m'
RESET = '\033[0m'
for bad1,bad2, f1, f2 in zip(data['badness'], data2['badness'], data['file'], data2['file']): for bad1,bad2, f1, f2 in zip(data['badness'], data2['badness'], data['file'], data2['file']):
assert f1==f2 assert f1==f2
if bad2>0 and bad2>1.1*bad1:
print(f"file {f1} got worse: {bad1} -> {bad2}") diff = f"{100*(bad2-bad1)/bad1:+.2f}%"
if bad2>0 and bad2<0.9*bad1: if bad2>0 and bad2>1.2*bad1:
print(f"file {f1} got better: {bad1} -> {bad2}") print(f"{RED}badness {f1} got worse: {bad1} -> {bad2}".ljust(w) + diff + RESET)
if bad2>0 and bad2<0.8*bad1:
print(f"{GREEN}badness {f1} got better: {bad1} -> {bad2}".ljust(w) + diff + RESET)
for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data2['file']): for bad1,bad2, f1, f2 in zip(data['#trigs'], data2['#trigs'], data['file'], data2['file']):
assert f1==f2 assert f1==f2
if bad2>0 and bad2>1.1*bad1: diff = f"{100*(bad2-bad1)/bad1:+.2f}%"
print(f"file {f1} got worse: {bad1} -> {bad2}") if bad2>0 and bad2>1.2*bad1:
if bad2>0 and bad2<0.9*bad1: print(f"{RED}ntrigs {f1} got worse: {bad1} -> {bad2}".ljust(w) + diff + RESET)
print(f"file {f1} got better: {bad1} -> {bad2}") if bad2>0 and bad2<0.8*bad1:
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=(10,7))

2524
tests/pytest/results.json
View File

File diff suppressed because it is too large Load Diff