Branch data     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