text:
// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "Framework/TestingFramework.hpp"

#include <array>
#include <boost/functional/hash.hpp>
#include <cstddef>
#include <unordered_map>
#include <utility>

#include "DataStructures/DataVector.hpp"
#include "DataStructures/Tags.hpp"
#include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "Domain/CoordinateMaps/Affine.hpp"
#include "Domain/CoordinateMaps/CoordinateMap.hpp"
#include "Domain/CoordinateMaps/CoordinateMap.tpp"
#include "Domain/CoordinateMaps/ProductMaps.hpp"
#include "Domain/CoordinateMaps/ProductMaps.tpp"
#include "Domain/ElementMap.hpp"
#include "Domain/LogicalCoordinates.hpp"
#include "Domain/SizeOfElement.hpp"
#include "Domain/Structure/Direction.hpp"
#include "Domain/Structure/Element.hpp"
#include "Domain/Structure/ElementId.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/Minmod.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/MinmodImpl.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp"
#include "Evolution/Systems/NewtonianEuler/Limiters/CharacteristicHelpers.hpp"
#include "Evolution/Systems/NewtonianEuler/Limiters/Minmod.hpp"
#include "Evolution/Systems/NewtonianEuler/Tags.hpp"
#include "Framework/TestCreation.hpp"
#include "Framework/TestHelpers.hpp"
#include "Helpers/Evolution/DiscontinuousGalerkin/Limiters/TestHelpers.hpp"
#include "NumericalAlgorithms/Spectral/Mesh.hpp"
#include "NumericalAlgorithms/Spectral/Spectral.hpp"
#include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeArray.hpp"
#include "Utilities/TMPL.hpp"
#include "Utilities/TaggedTuple.hpp"

namespace {

void test_neuler_minmod_option_parsing() noexcept {
  INFO("Testing option parsing");
  const auto lambda_pi1 =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: LambdaPi1\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");
  const auto lambda_pi1_cons =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: LambdaPi1\n"
          "VariablesToLimit: Conserved\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");
  const auto lambda_pi1_tvb =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: LambdaPi1\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 1.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");
  const auto lambda_pi1_noflat =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: LambdaPi1\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: False\n"
          "DisableForDebugging: False");
  const auto lambda_pi1_disabled =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: LambdaPi1\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: True");
  const auto muscl =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<1>>(
          "Type: Muscl\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");

  // Test operators == and !=
  CHECK(lambda_pi1 == lambda_pi1);
  CHECK(lambda_pi1 != lambda_pi1_cons);
  CHECK(lambda_pi1 != lambda_pi1_tvb);
  CHECK(lambda_pi1 != lambda_pi1_noflat);
  CHECK(lambda_pi1 != lambda_pi1_disabled);
  CHECK(lambda_pi1 != muscl);

  const auto lambda_pin_2d =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<2>>(
          "Type: LambdaPiN\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 0.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");
  const auto lambda_pin_3d =
      TestHelpers::test_creation<NewtonianEuler::Limiters::Minmod<3>>(
          "Type: LambdaPiN\n"
          "VariablesToLimit: Characteristic\n"
          "TvbConstant: 10.0\n"
          "ApplyFlattener: True\n"
          "DisableForDebugging: False");

  // Test that creation from options gives correct object
  const NewtonianEuler::Limiters::Minmod<1> expected_lambda_pi1(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 0., true);
  const NewtonianEuler::Limiters::Minmod<1> expected_lambda_pi1_cons(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Conserved, 0., true);
  const NewtonianEuler::Limiters::Minmod<1> expected_lambda_pi1_tvb(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 1., true);
  const NewtonianEuler::Limiters::Minmod<1> expected_lambda_pi1_noflat(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 0., false);
  const NewtonianEuler::Limiters::Minmod<1> expected_lambda_pi1_disabled(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 0., true,
      true);
  const NewtonianEuler::Limiters::Minmod<1> expected_muscl(
      Limiters::MinmodType::Muscl,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 0., true);
  const NewtonianEuler::Limiters::Minmod<2> expected_lambda_pin_2d(
      Limiters::MinmodType::LambdaPiN,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 0., true);
  const NewtonianEuler::Limiters::Minmod<3> expected_lambda_pin_3d(
      Limiters::MinmodType::LambdaPiN,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 10., true);

  CHECK(lambda_pi1 == expected_lambda_pi1);
  CHECK(lambda_pi1_cons == expected_lambda_pi1_cons);
  CHECK(lambda_pi1_tvb == expected_lambda_pi1_tvb);
  CHECK(lambda_pi1_noflat == expected_lambda_pi1_noflat);
  CHECK(lambda_pi1_disabled == expected_lambda_pi1_disabled);
  CHECK(muscl == expected_muscl);
  CHECK(lambda_pin_2d == expected_lambda_pin_2d);
  CHECK(lambda_pin_3d == expected_lambda_pin_3d);
}

void test_neuler_minmod_serialization() noexcept {
  INFO("Testing serialization");
  const NewtonianEuler::Limiters::Minmod<1> minmod(
      Limiters::MinmodType::LambdaPi1,
      NewtonianEuler::Limiters::VariablesToLimit::Characteristic, 1., true);
  test_serialization(minmod);
}

