Merge pull request #3642 from farscape-project/rt_opt

roystgnr · web-flow · commit bb136701c4cc · 2023-09-14T11:12:02.000-05:00
Optimise RT orientation code dealing w/ cyclic permutations
diff --git a/examples/vector_fe/vector_fe_ex6/vector_fe_ex6.C b/examples/vector_fe/vector_fe_ex6/vector_fe_ex6.C
@@ -31,9 +31,10 @@
 // The solver packages supported by libMesh.
 #include "libmesh/enum_solver_package.h"
 
-// The mesh object and mesh generation utilities.
+// The mesh object and mesh generation and modification utilities.
 #include "libmesh/mesh.h"
 #include "libmesh/mesh_generation.h"
+#include "libmesh/mesh_modification.h"
 
 // Matrix and vector types.
 #include "libmesh/dense_matrix.h"
@@ -131,6 +132,9 @@ int main (int argc, char ** argv)
                                        -1., 1.,
                                        Utility::string_to_enum<ElemType>(elem_str));
 
+  // Make sure the code is robust against nodal reorderings.
+  MeshTools::Modification::permute_elements(mesh);
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
diff --git a/src/fe/fe_raviart_shape_3D.C b/src/fe/fe_raviart_shape_3D.C
@@ -25,10 +25,11 @@ namespace libMesh
 {
 
 
-bool orientation(std::vector<Point> & arr)
+template<size_t N>
+bool orientation(std::array<Point, N> & arr)
 {
-  while (std::min_element(arr.begin(), arr.end()) != arr.begin())
-    std::rotate(arr.begin(), arr.begin() + 1, arr.end());
+  if (N % 2 == 0)
+    std::rotate(arr.begin(), std::min_element(arr.begin(), arr.end()), arr.end());
 
   size_t cnt = 0;
     for(size_t i = 0; i < arr.size(); i++)
@@ -76,47 +77,47 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,
                 {
                 case 0:
                   {
-                    std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
+                    std::array<Point, 4> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
                     if (orientation(arr))
                       return RealGradient( 0.0, 0.0,  0.125*(zeta-1.0) );
                     else
                       return RealGradient( 0.0, 0.0, -0.125*(zeta-1.0) );
                   }
                 case 1:
                   {
-                    std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
+                    std::array<Point, 4> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
                     if (orientation(arr))
                       return RealGradient( 0.0,  0.125*(eta-1.0), 0.0 );
                     else
                       return RealGradient( 0.0, -0.125*(eta-1.0), 0.0 );
                   }
                 case 2:
                   {
-                    std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
+                    std::array<Point, 4> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
                     if (orientation(arr))
                       return RealGradient(  0.125*(xi+1.0), 0.0, 0.0 );
                     else
                       return RealGradient( -0.125*(xi+1.0), 0.0, 0.0 );
                   }
                 case 3:
                   {
-                    std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
+                    std::array<Point, 4> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
                     if (orientation(arr))
                       return RealGradient( 0.0,  0.125*(1.0+eta), 0.0 );
                     else
                       return RealGradient( 0.0, -0.125*(1.0+eta), 0.0 );
                   }
                 case 4:
                   {
-                    std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
+                    std::array<Point, 4> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
                     if (orientation(arr))
                       return RealGradient(  0.125*(xi-1.0), 0.0, 0.0 );
                     else
                       return RealGradient( -0.125*(xi-1.0), 0.0, 0.0 );
                   }
                 case 5:
                   {
-                    std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
+                    std::array<Point, 4> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
                     if (orientation(arr))
                       return RealGradient( 0.0, 0.0,  0.125*(1.0+zeta) );
                     else
@@ -141,31 +142,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,
                 {
                 case 0:
                   {
-                    std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
+                    std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};
                     if (orientation(arr))
                       return RealGradient(  2.0*xi,  2.0*eta,  2.0*zeta-2.0 );
                     else
                       return RealGradient( -2.0*xi, -2.0*eta, -2.0*zeta+2.0 );
                   }
                 case 1:
                   {
-                    std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
+                    std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};
                     if (orientation(arr))
                       return RealGradient(  2.0*xi,  2.0*eta-2.0,  2.0*zeta );
                     else
                       return RealGradient( -2.0*xi, -2.0*eta+2.0, -2.0*zeta );
                   }
                 case 2:
                   {
-                    std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
+                    std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};
                     if (orientation(arr))
                       return RealGradient(  2.0*xi,  2.0*eta,  2.0*zeta );
                     else
                       return RealGradient( -2.0*xi, -2.0*eta, -2.0*zeta );
                   }
                 case 3:
                   {
-                    std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
+                    std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};
                     if (orientation(arr))
                       return RealGradient(  2.0*xi-2.0,  2.0*eta,  2.0*zeta );
                     else
@@ -257,15 +258,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                         return RealGradient();
                       case 2:
                         {
-                          std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
+                          std::array<Point, 4> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
                           if (orientation(arr))
                             return RealGradient(  0.125, 0.0, 0.0 );
                           else
                             return RealGradient( -0.125, 0.0, 0.0 );
                         }
                       case 4:
                         {
-                          std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
+                          std::array<Point, 4> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
                           if (orientation(arr))
                             return RealGradient(  0.125, 0.0, 0.0 );
                           else
@@ -289,15 +290,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                         return RealGradient();
                       case 1:
                         {
-                          std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
+                          std::array<Point, 4> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
                           if (orientation(arr))
                             return RealGradient( 0.0,  0.125, 0.0 );
                           else
                             return RealGradient( 0.0, -0.125, 0.0 );
                         }
                       case 3:
                         {
-                          std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
+                          std::array<Point, 4> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
                           if (orientation(arr))
                             return RealGradient( 0.0,  0.125, 0.0 );
                           else
@@ -321,15 +322,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                         return RealGradient();
                       case 0:
                         {
-                          std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
+                          std::array<Point, 4> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
                           if (orientation(arr))
                             return RealGradient( 0.0, 0.0,  0.125 );
                           else
                             return RealGradient( 0.0, 0.0, -0.125 );
                         }
                       case 5:
                         {
-                          std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
+                          std::array<Point, 4> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
                           if (orientation(arr))
                             return RealGradient( 0.0, 0.0,  0.125 );
                           else
@@ -361,31 +362,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                     {
                     case 0:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};
                         if (orientation(arr))
                           return RealGradient(  2.0, 0.0, 0.0 );
                         else
                           return RealGradient( -2.0, 0.0, 0.0 );
                       }
                     case 1:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};
                         if (orientation(arr))
                           return RealGradient(  2.0, 0.0, 0.0 );
                         else
                           return RealGradient( -2.0, 0.0, 0.0 );
                       }
                     case 2:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};
                         if (orientation(arr))
                           return RealGradient(  2.0, 0.0, 0.0 );
                         else
                           return RealGradient( -2.0, 0.0, 0.0 );
                       }
                     case 3:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};
                         if (orientation(arr))
                           return RealGradient(  2.0, 0.0, 0.0 );
                         else
