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

AVX512 for element-trafos

Browse Source
This commit is contained in:
Joachim Schöberl 2017-11-23 21:26:23 +01:00
parent b1cb2883e7
commit 7cf05d84a8
9 changed files with 212 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
libsrc/general/ngsimd.hpp
View File

@ -58,6 +58,16 @@ namespace ngsimd
#endif #endif
} }
#if defined __AVX512F__
typedef __m512 tAVX;
typedef __m512d tAVXd;
#elif defined __AVX__
typedef __m256 tAVX;
typedef __m256d tAVXd;
#endif
template <typename T, int N=GetDefaultSIMDSize()> class SIMD; template <typename T, int N=GetDefaultSIMDSize()> class SIMD;
template <typename T> template <typename T>

2
libsrc/general/profiler.hpp
View File

@ -18,7 +18,7 @@
#endif #endif
// #define USE_TSC #define USE_TSC
#ifdef USE_TSC #ifdef USE_TSC
#include <x86intrin.h> // for __rdtsc() CPU time step counter #include <x86intrin.h> // for __rdtsc() CPU time step counter
#endif #endif

44
libsrc/gprim/geomops.hpp
View File

@ -18,10 +18,10 @@ namespace netgen
*/ */
template <int D> template <int D, typename T>
inline Vec<D> operator+ (const Vec<D> & a, const Vec<D> & b) inline Vec<D,T> operator+ (Vec<D,T> a, Vec<D,T> b)
{ {
Vec<D> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) + b(i); res(i) = a(i) + b(i);
return res; return res;
@ -30,7 +30,7 @@ namespace netgen
template <int D, typename T> template <int D, typename T>
inline Point<D,T> operator+ (const Point<D,T> & a, const Vec<D,T> & b) inline Point<D,T> operator+ (Point<D,T> a, Vec<D,T> b)
{ {
Point<D,T> res; Point<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
@ -41,7 +41,7 @@ namespace netgen
template <int D, typename T> template <int D, typename T>
inline Vec<D,T> operator- (const Point<D,T> & a, const Point<D,T> & b) inline Vec<D,T> operator- (Point<D,T> a, Point<D,T> b)
{ {
Vec<D,T> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
@ -49,19 +49,19 @@ namespace netgen
return res; return res;
} }
template <int D> template <int D, typename T>
inline Point<D> operator- (const Point<D> & a, const Vec<D> & b) inline Point<D,T> operator- (Point<D,T> a, Vec<D,T> b)
{ {
Point<D> res; Point<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
} }
template <int D> template <int D, typename T>
inline Vec<D> operator- (const Vec<D> & a, const Vec<D> & b) inline Vec<D,T> operator- (Vec<D,T> a, Vec<D,T> b)
{ {
Vec<D> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
@ -70,7 +70,7 @@ namespace netgen
template <int D, typename T> template <int D, typename T>
inline Vec<D,T> operator* (T s, const Vec<D,T> & b) inline Vec<D,T> operator* (T s, Vec<D,T> b)
{ {
Vec<D,T> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
@ -80,7 +80,7 @@ namespace netgen
template <int D> template <int D>
inline double operator* (const Vec<D> & a, const Vec<D> & b) inline double operator* (Vec<D> a, Vec<D> b)
{ {
double sum = 0; double sum = 0;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
@ -90,26 +90,26 @@ namespace netgen
template <int D> template <int D, typename T>
inline Vec<D> operator- (const Vec<D> & b) inline Vec<D,T> operator- (Vec<D,T> b)
{ {
Vec<D> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = -b(i); res(i) = -b(i);
return res; return res;
} }
template <int D> template <int D, typename T>
inline Point<D> & operator+= (Point<D> & a, const Vec<D> & b) inline Point<D,T> & operator+= (Point<D,T> & a, Vec<D,T> b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) += b(i); a(i) += b(i);
return a; return a;
} }
template <int D> template <int D, typename T>
inline Vec<D> & operator+= (Vec<D> & a, const Vec<D> & b) inline Vec<D,T> & operator+= (Vec<D,T> & a, Vec<D> b)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) += b(i); a(i) += b(i);
@ -135,8 +135,8 @@ namespace netgen
template <int D> template <int D, typename T1, typename T2>
inline Vec<D> & operator*= (Vec<D> & a, double s) inline Vec<D,T1> & operator*= (Vec<D,T1> & a, T2 s)
{ {
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
a(i) *= s; a(i) *= s;

69
libsrc/interface/nginterface_v2.cpp
View File

@ -676,18 +676,18 @@ namespace netgen
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,1> (int elnr, int npts, MultiElementTransformation<1,1> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
cout << "multi-eltrafo simd called, 1,1,simd" << endl; cout << "multi-eltrafo simd called, 1,1,simd" << endl;
} }
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,2> (int elnr, int npts, MultiElementTransformation<2,2> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2>
(elnr, npts, (elnr, npts,
@ -720,9 +720,9 @@ namespace netgen
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<3,3> (int elnr, int npts, MultiElementTransformation<3,3> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
mesh->GetCurvedElements().CalcMultiPointElementTransformation mesh->GetCurvedElements().CalcMultiPointElementTransformation
(elnr, npts, (elnr, npts,
@ -755,28 +755,34 @@ namespace netgen
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,2> (int elnr, int npts, MultiElementTransformation<0,2> (int elnr, int npts,
const __m256d *xi, size_t sxi, const tAVXd *xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
cout << "MultiElementtransformation<0,2> simd not implemented" << endl; cout << "MultiElementtransformation<0,2> simd not implemented" << endl;
} }
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,1> (int elnr, int npts, MultiElementTransformation<0,1> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
cout << "multi-eltrafo simd called, 0,1,simd" << endl; cout << "multi-eltrafo simd called, 0,1,simd" << endl;
} }
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,3> (int elnr, int npts, MultiElementTransformation<1,3> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
/*
double hxi[4][1]; double hxi[4][1];
double hx[4][3]; double hx[4][3];
double hdxdxi[4][3]; double hdxdxi[4][3];
@ -793,14 +799,21 @@ namespace netgen
xi += sxi; xi += sxi;
x += sx; x += sx;
dxdxi += sdxdxi; dxdxi += sdxdxi;
*/
} }
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<1,2> (int elnr, int npts, MultiElementTransformation<1,2> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2>
(elnr, npts,
reinterpret_cast<const SIMD<double>*> (xi), sxi,
reinterpret_cast<SIMD<double>*> (x), sx,
reinterpret_cast<SIMD<double>*> (dxdxi), sdxdxi);
/*
for (int i = 0; i < npts; i++) for (int i = 0; i < npts; i++)
{ {
double hxi[4][1]; double hxi[4][1];
@ -821,14 +834,14 @@ namespace netgen
x += sx; x += sx;
dxdxi += sdxdxi; dxdxi += sdxdxi;
} }
*/
} }
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<2,3> (int elnr, int npts, MultiElementTransformation<2,3> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3>
(elnr, npts, (elnr, npts,
@ -861,9 +874,9 @@ namespace netgen
template<> DLL_HEADER void Ngx_Mesh :: template<> DLL_HEADER void Ngx_Mesh ::
MultiElementTransformation<0,3> (int elnr, int npts, MultiElementTransformation<0,3> (int elnr, int npts,
const __m256d * xi, size_t sxi, const tAVXd * xi, size_t sxi,
__m256d * x, size_t sx, tAVXd * x, size_t sx,
__m256d * dxdxi, size_t sdxdxi) const tAVXd * dxdxi, size_t sdxdxi) const
{ {
for (int i = 0; i < npts; i++) for (int i = 0; i < npts; i++)
{ {

90
libsrc/meshing/curvedelems.cpp
View File

@ -71,11 +71,12 @@ namespace netgen
} }
} }
static void CalcEdgeDx (int n, double x, double * dshape) template <typename T>
static void CalcEdgeDx (int n, T x, T * dshape)
{ {
double p1 = x, p2 = -1, p3 = 0; T p1 = x, p2 = -1, p3 = 0;
double p1dx = 1, p2dx = 0, p3dx = 0; T p1dx = 1, p2dx = 0, p3dx = 0;
for (int j=2; j<=n; j++) for (int j=2; j<=n; j++)
{ {
@ -1408,10 +1409,10 @@ namespace netgen
template <typename T>
void CurvedElements :: void CurvedElements ::
CalcSegmentTransformation (double xi, SegmentIndex elnr, CalcSegmentTransformation (T xi, SegmentIndex elnr,
Point<3> * x, Vec<3> * dxdxi, bool * curved) Point<3,T> * x, Vec<3,T> * dxdxi, bool * curved)
{ {
if (mesh.coarsemesh) if (mesh.coarsemesh)
{ {
@ -1419,11 +1420,11 @@ namespace netgen
(*mesh.hpelements) [mesh[elnr].hp_elnr]; (*mesh.hpelements) [mesh[elnr].hp_elnr];
// xi umrechnen // xi umrechnen
double lami[2] = { xi, 1-xi }; T lami[2] = { xi, 1-xi };
double dlami[2] = { 1, -1 }; T dlami[2] = { 1, -1 };
double coarse_xi = 0; T coarse_xi = 0;
double trans = 0; T trans = 0;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
{ {
coarse_xi += hpref_el.param[i][0] * lami[i]; coarse_xi += hpref_el.param[i][0] * lami[i];
@ -1437,8 +1438,9 @@ namespace netgen
} }
Vector shapes, dshapes;
Array<Vec<3> > coefs; // TVector<T> shapes, dshapes;
// Array<Vec<3> > coefs;
SegmentInfo info; SegmentInfo info;
info.elnr = elnr; info.elnr = elnr;
@ -1452,12 +1454,21 @@ namespace netgen
info.ndof += edgeorder[info.edgenr]-1; info.ndof += edgeorder[info.edgenr]-1;
} }
ArrayMem<Vec<3>,100> coefs(info.ndof);
ArrayMem<T, 100> shapes_mem(info.ndof);
TFlatVector<T> shapes(info.ndof, &shapes_mem[0]);
ArrayMem<T, 200> dshapes_mem(info.ndof);
TFlatVector<T> dshapes(info.ndof, &dshapes_mem[0]);
CalcElementShapes (info, xi, shapes); CalcElementShapes (info, xi, shapes);
GetCoefficients (info, coefs); GetCoefficients (info, coefs);
*x = 0; *x = 0;
for (int i = 0; i < shapes.Size(); i++) for (int i = 0; i < shapes.Size(); i++)
*x += shapes(i) * coefs[i]; // *x += shapes(i) * coefs[i];
for (int j = 0; j < 3; j++)
(*x)(j) += shapes(i) * coefs[i](j);
if (dxdxi) if (dxdxi)
@ -1477,10 +1488,11 @@ namespace netgen
} }
template <typename T>
void CurvedElements :: void CurvedElements ::
CalcElementShapes (SegmentInfo & info, double xi, Vector & shapes) const CalcElementShapes (SegmentInfo & info, T xi, TFlatVector<T> shapes) const
{ {
/*
if (rational && info.order == 2) if (rational && info.order == 2)
{ {
shapes.SetSize(3); shapes.SetSize(3);
@ -1491,9 +1503,9 @@ namespace netgen
shapes *= 1.0 / (1 + (w-1) *2*xi*(1-xi)); shapes *= 1.0 / (1 + (w-1) *2*xi*(1-xi));
return; return;
} }
*/
// shapes.SetSize(info.ndof);
shapes.SetSize(info.ndof);
shapes(0) = xi; shapes(0) = xi;
shapes(1) = 1-xi; shapes(1) = 1-xi;
@ -1505,9 +1517,11 @@ namespace netgen
} }
} }
template <typename T>
void CurvedElements :: void CurvedElements ::
CalcElementDShapes (SegmentInfo & info, double xi, Vector & dshapes) const CalcElementDShapes (SegmentInfo & info, T xi, TFlatVector<T> dshapes) const
{ {
/*
if (rational && info.order == 2) if (rational && info.order == 2)
{ {
dshapes.SetSize(3); dshapes.SetSize(3);
@ -1527,13 +1541,11 @@ namespace netgen
dshapes(j) = dshapes(j) / w - shapes[j] * dw / (w*w); dshapes(j) = dshapes(j) / w - shapes[j] * dw / (w*w);
return; return;
} }
*/
// dshapes.SetSize(info.ndof);
dshapes.SetSize(info.ndof);
dshapes = 0; dshapes = 0;
dshapes(0) = 1; dshapes(0) = 1;
dshapes(1) = -1; dshapes(1) = -1;
@ -1542,7 +1554,7 @@ namespace netgen
if (info.order >= 2) if (info.order >= 2)
{ {
double fac = 2; T fac = 2;
if (mesh[info.elnr][0] > mesh[info.elnr][1]) if (mesh[info.elnr][0] > mesh[info.elnr][1])
{ {
xi = 1-xi; xi = 1-xi;
@ -3651,7 +3663,7 @@ namespace netgen
/*
void CurvedElements :: void CurvedElements ::
CalcMultiPointSegmentTransformation (Array<double> * xi, SegmentIndex segnr, CalcMultiPointSegmentTransformation (Array<double> * xi, SegmentIndex segnr,
Array<Point<3> > * x, Array<Point<3> > * x,
@ -3659,22 +3671,22 @@ namespace netgen
{ {
; ;
} }
*/
template <int DIM_SPACE, typename T>
template <int DIM_SPACE>
void CurvedElements :: void CurvedElements ::
CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n, CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n,
const double * xi, size_t sxi, const T * xi, size_t sxi,
double * x, size_t sx, T * x, size_t sx,
double * dxdxi, size_t sdxdxi) T * dxdxi, size_t sdxdxi)
{ {
for (int ip = 0; ip < n; ip++) for (int ip = 0; ip < n; ip++)
{ {
Point<3> xg; Point<3,T> xg;
Vec<3> dx; Vec<3,T> dx;
// mesh->GetCurvedElements(). // mesh->GetCurvedElements().
CalcSegmentTransformation (xi[ip*sxi], elnr, xg, dx); CalcSegmentTransformation<T> (xi[ip*sxi], elnr, &xg, &dx);
if (x) if (x)
for (int i = 0; i < DIM_SPACE; i++) for (int i = 0; i < DIM_SPACE; i++)
@ -3699,6 +3711,18 @@ namespace netgen
double * x, size_t sx, double * x, size_t sx,
double * dxdxi, size_t sdxdxi); double * dxdxi, size_t sdxdxi);
template void CurvedElements ::
CalcMultiPointSegmentTransformation<2> (SegmentIndex elnr, int npts,
const SIMD<double> * xi, size_t sxi,
SIMD<double> * x, size_t sx,
SIMD<double> * dxdxi, size_t sdxdxi);
template void CurvedElements ::
CalcMultiPointSegmentTransformation<3> (SegmentIndex elnr, int npts,
const SIMD<double> * xi, size_t sxi,
SIMD<double> * x, size_t sx,
SIMD<double> * dxdxi, size_t sdxdxi);
void CurvedElements :: void CurvedElements ::

32
libsrc/meshing/curvedelems.hpp
View File

@ -57,15 +57,15 @@ public:
void CalcSegmentTransformation (double xi, SegmentIndex segnr, void CalcSegmentTransformation (double xi, SegmentIndex segnr,
Point<3> & x) Point<3> & x)
{ CalcSegmentTransformation (xi, segnr, &x, NULL); }; { CalcSegmentTransformation<double> (xi, segnr, &x, NULL); };
void CalcSegmentTransformation (double xi, SegmentIndex segnr, void CalcSegmentTransformation (double xi, SegmentIndex segnr,
Vec<3> & dxdxi) Vec<3> & dxdxi)
{ CalcSegmentTransformation (xi, segnr, NULL, &dxdxi); }; { CalcSegmentTransformation<double> (xi, segnr, NULL, &dxdxi); };
void CalcSegmentTransformation (double xi, SegmentIndex segnr, void CalcSegmentTransformation (double xi, SegmentIndex segnr,
Point<3> & x, Vec<3> & dxdxi) Point<3> & x, Vec<3> & dxdxi)
{ CalcSegmentTransformation (xi, segnr, &x, &dxdxi, NULL); }; { CalcSegmentTransformation<double> (xi, segnr, &x, &dxdxi, NULL); };
void CalcSegmentTransformation (double xi, SegmentIndex segnr, void CalcSegmentTransformation (double xi, SegmentIndex segnr,
Point<3> & x, Vec<3> & dxdxi, bool & curved) Point<3> & x, Vec<3> & dxdxi, bool & curved)
@ -115,16 +115,17 @@ public:
// { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , &curved * ); } // { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , &curved * ); }
/*
void CalcMultiPointSegmentTransformation (Array<double> * xi, SegmentIndex segnr, void CalcMultiPointSegmentTransformation (Array<double> * xi, SegmentIndex segnr,
Array<Point<3> > * x, Array<Point<3> > * x,
Array<Vec<3> > * dxdxi); Array<Vec<3> > * dxdxi);
*/
template <int DIM_SPACE>
template <int DIM_SPACE, typename T>
void CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n, void CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n,
const double * xi, size_t sxi, const T * xi, size_t sxi,
double * x, size_t sx, T * x, size_t sx,
double * dxdxi, size_t sdxdxi); T * dxdxi, size_t sdxdxi);
void CalcMultiPointSurfaceTransformation (Array< Point<2> > * xi, SurfaceElementIndex elnr, void CalcMultiPointSurfaceTransformation (Array< Point<2> > * xi, SurfaceElementIndex elnr,
Array< Point<3> > * x, Array< Point<3> > * x,
@ -150,9 +151,10 @@ public:
private: private:
void CalcSegmentTransformation (double xi, SegmentIndex segnr, template <typename T>
Point<3> * x = NULL, Vec<3> * dxdxi = NULL, bool * curved = NULL); void CalcSegmentTransformation (T xi, SegmentIndex segnr,
Point<3,T> * x = NULL, Vec<3,T> * dxdxi = NULL, bool * curved = NULL);
void CalcSurfaceTransformation (Point<2> xi, SurfaceElementIndex elnr, void CalcSurfaceTransformation (Point<2> xi, SurfaceElementIndex elnr,
Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL, bool * curved = NULL); Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL, bool * curved = NULL);
@ -176,9 +178,11 @@ private:
int edgenr; int edgenr;
}; };
void CalcElementShapes (SegmentInfo & elnr, double xi, Vector & shapes) const; template <typename T>
void CalcElementShapes (SegmentInfo & elnr, T xi, TFlatVector<T> shapes) const;
void GetCoefficients (SegmentInfo & elnr, Array<Vec<3> > & coefs) const; void GetCoefficients (SegmentInfo & elnr, Array<Vec<3> > & coefs) const;
void CalcElementDShapes (SegmentInfo & elnr, double xi, Vector & dshapes) const; template <typename T>
void CalcElementDShapes (SegmentInfo & elnr, T xi, TFlatVector<T> dshapes) const;
class ElementInfo class ElementInfo

65
libsrc/stlgeom/stlgeomchart.cpp
View File

@ -20,16 +20,24 @@ int chartdebug = 0;
void STLGeometry :: MakeAtlas(Mesh & mesh) void STLGeometry :: MakeAtlas(Mesh & mesh)
{ {
int timer1 = NgProfiler::CreateTimer ("makeatlas"); int timer1 = NgProfiler::CreateTimer ("makeatlas");
int timerb = NgProfiler::CreateTimer ("makeatlas - begin");
int timere = NgProfiler::CreateTimer ("makeatlas - end");
int timere1 = NgProfiler::CreateTimer ("makeatlas - end1");
int timere2 = NgProfiler::CreateTimer ("makeatlas - end2");
int timer2 = NgProfiler::CreateTimer ("makeatlas - part 2"); int timer2 = NgProfiler::CreateTimer ("makeatlas - part 2");
int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3"); int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3");
/*
int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4"); int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4");
int timer4a = NgProfiler::CreateTimer ("makeatlas - part 4a");
int timer4b = NgProfiler::CreateTimer ("makeatlas - part 4b");
int timer4c = NgProfiler::CreateTimer ("makeatlas - part 4c");
int timer4d = NgProfiler::CreateTimer ("makeatlas - part 4d");
int timer4e = NgProfiler::CreateTimer ("makeatlas - part 4e");
int timer5 = NgProfiler::CreateTimer ("makeatlas - part 5"); int timer5 = NgProfiler::CreateTimer ("makeatlas - part 5");
int timer5a = NgProfiler::CreateTimer ("makeatlas - part 5a"); int timer5a = NgProfiler::CreateTimer ("makeatlas - part 5a");
int timer5b = NgProfiler::CreateTimer ("makeatlas - part 5b"); int timer5b = NgProfiler::CreateTimer ("makeatlas - part 5b");
int timer5cs = NgProfiler::CreateTimer ("makeatlas - part 5cs"); int timer5cs = NgProfiler::CreateTimer ("makeatlas - part 5cs");
int timer5cl = NgProfiler::CreateTimer ("makeatlas - part 5cl"); int timer5cl = NgProfiler::CreateTimer ("makeatlas - part 5cl");
*/
PushStatusF("Make Atlas"); PushStatusF("Make Atlas");
double h = mparam.maxh; double h = mparam.maxh;
@ -110,9 +118,10 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
int markedtrigcnt = 0; int markedtrigcnt = 0;
while (markedtrigcnt < GetNT()) while (markedtrigcnt < GetNT())
{ {
if (multithread.terminate) { PopStatus(); return; } if (multithread.terminate) { PopStatus(); return; }
NgProfiler::StartTimer (timerb);
if (workedarea / atlasarea*100. >= nextshow) if (workedarea / atlasarea*100. >= nextshow)
{PrintDot(); nextshow+=showinc;} {PrintDot(); nextshow+=showinc;}
@ -177,7 +186,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
bool changed = true; bool changed = true;
int oldstartic = 1; int oldstartic = 1;
int oldstartic2; int oldstartic2;
NgProfiler::StopTimer (timerb);
NgProfiler::StartTimer (timer2); NgProfiler::StartTimer (timer2);
@ -342,8 +351,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
{ {
accepted = 1; accepted = 1;
// NgProfiler::StartTimer (timer4); NgProfiler::StartTimer (timer4);
bool isdirtytrig = false; bool isdirtytrig = false;
Vec<3> gn = GetTriangle(nt).GeomNormal(points); Vec<3> gn = GetTriangle(nt).GeomNormal(points);
double gnlen = gn.Length(); double gnlen = gn.Length();
@ -355,36 +363,45 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
//find overlapping charts exacter: //find overlapping charts exacter:
//do not check dirty trigs! //do not check dirty trigs!
NgProfiler::StartTimer (timer4a);
if (spiralcheckon && !isdirtytrig) if (spiralcheckon && !isdirtytrig)
for (int k = 1; k <= 3; k++) for (int k = 1; k <= 3; k++)
{ {
// NgProfiler::StartTimer (timer4b);
int nnt = NeighbourTrig(nt,k); int nnt = NeighbourTrig(nt,k);
if (outermark.Elem(nnt) != chartnum) if (outermark.Elem(nnt) != chartnum)
{ {
NgProfiler::StartTimer (timer4c);
int nnp1, nnp2; int nnp1, nnp2;
GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2);
NgProfiler::StopTimer (timer4c);
// NgProfiler::StartTimer (timer4a); NgProfiler::StartTimer (timer4d);
accepted = accepted =
chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2), chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2),
sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps); sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps);
// NgProfiler::StopTimer (timer4a); NgProfiler::StopTimer (timer4d);
NgProfiler::StartTimer (timer4e);
Vec<3> n3 = GetTriangle(nnt).Normal(); Vec<3> n3 = GetTriangle(nnt).Normal();
if ( (n3 * sn) >= cosouterchartangle && if ( (n3 * sn) >= cosouterchartangle &&
IsSmoothEdge (nnp1, nnp2) ) IsSmoothEdge (nnp1, nnp2) )
accepted = 1; accepted = 1;
NgProfiler::StopTimer (timer4e);
} }
// NgProfiler::StopTimer (timer4b);
if (!accepted) break; if (!accepted) break;
} }
NgProfiler::StopTimer (timer4a);
// NgProfiler::StopTimer (timer4); NgProfiler::StopTimer (timer4);
// NgProfiler::RegionTimer reg5(timer5); NgProfiler::RegionTimer reg5(timer5);
// outer chart is only small environment of // outer chart is only small environment of
@ -392,7 +409,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
if (accepted) if (accepted)
{ {
// NgProfiler::StartTimer (timer5a); NgProfiler::StartTimer (timer5a);
accepted = 0; accepted = 0;
for (int k = 1; k <= 3; k++) for (int k = 1; k <= 3; k++)
@ -402,9 +419,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
break; break;
} }
// NgProfiler::StopTimer (timer5a); NgProfiler::StopTimer (timer5a);
// int timer5csl = (innerchartpts.Size() < 100) ? timer5cs : timer5cl; int timer5csl = (innerchartpts.Size() < 100) ? timer5cs : timer5cl;
// NgProfiler::StartTimer (timer5csl); NgProfiler::StartTimer (timer5csl);
if (!accepted) if (!accepted)
for (int k = 1; k <= 3; k++) for (int k = 1; k <= 3; k++)
@ -434,9 +451,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
if (accepted) break; if (accepted) break;
} }
// NgProfiler::StopTimer (timer5csl); NgProfiler::StopTimer (timer5csl);
} }
// NgProfiler::StartTimer (timer5b); NgProfiler::StartTimer (timer5b);
if (accepted) if (accepted)
{ {
@ -458,14 +475,15 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
} }
} }
} }
// NgProfiler::StopTimer (timer5b); NgProfiler::StopTimer (timer5b);
} }
} }
} }
} }
NgProfiler::StopTimer (timer3); NgProfiler::StopTimer (timer3);
NgProfiler::StartTimer (timere);
NgProfiler::StartTimer (timere1);
//end of while loop for outer chart //end of while loop for outer chart
GetDirtyChartTrigs(chartnum, *chart, outermark, chartpointchecked, dirtycharttrigs); GetDirtyChartTrigs(chartnum, *chart, outermark, chartpointchecked, dirtycharttrigs);
//dirtycharttrigs are local (chart) point numbers!!!!!!!!!!!!!!!! //dirtycharttrigs are local (chart) point numbers!!!!!!!!!!!!!!!!
@ -492,12 +510,21 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
} }
chartbound.DeleteSearchTree(); chartbound.DeleteSearchTree();
NgProfiler::StopTimer (timere1);
NgProfiler::StartTimer (timere2);
// cout << "key" << endl; // cout << "key" << endl;
// char key; // char key;
// cin >> key; // cin >> key;
//calculate an estimate meshsize, not to produce too large outercharts, with factor 2 larger! //calculate an estimate meshsize, not to produce too large outercharts, with factor 2 larger!
RestrictHChartDistOneChart(chartnum, chartdistacttrigs, mesh, h, 0.5, atlasminh); RestrictHChartDistOneChart(chartnum, chartdistacttrigs, mesh, h, 0.5, atlasminh);
// NgProfiler::Print(stdout); // NgProfiler::Print(stdout);
NgProfiler::StopTimer (timere2);
NgProfiler::StopTimer (timere);
} }
NgProfiler::StopTimer (timer1); NgProfiler::StopTimer (timer1);