// Compare the specialized limiter to the generic limiter. The limiters are
// compared for limiting conserved/evolved and limiting characteristic
// variables.
//
// Note: comparing the specialized and generic limiters when limiting the
// characteristic variables is tedious to do thoroughly, because we need to
// call the generic limiter multiple times using characteristic fields w.r.t.
// multiple different unit vectors. In essence, the test would need to fully
// reimplement the specialized limiter. To avoid the test becoming so tedious,
// we test a simpler restricted problem, where in 2D and 3D we require
// input_momentum to be 0. With this simplification, the characteristic
// decomposition becomes independent of the choice of unit vector, and therefore
// the generic limiter can match the specialized limiter even when only one set
// of characteristic (e.g., w.r.t. \hat{x}) is used.
template <size_t VolumeDim>
void test_neuler_vs_generic_minmod_work(
    const Scalar<DataVector>& input_density,
    const tnsr::I<DataVector, VolumeDim>& input_momentum,
    const Scalar<DataVector>& input_energy, const Mesh<VolumeDim>& mesh,
    const tnsr::I<DataVector, VolumeDim, Frame::Logical>& logical_coords,
    const std::array<double, VolumeDim>& element_size,
    const Scalar<DataVector>& det_inv_logical_to_inertial_jacobian,
    const std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>,
        typename NewtonianEuler::Limiters::Minmod<VolumeDim>::PackagedData,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
        neighbor_data) noexcept {
  // Sanity check that input momentum satisfies simplifying assumptions
  if constexpr (VolumeDim > 1) {
    const auto zero_momentum = make_with_value<tnsr::I<DataVector, VolumeDim>>(
        get<0>(input_momentum), 0.);
    REQUIRE(input_momentum == zero_momentum);
  }

  const double tvb_constant = 0.;
  const auto element = TestHelpers::Limiters::make_element<VolumeDim>();
  const EquationsOfState::IdealFluid<false> equation_of_state{5. / 3.};

  {
    INFO("Compare specialized vs generic Minmod on conserved variables");
    auto density_generic = input_density;
    auto momentum_generic = input_momentum;
    auto energy_generic = input_energy;
    const Limiters::Minmod<
        VolumeDim, tmpl::list<NewtonianEuler::Tags::MassDensityCons,
                              NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
                              NewtonianEuler::Tags::EnergyDensity>>
        minmod_generic(Limiters::MinmodType::LambdaPi1, tvb_constant);
    const bool activated_generic = minmod_generic(
        make_not_null(&density_generic), make_not_null(&momentum_generic),
        make_not_null(&energy_generic), mesh, element, logical_coords,
        element_size, neighbor_data);

    auto density_specialized = input_density;
    auto momentum_specialized = input_momentum;
    auto energy_specialized = input_energy;
    const auto vars_to_limit =
        NewtonianEuler::Limiters::VariablesToLimit::Conserved;
    const bool apply_flattener = false;
    const NewtonianEuler::Limiters::Minmod<VolumeDim> minmod_specialized(
        Limiters::MinmodType::LambdaPi1, vars_to_limit, tvb_constant,
        apply_flattener);
    const bool activated_specialized = minmod_specialized(
        make_not_null(&density_specialized),
        make_not_null(&momentum_specialized),
        make_not_null(&energy_specialized), mesh, element, logical_coords,
        element_size, det_inv_logical_to_inertial_jacobian, equation_of_state,
        neighbor_data);

    CHECK(activated_generic);
    CHECK(activated_generic == activated_specialized);
    CHECK_ITERABLE_APPROX(density_generic, density_specialized);
    CHECK_ITERABLE_APPROX(momentum_generic, momentum_specialized);
    CHECK_ITERABLE_APPROX(energy_generic, energy_specialized);
  }

  {
    INFO("Compare specialized vs generic Minmod on characteristic variables");
    // Cellwise means, for cons/char transforms
    const auto mean_density =
        Scalar<double>{mean_value(get(input_density), mesh)};
    const auto mean_momentum = [&input_momentum, &mesh]() noexcept {
      tnsr::I<double, VolumeDim> result{};
      for (size_t i = 0; i < VolumeDim; ++i) {
        result.get(i) = mean_value(input_momentum.get(i), mesh);
      }
      return result;
    }();
    const auto mean_energy =
        Scalar<double>{mean_value(get(input_energy), mesh)};

    // Compute characteristic transformation using x-direction...
    // Note that in 2D and 3D we know the fluid velocity is 0, so this choice
    // is arbitrary in these higher dimensions.
    const auto unit_vector = []() noexcept {
      auto components = make_array<VolumeDim>(0.);
      components[0] = 1.;
      return tnsr::i<double, VolumeDim>(components);
    }();
    const auto right_and_left =
        NewtonianEuler::Limiters::right_and_left_eigenvectors(
            mean_density, mean_momentum, mean_energy, equation_of_state,
            unit_vector);
    const auto& right = right_and_left.first;
    const auto& left = right_and_left.second;

    // Transform tensors to characteristics
    auto char_v_minus = input_density;
    auto char_v_momentum = input_momentum;
    auto char_v_plus = input_energy;
    NewtonianEuler::Limiters::characteristic_fields(
        make_not_null(&char_v_minus), make_not_null(&char_v_momentum),
        make_not_null(&char_v_plus), input_density, input_momentum,
        input_energy, left);

    using GenericMinmod =
        Limiters::Minmod<VolumeDim,
                         tmpl::list<NewtonianEuler::Tags::VMinus,
                                    NewtonianEuler::Tags::VMomentum<VolumeDim>,
                                    NewtonianEuler::Tags::VPlus>>;
    std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>,
        typename GenericMinmod::PackagedData,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>
        neighbor_char_data{};
    for (const auto& [key, data] : neighbor_data) {
      neighbor_char_data[key].element_size = data.element_size;
      NewtonianEuler::Limiters::characteristic_fields(
          make_not_null(&(neighbor_char_data[key].means)), data.means, left);
    }

    const GenericMinmod minmod_generic(Limiters::MinmodType::LambdaPi1,
                                       tvb_constant);
    const bool activated_generic = minmod_generic(
        make_not_null(&char_v_minus), make_not_null(&char_v_momentum),
        make_not_null(&char_v_plus), mesh, element, logical_coords,
        element_size, neighbor_char_data);

    // Transform back to evolved variables
    // Note that because the fluid velocity is 0 and all sets of characteristics
    // are identical, we can skip the step of combining different limiter
    // results.
    auto density_generic = input_density;
    auto momentum_generic = input_momentum;
    auto energy_generic = input_energy;
    NewtonianEuler::Limiters::conserved_fields_from_characteristic_fields(
        make_not_null(&density_generic), make_not_null(&momentum_generic),
        make_not_null(&energy_generic), char_v_minus, char_v_momentum,
        char_v_plus, right);

    auto density_specialized = input_density;
    auto momentum_specialized = input_momentum;
    auto energy_specialized = input_energy;
    const auto vars_to_limit =
        NewtonianEuler::Limiters::VariablesToLimit::Characteristic;
    const bool apply_flattener = false;
    const NewtonianEuler::Limiters::Minmod<VolumeDim> minmod_specialized(
        Limiters::MinmodType::LambdaPi1, vars_to_limit, tvb_constant,
        apply_flattener);
    const bool activated_specialized = minmod_specialized(
        make_not_null(&density_specialized),
        make_not_null(&momentum_specialized),
        make_not_null(&energy_specialized), mesh, element, logical_coords,
        element_size, det_inv_logical_to_inertial_jacobian, equation_of_state,
        neighbor_data);

    CHECK(activated_generic);
    CHECK(activated_generic == activated_specialized);
    CHECK_ITERABLE_APPROX(density_generic, density_specialized);
    CHECK_ITERABLE_APPROX(momentum_generic, momentum_specialized);
    CHECK_ITERABLE_APPROX(energy_generic, energy_specialized);
  }
}

void test_neuler_vs_generic_minmod_1d() noexcept {
  INFO("Testing NewtonianEuler::Limiters::Minmod limiter in 1D");
  const auto mesh =
      Mesh<1>(3, Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto);
  const auto element = TestHelpers::Limiters::make_element<1>();
  const auto logical_coords = logical_coordinates(mesh);
  using Affine = domain::CoordinateMaps::Affine;
  const Affine xi_map{-1., 1., 3.7, 4.2};
  const auto coordmap =
      domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
          xi_map);
  const ElementMap<1, Frame::Inertial> element_map(element.id(),
                                                   coordmap->get_clone());
  const auto element_size = size_of_element(element_map);
  const auto det_inv_logical_to_inertial_jacobian =
      determinant(element_map.inv_jacobian(logical_coords));

  const auto& x = get<0>(logical_coords);
  const auto mass_density_cons = [&x]() noexcept {
    return Scalar<DataVector>{{{1. + 0.2 * x + 0.05 * square(x)}}};
  }();
  const auto momentum_density = [&x]() noexcept {
    return tnsr::I<DataVector, 1>{{{0.2 - 0.3 * x}}};
  }();
  const auto energy_density = [&mesh]() noexcept {
    return Scalar<DataVector>{DataVector(mesh.number_of_grid_points(), 1.)};
  }();

  std::unordered_map<std::pair<Direction<1>, ElementId<1>>,
                     NewtonianEuler::Limiters::Minmod<1>::PackagedData,
                     boost::hash<std::pair<Direction<1>, ElementId<1>>>>
      neighbor_data{};
  const std::array<std::pair<Direction<1>, ElementId<1>>, 2> dir_keys = {
      {{Direction<1>::lower_xi(), ElementId<1>(1)},
       {Direction<1>::upper_xi(), ElementId<1>(2)}}};
  for (const auto& dir_key : dir_keys) {
    neighbor_data[dir_key].element_size = element_size;
  }

  using mean_rho = Tags::Mean<NewtonianEuler::Tags::MassDensityCons>;
  get(get<mean_rho>(neighbor_data[dir_keys[0]].means)) = 0.4;
  get(get<mean_rho>(neighbor_data[dir_keys[1]].means)) = 1.1;

  using mean_rhou = Tags::Mean<NewtonianEuler::Tags::MomentumDensity<1>>;
  get<0>(get<mean_rhou>(neighbor_data[dir_keys[0]].means)) = 0.2;
  get<0>(get<mean_rhou>(neighbor_data[dir_keys[1]].means)) = -0.1;

  using mean_eps = Tags::Mean<NewtonianEuler::Tags::EnergyDensity>;
  get(get<mean_eps>(neighbor_data[dir_keys[0]].means)) = 1.;
  get(get<mean_eps>(neighbor_data[dir_keys[1]].means)) = 1.;

  test_neuler_vs_generic_minmod_work(
      mass_density_cons, momentum_density, energy_density, mesh, logical_coords,
      element_size, det_inv_logical_to_inertial_jacobian, neighbor_data);
}