@@ -404,31 +405,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                     {
                     case 0:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};
                         if (orientation(arr))
                           return RealGradient( 0.0,  2.0, 0.0 );
                         else
                           return RealGradient( 0.0, -2.0, 0.0 );
                       }
                     case 1:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};
                         if (orientation(arr))
                           return RealGradient( 0.0,  2.0, 0.0 );
                         else
                           return RealGradient( 0.0, -2.0, 0.0 );
                       }
                     case 2:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};
                         if (orientation(arr))
                           return RealGradient( 0.0,  2.0, 0.0 );
                         else
                           return RealGradient( 0.0, -2.0, 0.0 );
                       }
                     case 3:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};
                         if (orientation(arr))
                           return RealGradient( 0.0,  2.0, 0.0 );
                         else
@@ -447,31 +448,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
                     {
                     case 0:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};
                         if (orientation(arr))
                           return RealGradient( 0.0, 0.0,  2.0 );
                         else
                           return RealGradient( 0.0, 0.0, -2.0 );
                       }
                     case 1:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};
                         if (orientation(arr))
                           return RealGradient( 0.0, 0.0,  2.0 );
                         else
                           return RealGradient( 0.0, 0.0, -2.0 );
                       }
                     case 2:
                       {
-                        std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
+                        std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};
                         if (orientation(arr))
                           return RealGradient( 0.0, 0.0,  2.0 );
                         else
                           return RealGradient( 0.0, 0.0, -2.0 );
                       }
                     case 3:
                       {
-                        std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
+                        std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};
                         if (orientation(arr))
                           return RealGradient( 0.0, 0.0,  2.0 );
                         else

