diff --git a/libsrc/general/ngsimd.hpp b/libsrc/general/ngsimd.hpp index e38b4e3c..f78a33ad 100644 --- a/libsrc/general/ngsimd.hpp +++ b/libsrc/general/ngsimd.hpp @@ -58,6 +58,16 @@ namespace ngsimd #endif } +#if defined __AVX512F__ + typedef __m512 tAVX; + typedef __m512d tAVXd; +#elif defined __AVX__ + typedef __m256 tAVX; + typedef __m256d tAVXd; +#endif + + + template class SIMD; template diff --git a/libsrc/general/profiler.hpp b/libsrc/general/profiler.hpp index 98039114..5786ce98 100644 --- a/libsrc/general/profiler.hpp +++ b/libsrc/general/profiler.hpp @@ -18,7 +18,7 @@ #endif -// #define USE_TSC +#define USE_TSC #ifdef USE_TSC #include // for __rdtsc() CPU time step counter #endif diff --git a/libsrc/gprim/geomops.hpp b/libsrc/gprim/geomops.hpp index 217c539b..6b06cb6e 100644 --- a/libsrc/gprim/geomops.hpp +++ b/libsrc/gprim/geomops.hpp @@ -18,10 +18,10 @@ namespace netgen */ - template - inline Vec operator+ (const Vec & a, const Vec & b) + template + inline Vec operator+ (Vec a, Vec b) { - Vec res; + Vec res; for (int i = 0; i < D; i++) res(i) = a(i) + b(i); return res; @@ -30,7 +30,7 @@ namespace netgen template - inline Point operator+ (const Point & a, const Vec & b) + inline Point operator+ (Point a, Vec b) { Point res; for (int i = 0; i < D; i++) @@ -41,7 +41,7 @@ namespace netgen template - inline Vec operator- (const Point & a, const Point & b) + inline Vec operator- (Point a, Point b) { Vec res; for (int i = 0; i < D; i++) @@ -49,19 +49,19 @@ namespace netgen return res; } - template - inline Point operator- (const Point & a, const Vec & b) + template + inline Point operator- (Point a, Vec b) { - Point res; + Point res; for (int i = 0; i < D; i++) res(i) = a(i) - b(i); return res; } - template - inline Vec operator- (const Vec & a, const Vec & b) + template + inline Vec operator- (Vec a, Vec b) { - Vec res; + Vec res; for (int i = 0; i < D; i++) res(i) = a(i) - b(i); return res; @@ -70,7 +70,7 @@ namespace netgen template - inline Vec operator* (T s, const Vec & b) + inline Vec operator* (T s, Vec b) { Vec res; for (int i = 0; i < D; i++) @@ -80,7 +80,7 @@ namespace netgen template - inline double operator* (const Vec & a, const Vec & b) + inline double operator* (Vec a, Vec b) { double sum = 0; for (int i = 0; i < D; i++) @@ -90,26 +90,26 @@ namespace netgen - template - inline Vec operator- (const Vec & b) + template + inline Vec operator- (Vec b) { - Vec res; + Vec res; for (int i = 0; i < D; i++) res(i) = -b(i); return res; } - template - inline Point & operator+= (Point & a, const Vec & b) + template + inline Point & operator+= (Point & a, Vec b) { for (int i = 0; i < D; i++) a(i) += b(i); return a; } - template - inline Vec & operator+= (Vec & a, const Vec & b) + template + inline Vec & operator+= (Vec & a, Vec b) { for (int i = 0; i < D; i++) a(i) += b(i); @@ -135,8 +135,8 @@ namespace netgen - template - inline Vec & operator*= (Vec & a, double s) + template + inline Vec & operator*= (Vec & a, T2 s) { for (int i = 0; i < D; i++) a(i) *= s; diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 503b8577..64492434 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -676,18 +676,18 @@ namespace netgen template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<1,1> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { cout << "multi-eltrafo simd called, 1,1,simd" << endl; } template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<2,2> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> (elnr, npts, @@ -720,9 +720,9 @@ namespace netgen template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<3,3> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointElementTransformation (elnr, npts, @@ -755,28 +755,34 @@ namespace netgen template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<0,2> (int elnr, int npts, - const __m256d *xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd *xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { cout << "MultiElementtransformation<0,2> simd not implemented" << endl; } template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<0,1> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { cout << "multi-eltrafo simd called, 0,1,simd" << endl; } template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<1,3> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { + mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3> + (elnr, npts, + reinterpret_cast*> (xi), sxi, + reinterpret_cast*> (x), sx, + reinterpret_cast*> (dxdxi), sdxdxi); + /* double hxi[4][1]; double hx[4][3]; double hdxdxi[4][3]; @@ -793,14 +799,21 @@ namespace netgen xi += sxi; x += sx; dxdxi += sdxdxi; + */ } template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<1,2> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { + mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2> + (elnr, npts, + reinterpret_cast*> (xi), sxi, + reinterpret_cast*> (x), sx, + reinterpret_cast*> (dxdxi), sdxdxi); + /* for (int i = 0; i < npts; i++) { double hxi[4][1]; @@ -821,14 +834,14 @@ namespace netgen x += sx; dxdxi += sdxdxi; } - + */ } template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<2,3> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, @@ -861,9 +874,9 @@ namespace netgen template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<0,3> (int elnr, int npts, - const __m256d * xi, size_t sxi, - __m256d * x, size_t sx, - __m256d * dxdxi, size_t sdxdxi) const + const tAVXd * xi, size_t sxi, + tAVXd * x, size_t sx, + tAVXd * dxdxi, size_t sdxdxi) const { for (int i = 0; i < npts; i++) { diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 554f07f9..b0513ac3 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -71,11 +71,12 @@ namespace netgen } } - - static void CalcEdgeDx (int n, double x, double * dshape) + + template + static void CalcEdgeDx (int n, T x, T * dshape) { - double p1 = x, p2 = -1, p3 = 0; - double p1dx = 1, p2dx = 0, p3dx = 0; + T p1 = x, p2 = -1, p3 = 0; + T p1dx = 1, p2dx = 0, p3dx = 0; for (int j=2; j<=n; j++) { @@ -1408,10 +1409,10 @@ namespace netgen - + template void CurvedElements :: - CalcSegmentTransformation (double xi, SegmentIndex elnr, - Point<3> * x, Vec<3> * dxdxi, bool * curved) + CalcSegmentTransformation (T xi, SegmentIndex elnr, + Point<3,T> * x, Vec<3,T> * dxdxi, bool * curved) { if (mesh.coarsemesh) { @@ -1419,11 +1420,11 @@ namespace netgen (*mesh.hpelements) [mesh[elnr].hp_elnr]; // xi umrechnen - double lami[2] = { xi, 1-xi }; - double dlami[2] = { 1, -1 }; + T lami[2] = { xi, 1-xi }; + T dlami[2] = { 1, -1 }; - double coarse_xi = 0; - double trans = 0; + T coarse_xi = 0; + T trans = 0; for (int i = 0; i < 2; i++) { coarse_xi += hpref_el.param[i][0] * lami[i]; @@ -1437,8 +1438,9 @@ namespace netgen } - Vector shapes, dshapes; - Array > coefs; + + // TVector shapes, dshapes; + // Array > coefs; SegmentInfo info; info.elnr = elnr; @@ -1452,12 +1454,21 @@ namespace netgen info.ndof += edgeorder[info.edgenr]-1; } + ArrayMem,100> coefs(info.ndof); + ArrayMem shapes_mem(info.ndof); + TFlatVector shapes(info.ndof, &shapes_mem[0]); + ArrayMem dshapes_mem(info.ndof); + TFlatVector dshapes(info.ndof, &dshapes_mem[0]); + + CalcElementShapes (info, xi, shapes); GetCoefficients (info, coefs); *x = 0; 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) @@ -1477,10 +1488,11 @@ namespace netgen } - + template void CurvedElements :: - CalcElementShapes (SegmentInfo & info, double xi, Vector & shapes) const + CalcElementShapes (SegmentInfo & info, T xi, TFlatVector shapes) const { + /* if (rational && info.order == 2) { shapes.SetSize(3); @@ -1491,9 +1503,9 @@ namespace netgen shapes *= 1.0 / (1 + (w-1) *2*xi*(1-xi)); return; } + */ - - shapes.SetSize(info.ndof); + // shapes.SetSize(info.ndof); shapes(0) = xi; shapes(1) = 1-xi; @@ -1505,9 +1517,11 @@ namespace netgen } } + template void CurvedElements :: - CalcElementDShapes (SegmentInfo & info, double xi, Vector & dshapes) const + CalcElementDShapes (SegmentInfo & info, T xi, TFlatVector dshapes) const { + /* if (rational && info.order == 2) { dshapes.SetSize(3); @@ -1527,13 +1541,11 @@ namespace netgen dshapes(j) = dshapes(j) / w - shapes[j] * dw / (w*w); return; } + */ - - - - dshapes.SetSize(info.ndof); + // dshapes.SetSize(info.ndof); dshapes = 0; dshapes(0) = 1; dshapes(1) = -1; @@ -1542,7 +1554,7 @@ namespace netgen if (info.order >= 2) { - double fac = 2; + T fac = 2; if (mesh[info.elnr][0] > mesh[info.elnr][1]) { xi = 1-xi; @@ -3651,7 +3663,7 @@ namespace netgen - + /* void CurvedElements :: CalcMultiPointSegmentTransformation (Array * xi, SegmentIndex segnr, Array > * x, @@ -3659,22 +3671,22 @@ namespace netgen { ; } + */ - - template + template void CurvedElements :: CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + const T * xi, size_t sxi, + T * x, size_t sx, + T * dxdxi, size_t sdxdxi) { for (int ip = 0; ip < n; ip++) { - Point<3> xg; - Vec<3> dx; + Point<3,T> xg; + Vec<3,T> dx; // mesh->GetCurvedElements(). - CalcSegmentTransformation (xi[ip*sxi], elnr, xg, dx); + CalcSegmentTransformation (xi[ip*sxi], elnr, &xg, &dx); if (x) for (int i = 0; i < DIM_SPACE; i++) @@ -3699,6 +3711,18 @@ namespace netgen double * x, size_t sx, double * dxdxi, size_t sdxdxi); + template void CurvedElements :: + CalcMultiPointSegmentTransformation<2> (SegmentIndex elnr, int npts, + const SIMD * xi, size_t sxi, + SIMD * x, size_t sx, + SIMD * dxdxi, size_t sdxdxi); + + template void CurvedElements :: + CalcMultiPointSegmentTransformation<3> (SegmentIndex elnr, int npts, + const SIMD * xi, size_t sxi, + SIMD * x, size_t sx, + SIMD * dxdxi, size_t sdxdxi); + void CurvedElements :: diff --git a/libsrc/meshing/curvedelems.hpp b/libsrc/meshing/curvedelems.hpp index 90ae8d40..99747034 100644 --- a/libsrc/meshing/curvedelems.hpp +++ b/libsrc/meshing/curvedelems.hpp @@ -57,15 +57,15 @@ public: void CalcSegmentTransformation (double xi, SegmentIndex segnr, Point<3> & x) - { CalcSegmentTransformation (xi, segnr, &x, NULL); }; + { CalcSegmentTransformation (xi, segnr, &x, NULL); }; void CalcSegmentTransformation (double xi, SegmentIndex segnr, Vec<3> & dxdxi) - { CalcSegmentTransformation (xi, segnr, NULL, &dxdxi); }; + { CalcSegmentTransformation (xi, segnr, NULL, &dxdxi); }; void CalcSegmentTransformation (double xi, SegmentIndex segnr, Point<3> & x, Vec<3> & dxdxi) - { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, NULL); }; + { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, NULL); }; void CalcSegmentTransformation (double xi, SegmentIndex segnr, Point<3> & x, Vec<3> & dxdxi, bool & curved) @@ -115,16 +115,17 @@ public: // { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , &curved * ); } - + /* void CalcMultiPointSegmentTransformation (Array * xi, SegmentIndex segnr, Array > * x, Array > * dxdxi); - - template + */ + + template void CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi); + const T * xi, size_t sxi, + T * x, size_t sx, + T * dxdxi, size_t sdxdxi); void CalcMultiPointSurfaceTransformation (Array< Point<2> > * xi, SurfaceElementIndex elnr, Array< Point<3> > * x, @@ -150,9 +151,10 @@ public: private: - - void CalcSegmentTransformation (double xi, SegmentIndex segnr, - Point<3> * x = NULL, Vec<3> * dxdxi = NULL, bool * curved = NULL); + + template + 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, Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL, bool * curved = NULL); @@ -176,9 +178,11 @@ private: int edgenr; }; - void CalcElementShapes (SegmentInfo & elnr, double xi, Vector & shapes) const; + template + void CalcElementShapes (SegmentInfo & elnr, T xi, TFlatVector shapes) const; void GetCoefficients (SegmentInfo & elnr, Array > & coefs) const; - void CalcElementDShapes (SegmentInfo & elnr, double xi, Vector & dshapes) const; + template + void CalcElementDShapes (SegmentInfo & elnr, T xi, TFlatVector dshapes) const; class ElementInfo diff --git a/libsrc/stlgeom/stlgeomchart.cpp b/libsrc/stlgeom/stlgeomchart.cpp index 14a7da6b..ec610411 100644 --- a/libsrc/stlgeom/stlgeomchart.cpp +++ b/libsrc/stlgeom/stlgeomchart.cpp @@ -20,16 +20,24 @@ int chartdebug = 0; void STLGeometry :: MakeAtlas(Mesh & mesh) { 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 timer3 = NgProfiler::CreateTimer ("makeatlas - part 3"); - /* 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 timer5a = NgProfiler::CreateTimer ("makeatlas - part 5a"); int timer5b = NgProfiler::CreateTimer ("makeatlas - part 5b"); int timer5cs = NgProfiler::CreateTimer ("makeatlas - part 5cs"); int timer5cl = NgProfiler::CreateTimer ("makeatlas - part 5cl"); - */ + PushStatusF("Make Atlas"); double h = mparam.maxh; @@ -110,9 +118,10 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) int markedtrigcnt = 0; while (markedtrigcnt < GetNT()) - { + { if (multithread.terminate) { PopStatus(); return; } + NgProfiler::StartTimer (timerb); if (workedarea / atlasarea*100. >= nextshow) {PrintDot(); nextshow+=showinc;} @@ -177,7 +186,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) bool changed = true; int oldstartic = 1; int oldstartic2; - + NgProfiler::StopTimer (timerb); NgProfiler::StartTimer (timer2); @@ -342,8 +351,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) { accepted = 1; - // NgProfiler::StartTimer (timer4); - + NgProfiler::StartTimer (timer4); bool isdirtytrig = false; Vec<3> gn = GetTriangle(nt).GeomNormal(points); double gnlen = gn.Length(); @@ -355,36 +363,45 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) //find overlapping charts exacter: //do not check dirty trigs! + NgProfiler::StartTimer (timer4a); + if (spiralcheckon && !isdirtytrig) for (int k = 1; k <= 3; k++) - { + { + // NgProfiler::StartTimer (timer4b); int nnt = NeighbourTrig(nt,k); if (outermark.Elem(nnt) != chartnum) { + NgProfiler::StartTimer (timer4c); int nnp1, nnp2; GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2); + NgProfiler::StopTimer (timer4c); - // NgProfiler::StartTimer (timer4a); + NgProfiler::StartTimer (timer4d); accepted = chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2), sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps); - // NgProfiler::StopTimer (timer4a); + NgProfiler::StopTimer (timer4d); + NgProfiler::StartTimer (timer4e); Vec<3> n3 = GetTriangle(nnt).Normal(); if ( (n3 * sn) >= cosouterchartangle && IsSmoothEdge (nnp1, nnp2) ) accepted = 1; + NgProfiler::StopTimer (timer4e); } + // NgProfiler::StopTimer (timer4b); 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 @@ -392,7 +409,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) if (accepted) { - // NgProfiler::StartTimer (timer5a); + NgProfiler::StartTimer (timer5a); accepted = 0; for (int k = 1; k <= 3; k++) @@ -402,9 +419,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) break; } - // NgProfiler::StopTimer (timer5a); - // int timer5csl = (innerchartpts.Size() < 100) ? timer5cs : timer5cl; - // NgProfiler::StartTimer (timer5csl); + NgProfiler::StopTimer (timer5a); + int timer5csl = (innerchartpts.Size() < 100) ? timer5cs : timer5cl; + NgProfiler::StartTimer (timer5csl); if (!accepted) for (int k = 1; k <= 3; k++) @@ -434,9 +451,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) if (accepted) break; } - // NgProfiler::StopTimer (timer5csl); + NgProfiler::StopTimer (timer5csl); } - // NgProfiler::StartTimer (timer5b); + NgProfiler::StartTimer (timer5b); if (accepted) { @@ -458,14 +475,15 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) } } } - // NgProfiler::StopTimer (timer5b); + NgProfiler::StopTimer (timer5b); } } } } NgProfiler::StopTimer (timer3); - + NgProfiler::StartTimer (timere); + NgProfiler::StartTimer (timere1); //end of while loop for outer chart GetDirtyChartTrigs(chartnum, *chart, outermark, chartpointchecked, dirtycharttrigs); //dirtycharttrigs are local (chart) point numbers!!!!!!!!!!!!!!!! @@ -492,12 +510,21 @@ void STLGeometry :: MakeAtlas(Mesh & mesh) } chartbound.DeleteSearchTree(); + + NgProfiler::StopTimer (timere1); + NgProfiler::StartTimer (timere2); + + // cout << "key" << endl; // char key; // cin >> key; //calculate an estimate meshsize, not to produce too large outercharts, with factor 2 larger! RestrictHChartDistOneChart(chartnum, chartdistacttrigs, mesh, h, 0.5, atlasminh); // NgProfiler::Print(stdout); + NgProfiler::StopTimer (timere2); + + NgProfiler::StopTimer (timere); + } NgProfiler::StopTimer (timer1); diff --git a/libsrc/stlgeom/stltopology.cpp b/libsrc/stlgeom/stltopology.cpp index 53d3d9f0..54474130 100644 --- a/libsrc/stlgeom/stltopology.cpp +++ b/libsrc/stlgeom/stltopology.cpp @@ -7,6 +7,7 @@ #include #include "stlgeom.hpp" +#include namespace netgen { @@ -348,7 +349,7 @@ STLGeometry * STLTopology ::Load (istream & ist) int cntface = 0; int vertex = 0; bool badnormals = false; - + while (ist.good()) { ist >> buf; @@ -369,7 +370,7 @@ STLGeometry * STLTopology ::Load (istream & ist) >> normal(2); normal.Normalize(); } - + if (strcmp (buf, "vertex") == 0) { ist >> pts[vertex](0) @@ -404,7 +405,7 @@ STLGeometry * STLTopology ::Load (istream & ist) (Dist2 (pts[1], pts[2]) > 1e-16) ) { - readtrigs.Append (STLReadTriangle (pts, normal)); + readtrigs.Append (STLReadTriangle (pts, normal)); if (readtrigs.Size() % 100000 == 0) PrintMessageCR (3, readtrigs.Size(), " triangles loaded\r"); diff --git a/libsrc/visualization/soldata.hpp b/libsrc/visualization/soldata.hpp index 9778c5ff..d852de9b 100644 --- a/libsrc/visualization/soldata.hpp +++ b/libsrc/visualization/soldata.hpp @@ -7,6 +7,15 @@ namespace netgen using namespace std; +#if defined __AVX512F__ + typedef __m512 tAVX; + typedef __m512d tAVXd; +#elif defined __AVX__ + typedef __m256 tAVX; + typedef __m256d tAVXd; +#endif + + class SolutionData { protected: @@ -94,10 +103,10 @@ namespace netgen #ifdef __AVX__ virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts, - const __m256d * xref, - const __m256d * x, - const __m256d * dxdxref, - __m256d * values) + const tAVXd * xref, + const tAVXd * x, + const tAVXd * dxdxref, + tAVXd * values) { cerr << "GetMultiSurfVaue not overloaded for SIMD" << endl; return false;