void test_neuler_vs_generic_minmod_2d() noexcept {
  INFO("Testing NewtonianEuler::Limiters::Minmod limiter in 2D");
  const auto mesh =
      Mesh<2>(3, Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto);
  const auto element = TestHelpers::Limiters::make_element<2>();
  const auto logical_coords = logical_coordinates(mesh);
  using Affine = domain::CoordinateMaps::Affine;
  using Affine2D = domain::CoordinateMaps::ProductOf2Maps<Affine, Affine>;
  const Affine xi_map{-1., 1., 3.7, 4.2};
  const Affine eta_map{-1., 1., 3.2, 4.2};
  const Affine2D map2d{xi_map, eta_map};
  const auto coordmap =
      domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
          map2d);
  const ElementMap<2, Frame::Inertial> element_map(element.id(),
                                                   coordmap->get_clone());
  const auto element_size = size_of_element(element_map);
  const auto det_inv_logical_to_inertial_jacobian =
      determinant(element_map.inv_jacobian(logical_coords));

  const auto& x = get<0>(logical_coords);
  const auto& y = get<1>(logical_coords);
  const auto mass_density_cons = [&x, &y]() noexcept {
    return Scalar<DataVector>{
        {{1. + 0.2 * x - 0.1 * y + 0.05 * x * square(y)}}};
  }();
  const auto momentum_density = [&mesh]() noexcept {
    return tnsr::I<DataVector, 2>{DataVector(mesh.number_of_grid_points(), 0.)};
  }();
  const auto energy_density = [&mesh]() noexcept {
    return Scalar<DataVector>{DataVector(mesh.number_of_grid_points(), 1.)};
  }();

  std::unordered_map<std::pair<Direction<2>, ElementId<2>>,
                     NewtonianEuler::Limiters::Minmod<2>::PackagedData,
                     boost::hash<std::pair<Direction<2>, ElementId<2>>>>
      neighbor_data{};
  const std::array<std::pair<Direction<2>, ElementId<2>>, 4> dir_keys = {
      {{Direction<2>::lower_xi(), ElementId<2>(1)},
       {Direction<2>::upper_xi(), ElementId<2>(2)},
       {Direction<2>::lower_eta(), ElementId<2>(3)},
       {Direction<2>::upper_eta(), ElementId<2>(4)}}};
  for (const auto& dir_key : dir_keys) {
    neighbor_data[dir_key].element_size = element_size;
  }

  using mean_rho = Tags::Mean<NewtonianEuler::Tags::MassDensityCons>;
  get(get<mean_rho>(neighbor_data[dir_keys[0]].means)) = 0.4;
  get(get<mean_rho>(neighbor_data[dir_keys[1]].means)) = 1.1;
  get(get<mean_rho>(neighbor_data[dir_keys[2]].means)) = 2.1;
  get(get<mean_rho>(neighbor_data[dir_keys[3]].means)) = 0.9;

  for (const auto& dir_key : dir_keys) {
    using mean_rhou = Tags::Mean<NewtonianEuler::Tags::MomentumDensity<2>>;
    get<0>(get<mean_rhou>(neighbor_data[dir_key].means)) = 0.;
    get<1>(get<mean_rhou>(neighbor_data[dir_key].means)) = 0.;
    using mean_eps = Tags::Mean<NewtonianEuler::Tags::EnergyDensity>;
    get(get<mean_eps>(neighbor_data[dir_key].means)) = 1.;
  }

  test_neuler_vs_generic_minmod_work(
      mass_density_cons, momentum_density, energy_density, mesh, logical_coords,
      element_size, det_inv_logical_to_inertial_jacobian, neighbor_data);
}

void test_neuler_vs_generic_minmod_3d() noexcept {
  INFO("Testing NewtonianEuler::Limiters::Minmod limiter in 3D");
  const auto mesh =
      Mesh<3>(std::array<size_t, 3>{{3, 3, 4}}, Spectral::Basis::Legendre,
              Spectral::Quadrature::GaussLobatto);
  const auto element = TestHelpers::Limiters::make_element<3>();
  const auto logical_coords = logical_coordinates(mesh);
  using Affine = domain::CoordinateMaps::Affine;
  using Affine3D =
      domain::CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
  const Affine xi_map{-1., 1., 3.7, 4.2};
  const Affine eta_map{-1., 1., 3.2, 4.2};
  const Affine zeta_map{-1., 1., 1.3, 2.1};
  const Affine3D map3d{xi_map, eta_map, zeta_map};
  const auto coordmap =
      domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
          map3d);
  const ElementMap<3, Frame::Inertial> element_map(element.id(),
                                                   coordmap->get_clone());
  const auto element_size = size_of_element(element_map);
  const auto det_inv_logical_to_inertial_jacobian =
      determinant(element_map.inv_jacobian(logical_coords));

  const auto& x = get<0>(logical_coords);
  const auto& y = get<1>(logical_coords);
  const auto& z = get<2>(logical_coords);
  const auto mass_density_cons = [&x, &y, &z]() noexcept {
    return Scalar<DataVector>{{{1. + 0.2 * x - 0.1 * y + 0.4 * z}}};
  }();
  const auto momentum_density = [&mesh]() noexcept {
    return tnsr::I<DataVector, 3>{DataVector(mesh.number_of_grid_points(), 0.)};
  }();
  const auto energy_density = [&mesh]() noexcept {
    return Scalar<DataVector>{DataVector(mesh.number_of_grid_points(), 1.)};
  }();

  std::unordered_map<std::pair<Direction<3>, ElementId<3>>,
                     NewtonianEuler::Limiters::Minmod<3>::PackagedData,
                     boost::hash<std::pair<Direction<3>, ElementId<3>>>>
      neighbor_data{};
  const std::array<std::pair<Direction<3>, ElementId<3>>, 6> dir_keys = {
      {{Direction<3>::lower_xi(), ElementId<3>(1)},
       {Direction<3>::upper_xi(), ElementId<3>(2)},
       {Direction<3>::lower_eta(), ElementId<3>(3)},
       {Direction<3>::upper_eta(), ElementId<3>(4)},
       {Direction<3>::lower_zeta(), ElementId<3>(5)},
       {Direction<3>::upper_zeta(), ElementId<3>(6)}}};
  for (const auto& dir_key : dir_keys) {
    neighbor_data[dir_key].element_size = element_size;
  }

  using mean_rho = Tags::Mean<NewtonianEuler::Tags::MassDensityCons>;
  get(get<mean_rho>(neighbor_data[dir_keys[0]].means)) = 0.;
  get(get<mean_rho>(neighbor_data[dir_keys[1]].means)) = 0.;
  get(get<mean_rho>(neighbor_data[dir_keys[2]].means)) = 0.1;
  get(get<mean_rho>(neighbor_data[dir_keys[3]].means)) = 1.2;
  get(get<mean_rho>(neighbor_data[dir_keys[4]].means)) = 1.1;
  get(get<mean_rho>(neighbor_data[dir_keys[5]].means)) = 0.9;

  for (const auto& dir_key : dir_keys) {
    using mean_rhou = Tags::Mean<NewtonianEuler::Tags::MomentumDensity<3>>;
    get<0>(get<mean_rhou>(neighbor_data[dir_key].means)) = 0.;
    get<1>(get<mean_rhou>(neighbor_data[dir_key].means)) = 0.;
    get<2>(get<mean_rhou>(neighbor_data[dir_key].means)) = 0.;
    using mean_eps = Tags::Mean<NewtonianEuler::Tags::EnergyDensity>;
    get(get<mean_eps>(neighbor_data[dir_key].means)) = 1.;
  }

  test_neuler_vs_generic_minmod_work(
      mass_density_cons, momentum_density, energy_density, mesh, logical_coords,
      element_size, det_inv_logical_to_inertial_jacobian, neighbor_data);
}