Original file line number	Diff line number	Diff line change
`@@ -25,10 +25,11 @@ namespace libMesh`
`25`	`25`	`{`
`26`	`26`
`27`	`27`
`28`		`-bool orientation(std::vector<Point> & arr)`
	`28`	`+template<size_t N>`
	`29`	`+bool orientation(std::array<Point, N> & arr)`
`29`	`30`	`{`
`30`		`- while (std::min_element(arr.begin(), arr.end()) != arr.begin())`
`31`		`- std::rotate(arr.begin(), arr.begin() + 1, arr.end());`
	`31`	`+ if (N % 2 == 0)`
	`32`	`+ std::rotate(arr.begin(), std::min_element(arr.begin(), arr.end()), arr.end());`
`32`	`33`
`33`	`34`	`size_t cnt = 0;`
`34`	`35`	`for(size_t i = 0; i < arr.size(); i++)`
`@@ -76,47 +77,47 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,`
`76`	`77`	`{`
`77`	`78`	`case 0:`
`78`	`79`	`{`
`79`		`- std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};`
	`80`	`+ std::array<Point, 4> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};`
`80`	`81`	`if (orientation(arr))`
`81`	`82`	`return RealGradient( 0.0, 0.0, 0.125*(zeta-1.0) );`
`82`	`83`	`else`
`83`	`84`	`return RealGradient( 0.0, 0.0, -0.125*(zeta-1.0) );`
`84`	`85`	`}`
`85`	`86`	`case 1:`
`86`	`87`	`{`
`87`		`- std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};`
	`88`	`+ std::array<Point, 4> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};`
`88`	`89`	`if (orientation(arr))`
`89`	`90`	`return RealGradient( 0.0, 0.125*(eta-1.0), 0.0 );`
`90`	`91`	`else`
`91`	`92`	`return RealGradient( 0.0, -0.125*(eta-1.0), 0.0 );`
`92`	`93`	`}`
`93`	`94`	`case 2:`
`94`	`95`	`{`
`95`		`- std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};`
	`96`	`+ std::array<Point, 4> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};`
`96`	`97`	`if (orientation(arr))`
`97`	`98`	`return RealGradient( 0.125*(xi+1.0), 0.0, 0.0 );`
`98`	`99`	`else`
`99`	`100`	`return RealGradient( -0.125*(xi+1.0), 0.0, 0.0 );`
`100`	`101`	`}`
`101`	`102`	`case 3:`
`102`	`103`	`{`
`103`		`- std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};`
	`104`	`+ std::array<Point, 4> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};`
`104`	`105`	`if (orientation(arr))`
`105`	`106`	`return RealGradient( 0.0, 0.125*(1.0+eta), 0.0 );`
`106`	`107`	`else`
`107`	`108`	`return RealGradient( 0.0, -0.125*(1.0+eta), 0.0 );`
`108`	`109`	`}`
`109`	`110`	`case 4:`
`110`	`111`	`{`
`111`		`- std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};`
	`112`	`+ std::array<Point, 4> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};`
`112`	`113`	`if (orientation(arr))`
`113`	`114`	`return RealGradient( 0.125*(xi-1.0), 0.0, 0.0 );`
`114`	`115`	`else`
`115`	`116`	`return RealGradient( -0.125*(xi-1.0), 0.0, 0.0 );`
`116`	`117`	`}`
`117`	`118`	`case 5:`
`118`	`119`	`{`
`119`		`- std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};`
	`120`	`+ std::array<Point, 4> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};`
`120`	`121`	`if (orientation(arr))`
`121`	`122`	`return RealGradient( 0.0, 0.0, 0.125*(1.0+zeta) );`
`122`	`123`	`else`
`@@ -141,31 +142,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,`
`141`	`142`	`{`
`142`	`143`	`case 0:`
`143`	`144`	`{`
`144`		`- std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};`
	`145`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};`
`145`	`146`	`if (orientation(arr))`
`146`	`147`	`return RealGradient( 2.0xi, 2.0eta, 2.0*zeta-2.0 );`
`147`	`148`	`else`
`148`	`149`	`return RealGradient( -2.0xi, -2.0eta, -2.0*zeta+2.0 );`
`149`	`150`	`}`
`150`	`151`	`case 1:`
`151`	`152`	`{`
`152`		`- std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};`
	`153`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};`