7
libsrc/stlgeom/stltopology.cpp
View File

@ -7,6 +7,7 @@
#include <meshing.hpp> #include <meshing.hpp>
#include "stlgeom.hpp" #include "stlgeom.hpp"
#include <vector>
namespace netgen namespace netgen
{ {
@ -348,7 +349,7 @@ STLGeometry * STLTopology ::Load (istream & ist)
int cntface = 0; int cntface = 0;
int vertex = 0; int vertex = 0;
bool badnormals = false; bool badnormals = false;
while (ist.good()) while (ist.good())
{ {
ist >> buf; ist >> buf;
@ -369,7 +370,7 @@ STLGeometry * STLTopology ::Load (istream & ist)
>> normal(2); >> normal(2);
normal.Normalize(); normal.Normalize();
} }
if (strcmp (buf, "vertex") == 0) if (strcmp (buf, "vertex") == 0)
{ {
ist >> pts[vertex](0) ist >> pts[vertex](0)
@ -404,7 +405,7 @@ STLGeometry * STLTopology ::Load (istream & ist)
(Dist2 (pts[1], pts[2]) > 1e-16) ) (Dist2 (pts[1], pts[2]) > 1e-16) )
{ {
readtrigs.Append (STLReadTriangle (pts, normal)); readtrigs.Append (STLReadTriangle (pts, normal));
if (readtrigs.Size() % 100000 == 0) if (readtrigs.Size() % 100000 == 0)
PrintMessageCR (3, readtrigs.Size(), " triangles loaded\r"); PrintMessageCR (3, readtrigs.Size(), " triangles loaded\r");

17
libsrc/visualization/soldata.hpp
View File

@ -7,6 +7,15 @@ namespace netgen
using namespace std; using namespace std;
#if defined __AVX512F__
typedef __m512 tAVX;
typedef __m512d tAVXd;
#elif defined __AVX__
typedef __m256 tAVX;
typedef __m256d tAVXd;
#endif
class SolutionData class SolutionData
{ {
protected: protected:
@ -94,10 +103,10 @@ namespace netgen
#ifdef __AVX__ #ifdef __AVX__
virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts, virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts,
const __m256d * xref, const tAVXd * xref,
const __m256d * x, const tAVXd * x,
const __m256d * dxdxref, const tAVXd * dxdxref,
__m256d * values) tAVXd * values)
{ {
cerr << "GetMultiSurfVaue not overloaded for SIMD<double>" << endl; cerr << "GetMultiSurfVaue not overloaded for SIMD<double>" << endl;
return false; return false;