Walker test code coverage report
Current view: top level - Base - Vector.hpp (source / functions) Hit Total Coverage
Commit: test_coverage.info Lines: 9 9 100.0 %
Date: 2022-09-21 13:52:12 Functions: 1 1 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 0 0 -

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Base/Vector.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Vector algebra
       9                 :            :   \details   Vector algebra.
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : #ifndef Vector_h
      13                 :            : #define Vector_h
      14                 :            : 
      15                 :            : #include <array>
      16                 :            : #include <cmath>
      17                 :            : #include <vector>
      18                 :            : 
      19                 :            : #include "Types.hpp"
      20                 :            : #include "Exception.hpp"
      21                 :            : 
      22                 :            : namespace tk {
      23                 :            : 
      24                 :            : //! Flip sign of vector components
      25                 :            : //! \param[in] v Vector whose components to multiply by -1.0
      26                 :            : inline void
      27                 :            : flip( std::array< real, 3 >& v )
      28                 :            : {
      29                 :            :   v[0] = -v[0];
      30                 :            :   v[1] = -v[1];
      31                 :            :   v[2] = -v[2];
      32                 :            : }
      33                 :            : 
      34                 :            : //! Compute the cross-product of two vectors
      35                 :            : //! \param[in] v1x x coordinate of the 1st vector
      36                 :            : //! \param[in] v1y y coordinate of the 1st vector
      37                 :            : //! \param[in] v1z z coordinate of the 1st vector
      38                 :            : //! \param[in] v2x x coordinate of the 2nd vector
      39                 :            : //! \param[in] v2y y coordinate of the 2nd vector
      40                 :            : //! \param[in] v2z z coordinate of the 2nd vector
      41                 :            : //! \param[out] rx x coordinate of the product vector
      42                 :            : //! \param[out] ry y coordinate of the product vector
      43                 :            : //! \param[out] rz z coordinate of the product vector
      44                 :            : #pragma omp declare simd
      45                 :            : inline void
      46                 :            : cross( real v1x, real v1y, real v1z,
      47                 :            :        real v2x, real v2y, real v2z,
      48                 :            :        real& rx, real& ry, real& rz )
      49                 :            : {
      50                 :            :   rx = v1y*v2z - v2y*v1z;
      51                 :            :   ry = v1z*v2x - v2z*v1x;
      52                 :            :   rz = v1x*v2y - v2x*v1y;
      53                 :            : }
      54                 :            : 
      55                 :            : //! Compute the cross-product of two vectors
      56                 :            : //! \param[in] v1 1st vector
      57                 :            : //! \param[in] v2 2nd vector
      58                 :            : //! \return Cross-product
      59                 :            : inline std::array< real, 3 >
      60                 :            : cross( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
      61                 :            : {
      62                 :            :   real rx, ry, rz;
      63                 :            :   cross( v1[0], v1[1], v1[2], v2[0], v2[1], v2[2], rx, ry, rz );
      64                 :          1 :   return { std::move(rx), std::move(ry), std::move(rz) };
      65                 :            : }
      66                 :            : 
      67                 :            : //! Compute the cross-product of two vectors divided by a scalar
      68                 :            : //! \param[in] v1x x coordinate of the 1st vector
      69                 :            : //! \param[in] v1y y coordinate of the 1st vector
      70                 :            : //! \param[in] v1z z coordinate of the 1st vector
      71                 :            : //! \param[in] v2x x coordinate of the 2nd vector
      72                 :            : //! \param[in] v2y y coordinate of the 2nd vector
      73                 :            : //! \param[in] v2z z coordinate of the 2nd vector
      74                 :            : //! \param[in] j The scalar to divide the product with
      75                 :            : //! \param[out] rx x coordinate of the product vector
      76                 :            : //! \param[out] ry y coordinate of the product vector
      77                 :            : //! \param[out] rz z coordinate of the product vector
      78                 :            : #pragma omp declare simd uniform(j)
      79                 :            : inline void
      80                 :            : crossdiv( real v1x, real v1y, real v1z,
      81                 :            :           real v2x, real v2y, real v2z,
      82                 :            :           real j,
      83                 :            :           real& rx, real& ry, real& rz )
      84                 :            : {
      85                 :            :   cross( v1x, v1y, v1z, v2x, v2y, v2z, rx, ry, rz );
      86                 :            :   rx /= j;
      87                 :            :   ry /= j;
      88                 :            :   rz /= j;
      89                 :            : }
      90                 :            : 
      91                 :            : //! Compute the cross-product of two vectors divided by a scalar
      92                 :            : //! \param[in] v1 1st vector
      93                 :            : //! \param[in] v2 2nd vector
      94                 :            : //! \param[in] j Scalar to divide each component by
      95                 :            : //! \return Cross-product divided by scalar
      96                 :            : inline std::array< real, 3 >
      97                 :            : crossdiv( const std::array< real, 3 >& v1,
      98                 :            :           const std::array< real, 3 >& v2,
      99                 :            :           real j )
     100                 :            : {
     101                 :            :   return {{ (v1[1]*v2[2] - v2[1]*v1[2]) / j,
     102                 :            :             (v1[2]*v2[0] - v2[2]*v1[0]) / j,
     103                 :          1 :             (v1[0]*v2[1] - v2[0]*v1[1]) / j }};
     104                 :            : }
     105                 :            : 
     106                 :            : //! Compute the dot-product of two vectors
     107                 :            : //! \param[in] v1 1st vector
     108                 :            : //! \param[in] v2 2nd vector
     109                 :            : //! \return Dot-product
     110                 :            : inline real
     111                 :            : dot( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
     112                 :            : {
     113                 :          4 :   return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
     114                 :            : }
     115                 :            : 
     116                 :            : //! Compute length of a vector
     117                 :            : //! \param[in] x X coordinate of vector
     118                 :            : //! \param[in] y Y coordinate of vector
     119                 :            : //! \param[in] z Z coordinate of vector
     120                 :            : //! \return length
     121                 :            : #pragma omp declare simd
     122                 :            : inline real
     123                 :            : length( real x, real y, real z )
     124                 :            : {
     125                 :            :   return std::sqrt( x*x + y*y + z*z );
     126                 :            : }
     127                 :            : 
     128                 :            : //! Compute length of a vector
     129                 :            : //! \param[in] v vector
     130                 :            : //! \return length
     131                 :            : inline real
     132                 :            : length( const std::array< real, 3 >& v )
     133                 :            : {
     134                 :          1 :   return std::sqrt( dot(v,v) );
     135                 :            : }
     136                 :            : 
     137                 :            : //! Scale vector to unit length
     138                 :            : //! \param[in,out] v Vector to normalize
     139                 :            : inline void
     140                 :          2 : unit( std::array< real, 3 >& v ) noexcept(ndebug)
     141                 :            : {
     142                 :            :   auto len = length( v );
     143                 :            :   Assert( len > std::numeric_limits< tk::real >::epsilon(), "div by zero" );
     144                 :          2 :   v[0] /= len;
     145                 :          2 :   v[1] /= len;
     146                 :          2 :   v[2] /= len;
     147                 :          2 : }
     148                 :            : 
     149                 :            : //! Compute the triple-product of three vectors
     150                 :            : //! \param[in] v1x x coordinate of the 1st vector
     151                 :            : //! \param[in] v1y y coordinate of the 1st vector
     152                 :            : //! \param[in] v1z z coordinate of the 1st vector
     153                 :            : //! \param[in] v2x x coordinate of the 2nd vector
     154                 :            : //! \param[in] v2y y coordinate of the 2nd vector
     155                 :            : //! \param[in] v2z z coordinate of the 2nd vector
     156                 :            : //! \param[in] v3x x coordinate of the 3rd vector
     157                 :            : //! \param[in] v3y y coordinate of the 3rd vector
     158                 :            : //! \param[in] v3z z coordinate of the 3rd vector
     159                 :            : //! \return Scalar value of the triple product
     160                 :            : #pragma omp declare simd
     161                 :            : inline tk::real
     162                 :            : triple( real v1x, real v1y, real v1z,
     163                 :            :         real v2x, real v2y, real v2z,
     164                 :            :         real v3x, real v3y, real v3z )
     165                 :            : {
     166                 :            :   real cx, cy, cz;
     167                 :            :   cross( v2x, v2y, v2z, v3x, v3y, v3z, cx, cy, cz );
     168                 :            :   return v1x*cx + v1y*cy + v1z*cz;
     169                 :            : }
     170                 :            : 
     171                 :            : //! Compute the triple-product of three vectors
     172                 :            : //! \param[in] v1 1st vector
     173                 :            : //! \param[in] v2 2nd vector
     174                 :            : //! \param[in] v3 3rd vector
     175                 :            : //! \return Triple-product
     176                 :            : inline real
     177                 :            : triple( const std::array< real, 3 >& v1,
     178                 :            :         const std::array< real, 3 >& v2,
     179                 :            :         const std::array< real, 3 >& v3 )
     180                 :            : {
     181                 :            :   return dot( v1, cross(v2,v3) );
     182                 :            : }
     183                 :            : 
     184                 :            : //! Rotate vector about X axis
     185                 :            : //! \param[in] v Vector to rotate
     186                 :            : //! \param[in] angle Angle to use to rotate with
     187                 :            : //! \return Rotated vector
     188                 :            : inline std::array< real, 3 >
     189                 :            : rotateX( const std::array< real, 3 >& v, real angle )
     190                 :            : {
     191                 :            :   using std::cos;  using std::sin;
     192                 :            : 
     193                 :            :   std::array< std::array< real, 3 >, 3 >
     194                 :            :     R{{ {{ 1.0,         0.0,          0.0 }},
     195                 :            :         {{ 0.0,   cos(angle), -sin(angle) }},
     196                 :            :         {{ 0.0,   sin(angle),  cos(angle) }} }};
     197                 :            : 
     198                 :            :   return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
     199                 :            : }
     200                 :            : 
     201                 :            : //! Rotate vector about Y axis
     202                 :            : //! \param[in] v Vector to rotate
     203                 :            : //! \param[in] angle Angle to use to rotate with
     204                 :            : //! \return Rotated vector
     205                 :            : inline std::array< real, 3 >
     206                 :            : rotateY( const std::array< real, 3 >& v, real angle )
     207                 :            : {
     208                 :            :   using std::cos;  using std::sin;
     209                 :            : 
     210                 :            :   std::array< std::array< real, 3 >, 3 >
     211                 :            :     R{{ {{ cos(angle),  0.0, sin(angle) }},
     212                 :            :         {{ 0.0,         1.0,        0.0 }},
     213                 :            :         {{ -sin(angle), 0.0, cos(angle) }} }};
     214                 :            : 
     215                 :            :   return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
     216                 :            : }
     217                 :            : 
     218                 :            : //! Rotate vector about Z axis
     219                 :            : //! \param[in] v Vector to rotate
     220                 :            : //! \param[in] angle Angle to use to rotate with
     221                 :            : //! \return Rotated vector
     222                 :            : inline std::array< real, 3 >
     223                 :            : rotateZ( const std::array< real, 3 >& v, real angle )
     224                 :            : {
     225                 :            :   using std::cos;  using std::sin;
     226                 :            : 
     227                 :            :   std::array< std::array< real, 3 >, 3 >
     228                 :            :     R{{ {{ cos(angle), -sin(angle), 0.0 }},
     229                 :            :         {{ sin(angle),  cos(angle), 0.0 }},
     230                 :            :         {{ 0.0,         0.0,        1.0 }} }};
     231                 :            : 
     232                 :            :   return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
     233                 :            : }
     234                 :            : 
     235                 :            : //! \brief Compute the determinant of the Jacobian of a coordinate
     236                 :            : //!  transformation over a tetrahedron
     237                 :            : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
     238                 :            : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
     239                 :            : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
     240                 :            : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
     241                 :            : //! \return Determinant of the Jacobian of transformation of physical
     242                 :            : //!   tetrahedron to reference (xi, eta, zeta) space
     243                 :            : inline real
     244                 :            : Jacobian( const std::array< real, 3 >& v1,
     245                 :            :           const std::array< real, 3 >& v2,
     246                 :            :           const std::array< real, 3 >& v3,
     247                 :            :           const std::array< real, 3 >& v4 )
     248                 :            : {
     249                 :            :   std::array< real, 3 > ba{{ v2[0]-v1[0], v2[1]-v1[1], v2[2]-v1[2] }},
     250                 :            :                         ca{{ v3[0]-v1[0], v3[1]-v1[1], v3[2]-v1[2] }},
     251                 :            :                         da{{ v4[0]-v1[0], v4[1]-v1[1], v4[2]-v1[2] }};
     252                 :            :   return triple( ba, ca, da );
     253                 :            : }
     254                 :            : 
     255                 :            : //! \brief Compute the inverse of the Jacobian of a coordinate transformation
     256                 :            : //!   over a tetrahedron
     257                 :            : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
     258                 :            : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
     259                 :            : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
     260                 :            : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
     261                 :            : //! \return Inverse of the Jacobian of transformation of physical
     262                 :            : //!   tetrahedron to reference (xi, eta, zeta) space
     263                 :            : inline std::array< std::array< real, 3 >, 3 >
     264                 :            : inverseJacobian( const std::array< real, 3 >& v1,
     265                 :            :                  const std::array< real, 3 >& v2,
     266                 :            :                  const std::array< real, 3 >& v3,
     267                 :            :                  const std::array< real, 3 >& v4 )
     268                 :            : {
     269                 :            :   std::array< std::array< real, 3 >, 3 > jacInv;
     270                 :            : 
     271                 :            :   auto detJ = Jacobian( v1, v2, v3, v4 );
     272                 :            : 
     273                 :            :   jacInv[0][0] =  (  (v3[1]-v1[1])*(v4[2]-v1[2])
     274                 :            :                    - (v4[1]-v1[1])*(v3[2]-v1[2])) / detJ;
     275                 :            :   jacInv[1][0] = -(  (v2[1]-v1[1])*(v4[2]-v1[2])
     276                 :            :                    - (v4[1]-v1[1])*(v2[2]-v1[2])) / detJ;
     277                 :            :   jacInv[2][0] =  (  (v2[1]-v1[1])*(v3[2]-v1[2])
     278                 :            :                    - (v3[1]-v1[1])*(v2[2]-v1[2])) / detJ;
     279                 :            : 
     280                 :            :   jacInv[0][1] = -(  (v3[0]-v1[0])*(v4[2]-v1[2])
     281                 :            :                    - (v4[0]-v1[0])*(v3[2]-v1[2])) / detJ;
     282                 :            :   jacInv[1][1] =  (  (v2[0]-v1[0])*(v4[2]-v1[2])
     283                 :            :                    - (v4[0]-v1[0])*(v2[2]-v1[2])) / detJ;
     284                 :            :   jacInv[2][1] = -(  (v2[0]-v1[0])*(v3[2]-v1[2])
     285                 :            :                    - (v3[0]-v1[0])*(v2[2]-v1[2])) / detJ;
     286                 :            : 
     287                 :            :   jacInv[0][2] =  (  (v3[0]-v1[0])*(v4[1]-v1[1])
     288                 :            :                    - (v4[0]-v1[0])*(v3[1]-v1[1])) / detJ;
     289                 :            :   jacInv[1][2] = -(  (v2[0]-v1[0])*(v4[1]-v1[1])
     290                 :            :                    - (v4[0]-v1[0])*(v2[1]-v1[1])) / detJ;
     291                 :            :   jacInv[2][2] =  (  (v2[0]-v1[0])*(v3[1]-v1[1])
     292                 :            :                    - (v3[0]-v1[0])*(v2[1]-v1[1])) / detJ;
     293                 :            : 
     294                 :            :   return jacInv;
     295                 :            : }
     296                 :            : 
     297                 :            : //! Compute the determinant of 3x3 matrix
     298                 :            : //!  \param[in] a 3x3 matrix
     299                 :            : //!  \return Determinant of the 3x3 matrix
     300                 :            : inline tk::real
     301                 :            : determinant( const std::array< std::array< tk::real, 3 >, 3 >& a )
     302                 :            : {
     303                 :            :   return ( a[0][0] * (a[1][1]*a[2][2]-a[1][2]*a[2][1])
     304                 :            :          - a[0][1] * (a[1][0]*a[2][2]-a[1][2]*a[2][0])
     305                 :            :          + a[0][2] * (a[1][0]*a[2][1]-a[1][1]*a[2][0]) );
     306                 :            : }
     307                 :            : 
     308                 :            : //! Solve a 3x3 system of equations using Cramer's rule
     309                 :            : //!  \param[in] a 3x3 lhs matrix
     310                 :            : //!  \param[in] b 3x1 rhs matrix
     311                 :            : //!  \return Array of solutions of the 3x3 system
     312                 :            : inline std::array < tk::real, 3 >
     313                 :            : cramer( const std::array< std::array< tk::real, 3 >, 3>& a,
     314                 :            :         const std::array< tk::real, 3 >& b )
     315                 :            : {
     316                 :            :   auto de = determinant( a );
     317                 :            : 
     318                 :            :   auto nu(0.0);
     319                 :            :   std::array < real, 3 > x;
     320                 :            : 
     321                 :            :   nu = determinant( {{{{b[0], a[0][1], a[0][2]}},
     322                 :            :                       {{b[1], a[1][1], a[1][2]}},
     323                 :            :                       {{b[2], a[2][1], a[2][2]}}}} );
     324                 :            :   x[0] = nu/de;
     325                 :            : 
     326                 :            :   nu = determinant( {{{{a[0][0], b[0], a[0][2]}},
     327                 :            :                       {{a[1][0], b[1], a[1][2]}},
     328                 :            :                       {{a[2][0], b[2], a[2][2]}}}} );
     329                 :            :   x[1] = nu/de;
     330                 :            : 
     331                 :            :   nu = determinant( {{{{a[0][0], a[0][1], b[0]}},
     332                 :            :                       {{a[1][0], a[1][1], b[1]}},
     333                 :            :                       {{a[2][0], a[2][1], b[2]}}}} );
     334                 :            :   x[2] = nu/de;
     335                 :            : 
     336                 :            :   return x;
     337                 :            : }
     338                 :            : 
     339                 :            : } // tk::
     340                 :            : 
     341                 :            : #endif // Vector_h

Generated by: LCOV version 1.14