`153`	`154`	`if (orientation(arr))`
`154`	`155`	`return RealGradient( 2.0xi, 2.0eta-2.0, 2.0*zeta );`
`155`	`156`	`else`
`156`	`157`	`return RealGradient( -2.0xi, -2.0eta+2.0, -2.0*zeta );`
`157`	`158`	`}`
`158`	`159`	`case 2:`
`159`	`160`	`{`
`160`		`- std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};`
	`161`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};`
`161`	`162`	`if (orientation(arr))`
`162`	`163`	`return RealGradient( 2.0xi, 2.0eta, 2.0*zeta );`
`163`	`164`	`else`
`164`	`165`	`return RealGradient( -2.0xi, -2.0eta, -2.0*zeta );`
`165`	`166`	`}`
`166`	`167`	`case 3:`
`167`	`168`	`{`
`168`		`- std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};`
	`169`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};`
`169`	`170`	`if (orientation(arr))`
`170`	`171`	`return RealGradient( 2.0xi-2.0, 2.0eta, 2.0*zeta );`
`171`	`172`	`else`
`@@ -257,15 +258,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`257`	`258`	`return RealGradient();`
`258`	`259`	`case 2:`
`259`	`260`	`{`
`260`		`- std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};`
	`261`	`+ std::array<Point, 4> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};`
`261`	`262`	`if (orientation(arr))`
`262`	`263`	`return RealGradient( 0.125, 0.0, 0.0 );`
`263`	`264`	`else`
`264`	`265`	`return RealGradient( -0.125, 0.0, 0.0 );`
`265`	`266`	`}`
`266`	`267`	`case 4:`
`267`	`268`	`{`
`268`		`- std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};`
	`269`	`+ std::array<Point, 4> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};`
`269`	`270`	`if (orientation(arr))`
`270`	`271`	`return RealGradient( 0.125, 0.0, 0.0 );`
`271`	`272`	`else`
`@@ -289,15 +290,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`289`	`290`	`return RealGradient();`
`290`	`291`	`case 1:`
`291`	`292`	`{`
`292`		`- std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};`
	`293`	`+ std::array<Point, 4> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};`
`293`	`294`	`if (orientation(arr))`
`294`	`295`	`return RealGradient( 0.0, 0.125, 0.0 );`
`295`	`296`	`else`
`296`	`297`	`return RealGradient( 0.0, -0.125, 0.0 );`
`297`	`298`	`}`
`298`	`299`	`case 3:`
`299`	`300`	`{`
`300`		`- std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};`
	`301`	`+ std::array<Point, 4> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};`
`301`	`302`	`if (orientation(arr))`
`302`	`303`	`return RealGradient( 0.0, 0.125, 0.0 );`
`303`	`304`	`else`
`@@ -321,15 +322,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`321`	`322`	`return RealGradient();`
`322`	`323`	`case 0:`
`323`	`324`	`{`
`324`		`- std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};`
	`325`	`+ std::array<Point, 4> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};`
`325`	`326`	`if (orientation(arr))`
`326`	`327`	`return RealGradient( 0.0, 0.0, 0.125 );`
`327`	`328`	`else`
`328`	`329`	`return RealGradient( 0.0, 0.0, -0.125 );`
`329`	`330`	`}`
`330`	`331`	`case 5:`
`331`	`332`	`{`
`332`		`- std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};`
	`333`	`+ std::array<Point, 4> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};`
`333`	`334`	`if (orientation(arr))`
`334`	`335`	`return RealGradient( 0.0, 0.0, 0.125 );`
`335`	`336`	`else`
`@@ -361,31 +362,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`361`	`362`	`{`
`362`	`363`	`case 0:`
`363`	`364`	`{`
`364`		`- std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};`
	`365`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};`
`365`	`366`	`if (orientation(arr))`
`366`	`367`	`return RealGradient( 2.0, 0.0, 0.0 );`
`367`	`368`	`else`
`368`	`369`	`return RealGradient( -2.0, 0.0, 0.0 );`
`369`	`370`	`}`
`370`	`371`	`case 1:`
`371`	`372`	`{`
`372`		`- std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};`
	`373`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};`
