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