template <size_t VolumeDim>
void test_neuler_minmod_flattener() noexcept {
  INFO("Testing flattener use in NewtonianEuler::Limiters::Minmod limiter");
  CAPTURE(VolumeDim);

  const auto mesh = Mesh<VolumeDim>(2, Spectral::Basis::Legendre,
                                    Spectral::Quadrature::GaussLobatto);
  // We use an element with no neighbors so the limiter does nothing
  const auto element = Element<VolumeDim>{ElementId<VolumeDim>(0), {}};
  const auto logical_coords = logical_coordinates(mesh);
  const auto element_size = make_array<VolumeDim>(0.6);
  // inv_jac = logical volume / inertial volume
  const Scalar<DataVector> det_inv_logical_to_inertial_jacobian{DataVector(
      mesh.number_of_grid_points(), pow<VolumeDim>(2.) / pow<VolumeDim>(0.6))};
  const EquationsOfState::IdealFluid<false> equation_of_state{5. / 3.};

  const size_t num_points = mesh.number_of_grid_points();
  const auto input_density = [&num_points]() noexcept {
    // One negative value to trigger flattener
    Scalar<DataVector> density{DataVector(num_points, 0.8)};
    get(density)[0] = -0.2;
    return density;
  }();
  const auto input_momentum = [&num_points]() noexcept {
    return tnsr::I<DataVector, VolumeDim>{DataVector(num_points, 0.1)};
  }();
  const auto input_energy = [&num_points]() noexcept {
    return Scalar<DataVector>{DataVector(num_points, 1.4)};
  }();

  // Empty map because no neighbors
  const std::unordered_map<
      std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>,
      typename NewtonianEuler::Limiters::Minmod<VolumeDim>::PackagedData,
      boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>
      neighbor_data{};

  // First a sanity check: there should be a negative density for this test to
  // make sense... otherwise we aren't testing anything useful
  const bool has_negative_density = min(get(input_density)) < 0.;
  REQUIRE(has_negative_density);
  CHECK(true);

  // With flattener off, check the limiter does nothing (because no neighbors).
  //
  // Note that the post-limiter result is unphysical: the density is still
  // negative at one point. Because of this, we have to use the generic limiter
  // when testing the limiter does nothing: the specialized limiter has an
  // ASSERT that the post-limiting solution is valid which would fail here.
  auto density = input_density;
  auto momentum = input_momentum;
  auto energy = input_energy;
  const double tvb_constant = 0.;
  const Limiters::Minmod<
      VolumeDim, tmpl::list<NewtonianEuler::Tags::MassDensityCons,
                            NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
                            NewtonianEuler::Tags::EnergyDensity>>
      minmod_noflat(Limiters::MinmodType::LambdaPi1, tvb_constant);
  const bool activated_noflat = minmod_noflat(
      make_not_null(&density), make_not_null(&momentum), make_not_null(&energy),
      mesh, element, logical_coords, element_size, neighbor_data);
  CHECK_FALSE(activated_noflat);

  // With flattener on, check limiter activates
  density = input_density;
  momentum = input_momentum;
  energy = input_energy;
  const auto vars_to_limit =
      NewtonianEuler::Limiters::VariablesToLimit::Conserved;
  const NewtonianEuler::Limiters::Minmod<VolumeDim> minmod(
      Limiters::MinmodType::LambdaPi1, vars_to_limit, tvb_constant, true);
  const bool activated = minmod(
      make_not_null(&density), make_not_null(&momentum), make_not_null(&energy),
      mesh, element, logical_coords, element_size,
      det_inv_logical_to_inertial_jacobian, equation_of_state, neighbor_data);
  CHECK(activated);

  // Finally, check the negative density is brought back to positive
  const bool density_is_positive_everywhere = min(get(density)) > 0.;
  CHECK(density_is_positive_everywhere);
}

}  // namespace

SPECTRE_TEST_CASE("Unit.Evolution.Systems.NewtonianEuler.Limiters.Minmod",
                  "[Limiters][Unit]") {
  test_neuler_minmod_option_parsing();
  test_neuler_minmod_serialization();

  // The package_data function for NewtonianEuler::Limiters::Minmod is just a
  // direct call to the generic Minmod package_data. We do not test it.

  // We test the action of NewtonianEuler::Limiters::Minmod by comparing it to
  // the action of the generic Minmod, for certain inputs where the comparison
  // is not too difficult. To fully test the specialized limiter, would need to
  // completely reimplement the logic within the test.
  test_neuler_vs_generic_minmod_1d();
  test_neuler_vs_generic_minmod_2d();
  test_neuler_vs_generic_minmod_3d();

  // Test the use of the post-processing flattener
  test_neuler_minmod_flattener<1>();
  test_neuler_minmod_flattener<2>();
  test_neuler_minmod_flattener<3>();
}

max_stars_repo_name: fmahebert/spectre

max_issues_repo_name: fmahebert/spectre

max_forks_repo_name: fmahebert/spectre


text:
#include <boost/tokenizer.hpp>
#include <fstream>
#include <stdlib.h>
#include "test/TestUtils.hpp"

namespace tsdb{
namespace test{

std::deque<std::deque<double>> load_sample_data(const std::string & filename){
    std::deque<std::deque<double>> r;
    std::ifstream file(filename);
    std::string line;
    boost::char_separator<char> tokenSep(",");
    while(getline(file, line)){
        r.push_back(std::deque<double>());
        boost::tokenizer<boost::char_separator<char> > tokens(line, tokenSep);
        for(auto const& token: tokens){
            r.back().push_back(atof(token.c_str()));
        }
    }
    return r;
}

}
}

ext: cpp

max_stars_repo_name: Jimx-/tsdb-fork

max_issues_repo_name: Jimx-/tsdb-fork

max_forks_repo_name: Jimx-/tsdb-fork


text:
#include "socks5/tcp_socket.hpp"

#include <boost/asio/io_context.hpp>
#include <boost/asio/ip/tcp.hpp>
#include <boost/asio/read_until.hpp>
#include <boost/asio/streambuf.hpp>
#include <boost/system/error_code.hpp>

#include <cstdint>
#include <iostream>
#include <regex>
#include <sstream>
#include <stdexcept>
#include <string_view>
#include <tuple>