`373`	`374`	`if (orientation(arr))`
`374`	`375`	`return RealGradient( 2.0, 0.0, 0.0 );`
`375`	`376`	`else`
`376`	`377`	`return RealGradient( -2.0, 0.0, 0.0 );`
`377`	`378`	`}`
`378`	`379`	`case 2:`
`379`	`380`	`{`
`380`		`- std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};`
	`381`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};`
`381`	`382`	`if (orientation(arr))`
`382`	`383`	`return RealGradient( 2.0, 0.0, 0.0 );`
`383`	`384`	`else`
`384`	`385`	`return RealGradient( -2.0, 0.0, 0.0 );`
`385`	`386`	`}`
`386`	`387`	`case 3:`
`387`	`388`	`{`
`388`		`- std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};`
	`389`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};`
`389`	`390`	`if (orientation(arr))`
`390`	`391`	`return RealGradient( 2.0, 0.0, 0.0 );`
`391`	`392`	`else`
`@@ -404,31 +405,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`404`	`405`	`{`
`405`	`406`	`case 0:`
`406`	`407`	`{`
`407`		`- std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};`
	`408`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};`
`408`	`409`	`if (orientation(arr))`
`409`	`410`	`return RealGradient( 0.0, 2.0, 0.0 );`
`410`	`411`	`else`
`411`	`412`	`return RealGradient( 0.0, -2.0, 0.0 );`
`412`	`413`	`}`
`413`	`414`	`case 1:`
`414`	`415`	`{`
`415`		`- std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};`
	`416`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};`
`416`	`417`	`if (orientation(arr))`
`417`	`418`	`return RealGradient( 0.0, 2.0, 0.0 );`
`418`	`419`	`else`
`419`	`420`	`return RealGradient( 0.0, -2.0, 0.0 );`
`420`	`421`	`}`
`421`	`422`	`case 2:`
`422`	`423`	`{`
`423`		`- std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};`
	`424`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};`
`424`	`425`	`if (orientation(arr))`
`425`	`426`	`return RealGradient( 0.0, 2.0, 0.0 );`
`426`	`427`	`else`
`427`	`428`	`return RealGradient( 0.0, -2.0, 0.0 );`
`428`	`429`	`}`
`429`	`430`	`case 3:`
`430`	`431`	`{`
`431`		`- std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};`
	`432`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};`
`432`	`433`	`if (orientation(arr))`
`433`	`434`	`return RealGradient( 0.0, 2.0, 0.0 );`
`434`	`435`	`else`
`@@ -447,31 +448,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,`
`447`	`448`	`{`
`448`	`449`	`case 0:`
`449`	`450`	`{`
`450`		`- std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};`
	`451`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(2), elem->point(1)};`
`451`	`452`	`if (orientation(arr))`
`452`	`453`	`return RealGradient( 0.0, 0.0, 2.0 );`
`453`	`454`	`else`
`454`	`455`	`return RealGradient( 0.0, 0.0, -2.0 );`
`455`	`456`	`}`
`456`	`457`	`case 1:`
`457`	`458`	`{`
`458`		`- std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};`
	`459`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(3), elem->point(0)};`
`459`	`460`	`if (orientation(arr))`
`460`	`461`	`return RealGradient( 0.0, 0.0, 2.0 );`
`461`	`462`	`else`
`462`	`463`	`return RealGradient( 0.0, 0.0, -2.0 );`
`463`	`464`	`}`
`464`	`465`	`case 2:`
`465`	`466`	`{`
`466`		`- std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};`
	`467`	`+ std::array<Point, 3> arr = {elem->point(1), elem->point(2), elem->point(3)};`
`467`	`468`	`if (orientation(arr))`
`468`	`469`	`return RealGradient( 0.0, 0.0, 2.0 );`
`469`	`470`	`else`
`470`	`471`	`return RealGradient( 0.0, 0.0, -2.0 );`
`471`	`472`	`}`
`472`	`473`	`case 3:`
`473`	`474`	`{`
`474`		`- std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};`
	`475`	`+ std::array<Point, 3> arr = {elem->point(0), elem->point(3), elem->point(2)};`
`475`	`476`	`if (orientation(arr))`
`476`	`477`	`return RealGradient( 0.0, 0.0, 2.0 );`
`477`	`478`	`else`