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