Line data    Source code

       1              : //
       2              : // File FieldOperations
       3              : //   Differential operators for fields
       4              : //
       5              : 
       6              : namespace ippl {
       7              :     /*!
       8              :      * User interface of gradient
       9              :      * @param u field
      10              :      */
      11              :     template <typename Field>
      12           48 :     detail::meta_grad<Field> grad(Field& u) {
      13           48 :         constexpr unsigned Dim = Field::dim;
      14              : 
      15           48 :         u.fillHalo();
      16           48 :         BConds<Field, Dim>& bcField = u.getFieldBC();
      17           48 :         bcField.apply(u);
      18              : 
      19              :         using mesh_type   = typename Field::Mesh_t;
      20              :         using vector_type = typename mesh_type::vector_type;
      21              : 
      22           48 :         mesh_type& mesh = u.get_mesh();
      23          432 :         vector_type vectors[Dim];
      24          216 :         for (unsigned d = 0; d < Dim; d++) {
      25          168 :             vectors[d]    = 0;
      26          168 :             vectors[d][d] = 0.5 / mesh.getMeshSpacing(d);
      27              :         }
      28           96 :         return detail::meta_grad<Field>(u, vectors);
      29          216 :     }
      30              : 
      31              :     /*!
      32              :      * User interface of divergence in three dimensions.
      33              :      * @param u field
      34              :      */
      35              :     template <typename Field>
      36           48 :     detail::meta_div<Field> div(Field& u) {
      37           48 :         constexpr unsigned Dim = Field::dim;
      38              : 
      39           48 :         u.fillHalo();
      40           48 :         BConds<Field, Dim>& bcField = u.getFieldBC();
      41           48 :         bcField.apply(u);
      42              : 
      43              :         using mesh_type   = typename Field::Mesh_t;
      44              :         using vector_type = typename mesh_type::vector_type;
      45              : 
      46           48 :         mesh_type& mesh = u.get_mesh();
      47          432 :         vector_type vectors[Dim];
      48          216 :         for (unsigned d = 0; d < Dim; d++) {
      49          168 :             vectors[d]    = 0;
      50          168 :             vectors[d][d] = 0.5 / mesh.getMeshSpacing(d);
      51              :         }
      52           96 :         return detail::meta_div<Field>(u, vectors);
      53          216 :     }
      54              : 
      55              :     /*!
      56              :      * User interface of Laplacian
      57              :      * @param u field
      58              :      */
      59              :     template <typename Field>
      60           48 :     detail::meta_laplace<Field> laplace(Field& u) {
      61           48 :         constexpr unsigned Dim = Field::dim;
      62              : 
      63           48 :         u.fillHalo();
      64           48 :         BConds<Field, Dim>& bcField = u.getFieldBC();
      65           48 :         bcField.apply(u);
      66              : 
      67              :         using mesh_type = typename Field::Mesh_t;
      68           48 :         mesh_type& mesh = u.get_mesh();
      69           48 :         typename mesh_type::vector_type hvector(0);
      70          216 :         for (unsigned d = 0; d < Dim; d++) {
      71          168 :             hvector[d] = 1.0 / std::pow(mesh.getMeshSpacing(d), 2);
      72              :         }
      73           96 :         return detail::meta_laplace<Field>(u, hvector);
      74           48 :     }
      75              : 
      76              :     /*!
      77              :      * User interface of curl in three dimensions.
      78              :      * @param u field
      79              :      */
      80              :     template <typename Field>
      81            8 :     detail::meta_curl<Field> curl(Field& u) {
      82            8 :         constexpr unsigned Dim = Field::dim;
      83              : 
      84            8 :         u.fillHalo();
      85            8 :         BConds<Field, Dim>& bcField = u.getFieldBC();
      86            8 :         bcField.apply(u);
      87              : 
      88              :         using mesh_type = typename Field::Mesh_t;
      89            8 :         mesh_type& mesh = u.get_mesh();
      90            8 :         typename mesh_type::vector_type xvector(0);
      91            8 :         xvector[0] = 1.0;
      92            8 :         typename mesh_type::vector_type yvector(0);
      93            8 :         yvector[1] = 1.0;
      94            8 :         typename mesh_type::vector_type zvector(0);
      95            8 :         zvector[2] = 1.0;
      96            8 :         typename mesh_type::vector_type hvector(0);
      97            8 :         hvector = mesh.getMeshSpacing();
      98           16 :         return detail::meta_curl<Field>(u, xvector, yvector, zvector, hvector);
      99            8 :     }
     100              : 
     101              :     /*!
     102              :      * User interface of Hessian in three dimensions.
     103              :      * @param u field
     104              :      */
     105              :     template <typename Field>
     106           48 :     detail::meta_hess<Field> hess(Field& u) {
     107           48 :         constexpr unsigned Dim = Field::dim;
     108              : 
     109           48 :         u.fillHalo();
     110           48 :         BConds<Field, Dim>& bcField = u.getFieldBC();
     111           48 :         bcField.apply(u);
     112              : 
     113              :         using mesh_type   = typename Field::Mesh_t;
     114              :         using vector_type = typename mesh_type::vector_type;
     115              : 
     116           48 :         mesh_type& mesh = u.get_mesh();
     117          432 :         vector_type vectors[Dim];
     118          216 :         for (unsigned d = 0; d < Dim; d++) {
     119          168 :             vectors[d]    = 0;
     120          168 :             vectors[d][d] = 1;
     121              :         }
     122           48 :         auto hvector = mesh.getMeshSpacing();
     123              : 
     124           96 :         return detail::meta_hess<Field>(u, vectors, hvector);
     125          264 :     }
     126              : }  // namespace ippl