namespace {

using namespace std::literals;

inline auto parse_url(std::string_view url)
    -> std::tuple<std::string, std::string> {
    static std::regex r {
        R"(^(http://)?([a-zA-Z0-9\.\-]+[a-zA-Z0-9]\.[a-zA-Z0-9]+)(/.*)?$)"};

    std::string urlstr {url};
    std::smatch m;

    if (std::regex_search(urlstr, m, r)) {
        return std::make_tuple(m[2], m[3].length() > 0 ? m[3] : "/"s);
    } else {
        throw std::runtime_error {"Invalid url"};
    }
}

inline auto create_tor_socket(boost::asio::io_context &io)
    -> socks5::tcp_stream {
    socks5::tcp_stream sock {io};

    std::clog << "[*] Opening SOCKS5 connection to 127.0.0.1:9050" << std::endl;
    sock.socks5_open(
        boost::asio::ip::tcp::endpoint {boost::asio::ip::tcp::v4(), 9050});

    std::clog << "[*] Authenticating SOCKS5 with none" << std::endl;
    sock.socks5_auth();

    return sock;
}

inline auto make_get_request(std::string_view domain, std::string_view path)
    -> std::string {
    std::stringstream req {};

    req << "GET " << (path.size() > 0 ? path : "/") << " HTTP/1.1\r\n";
    req << "Host: " << domain << "\r\n";
    req << "Connection: close\r\n";
    req << "\r\n";

    return req.str();
}

inline auto read_response(socks5::tcp_stream &sock) -> std::string {
    std::stringstream              resp {};
    std::array<std::uint8_t, 0xFF> buf {};
    std::size_t                    tx;

    std::cout << resp.str() << std::endl;
    boost::system::error_code ec;

    while ((tx = boost::asio::read(sock, boost::asio::buffer(buf), ec))) {
        std::string str {buf.cbegin(), buf.cbegin() + tx};
        resp << str;
    }

    if (ec != boost::asio::error::eof) {
        throw boost::system::system_error {ec};
    }

    return resp.str();
}

inline void tor_http_get(std::string_view url) {
    boost::asio::io_context io;

    auto sock                 = create_tor_socket(io);
    auto [domain, path] = parse_url(url);

    std::clog << "[*] Connecting to: " << domain << std::endl;
    sock.connect(domain, 80);

    auto req = make_get_request(domain, path);

    std::clog << "[*] Sending GET request" << std::endl;
    boost::asio::write(sock, boost::asio::buffer(req));

    std::clog << "[*] Reading response" << std::endl;
    std::cout << read_response(sock) << std::endl;
}

} // namespace

int main(int argc, char *argv[]) {
    if (argc != 2) {
        std::cerr << "Usage: " << argv[0] << " url" << std::endl;
        return 1;
    }

    try {
        tor_http_get(argv[1]);
        return 0;
    } catch (const std::exception &ex) {
        std::cerr << "[!] " << ex.what() << std::endl;
        return 1;
    }
}


text:
//  process.hpp  --------------------------------------------------------------//

//  Copyright 2010 Vicente J. Botet Escriba

//  Distributed under the Boost Software License, Version 1.0.
//  See http://www.boost.org/LICENSE_1_0.txt

#ifndef BOOST_DETAIL_WIN_PROCESS_HPP
#define BOOST_DETAIL_WIN_PROCESS_HPP

#include <boost/detail/win/basic_types.hpp>
#include <boost/detail/win/GetCurrentProcess.hpp>

namespace boost {
namespace detail {
namespace win32 {
#if defined( BOOST_USE_WINDOWS_H )
    using ::GetCurrentProcessId;
#else
# ifndef UNDER_CE
extern "C" {
    __declspec(dllimport) unsigned long __stdcall
        GetCurrentProcessId(void);
}
# else
    using ::GetCurrentProcessId;
# endif
#endif
}
}
}
#endif // BOOST_DETAIL_WIN_PROCESS_HPP

ext: hpp

max_stars_repo_name: PXLVision/opengv

max_issues_repo_name: isuhao/klib

max_forks_repo_name: isuhao/klib


text:
#include <iostream>
#include <armadillo>
#include <cmath>
#include <time.h>
#include <fstream>

#include "../code-fredrik/comp_eig.hh"

using namespace std;
using namespace arma;

void RHO_A_FILL(vec &rho, mat &A, int N,double rhoN); //rho, kind of like a linespace
                                                 //A, Tridiagonal matrix
void Maxoff(mat &A, int N,int &k, int &l, double &max);        //Finds the max element of a matrix

int main(){
    int values[] = { 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100 };
    fstream Outfile;
    Outfile.open("steps.dat",ios::out);
    Outfile << "N     epsilon    steps  step_time  error"<<endl;
    for(int iter = 0; iter < 20; iter++){
        int N = values[iter]; //matrix size; N x N
        double rhoN = 6;
        double eps = 1E-10;
        double max;
        int k,l;

        double total_error = 0;
        double total_time = 0;
        int total_iterations = 0;

        // run 10 times and average
        int zcount = 10;
        for(int z = 0; z < zcount; z++) {
          std::cout << N << std::endl;
          mat A(N,N,fill::zeros); //indexes go from (0) to (N-1)
          mat S(N,N,fill::eye);
          vec rho(N,fill::zeros);
          vec eigen(N);

          RHO_A_FILL(rho,A,N,rhoN);
          Maxoff(A,N,k,l,max);
          //A.print();
          int iterations = 0;

          // solve with Armadillo's eig_sym
          vec arma_eigenvalues;
          mat arma_eigenvectors;
          eig_sym(arma_eigenvalues, arma_eigenvectors, A);

          double tau, t, s, c, il, ik, kk, ll, s_ik, s_il;
          double start, finish;

          start = clock();                      //clock value before eigen solve
          while(max > eps){
              tau = (A(l,l)-A(k,k))/(2.*A(k,l));
              if(tau>0){
                  t = 1.0/(tau + sqrt(1.0 + tau*tau));
              } else{
                  t = -1.0/( -tau + sqrt(1.0 + tau*tau));
              }

              //cosine and sine
              c = 1./sqrt(1.+t*t);
              s = t*c;

              //Jacobi rotating A round theta in N-dim space
              for(int i = 0; i<N; i++){
                  if ((i != k) && (i !=l)){

                      ik = A(i,k)*c - A(i,l)*s;
                      il = A(i,l)*c + A(i,k)*s;
                      A(i,k) = ik;
                      A(i,l) = il;
                      A(k,i) = ik;
                      A(l,i) = il;
                  }
                  s_ik = S(i,k);
                  s_il = S(i,l);
                  S(i,k) = c*s_ik - s*s_il;
                  S(i,l) = c*s_il + s*s_ik;
              }

              kk = A(k,k)*c*c - 2.*A(k,l)*c*s + A(l,l)*s*s;
              ll = A(l,l)*c*c + 2.*A(k,l)*c*s + A(k,k)*s*s;
              A(k,k) = kk;
              A(l,l) = ll ;
              A(k,l) = 0;
              A(l,k) = 0;

              iterations++;
              Maxoff(A,N,k,l,max);

          } //end of while
          finish = clock();                    //clock value after eigen solve

          //[eigen] is now eigenvalues
          for(int i = 0; i<N;i++){
              eigen(i) = A(i,i);
          }

          // eigen.print();

          // compare eigenvectors from own method with armadillo
          total_error += comp_eig(eigen, S, arma_eigenvalues, arma_eigenvectors);

          total_time += (finish -start)/CLOCKS_PER_SEC/iterations;
          total_iterations += iterations;

          /*
          fstream outfile;
          outfile.open("eigenvectors.dat",ios::out);

          for (int i = 0; i<N; i++){
              for (int j = 0; j<N; j++){
                  if(j%N == 0){outfile<<endl;}
                  outfile << S(i,j)<<" ";
                  }
              }
          outfile.close();
          */
        }

        Outfile <<N <<" "<<eps<<" "<<((double)total_iterations / zcount)<<" "<<(total_time / zcount) <<" "<<(total_error / zcount)<<endl;

    //    Outfile.close();
    }
} //end of main

void RHO_A_FILL(vec &rho, mat &A, int N,double rhoN){
    double h = rhoN/(N);
    for(int i = 0; i < N;i++){
        rho(i) = i*h;
        A(i,i) = 2./(h*h)+(rho(i)*rho(i));
        if(i<N-1){
            A(i+1,i) = -1./(h*h);
            A(i,i+1) = -1./(h*h);
            }
        }
    }

void Maxoff(mat &A, int N,int &k, int &l, double &max){
    max = 0;
    for(int i = 0; i < N ; i++){
        for(int j = i+1; j < N ; j++){
            if ( A(i,j)*A(i,j) > max){
                max =A(i,j)*A(i,j);
                k = i;
                l = j;
                }
            }
        }

    }


