solution visualization using AVX

libsrc/csg/python_csg.cpp

@@ -167,6 +167,9 @@ namespace netgen
   extern CSGeometry * ParseCSG (istream & istr);
 }
 
+
+static Transformation<3> global_trafo(Vec<3> (0,0,0));
+
 DLL_HEADER void ExportCSG() 
 {
   ModuleScope module("csg");
@@ -186,10 +189,22 @@ DLL_HEADER void ExportCSG()
     ;
 
   bp::def ("Pnt", FunctionPointer
-           ([](double x, double y, double z) { return Point<3>(x,y,z); }));
+           ([](double x, double y, double z) { return global_trafo(Point<3>(x,y,z)); }));
   bp::def ("Pnt", FunctionPointer
            ([](double x, double y) { return Point<2>(x,y); }));
 
+  bp::def ("SetTransformation", FunctionPointer
+           ([](int dir, double angle)
+            {
+              if (dir > 0)
+                global_trafo.SetAxisRotation (dir, angle*M_PI/180);
+              else
+                global_trafo = Transformation<3> (Vec<3>(0,0,0));
+            }),
+           (bp::arg("dir")=int(0), bp::arg("angle")=0));
+
+           
+  
   bp::class_<Vec<2>> ("Vec2d", bp::init<double,double>()) 
     .def ("__str__", &ToString<Vec<3>>)
     .def(bp::self+bp::self)
@@ -209,7 +224,7 @@ DLL_HEADER void ExportCSG()
     ;
 
   bp::def ("Vec", FunctionPointer
-           ([] (double x, double y, double z) { return Vec<3>(x,y,z); }));
+           ([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); }));
   bp::def ("Vec", FunctionPointer
            ([] (double x, double y) { return Vec<2>(x,y); }));
 

libsrc/general/mysimd.hpp

@@ -54,8 +54,18 @@ namespace netgen
 
 #ifdef __AVX__
 
+  template <typename T>
+  class AlignedAlloc
+  {
+  public:
+    void * operator new (size_t s) { return  _mm_malloc(s, alignof(T)); }
+    void * operator new[] (size_t s) { return  _mm_malloc(s, alignof(T)); }
+    void operator delete (void * p) { _mm_free(p); }
+    void operator delete[] (void * p) { _mm_free(p); }
+  };
+  
   template<>
