templated and instantiated CalcShape

libsrc/linalg/vector.hpp

@@ -8,6 +8,50 @@
 /* *************************************************************************/
 
 
+
+
+template <typename T>
+class TFlatVector
+{
+protected:
+  int s;
+  T * data;
+public:
+  TFlatVector () { ; }
+  TFlatVector (int as, T * adata)
+  { s = as; data = adata; }
+    
+  int Size () const
+  { return s; }
+
+  TFlatVector & operator= (const TFlatVector & v) 
+  { memcpy (data, v.data, s*sizeof(T)); return *this; }
+
+  TFlatVector & operator= (T scal) 
+  {
+    for (int i = 0; i < s; i++) data[i] = scal; 
+    return *this;
+  }
+
+  T & operator[] (int i) { return data[i]; }
+  const T & operator[] (int i) const { return data[i]; }
+  T & operator() (int i) { return data[i]; }
+  const T & operator() (int i) const { return data[i]; }
+
+  // double & Elem (int i) { return data[i-1]; }
+  // const double & Get (int i) const { return data[i-1]; }
+  // void Set (int i, double val) { data[i-1] = val; }
+
+  TFlatVector & operator*= (T scal)
+  {
+    for (int i = 0; i < s; i++) data[i] *= scal;
+    return *this;
+  }
+};
+
+
+
+
 class FlatVector
 {
 protected:
@@ -75,50 +119,11 @@ public:
     return sqrt (sum);
   }
 
+  operator TFlatVector<double> () const { return TFlatVector<double> (s, data); } 
   friend double operator* (const FlatVector & v1, const FlatVector & v2);
 };
 
 
-template <typename T>
-class TFlatVector
-{
-protected:
-  int s;
-  T * data;
-public:
-  TFlatVector () { ; }
-  TFlatVector (int as, T * adata)
-  { s = as; data = adata; }
-
-  int Size () const
-  { return s; }
-
-  TFlatVector & operator= (const TFlatVector & v) 
-  { memcpy (data, v.data, s*sizeof(T)); return *this; }
-
-  TFlatVector & operator= (T scal) 
-  {
-    for (int i = 0; i < s; i++) data[i] = scal; 
-    return *this;
-  }
-
-  T & operator[] (int i) { return data[i]; }
-  const T & operator[] (int i) const { return data[i]; }
-  T & operator() (int i) { return data[i]; }
-  const T & operator() (int i) const { return data[i]; }
-
-  // double & Elem (int i) { return data[i-1]; }
-  // const double & Get (int i) const { return data[i-1]; }
-  // void Set (int i, double val) { data[i-1] = val; }
-
-  TFlatVector & operator*= (T scal)
-  {
-    for (int i = 0; i < s; i++) data[i] *= scal;
-    return *this;
-  }
-};
-
-
 
 
 class Vector : public FlatVector
@@ -154,6 +159,7 @@ public:
       }
   }
 
+  operator TFlatVector<double> () const { return TFlatVector<double> (s, data); } 
 };
 
 template <int S>

libsrc/meshing/curvedelems.cpp

@@ -2313,7 +2313,7 @@ namespace netgen
 	double lami[8];
 	FlatVector vlami(8, lami);
 	vlami = 0;
-	mesh[elnr].GetShapeNew (xi, vlami);
+	mesh[elnr].GetShapeNew<double> (xi, vlami);
 
 	Mat<3,3> trans, dxdxic;
 	if (dxdxi)

libsrc/meshing/meshtype.cpp

@@ -1850,8 +1850,8 @@ namespace netgen
   }
 
 
-
+  template <typename T>
-  void Element :: GetShapeNew (const Point<3> & p, FlatVector & shape) const
+  void Element :: GetShapeNew (const Point<3,T> & p, TFlatVector<T> shape) const
   {
     /*
       if (shape.Size() < GetNP())
@@ -1874,10 +1874,10 @@ namespace netgen
 
       case TET10:
         {
-          double lam1 = p(0);
+          T lam1 = p(0);
-          double lam2 = p(1);
+          T lam2 = p(1);
-          double lam3 = p(2);
+          T lam3 = p(2);
-          double lam4 = 1-p(0)-p(1)-p(2);
+          T lam4 = 1-p(0)-p(1)-p(2);
 	
           shape(0) = 2 * lam1 * (lam1-0.5);
           shape(1) = 2 * lam2 * (lam2-0.5);
@@ -1897,11 +1897,12 @@ namespace netgen
 
       case PYRAMID:
         {
-          double noz = 1-p(2);
+          T noz = 1-p(2);
-          if (noz == 0.0) noz = 1e-10;
+          // if (noz == 0.0) noz = 1e-10;
+          noz += T(1e-12);
 
-          double xi  = p(0) / noz;
+          T xi  = p(0) / noz;
-          double eta = p(1) / noz;
+          T eta = p(1) / noz;
           shape(0) = (1-xi)*(1-eta) * (noz);
           shape(1) = (  xi)*(1-eta) * (noz);
           shape(2) = (  xi)*(  eta) * (noz);
@@ -2026,6 +2027,8 @@ namespace netgen
         }
       }
   }
+  template void Element :: GetShapeNew (const Point<3,double> & p, TFlatVector<double> shape) const;
+  template void Element :: GetShapeNew (const Point<3,SIMD<double>> & p, TFlatVector<SIMD<double>> shape) const;
   
   template void Element::GetDShapeNew<double> (const Point<3> &, MatrixFixWidth<3> &) const;
   template void Element::GetDShapeNew<SIMD<double>> (const Point<3,SIMD<double>> &, MatrixFixWidth<3,SIMD<double>> &) const;

libsrc/meshing/meshtype.hpp

@@ -716,7 +716,7 @@ namespace netgen
 			    class DenseMatrix & trans) const;
 
     void GetShape (const Point<3> & p, class Vector & shape) const;
-    void GetShapeNew (const Point<3> & p, class FlatVector & shape) const;
+    // void GetShapeNew (const Point<3> & p, class FlatVector & shape) const;
     template <typename T>
     void GetShapeNew (const Point<3,T> & p, TFlatVector<T> shape) const;    
     /// matrix 2 * np

libsrc/visualization/vssolution.cpp

@@ -1846,7 +1846,7 @@ namespace netgen
                 
                 for (int i = 0; i < cnt_valid; i++)
                   {
-                    el.GetShapeNew (locgrid[i], shape);
+                    el.GetShapeNew<double> (locgrid[i], shape);
                     Point<3> pglob;
                     for (int j = 0; j < 3; j++)
                       {
@@ -3913,7 +3913,7 @@ namespace netgen
                 
                 for (int i = 0; i < cnt_valid; i++)
                   {
-                    el.GetShapeNew (locgrid[i], shape);
+                    el.GetShapeNew<double> (locgrid[i], shape);
                     Point<3> pglob;
                     for (int j = 0; j < 3; j++)
                       {