text:
/*
*
*      Filename: proj1.cpp
*
*        Author: Haibo Hao
*        Email : haohaibo@ncic.ac.cn
*   Description: ---
*        Create: 2017-05-22 14:15:11
* Last Modified: 2017-05-22 17:31:18
**/
#include <sys/syscall.h>
#include <unistd.h>
#include <boost/thread.hpp>
#include <iostream>
using namespace std;

void show() { cout << "hello world " << syscall(SYS_gettid) << endl; }
int main() {
  // 这里也可以用t(&show)也能成功和pthread_create一样
  boost::thread t(show);
  cout << "main pid " << syscall(SYS_gettid) << " " << t.get_id() << endl;
  // t.detach();
  if (!t.joinable()) {
    cout << "thread unjoinable" << endl;
  } else {
    cout << "thread joinable" << endl;
    t.join();
  }
  return 0;
}

max_stars_repo_name: haohaibo/learn

max_issues_repo_name: haohaibo/learn

max_forks_repo_name: haohaibo/learn


text:
#include "io/detail/ArhReader.hpp"

#include <fstream>
#include <stdint.h>
#include <boost/filesystem.hpp>

#include "Mesh.hpp"
#include "MeshStructure.hpp"
#include "Traverse.hpp"
#include "DOFVector.hpp"
#include "MacroElement.hpp"

namespace AMDiS
{
  namespace io
  {

    using namespace std;

    namespace ArhReader
    {
      namespace detail
      {

        void setDofValues(int macroElIndex, Mesh* mesh,
                          vector<double>& values, DOFVector<double>* vec)
        {
          bool macroElement = false;
          int valuePos = 0;
          TraverseStack stack;
          ElInfo* elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
                           Mesh::CALL_EVERY_EL_PREORDER);
          while (elInfo)
          {
            if (!macroElement)
            {
              Element* mEl = elInfo->getMacroElement()->getElement();
              for (int i = 0; i < mesh->getGeo(VERTEX); i++)
                (*vec)[mEl->getDof(i, vec->getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX))] = values[valuePos++];
              macroElement = true;
            }

            Element* el = elInfo->getElement();
            if (!el->isLeaf())
              (*vec)[el->getChild(0)->getDof(mesh->getDim(),
                                             vec->getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX))] = values[valuePos++];

            elInfo = stack.traverseNext(elInfo);
          }
        }

        void read(string filename,
                  Mesh* mesh,
                  vector<DOFVector<double>*> vecs)
        {
          FUNCNAME("ArhReader::detail::read()");

          // === Get set of all macro elements in mesh. ===
          std::set<int> macroInMesh;
          for (std::deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
               it != mesh->getMacroElements().end(); ++it)
            macroInMesh.insert((*it)->getIndex());

          RefinementManager* refManager = NULL;
          switch (mesh->getDim())
          {
          case 1:
            refManager = new RefinementManager1d();
            break;
          case 2:
            refManager = new RefinementManager2d();
            break;
          case 3:
            refManager = new RefinementManager3d();
            break;
          default:
            ERROR_EXIT("Should not happen!\n");
          }

          ifstream file;
          file.open(filename.c_str(), ios::in | ios::binary);
          TEST_EXIT(file.is_open())
          ("Cannot open file %s\n", filename.c_str());

          string typeId(4, ' ');
          uint32_t nMacroElements = 0;
          uint32_t nValueVectors = 0;
          uint32_t nAllValues = 0;

          file.read(const_cast<char*>(typeId.data()), 4);
          TEST_EXIT(typeId == "arhP" || typeId == "arhS")
          ("Cannot read Arh format: this file is not \"arhP\" or \"arhS\" format.\n");

          file.read(reinterpret_cast<char*>(&nMacroElements), 4);
          file.read(reinterpret_cast<char*>(&nValueVectors), 4);
          file.read(reinterpret_cast<char*>(&nAllValues), 4);

          TEST_EXIT(nValueVectors == vecs.size())
          ("File %s has %d vectors, but %d DOFVectors are provided!\n",
           filename.c_str(), nValueVectors, vecs.size());

          for (unsigned int i = 0; i < nMacroElements; i++)
          {
            uint32_t elIndex = 0;
            uint32_t nStructureCodes = 0;
            uint32_t codeSize = 0;

            file.read(reinterpret_cast<char*>(&elIndex), 4);
            file.read(reinterpret_cast<char*>(&nStructureCodes), 4);
            file.read(reinterpret_cast<char*>(&codeSize), 4);

            vector<uint64_t> structureCode(nStructureCodes);
            file.read(reinterpret_cast<char*>(&(structureCode[0])), 8 * nStructureCodes);

            MeshStructure elementStructure;
            elementStructure.init(structureCode, codeSize);
            if (macroInMesh.count(elIndex))
              elementStructure.fitMeshToStructure(mesh, refManager, false, elIndex);

            if (nValueVectors > 0)
            {
              uint32_t nValuesPerVector = 0;
              file.read(reinterpret_cast<char*>(&nValuesPerVector), 4);

              for (unsigned int j = 0; j < nValueVectors; j++)
              {
                vector<double> values(nValuesPerVector);
                file.read(reinterpret_cast<char*>(&(values[0])), 8 * nValuesPerVector);
                if (vecs[j] != NULL)
                {
                  if (macroInMesh.count(elIndex))
                    setDofValues(elIndex, mesh, values, vecs[j]);
                }
              }
            }
          }

          file.close();
          delete refManager;

          MSG("ARH file read from: %s\n", filename.c_str());
        }

        void readBlock(vector<char>& data,
                       Mesh* mesh,
                       vector<DOFVector<double>*> vecs)
        {
          FUNCNAME("ArhReader::detail::readBlock()");

          // === Get set of all macro elements in mesh. ===
          std::set<int> macroInMesh;
          for (std::deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
               it != mesh->getMacroElements().end(); ++it)
            macroInMesh.insert((*it)->getIndex());

          RefinementManager* refManager = NULL;
          switch (mesh->getDim())
          {
          case 2:
            refManager = new RefinementManager2d();
            break;
          case 3:
            refManager = new RefinementManager3d();
            break;
          default:
            ERROR_EXIT("Should not happen!\n");
          }

          string typeId = "";
          uint32_t nMacroElements = 0;
          uint32_t nValueVectors = 0;
          uint32_t nAllValues = 0;
          uint32_t mem_index = 0;

          memcpy(const_cast<char*>(typeId.data()), &data[mem_index], 4);
          mem_index += 4;
          memcpy(reinterpret_cast<char*>(&nMacroElements), &data[mem_index], 4);
          mem_index += 4;
          memcpy(reinterpret_cast<char*>(&nValueVectors), &data[mem_index], 4);
          mem_index += 4;
          memcpy(reinterpret_cast<char*>(&nAllValues), &data[mem_index], 4);
          mem_index += 4;

          TEST_EXIT(nValueVectors == vecs.size())
          ("data has %d vectors, but %d DOFVectors are provided!\n", nValueVectors, vecs.size());

          for (unsigned int i = 0; i < nMacroElements; i++)
          {
            uint32_t elIndex = 0;
            uint32_t nStructureCodes = 0;
            uint32_t codeSize = 0;

            memcpy(reinterpret_cast<char*>(&elIndex), &data[mem_index], 4);
            mem_index += 4;
            memcpy(reinterpret_cast<char*>(&nStructureCodes), &data[mem_index], 4);
            mem_index += 4;
            memcpy(reinterpret_cast<char*>(&codeSize), &data[mem_index], 4);
            mem_index += 4;

            vector<uint64_t> structureCode(nStructureCodes);
            memcpy(reinterpret_cast<char*>(&structureCode[0]), &data[mem_index], 8 * nStructureCodes);
            mem_index += 8 * nStructureCodes;

            MeshStructure elementStructure;
            elementStructure.init(structureCode, codeSize);
            if (macroInMesh.count(elIndex) == 1)
              elementStructure.fitMeshToStructure(mesh, refManager, false, elIndex);

            if(nValueVectors > 0)
            {
              uint32_t nValuesPerVector = 0;
              memcpy(reinterpret_cast<char*>(&nValuesPerVector), &data[mem_index], 4);
              mem_index += 4;

              for (unsigned int j = 0; j < nValueVectors; j++)
              {
                vector<double> values(nValuesPerVector);
                memcpy(reinterpret_cast<char*>(&values[0]), &data[mem_index], 8 * nValuesPerVector);
                mem_index += 8 * nValuesPerVector;
                if (vecs[j] != NULL)
                {
                  if (macroInMesh.count(elIndex) == 1)
                    setDofValues(elIndex, mesh, values, vecs[j]);
                }
              }
            }
          }

          delete refManager;
        }

