1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
// *****************************************************************************
/*!
  \file      src/Inciter/AMR/Error.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.,
             2022-2025 J. Bakosi
             All rights reserved. See the LICENSE file for details.
  \brief     Class for computing error estimates for mesh refinement
  \details   Class for computing error estimates for mesh refinement.
*/
// *****************************************************************************

#include <cmath><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.
#include <limits><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.

#include "Exception.hpp"
#include "Error.hpp"
#include "Vector.hpp"
#include "Gradients.hpp"<--- Include file: "Gradients.hpp" not found.

using AMR::Error;

tk::real
Error::scalar( const tk::Fields& u,
               const edge_t& edge,
               uint64_t c,
               const std::array< std::vector< tk::real >, 3 >& coord,
               const std::vector< std::size_t >& inpoel,
               const std::pair< std::vector< std::size_t >,
                                std::vector< std::size_t > >& esup,
               const std::string& err ) const
// *****************************************************************************
//  Estimate error for scalar quantity
//! \param[in] u Solution vector
//! \param[in] edge Edge defined by its two end-point IDs
//! \param[in] c Scalar component to compute error of
//! \param[in] coord Mesh node coordinates
//! \param[in] inpoel Mesh element connectivity
//! \param[in] esup Linked lists storing elements surrounding points, see
//!    tk::genEsup()
//! \param[in] err AMR Error indicator type
//! \return Error indicator: a real number between [0...1] inclusive
// *****************************************************************************
{
  if (err == "jump") {
    return error_jump( u, edge, c );
  } else if (err == "hessian") {
    return error_hessian( u, edge, c, coord, inpoel, esup );
  } else {
    Throw( "No such AMR error indicator type" );
  }
}

tk::real
Error::error_jump( const tk::Fields& u,
                   const edge_t& edge,
                   uint64_t c ) const
// *****************************************************************************
//  Estimate error for scalar quantity on edge based on jump in solution
//! \param[in] u Solution vector
//! \param[in] edge Edge defined by its two end-point IDs
//! \param[in] c Scalar component to compute error of
//! \return Error indicator: a real number between [0...1] inclusive
// *****************************************************************************
{
  auto a = edge.first();
  auto b = edge.second();

  // Normalization factor with a noise filter
  auto norm = std::abs(u(a,c)) + std::abs(u(b,c)) + 1e-3;

  // Normalized error
  return std::abs( u(a,c) - u(b,c) ) / norm;
}

tk::real
Error::error_hessian( const tk::Fields& u,
                      const edge_t& edge,
                      uint64_t c,
                      const std::array< std::vector< tk::real >, 3 >& coord,
                      const std::vector< std::size_t >& inpoel,
                      const std::pair< std::vector< std::size_t >,
                                       std::vector< std::size_t > >& esup )
const
// *****************************************************************************
//  Estimate error for scalar quantity on edge based on Hessian of solution
//! \param[in] u Solution vector
//! \param[in] edge Edge defined by its two end-point IDs
//! \param[in] c Scalar component to compute error of
//! \param[in] coord Mesh node coordinates
//! \param[in] inpoel Mesh element connectivity
//! \param[in] esup Linked lists storing elements surrounding points, see
//!    tk::genEsup()
//! \return Error indicator: a real number between [0...1] inclusive
// *****************************************************************************
{
  const tk::real small = std::numeric_limits< tk::real >::epsilon();

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];
  auto a = edge.first();
  auto b = edge.second();

  // Compute edge vector
  std::array< tk::real, 3 > h {{ x[a]-x[b], y[a]-y[b], z[a]-z[b] }};

  // Compute gradients at edge-end points
  auto ga = nodegrad( a, coord, inpoel, esup, u, c );
  auto gb = nodegrad( b, coord, inpoel, esup, u, c );

  // Compute dot products of gradients and edge vectors
  auto dua = tk::dot( ga, h );
  auto dub = tk::dot( gb, h );

  // If the normalization factor is zero, return zero error
  auto norm = std::abs(dua) + std::abs(dub);
  if (norm < small) return 0.0;

  return std::abs(dub-dua) / norm;
}