-  class alignas(32) SIMD<double>
+  class alignas(32) SIMD<double> : public AlignedAlloc<SIMD<double>>
   {
     __m256d data;
     

libsrc/gprim/geomfuncs.hpp

@@ -82,9 +82,10 @@ namespace netgen
    */
   // inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2)
 
-  inline Vec<3> Cross (Vec<3> v1, Vec<3> v2)
+  template <typename T>
+  inline Vec<3,T> Cross (Vec<3,T> v1, Vec<3,T> v2)
   {
-    return Vec<3> 
+    return Vec<3,T> 
       ( v1(1) * v2(2) - v1(2) * v2(1),
 	v1(2) * v2(0) - v1(0) * v2(2),
 	v1(0) * v2(1) - v1(1) * v2(0) );

libsrc/gprim/geomobjects.hpp

@@ -17,7 +17,7 @@ namespace netgen
 
 
   template <int D, typename T>
-  class Point
+  class Point : public AlignedAlloc<Point<D,T>>
   {
 
   protected:
@@ -39,8 +39,9 @@ namespace netgen
     Point (T ax, T ay, T az, T au)
     { x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;}
 
-    Point (const Point<D> & p2)
+    template <typename T2>
-    { for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
+    Point (const Point<D,T2> & p2)
+    { for (int i = 0; i < D; i++) x[i] = p2(i); }
 
     explicit Point (const Vec<D> & v)
     { for (int i = 0; i < D; i++) x[i] = v(i); }
@@ -65,7 +66,7 @@ namespace netgen
   };
 
   template <int D, typename T>
-  class Vec
+  class Vec : public AlignedAlloc<Vec<D,T>>
   {
 
   protected:
@@ -90,15 +91,15 @@ namespace netgen
     Vec (const Vec<D> & p2)
     { for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
 
-    explicit Vec (const Point<D> & p)
+    explicit Vec (const Point<D,T> & p)
     { for (int i = 0; i < D; i++) x[i] = p(i); }
 
-    Vec (const Vec<D> & p1, const Vec<D> & p2)
+    Vec (const Vec & p1, const Vec & p2)
     { for(int i=0; i<D; i++) x[i] = p2(i)-p1(1); }
   
 
 
-    Vec & operator= (const Vec<D> & p2)
+    Vec & operator= (const Vec & p2)
     {
       for (int i = 0; i < D; i++) x[i] = p2.x[i]; 
       return *this;
@@ -131,12 +132,12 @@ namespace netgen
       return l;
     }
 
-    const Vec<D> & Normalize ()
+    Vec & Normalize ()
     {
       T l = Length();
-      if (l != 0)
+      // if (l != 0)
       for (int i = 0; i < D; i++)
-	  x[i] /= l;
+        x[i] /= (l+1e-40);
       return *this;
     }
 
@@ -148,7 +149,7 @@ namespace netgen
 
 
   template <int H, int W=H, typename T = double>
-  class Mat
+  class Mat : public AlignedAlloc<Mat<H,W,T>>
   {
 
   protected:

libsrc/gprim/geomops.hpp

@@ -29,10 +29,10 @@ namespace netgen
 
 
 
-  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+ (const Point<D,T> & a, const Vec<D,T> & b)
   {
-    Point<D> res;
+    Point<D,T> res;
     for (int i = 0; i < D; i++)
       res(i) = a(i) + b(i);
     return res;
@@ -40,10 +40,10 @@ namespace netgen
 
 
 
-  template <int D>
+  template <int D, typename T>
-  inline Vec<D> operator- (const Point<D> & a, const Point<D> & b)
+  inline Vec<D,T> operator- (const Point<D,T> & a, const Point<D,T> & b)
   {
-    Vec<D> res;
+    Vec<D,T> res;
     for (int i = 0; i < D; i++)
       res(i) = a(i) - b(i);
     return res;
@@ -69,10 +69,10 @@ namespace netgen
 
 
 
-  template <int D>
+  template <int D, typename T>
-  inline Vec<D> operator* (double s, const Vec<D> & b)
+  inline Vec<D,T> operator* (T s, const Vec<D,T> & b)
   {
-    Vec<D> res;
+    Vec<D,T> res;
     for (int i = 0; i < D; i++)
       res(i) = s * b(i);
     return res;

libsrc/gprim/transform3d.hpp

@@ -182,6 +182,9 @@ public:
   {
     to = m * from;
   }
+
+  Point<D> operator() (Point<D> from) const { Point<D> to; Transform(from, to); return to; }
+  Vec<D> operator() (Vec<D> from) const { Vec<D> to; Transform(from, to); return to; }
 };
 
 template <int D>

libsrc/visualization/soldata.hpp

@@ -92,6 +92,18 @@ namespace netgen
       return res;
     }
 
+#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)
+    {
+      cerr << "GetMultiSurfVaue not overloaded" << endl;
+      return false;
+    }
+#endif
+    
     virtual bool GetSegmentValue (int segnr, double xref, double * values)
     { return false; }
     

libsrc/visualization/vssolution.cpp

@@ -1256,6 +1256,14 @@ namespace netgen
     Array<double> values(npt);
 
     Array<double> mvalues(npt);
+#ifdef __AVX__
+    Array<Point<2,SIMD<double>> > simd_pref ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
+    Array<Point<3,SIMD<double>> > simd_points ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
+    Array<Mat<3,2,SIMD<double>> > simd_dxdxis ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
+    Array<Vec<3,SIMD<double>> > simd_nvs( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
+    Array<SIMD<double>> simd_values( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() * sol->components);
+#endif
+    
     if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components);
       
     Array<complex<double> > valuesc(npt);
@@ -1436,6 +1444,93 @@ namespace netgen
         
         if ( el.GetType() == TRIG || el.GetType() == TRIG6 )
           {
+#ifdef __AVX_try_it_out__
+	    bool curved = curv.IsSurfaceElementCurved(sei);
+            
+            for (int iy = 0, ii = 0; iy <= n; iy++)
+              for (int ix = 0; ix <= n-iy; ix++, ii++)
+                pref[ii] = Point<2> (ix*invn, iy*invn);
+
+            constexpr size_t simd_size = SIMD<double>::Size();
+            size_t simd_npt = (npt+simd_size-1)/simd_size;
+
+            for (size_t i = 0; i < simd_npt; i++)
+              {
+                simd_pref[i](0).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](0) : 0; }, std::true_type());
+                simd_pref[i](1).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](1) : 0; }, std::true_type());
+              }
+
+            if (curved)
+              {
+                mesh->GetCurvedElements().
+                  CalcMultiPointSurfaceTransformation<3> (sei, simd_npt,
+                                                          &simd_pref[0](0), 2,
+                                                          &simd_points[0](0), 3,
+                                                          &simd_dxdxis[0](0,0), 6);
+                
+                for (size_t ii = 0; ii < simd_npt; ii++)
+                  simd_nvs[ii] = Cross (simd_dxdxis[ii].Col(0), simd_dxdxis[ii].Col(1)).Normalize();
+              }
+            else
+              {
+		Point<3,SIMD<double>> p1 = mesh->Point (el[0]);
+		Point<3,SIMD<double>> p2 = mesh->Point (el[1]);
+		Point<3,SIMD<double>> p3 = mesh->Point (el[2]);
+
+                Vec<3,SIMD<double>> vx = p1-p3;
+                Vec<3,SIMD<double>> vy = p2-p3;
+                for (size_t ii = 0; ii < simd_npt; ii++)
+                  {
+                    simd_points[ii] = p3 + simd_pref[ii](0) * vx + simd_pref[ii](1) * vy;
+                    for (size_t j = 0; j < 3; j++)
+                      {
+                        simd_dxdxis[ii](j,0) = vx(j);
+                        simd_dxdxis[ii](j,1) = vy(j);
+                      }
+                  }
+
+                Vec<3,SIMD<double>> nv = Cross (vx, vy).Normalize();
+                for (size_t ii = 0; ii < simd_npt; ii++)
+                  simd_nvs[ii] = nv;
+              }
+
+            bool drawelem = false;
+            if (sol && sol->draw_surface) 
+              {
+		drawelem = sol->solclass->GetMultiSurfValue (sei, -1, simd_npt, 
+                                                             &simd_pref[0](0).Data(),
+                                                             &simd_points[0](0).Data(),
+                                                             &simd_dxdxis[0](0).Data(),
+                                                             &simd_values[0].Data());
+
+                for (size_t j = 0; j < sol->components; j++)
+                  for (size_t i = 0; i < npt; i++)
+                    mvalues[i*sol->components+j] = ((double*)&simd_values[j*simd_npt])[i];
+
+                if (usetexture == 2)
+		  for (int ii = 0; ii < npt; ii++)
+		    valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]);
+                else
+		  for (int ii = 0; ii < npt; ii++)
+		    values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]);
+              }
+
+            for (size_t i = 0; i < npt; i++)
+              {
+                size_t ii = i/4;
+                size_t r = i%4;
+                for (int j = 0; j < 2; j++)
+                  pref[i](j) = simd_pref[ii](j)[r];
+                for (int j = 0; j < 3; j++)
+                  points[i](j) = simd_points[ii](j)[r];
+                for (int j = 0; j < 3; j++)
+                  nvs[i](j) = simd_nvs[ii](j)[r];
+              }
+
+            if (deform)
+              for (int ii = 0; ii < npt; ii++)
+                points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
+#else
 	    bool curved = curv.IsSurfaceElementCurved(sei);
 
             for (int iy = 0, ii = 0; iy <= n; iy++)
@@ -1492,6 +1587,7 @@ namespace netgen
             if (deform)
               for (int ii = 0; ii < npt; ii++)
                 points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
+#endif
             
             int save_usetexture = usetexture;
             if (!drawelem)