        void readFile(string filename, Mesh* mesh,
                      vector<DOFVector<double>*> vecs,
                      bool writeParallel,
                      int nProcs)
        {
          FUNCNAME("ArhReader::detail::readFile()");

          for(size_t i = 0; i < vecs.size(); i++)
          {
            if(vecs[i])
            {
              if(!mesh)
                mesh = vecs[i]->getFeSpace()->getMesh();
              else
                TEST_EXIT(mesh == vecs[i]->getFeSpace()->getMesh())
                ("The mesh of the DOFVectors should be the same for Reader because in one file there is only one mesh.\n");
            }
          }
          if(!mesh)
          {
            WARNING("You haven't specified the target, no mesh or DOFVectors is given.\n");
            return;
          }

          if (writeParallel)
          {
            int sPos = filename.find(".arh");
            TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n");
            string name = filename.substr(0, sPos);

            if (nProcs == -1)
            {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
              string procFilename = name + "-p" + std::to_string(MPI::COMM_WORLD.Get_rank()) + "-.arh";
              read(procFilename, mesh, vecs);
#else
              ERROR_EXIT("Reading parallel ARH files in sequential computations requires to specify the number of nodes on which the ARH file was created!\n");
#endif
            }
            else
            {
              for (int i = 0; i < nProcs; i++)
              {
                string procFilename = name + "-p" + std::to_string(i) + "-.arh";
                read(procFilename, mesh, vecs);
              }
            }
          }
          else
          {
            read(filename, mesh, vecs);
          }
        }
      }//end namespace detail
    } // end namespace ArhReader
  }
} // end namespace io, AMDiS


text:
//  Distributed under the Boost Software License, Version 1.0.
//  (See accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)

#include <boost/core/lightweight_test.hpp>
#include <boost/outbuf.hpp>
#include "test_utils.hpp"

template <typename CharT>
void test_discarded_outbuf()
{
    boost::discarded_outbuf<CharT> dob;
    BOOST_TEST(!dob.good());

    auto tiny_str = test_utils::make_tiny_string<CharT>();
    auto double_str = test_utils::make_double_string<CharT>();

    write(dob, tiny_str.data(), tiny_str.size());
    write(dob, double_str.data(), double_str.size());

    dob.recycle();
    BOOST_TEST_GE(dob.size(), boost::min_size_after_recycle<CharT>());
    BOOST_TEST(dob.pos() == boost::outbuf_garbage_buf<CharT>());
    BOOST_TEST(dob.end() == boost::outbuf_garbage_buf_end<CharT>());
}

int main()
{
    test_discarded_outbuf<char>();
    test_discarded_outbuf<char16_t>();

    return boost::report_errors();
}


text:
#pragma once

//#define RIKFAN_DEBUG

#include <boost/asio/io_service.hpp>
#include <sdbusplus/asio/object_server.hpp>

#include "pid.hpp"
#include "util.hpp"

#include <vector>
#include <filesystem>
namespace fs = std::filesystem;

#ifdef RIKFAN_DEBUG
#include <unordered_map>
#endif

#ifdef RIKFAN_DEBUG

struct DBusReport
{
	int main_cnt;
    int nan_cnt;
    int val_cnt;
};

#endif // RIKFAN_DEBUG

class SensorBase
{
public:
    virtual double get_value() = 0;
    virtual void set_value(double in) = 0;
    virtual std::string repr() = 0;
    virtual bool init_complete() {return false;}
};

class SensorFS: public SensorBase
{
public:
    SensorFS(const std::string &p, double err_val);
    double get_value() override;
    void set_value(double in) override;
    std::string repr() override;
private:
    std::string real_path;
    int state;
    double error_value;
    std::string sensor_id;
    bool error_sensor_state;
};

class SensorDBus: public SensorBase
{
public:
    SensorDBus(const std::string &p, double err_val, std::shared_ptr<sdbusplus::asio::connection>& b);
    double get_value() override ;
    void set_value(double in) override {throw "Not implemented yet.";}
    std::string repr() override;

#ifdef RIKFAN_DEBUG
    void set_report(std::shared_ptr<DBusReport> r)
    {
    	report = r;
    }
#endif

private:
    std::string real_path;
    int state;
    double error_value;
    std::shared_ptr<sdbusplus::asio::connection> passive_bus;

#ifdef RIKFAN_DEBUG
    std::shared_ptr<DBusReport> report;
#endif

};

class Zone
{
private:
	static const constexpr long long loop_min_delay = 300;
	static const constexpr double stop_output_const = 130.0;
	static const constexpr double margin_error_read_temp = 100.0;

public:

	// Zone() = delete;
	Zone(const Zone &) = delete;
	Zone(Zone &&) = delete;
	void operator=(const Zone&) = delete;
	void operator=(Zone&&) = delete;

	Zone(std::string n,
	     std::string t,
	     ec::pidinfo &pidinfo_initial,
	     std::vector<std::string> &s,
	     std::vector<std::string> &f,
	     double sp,
	     long long ms,
	     std::shared_ptr<sdbusplus::asio::connection>& conn_
	    );

	~Zone();
	void start();
	void stop();
	void command(const char *cmd);

private:
	bool manualmode;
	bool stop_flag;
	long long millisec;
	std::unique_ptr<std::thread> pmainthread;

	std::string name;
	std::string type;
	ec::pid_info_t pid_info;
	double setpt;
	std::vector<std::unique_ptr<SensorBase>> sensors;
	std::vector<std::unique_ptr<SensorBase>> pwms;
	double error_read_temp;
	std::shared_ptr<sdbusplus::asio::connection> conn;

#ifdef RIKFAN_DEBUG
	decltype(std::chrono::system_clock::now()) sample_time;
#endif

	double processInputs();
	double processPID(double in);
	void processOutputs(double in);
	static void zone_control_loop(Zone *zone);
};

class ZoneManager
{
	/* async io context for operation */
	boost::asio::io_context io;

	std::shared_ptr<sdbusplus::asio::connection> conn;
	std::vector<std::unique_ptr<Zone>> zones;

public:
	// ZoneManager() = delete;
	ZoneManager(const ZoneManager &) = delete;
	ZoneManager(ZoneManager &&) = delete;
	void operator=(const ZoneManager&) = delete;
	void operator=(ZoneManager&&) = delete;

	ZoneManager(fs::path conf_fname, boost::asio::io_service& io_);

	void setFanMode(unsigned int mode);
	void start();

};

max_stars_repo_name: ShajTester/rikfan

max_issues_repo_name: ShajTester/rikfan

max_forks_repo_name: ShajTester/rikfan


text:
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    ArrayCasting.cxx

-------------------------------------------------------------------------
  Copyright 2008 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include <vtkDenseArray.h>
#include <vtkSmartPointer.h>
#include <vtkSparseArray.h>
#include <vtkTryDowncast.h>

#include <iostream>
#include <sstream>
#include <stdexcept>

#include <boost/algorithm/string.hpp>

#define VTK_CREATE(type, name) vtkSmartPointer<type> name = vtkSmartPointer<type>::New()

#define test_expression(expression)                                                                \
  {                                                                                \
    if (!(expression))                                                                             \
    {                                                                                \
      std::ostringstream buffer;                                                                   \
      buffer << "Expression failed at line " << __LINE__ << ": " << #expression;                   \
      throw std::runtime_error(buffer.str());                                                      \
    }                                                                                \
  }

class DowncastTest
{
public:
  DowncastTest(int& count)
    : Count(count)
  {
  }

  template <typename T>
  void operator()(T* vtkNotUsed(array)) const
  {
    ++Count;
  }

  int& Count;

private:
  DowncastTest& operator=(const DowncastTest&);
};

template <template <typename> class TargetT, typename TypesT>
void SuccessTest(vtkObject* source, int line)
{
  int count = 0;
  if (!vtkTryDowncast<TargetT, TypesT>(source, DowncastTest(count)))
  {
    std::ostringstream buffer;
    buffer << "Expression failed at line " << line;
    throw std::runtime_error(buffer.str());
  }

  if (count != 1)
  {
    std::ostringstream buffer;
    buffer << "Functor was called " << count << " times at line " << line;
    throw std::runtime_error(buffer.str());
  }
}

template <template <typename> class TargetT, typename TypesT>
void FailTest(vtkObject* source, int line)
{
  int count = 0;
  if (vtkTryDowncast<TargetT, TypesT>(source, DowncastTest(count)))
  {
    std::ostringstream buffer;
    buffer << "Expression failed at line " << line;
    throw std::runtime_error(buffer.str());
  }

  if (count != 0)
  {
    std::ostringstream buffer;
    buffer << "Functor was called " << count << " times at line " << line;
    throw std::runtime_error(buffer.str());
  }
}

/*
// This functor increments array values in-place using a parameter passed via the algorithm (instead
of a parameter
// stored in the functor).  It can work with any numeric array type.
struct IncrementValues
{
  template<typename T>
  void operator()(T* array, int amount) const
  {
    for(vtkIdType n = 0; n != array->GetNonNullSize(); ++n)
      array->SetValueN(n, array->GetValueN(n) + amount);
  }
};

// This functor converts strings in-place to a form suitable for case-insensitive comparison.  It's
an example of
// how you can write generic code while still specializing functionality on a case-by-case basis,
since
// in this situation we want to use some special functionality provided by vtkUnicodeString.
struct FoldCase
{
  template<typename ValueT>
  void operator()(vtkTypedArray<ValueT>* array) const
  {
    for(vtkIdType n = 0; n != array->GetNonNullSize(); ++n)
      {
      ValueT value = array->GetValueN(n);
      boost::algorithm::to_lower(value);
      array->SetValueN(n, value);
      }
  }

  void operator()(vtkTypedArray<vtkUnicodeString>* array) const
  {
    for(vtkIdType n = 0; n != array->GetNonNullSize(); ++n)
      array->SetValueN(n, array->GetValueN(n).fold_case());
  }
};

// This functor efficiently creates a transposed array.  It's one example of how you can create an
output array
// with the same type as an input array.
struct Transpose
{
  Transpose(vtkSmartPointer<vtkArray>& result_matrix) : ResultMatrix(result_matrix) {}

  template<typename ValueT>
  void operator()(vtkDenseArray<ValueT>* input) const
  {
    if(input->GetDimensions() != 2 || input->GetExtents()[0] != input->GetExtents()[1])
      throw std::runtime_error("A square matrix is required.");

    vtkDenseArray<ValueT>* output = vtkDenseArray<ValueT>::SafeDownCast(input->DeepCopy());
    for(vtkIdType i = 0; i != input->GetExtents()[0]; ++i)
      {
      for(vtkIdType j = i + 1; j != input->GetExtents()[1]; ++j)
        {
        output->SetValue(i, j, input->GetValue(j, i));
        output->SetValue(j, i, input->GetValue(i, j));
        }
      }

    this->ResultMatrix = output;
  }

  vtkSmartPointer<vtkArray>& ResultMatrix;
};
*/

//
//
// Here are some examples of how the algorithm might be called.
//
//

int TestArrayCasting(int vtkNotUsed(argc), char* vtkNotUsed(argv)[])
{
  try
  {
    VTK_CREATE(vtkDenseArray<int>, dense_int);
    VTK_CREATE(vtkDenseArray<double>, dense_double);
    VTK_CREATE(vtkDenseArray<vtkStdString>, dense_string);
    VTK_CREATE(vtkSparseArray<int>, sparse_int);
    VTK_CREATE(vtkSparseArray<double>, sparse_double);
    VTK_CREATE(vtkSparseArray<vtkStdString>, sparse_string);

    SuccessTest<vtkTypedArray, vtkIntegerTypes>(dense_int, __LINE__);
    FailTest<vtkTypedArray, vtkIntegerTypes>(dense_double, __LINE__);
    FailTest<vtkTypedArray, vtkIntegerTypes>(dense_string, __LINE__);
    SuccessTest<vtkTypedArray, vtkIntegerTypes>(sparse_int, __LINE__);
    FailTest<vtkTypedArray, vtkIntegerTypes>(sparse_double, __LINE__);
    FailTest<vtkTypedArray, vtkIntegerTypes>(sparse_string, __LINE__);

    FailTest<vtkTypedArray, vtkFloatingPointTypes>(dense_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkFloatingPointTypes>(dense_double, __LINE__);
    FailTest<vtkTypedArray, vtkFloatingPointTypes>(dense_string, __LINE__);
    FailTest<vtkTypedArray, vtkFloatingPointTypes>(sparse_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkFloatingPointTypes>(sparse_double, __LINE__);
    FailTest<vtkTypedArray, vtkFloatingPointTypes>(sparse_string, __LINE__);

    SuccessTest<vtkTypedArray, vtkNumericTypes>(dense_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkNumericTypes>(dense_double, __LINE__);
    FailTest<vtkTypedArray, vtkNumericTypes>(dense_string, __LINE__);
    SuccessTest<vtkTypedArray, vtkNumericTypes>(sparse_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkNumericTypes>(sparse_double, __LINE__);
    FailTest<vtkTypedArray, vtkNumericTypes>(sparse_string, __LINE__);

    FailTest<vtkTypedArray, vtkStringTypes>(dense_int, __LINE__);
    FailTest<vtkTypedArray, vtkStringTypes>(dense_double, __LINE__);
    SuccessTest<vtkTypedArray, vtkStringTypes>(dense_string, __LINE__);
    FailTest<vtkTypedArray, vtkStringTypes>(sparse_int, __LINE__);
    FailTest<vtkTypedArray, vtkStringTypes>(sparse_double, __LINE__);
    SuccessTest<vtkTypedArray, vtkStringTypes>(sparse_string, __LINE__);

    SuccessTest<vtkTypedArray, vtkAllTypes>(dense_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkAllTypes>(dense_double, __LINE__);
    SuccessTest<vtkTypedArray, vtkAllTypes>(dense_string, __LINE__);
    SuccessTest<vtkTypedArray, vtkAllTypes>(sparse_int, __LINE__);
    SuccessTest<vtkTypedArray, vtkAllTypes>(sparse_double, __LINE__);
    SuccessTest<vtkTypedArray, vtkAllTypes>(sparse_string, __LINE__);

    SuccessTest<vtkDenseArray, vtkAllTypes>(dense_int, __LINE__);
    FailTest<vtkDenseArray, vtkAllTypes>(sparse_int, __LINE__);
    FailTest<vtkSparseArray, vtkAllTypes>(dense_int, __LINE__);
    SuccessTest<vtkSparseArray, vtkAllTypes>(sparse_int, __LINE__);

    return 0;
  }
  catch (std::exception& e)
  {
    cerr << e.what() << endl;
    return 1;
  }
}

ext: cxx

max_stars_repo_name: txwhhny/vtk

max_issues_repo_name: txwhhny/vtk

max_forks_repo_name: txwhhny/vtk