text:
/* ---------------------------------------------------------------------
 *
 * Copyright (C) 2009 - 2021 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * ---------------------------------------------------------------------
 *
 * Author: Luca Heltai, Cataldo Manigrasso, 2009
 */

// @sect3{Include files}

// 程序一开始就包括了一堆include文件，我们将在程序的各个部分使用这些文件。其中大部分在以前的教程中已经讨论过了。

#include <deal.II/base/smartpointer.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/quadrature_selector.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/utilities.h>

#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>

#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>

// 这里有一些我们需要的C++标准头文件。

#include <cmath>
#include <iostream>
#include <fstream>
#include <string>

// 这个序言的最后部分是将dealii命名空间中的所有内容导入到这个程序中的所有内容中。

namespace Step34
{
  using namespace dealii;
// @sect3{Single and double layer operator kernels}

// 首先，让我们定义一下边界积分方程的机制。

// 以下两个函数是单层和双层势能核的实际计算，即  $G$  和  $\nabla G$  。只有当矢量 $R = \mathbf{y}-\mathbf{x}$ 不同于零时，它们才是定义良好的。

  namespace LaplaceKernel
  {
    template <int dim>
    double single_layer(const Tensor<1, dim> &R)
    {
      switch (dim)
        {
          case 2:
            return (-std::log(R.norm()) / (2 * numbers::PI));

          case 3:
            return (1. / (R.norm() * 4 * numbers::PI));

          default:
            Assert(false, ExcInternalError());
            return 0.;
        }
    }

    template <int dim>
    Tensor<1, dim> double_layer(const Tensor<1, dim> &R)
    {
      switch (dim)
        {
          case 2:
            return R / (-2 * numbers::PI * R.norm_square());
          case 3:
            return R / (-4 * numbers::PI * R.norm_square() * R.norm());

          default:
            Assert(false, ExcInternalError());
            return Tensor<1, dim>();
        }
    }
  } // namespace LaplaceKernel
// @sect3{The BEMProblem class}

// 边界元素方法代码的结构与有限元素代码的结构非常相似，所以这个类的成员函数与其他大多数教程程序的成员函数一样。特别是，现在你应该熟悉从外部文件中读取参数，以及将不同的任务分割成不同的模块。这同样适用于边界元素方法，我们不会对其进行过多的评论，只是对其中的差异进行评论。

  template <int dim>
  class BEMProblem
  {
  public:
    BEMProblem(const unsigned int fe_degree      = 1,
               const unsigned int mapping_degree = 1);

    void run();

  private:
    void read_parameters(const std::string &filename);

    void read_domain();

    void refine_and_resize();

// 我们在这里发现的唯一真正不同的函数是装配程序。我们以最通用的方式编写了这个函数，以便能够方便地推广到高阶方法和不同的基本解（例如斯托克斯或麦克斯韦）。

// 最明显的区别是，最终的矩阵是完整的，而且我们在通常的单元格循环内有一个嵌套的循环，访问所有自由度的支持点。 此外，当支持点位于我们所访问的单元内时，我们所执行的积分就会变成单数。

// 实际的结果是，我们有两套正交公式、有限元值和临时存储，一套用于标准积分，另一套用于奇异积分，在必要时使用。

    void assemble_system();

// 对于这个问题的解决有两种选择。第一个是使用直接求解器，第二个是使用迭代求解器。我们选择了第二种方案。

// 我们组装的矩阵不是对称的，我们选择使用GMRES方法；然而为边界元素方法构建一个有效的预处理程序并不是一个简单的问题。这里我们使用一个非预处理的GMRES求解器。迭代求解器的选项，如公差、最大迭代次数等，都是通过参数文件选择的。

    void solve_system();

// 一旦我们得到了解决方案，我们将计算计算出的势的 $L^2$ 误差，以及实体角的近似值的 $L^\infty$ 误差。我们使用的网格是平滑曲线的近似值，因此计算出的角的分量或实体角的对角线矩阵  $\alpha(\mathbf{x})$  应该一直等于  $\frac 12$  。在这个例程中，我们输出势的误差和计算角度的近似值的误差。注意，后者的误差实际上不是计算角度的误差，而是衡量我们对球体和圆的近似程度。

// 对角度的计算做一些实验，对于较简单的几何形状，可以得到非常准确的结果。为了验证这一点，你可以在read_domain()方法中注释掉tria.set_manifold(1, manifold)一行，并检查程序生成的alpha。通过删除这个调用，每当细化网格时，新的节点将沿着构成粗略网格的直线放置，而不是被拉到我们真正想要近似的表面。在三维案例中，球体的粗网格是从一个立方体开始得到的，得到的字母值正好是面的节点上的 $\frac 12$ ，边的节点上的 $\frac 34$ 和顶点的8个节点上的 $\frac 78$ 。

    void compute_errors(const unsigned int cycle);

// 一旦我们在一维领域得到了一个解决方案，我们就想把它插值到空间的其他部分。这可以通过在compute_exterior_solution()函数中再次进行解与核的卷积来实现。

// 我们想绘制速度变量，也就是势解的梯度。势解只在边界上是已知的，但我们使用与基本解的卷积在标准的二维连续有限元空间上进行插值。外推解的梯度图将给我们提供我们想要的速度。

// 除了外域上的解，我们还在output_results()函数中输出域的边界上的解，当然了。

    void compute_exterior_solution();

    void output_results(const unsigned int cycle);

// 为了实现不受维度限制的编程，我们对这个单一的函数进行了专业化处理，以提取整合单元内部的奇异核所需的奇异正交公式。

    const Quadrature<dim - 1> &get_singular_quadrature(
      const typename DoFHandler<dim - 1, dim>::active_cell_iterator &cell,
      const unsigned int index) const;

// 通常的deal.II类可以通过指定问题的 "二维 "来用于边界元素方法。这是通过将Triangulation, FiniteElement和DoFHandler的可选第二模板参数设置为嵌入空间的维度来实现的。在我们的例子中，我们生成了嵌入在二维或三维空间的一维或二维网格。

// 可选参数默认等于第一个参数，并产生我们在之前所有例子中看到的通常的有限元类。

// 该类的构造方式是允许任意的域（通过高阶映射）和有限元空间的逼近顺序。有限元空间和映射的顺序可以在该类的构造函数中选择。

    Triangulation<dim - 1, dim> tria;
    FE_Q<dim - 1, dim>          fe;
    DoFHandler<dim - 1, dim>    dof_handler;
    MappingQ<dim - 1, dim>      mapping;

// 在BEM方法中，生成的矩阵是密集的。根据问题的大小，最终的系统可能通过直接的LU分解来解决，或者通过迭代方法来解决。在这个例子中，我们使用了一个无条件的GMRES方法。为BEM方法建立一个预处理程序是不容易的，我们在此不做处理。

    FullMatrix<double> system_matrix;
    Vector<double>     system_rhs;

// 接下来的两个变量将表示解决方案 $\phi$ 以及一个向量，它将保存 $\alpha(\mathbf x)$ 的值（从一个点 $\mathbf x$ 可见的 $\Omega$ 的部分）在我们形状函数的支持点。

    Vector<double> phi;
    Vector<double> alpha;

// 收敛表是用来输出精确解和计算出的字母的误差的。

    ConvergenceTable convergence_table;

// 下面的变量是我们通过参数文件来填充的。 本例中我们使用的新对象是 Functions::ParsedFunction 对象和QuadratureSelector对象。

//  Functions::ParsedFunction 类允许我们通过参数文件方便快捷地定义新的函数对象，自定义的定义可以非常复杂（关于所有可用的选项，见该类的文档）。

// 我们将使用QuadratureSelector类来分配正交对象，该类允许我们根据一个识别字符串和公式本身的可能程度来生成正交公式。我们用它来允许自定义选择标准积分的正交公式，并定义奇异正交规则的顺序。

// 我们还定义了几个参数，这些参数是在我们想把解决方案扩展到整个领域的情况下使用的。

    Functions::ParsedFunction<dim> wind;
    Functions::ParsedFunction<dim> exact_solution;

    unsigned int                         singular_quadrature_order;
    std::shared_ptr<Quadrature<dim - 1>> quadrature;

    SolverControl solver_control;

    unsigned int n_cycles;
    unsigned int external_refinement;

    bool run_in_this_dimension;
    bool extend_solution;
  };
// @sect4{BEMProblem::BEMProblem and BEMProblem::read_parameters}

//构造函数初始化各种对象的方式与有限元程序（如  step-4  或  step-6  ）中的方式基本相同。这里唯一的新成分是ParsedFunction对象，它在构造时需要说明组件的数量。

// 对于精确解来说，向量分量的数量是1，而且不需要任何操作，因为1是ParsedFunction对象的默认值。然而，风需要指定dim组件。注意，在为 Functions::ParsedFunction, 的表达式声明参数文件中的条目时，我们需要明确指定分量的数量，因为函数 Functions::ParsedFunction::declare_parameters 是静态的，对分量的数量没有了解。

  template <int dim>
  BEMProblem<dim>::BEMProblem(const unsigned int fe_degree,
                              const unsigned int mapping_degree)
    : fe(fe_degree)
    , dof_handler(tria)
    , mapping(mapping_degree, true)
    , wind(dim)
    , singular_quadrature_order(5)
    , n_cycles(4)
    , external_refinement(5)
    , run_in_this_dimension(true)
    , extend_solution(true)
  {}

  template <int dim>
  void BEMProblem<dim>::read_parameters(const std::string &filename)
  {
    deallog << std::endl
            << "Parsing parameter file " << filename << std::endl
            << "for a " << dim << " dimensional simulation. " << std::endl;

    ParameterHandler prm;

    prm.declare_entry("Number of cycles", "4", Patterns::Integer());
    prm.declare_entry("External refinement", "5", Patterns::Integer());
    prm.declare_entry("Extend solution on the -2,2 box",
                      "true",
                      Patterns::Bool());
    prm.declare_entry("Run 2d simulation", "true", Patterns::Bool());
    prm.declare_entry("Run 3d simulation", "true", Patterns::Bool());

    prm.enter_subsection("Quadrature rules");
    {
      prm.declare_entry(
        "Quadrature type",
        "gauss",
        Patterns::Selection(
          QuadratureSelector<(dim - 1)>::get_quadrature_names()));
      prm.declare_entry("Quadrature order", "4", Patterns::Integer());
      prm.declare_entry("Singular quadrature order", "5", Patterns::Integer());
    }
    prm.leave_subsection();

// 对于二维和三维，我们将默认的输入数据设置为：解为  $x+y$  或  $x+y+z$  。实际计算出的解在无穷大时的数值为零。在这种情况下，这与精确解相吻合，不需要额外的修正，但是你应该注意，我们任意设置了 $\phi_\infty$ ，而我们传递给程序的精确解需要在无穷远处有相同的值，才能正确计算出误差。

//  Functions::ParsedFunction 对象的使用是非常直接的。 Functions::ParsedFunction::declare_parameters 函数需要一个额外的整数参数，指定给定函数的分量数量。它的默认值是1。当相应的 Functions::ParsedFunction::parse_parameters 方法被调用时，调用对象必须有与这里定义的相同数量的组件，否则会产生异常。

// 在声明条目时，我们同时声明了二维和三维的函数。然而只有二维的最终被解析。这使得我们对二维和三维问题都只需要一个参数文件。

// 注意，从数学的角度来看，边界上的风函数应该满足条件 $\int_{\partial\Omega} \mathbf{v}\cdot \mathbf{n} d \Gamma = 0$  ，这样问题才会有解。如果不满足这个条件，那么就找不到解，求解器也就不会收敛。

    prm.enter_subsection("Wind function 2d");
    {
      Functions::ParsedFunction<2>::declare_parameters(prm, 2);
      prm.set("Function expression", "1; 1");
    }
    prm.leave_subsection();

    prm.enter_subsection("Wind function 3d");
    {
      Functions::ParsedFunction<3>::declare_parameters(prm, 3);
      prm.set("Function expression", "1; 1; 1");
    }
    prm.leave_subsection();

    prm.enter_subsection("Exact solution 2d");
    {
      Functions::ParsedFunction<2>::declare_parameters(prm);
      prm.set("Function expression", "x+y");
    }
    prm.leave_subsection();

    prm.enter_subsection("Exact solution 3d");
    {
      Functions::ParsedFunction<3>::declare_parameters(prm);
      prm.set("Function expression", "x+y+z");
    }
    prm.leave_subsection();

// 在求解器部分，我们设置所有的SolverControl参数。然后，该对象将在solve_system()函数中被送入GMRES求解器。

    prm.enter_subsection("Solver");
    SolverControl::declare_parameters(prm);
    prm.leave_subsection();

// 在向ParameterHandler对象声明了所有这些参数后，让我们读取一个输入文件，该文件将为这些参数提供其值。然后我们继续从ParameterHandler对象中提取这些值。

    prm.parse_input(filename);

    n_cycles            = prm.get_integer("Number of cycles");
    external_refinement = prm.get_integer("External refinement");
    extend_solution     = prm.get_bool("Extend solution on the -2,2 box");

    prm.enter_subsection("Quadrature rules");
    {
      quadrature = std::shared_ptr<Quadrature<dim - 1>>(
        new QuadratureSelector<dim - 1>(prm.get("Quadrature type"),
                                        prm.get_integer("Quadrature order")));
      singular_quadrature_order = prm.get_integer("Singular quadrature order");
    }
    prm.leave_subsection();

    prm.enter_subsection("Wind function " + std::to_string(dim) + "d");
    {
      wind.parse_parameters(prm);
    }
    prm.leave_subsection();

    prm.enter_subsection("Exact solution " + std::to_string(dim) + "d");
    {
      exact_solution.parse_parameters(prm);
    }
    prm.leave_subsection();

    prm.enter_subsection("Solver");
    solver_control.parse_parameters(prm);
    prm.leave_subsection();

// 最后，这里是另一个如何在独立维度编程中使用参数文件的例子。 如果我们想关闭两个模拟中的一个，我们可以通过设置相应的 "运行2D模拟 "或 "运行3D模拟 "标志为假来实现。

    run_in_this_dimension =
      prm.get_bool("Run " + std::to_string(dim) + "d simulation");
  }
// @sect4{BEMProblem::read_domain}

// 边界元素法三角剖分与（dim-1）维三角剖分基本相同，不同之处在于顶点属于（dim）维空间。

// deal.II中支持的一些网格格式默认使用三维点来描述网格。这些格式与deal.II的边界元素方法功能兼容。特别是我们可以使用UCD或GMSH格式。在这两种情况下，我们必须特别注意网格的方向，因为与标准有限元的情况不同，这里没有进行重新排序或兼容性检查。 所有的网格都被认为是有方向性的，因为它们被嵌入到一个高维空间中。参见GridIn和Triangulation的文档，以进一步了解三角结构中单元的方向。在我们的例子中，网格的法线是外在于2D的圆或3D的球体。

// 对边界元素网格进行适当细化所需要的另一个细节是对网格所逼近的流形的准确描述。对于标准有限元网格的边界，我们已经多次看到了这一点（例如在 step-5 和 step-6 中），这里的原理和用法是一样的，只是SphericalManifold类需要一个额外的模板参数来指定嵌入空间维度。

  template <int dim>
  void BEMProblem<dim>::read_domain()
  {
    const Point<dim>                      center = Point<dim>();
    const SphericalManifold<dim - 1, dim> manifold(center);

    std::ifstream in;
    switch (dim)
      {
        case 2:
          in.open("coarse_circle.inp");
          break;

        case 3:
          in.open("coarse_sphere.inp");
          break;

        default:
          Assert(false, ExcNotImplemented());
      }

    GridIn<dim - 1, dim> gi;
    gi.attach_triangulation(tria);
    gi.read_ucd(in);

    tria.set_all_manifold_ids(1);

// 对  Triangulation::set_manifold  的调用复制了流形（通过  Manifold::clone()),  所以我们不需要担心对  manifold  的无效指针。

    tria.set_manifold(1, manifold);
  }
// @sect4{BEMProblem::refine_and_resize}

// 这个函数对网格进行全局细化，分配自由度，并调整矩阵和向量的大小。

  template <int dim>
  void BEMProblem<dim>::refine_and_resize()
  {
    tria.refine_global(1);

    dof_handler.distribute_dofs(fe);

    const unsigned int n_dofs = dof_handler.n_dofs();

    system_matrix.reinit(n_dofs, n_dofs);

    system_rhs.reinit(n_dofs);
    phi.reinit(n_dofs);
    alpha.reinit(n_dofs);
  }
// @sect4{BEMProblem::assemble_system}

// 下面是这个程序的主要功能，组装与边界积分方程相对应的矩阵。

  template <int dim>
  void BEMProblem<dim>::assemble_system()
  {

// 首先我们用正交公式初始化一个FEValues对象，用于在非奇异单元中进行内核积分。这个正交公式是通过参数文件选择的，并且需要相当精确，因为我们要积分的函数不是多项式函数。

    FEValues<dim - 1, dim> fe_v(mapping,
                                fe,
                                *quadrature,
                                update_values | update_normal_vectors |
                                  update_quadrature_points | update_JxW_values);

    const unsigned int n_q_points = fe_v.n_quadrature_points;

    std::vector<types::global_dof_index> local_dof_indices(
      fe.n_dofs_per_cell());

    std::vector<Vector<double>> cell_wind(n_q_points, Vector<double>(dim));
    double                      normal_wind;

// 与有限元方法不同的是，如果我们使用拼合边界元方法，那么在每个装配循环中，我们只装配与一个自由度（与支撑点 $i$ 相关的自由度）和当前单元之间的耦合信息。这是用fe.dofs_per_cell元素的向量完成的，然后将其分配到全局行的矩阵中  $i$  。以下对象将持有这些信息。

    Vector<double> local_matrix_row_i(fe.n_dofs_per_cell());

// 索引  $i$  运行在拼合点上，这是  $i$  第三个基函数的支持点，而  $j$  运行在内部积分点上。

// 我们构建一个支持点的向量，它将用于局部积分。

    std::vector<Point<dim>> support_points(dof_handler.n_dofs());
    DoFTools::map_dofs_to_support_points<dim - 1, dim>(mapping,
                                                       dof_handler,
                                                       support_points);

// 这样做之后，我们就可以开始对所有单元进行积分循环，首先初始化FEValues对象，得到正交点的 $\mathbf{\tilde v}$ 的值（这个向量场应该是常数，但更通用也无妨）。

    for (const auto &cell : dof_handler.active_cell_iterators())
      {
        fe_v.reinit(cell);
        cell->get_dof_indices(local_dof_indices);

        const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
        const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();
        wind.vector_value_list(q_points, cell_wind);

// 然后我们在当前单元上形成所有自由度的积分（注意，这包括不在当前单元上的自由度，这与通常的有限元积分有偏差）。如果其中一个局部自由度与支持点 $i$ 相同，我们需要执行的积分是单数。因此，在循环的开始，我们检查是否是这种情况，并存储哪一个是奇异指数。

        for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
          {
            local_matrix_row_i = 0;

            bool         is_singular    = false;
            unsigned int singular_index = numbers::invalid_unsigned_int;

            for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
              if (local_dof_indices[j] == i)
                {
                  singular_index = j;
                  is_singular    = true;
                  break;
                }

// 然后我们进行积分。如果指数 $i$ 不是局部自由度之一，我们只需将单层项加到右边，将双层项加到矩阵中。

            if (is_singular == false)
              {
                for (unsigned int q = 0; q < n_q_points; ++q)
                  {
                    normal_wind = 0;
                    for (unsigned int d = 0; d < dim; ++d)
                      normal_wind += normals[q][d] * cell_wind[q](d);

                    const Tensor<1, dim> R = q_points[q] - support_points[i];

                    system_rhs(i) += (LaplaceKernel::single_layer(R) *
                                      normal_wind * fe_v.JxW(q));

                    for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)

                      local_matrix_row_i(j) -=
                        ((LaplaceKernel::double_layer(R) * normals[q]) *
                         fe_v.shape_value(j, q) * fe_v.JxW(q));
                  }
              }
            else
              {

// 现在我们处理更微妙的情况。如果我们在这里，这意味着在 $j$ 索引上运行的单元包含support_point[i]。在这种情况下，单层和双层势都是单数，它们需要特殊处理。
//每当在给定单元内进行积分时，就会使用一个特殊的正交公式，允许人们对参考单元上的奇异权重进行任意函数的积分。
//正确的正交公式由get_singular_quadrature函数选择，下面将详细说明。

                Assert(singular_index != numbers::invalid_unsigned_int,
                       ExcInternalError());

                const Quadrature<dim - 1> &singular_quadrature =
                  get_singular_quadrature(cell, singular_index);

                FEValues<dim - 1, dim> fe_v_singular(
                  mapping,
                  fe,
                  singular_quadrature,
                  update_jacobians | update_values | update_normal_vectors |
                    update_quadrature_points);

                fe_v_singular.reinit(cell);

                std::vector<Vector<double>> singular_cell_wind(
                  singular_quadrature.size(), Vector<double>(dim));

                const std::vector<Tensor<1, dim>> &singular_normals =
                  fe_v_singular.get_normal_vectors();
                const std::vector<Point<dim>> &singular_q_points =
                  fe_v_singular.get_quadrature_points();

                wind.vector_value_list(singular_q_points, singular_cell_wind);

                for (unsigned int q = 0; q < singular_quadrature.size(); ++q)
                  {
                    const Tensor<1, dim> R =
                      singular_q_points[q] - support_points[i];
                    double normal_wind = 0;
                    for (unsigned int d = 0; d < dim; ++d)
                      normal_wind +=
                        (singular_cell_wind[q](d) * singular_normals[q][d]);

                    system_rhs(i) += (LaplaceKernel::single_layer(R) *
                                      normal_wind * fe_v_singular.JxW(q));

                    for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
                      {
                        local_matrix_row_i(j) -=
                          ((LaplaceKernel::double_layer(R) *
                            singular_normals[q]) *
                           fe_v_singular.shape_value(j, q) *
                           fe_v_singular.JxW(q));
                      }
                  }
              }

// 最后，我们需要将当前单元格的贡献添加到全局矩阵中。

            for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
              system_matrix(i, local_dof_indices[j]) += local_matrix_row_i(j);
          }
      }

// 积分运算符的第二部分是术语  $\alpha(\mathbf{x}_i) \phi_j(\mathbf{x}_i)$  。由于我们使用的是配位方案， $\phi_j(\mathbf{x}_i)=\delta_{ij}$  而相应的矩阵是一个对角线的矩阵，其条目等于 $\alpha(\mathbf{x}_i)$  。

// 计算这个实体角的对角矩阵的一个快速方法是使用诺伊曼矩阵本身。只需将该矩阵与一个元素都等于-1的向量相乘，就可以得到阿尔法角或实体角的对角线矩阵（见介绍中的公式）。然后将这个结果加回到系统矩阵对象上，得到矩阵的最终形式。

    Vector<double> ones(dof_handler.n_dofs());
    ones.add(-1.);

    system_matrix.vmult(alpha, ones);
    alpha.add(1);
    for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
      system_matrix(i, i) += alpha(i);
  }
// @sect4{BEMProblem::solve_system}

// 下一个函数简单地解决了线性系统。

  template <int dim>
  void BEMProblem<dim>::solve_system()
  {
    SolverGMRES<Vector<double>> solver(solver_control);
    solver.solve(system_matrix, phi, system_rhs, PreconditionIdentity());
  }
// @sect4{BEMProblem::compute_errors}

// 误差的计算在其他所有的例子程序中都是完全一样的，我们就不做过多的评论。请注意，在有限元方法中使用的方法在这里也可以使用。

  template <int dim>
  void BEMProblem<dim>::compute_errors(const unsigned int cycle)
  {
    Vector<float> difference_per_cell(tria.n_active_cells());
    VectorTools::integrate_difference(mapping,
                                      dof_handler,
                                      phi,
                                      exact_solution,
                                      difference_per_cell,
                                      QGauss<(dim - 1)>(2 * fe.degree + 1),
                                      VectorTools::L2_norm);
    const double L2_error =
      VectorTools::compute_global_error(tria,
                                        difference_per_cell,
                                        VectorTools::L2_norm);

//可以直接使用 Vector::linfty_norm() 函数来计算α向量的误差，因为在每个节点上，该值应该是 $\frac 12$  。然后，所有的误差都会被输出并附加到我们的ConvergenceTable对象中，以便以后计算收敛率。

    Vector<double> difference_per_node(alpha);
    difference_per_node.add(-.5);

    const double       alpha_error    = difference_per_node.linfty_norm();
    const unsigned int n_active_cells = tria.n_active_cells();
    const unsigned int n_dofs         = dof_handler.n_dofs();

    deallog << "Cycle " << cycle << ':' << std::endl
            << "   Number of active cells:       " << n_active_cells
            << std::endl
            << "   Number of degrees of freedom: " << n_dofs << std::endl;

    convergence_table.add_value("cycle", cycle);
    convergence_table.add_value("cells", n_active_cells);
    convergence_table.add_value("dofs", n_dofs);
    convergence_table.add_value("L2(phi)", L2_error);
    convergence_table.add_value("Linfty(alpha)", alpha_error);
  }

// 奇异积分需要仔细选择正交规则。特别是deal.II库提供了为对数奇异性（QGaussLog, QGaussLogR）以及1/R奇异性（QGaussOneOverR）量身定制的正交规则。

// 奇异积分通常是通过构建具有奇异权重的加权正交公式得到的，因此可以写成

// \f[ \int_K f(x) s(x) dx = \sum_{i=1}^N w_i f(q_i) \f]

// 其中 $s(x)$ 是一个给定的奇点，权重和正交点 $w_i,q_i$ 是精心选择的，以使上述公式对某类函数 $f(x)$ 是一个等式。

// 在我们目前看到的所有有限元例子中，正交点本身的权重（即函数  $s(x)$  ），总是不断地等于1。 对于奇异积分，我们有两个选择：我们可以使用上面的定义，从积分中剔除奇异性（即用特殊的正交规则对 $f(x)$ 进行积分），或者我们可以要求正交规则用 $s(q_i)$ 对权重 $w_i$ 进行 "标准化"。

// \f[ \int_K f(x) s(x) dx = \int_K g(x) dx = \sum_{i=1}^N \frac{w_i}{s(q_i)} g(q_i) \f]

// 我们通过QGaussLogR和QGaussOneOverR的 @p factor_out_singularity 参数来使用这第二种选择。

// 这些积分有些微妙，特别是在二维空间，由于从实数到参考单元的转换，积分的变量是以转换的行列式为尺度的。

// 在二维空间中，这个过程不仅会导致一个因子作为常数出现在整个积分上，而且还会导致一个需要评估的额外积分。

// \f[ \int_0^1 f(x)\ln(x/\alpha) dx = \int_0^1 f(x)\ln(x) dx - \int_0^1  f(x) \ln(\alpha) dx.  \f]

// 这个过程由QGaussLogR类的构造函数来处理，它增加了额外的正交点和权重，以考虑到积分的第二部分。

// 类似的推理应该在三维情况下进行，因为奇异正交是在参考单元的半径 $r$ 的逆上定制的，而我们的奇异函数生活在实空间，然而在三维情况下一切都更简单，因为奇异性与变换的行列式成线性比例。这使我们可以只建立一次奇异的二维正交规则，并在所有单元格中重复使用。

// 在一维的奇异积分中，这是不可能的，因为我们需要知道正交的缩放参数，而这个参数并不是先验的。这里，正交规则本身也取决于当前单元格的大小。出于这个原因，有必要为每个单数积分创建一个新的正交。

// 不同的正交规则是在get_singular_quadrature中建立的，它专门用于dim=2和dim=3，它们在assemble_system函数中被检索。作为参数给出的索引是奇异点所在的单位支持点的索引。

  template <>
  const Quadrature<2> &BEMProblem<3>::get_singular_quadrature(
    const DoFHandler<2, 3>::active_cell_iterator &,
    const unsigned int index) const
  {
    Assert(index < fe.n_dofs_per_cell(),
           ExcIndexRange(0, fe.n_dofs_per_cell(), index));

    static std::vector<QGaussOneOverR<2>> quadratures;
    if (quadratures.size() == 0)
      for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
        quadratures.emplace_back(singular_quadrature_order,
                                 fe.get_unit_support_points()[i],
                                 true);
    return quadratures[index];
  }

  template <>
  const Quadrature<1> &BEMProblem<2>::get_singular_quadrature(
    const DoFHandler<1, 2>::active_cell_iterator &cell,
    const unsigned int                            index) const
  {
    Assert(index < fe.n_dofs_per_cell(),
           ExcIndexRange(0, fe.n_dofs_per_cell(), index));

    static Quadrature<1> *q_pointer = nullptr;
    if (q_pointer)
      delete q_pointer;

    q_pointer = new QGaussLogR<1>(singular_quadrature_order,
                                  fe.get_unit_support_points()[index],
                                  1. / cell->measure(),
                                  true);
    return (*q_pointer);
  }

//  @sect4{BEMProblem::compute_exterior_solution}

// 我们还想知道一些关于外域中电势 $\phi$ 的值：毕竟我们考虑边界积分问题的动机是想知道外域中的速度!

// 为此，我们在此假设边界元素域包含在盒子 $[-2,2]^{\text{dim}}$ 中，我们用与基本解的卷积来推算这个盒子内的实际解。这方面的公式在介绍中已经给出。

// 整个空间的解的重建是在一个连续的、尺寸为dim的有限元网格上完成的。这些都是常用的，我们不做进一步评论。在函数的最后，我们再次以通常的方式输出这个外部解。

  template <int dim>
  void BEMProblem<dim>::compute_exterior_solution()
  {
    Triangulation<dim> external_tria;
    GridGenerator::hyper_cube(external_tria, -2, 2);

    FE_Q<dim>       external_fe(1);
    DoFHandler<dim> external_dh(external_tria);
    Vector<double>  external_phi;

    external_tria.refine_global(external_refinement);
    external_dh.distribute_dofs(external_fe);
    external_phi.reinit(external_dh.n_dofs());

    FEValues<dim - 1, dim> fe_v(mapping,
                                fe,
                                *quadrature,
                                update_values | update_normal_vectors |
                                  update_quadrature_points | update_JxW_values);

    const unsigned int n_q_points = fe_v.n_quadrature_points;

    std::vector<types::global_dof_index> dofs(fe.n_dofs_per_cell());

    std::vector<double>         local_phi(n_q_points);
    std::vector<double>         normal_wind(n_q_points);
    std::vector<Vector<double>> local_wind(n_q_points, Vector<double>(dim));

    std::vector<Point<dim>> external_support_points(external_dh.n_dofs());
    DoFTools::map_dofs_to_support_points<dim>(StaticMappingQ1<dim>::mapping,
                                              external_dh,
                                              external_support_points);

    for (const auto &cell : dof_handler.active_cell_iterators())
      {
        fe_v.reinit(cell);

        const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
        const std::vector<Tensor<1, dim>> &normals = fe_v.get_normal_vectors();

        cell->get_dof_indices(dofs);
        fe_v.get_function_values(phi, local_phi);

        wind.vector_value_list(q_points, local_wind);

        for (unsigned int q = 0; q < n_q_points; ++q)
          {
            normal_wind[q] = 0;
            for (unsigned int d = 0; d < dim; ++d)
              normal_wind[q] += normals[q][d] * local_wind[q](d);
          }

        for (unsigned int i = 0; i < external_dh.n_dofs(); ++i)
          for (unsigned int q = 0; q < n_q_points; ++q)
            {
              const Tensor<1, dim> R = q_points[q] - external_support_points[i];

              external_phi(i) +=
                ((LaplaceKernel::single_layer(R) * normal_wind[q] +
                  (LaplaceKernel::double_layer(R) * normals[q]) *
                    local_phi[q]) *
                 fe_v.JxW(q));
            }
      }

    DataOut<dim> data_out;

    data_out.attach_dof_handler(external_dh);
    data_out.add_data_vector(external_phi, "external_phi");
    data_out.build_patches();

    const std::string filename = std::to_string(dim) + "d_external.vtk";
    std::ofstream     file(filename);

  }
// @sect4{BEMProblem::output_results}

// 输出我们的计算结果是一个相当机械的任务。这个函数的所有组成部分之前已经讨论过了。

  template <int dim>
  void BEMProblem<dim>::output_results(const unsigned int cycle)
  {
    DataOut<dim - 1, dim> dataout;

    dataout.attach_dof_handler(dof_handler);
    dataout.add_data_vector(phi, "phi", DataOut<dim - 1, dim>::type_dof_data);
    dataout.add_data_vector(alpha,
                            "alpha",
                            DataOut<dim - 1, dim>::type_dof_data);
    dataout.build_patches(mapping,
                          mapping.get_degree(),
                          DataOut<dim - 1, dim>::curved_inner_cells);

    const std::string filename = std::to_string(dim) + "d_boundary_solution_" +
                                 std::to_string(cycle) + ".vtk";
    std::ofstream file(filename);

    dataout.write_vtk(file);

    if (cycle == n_cycles - 1)
      {
        convergence_table.set_precision("L2(phi)", 3);
        convergence_table.set_precision("Linfty(alpha)", 3);

        convergence_table.set_scientific("L2(phi)", true);
        convergence_table.set_scientific("Linfty(alpha)", true);

        convergence_table.evaluate_convergence_rates(
          "L2(phi)", ConvergenceTable::reduction_rate_log2);
        convergence_table.evaluate_convergence_rates(
          "Linfty(alpha)", ConvergenceTable::reduction_rate_log2);
        deallog << std::endl;
        convergence_table.write_text(std::cout);
      }
  }
// @sect4{BEMProblem::run}

// 这是最主要的功能。它应该是不言自明的。

  template <int dim>
  void BEMProblem<dim>::run()
  {
    read_parameters("parameters.prm");

    if (run_in_this_dimension == false)
      {
        deallog << "Run in dimension " << dim
                << " explicitly disabled in parameter file. " << std::endl;
        return;
      }

    read_domain();

    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
      {
        refine_and_resize();
        assemble_system();
        solve_system();
        compute_errors(cycle);
        output_results(cycle);
      }

    if (extend_solution == true)
      compute_exterior_solution();
  }
} // namespace Step34
// @sect3{The main() function}

// 这是本程序的主要功能。它与以前所有的教程程序完全一样。

int main()
{
  try
    {
      using namespace Step34;

      const unsigned int degree         = 1;
      const unsigned int mapping_degree = 1;

      deallog.depth_console(3);
      BEMProblem<2> laplace_problem_2d(degree, mapping_degree);
      laplace_problem_2d.run();

      BEMProblem<3> laplace_problem_3d(degree, mapping_degree);
      laplace_problem_3d.run();
    }
  catch (std::exception &exc)
    {
      std::cerr << std::endl
                << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Exception on processing: " << std::endl
                << exc.what() << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;

      return 1;
    }
  catch (...)
    {
      std::cerr << std::endl
                << std::endl
                << "----------------------------------------------------"
                << std::endl;
      std::cerr << "Unknown exception!" << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
      return 1;
    }

  return 0;
}

ext: cc


text:
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing,
// software distributed under the License is distributed on an
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.  See the License for the
// specific language governing permissions and limitations
// under the License.

#include <algorithm>
#include <atomic>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <initializer_list>
#include <memory>
#include <numeric>
#include <ostream>
#include <string>
#include <thread>
#include <utility>
#include <vector>

#include <boost/optional/optional.hpp>
#include <gflags/gflags.h>
#include <glog/logging.h>
#include <google/protobuf/arena.h>
#include <gtest/gtest.h>

#include "kudu/client/client.h"
#include "kudu/client/schema.h"
#include "kudu/common/common.pb.h"
#include "kudu/common/partial_row.h"
#include "kudu/common/row_operations.h"
#include "kudu/common/row_operations.pb.h"
#include "kudu/common/schema.h"
#include "kudu/common/wire_protocol.h"
#include "kudu/gutil/ref_counted.h"
#include "kudu/gutil/strings/substitute.h"
#include "kudu/integration-tests/cluster_itest_util.h"
#include "kudu/integration-tests/test_workload.h"
#include "kudu/master/catalog_manager.h"
#include "kudu/master/master.h"
#include "kudu/master/master.pb.h"
#include "kudu/master/master.proxy.h"
#include "kudu/master/mini_master.h"
#include "kudu/mini-cluster/external_mini_cluster.h"
#include "kudu/mini-cluster/internal_mini_cluster.h"
#include "kudu/rpc/messenger.h"
#include "kudu/rpc/rpc_controller.h"
#include "kudu/tserver/mini_tablet_server.h"
#include "kudu/util/hdr_histogram.h"
#include "kudu/util/metrics.h"
#include "kudu/util/monotime.h"
#include "kudu/util/net/sockaddr.h"
#include "kudu/util/path_util.h"
#include "kudu/util/pb_util.h"
#include "kudu/util/status.h"
#include "kudu/util/test_macros.h"
#include "kudu/util/test_util.h"

using kudu::client::KuduClient;
using kudu::client::KuduSchema;
using kudu::client::KuduTableAlterer;
using kudu::client::KuduTableCreator;
using kudu::cluster::ExternalMiniCluster;
using kudu::cluster::ExternalMiniClusterOptions;
using kudu::cluster::InternalMiniCluster;
using kudu::cluster::InternalMiniClusterOptions;
using kudu::master::ReplicaTypeFilter;
using kudu::pb_util::SecureDebugString;
using kudu::rpc::Messenger;
using kudu::rpc::MessengerBuilder;
using kudu::rpc::RpcController;
using std::ostringstream;
using std::pair;
using std::shared_ptr;
using std::string;
using std::thread;
using std::unique_ptr;
using std::vector;
using strings::Substitute;

DEFINE_int32(benchmark_runtime_secs, 5, "Number of seconds to run the benchmark");
DEFINE_int32(benchmark_num_threads, 16, "Number of threads to run the benchmark");
DEFINE_int32(benchmark_num_tablets, 60, "Number of tablets to create");

DECLARE_bool(enable_per_range_hash_schemas);
DECLARE_double(leader_failure_max_missed_heartbeat_periods);
DECLARE_int32(follower_unavailable_considered_failed_sec);
DECLARE_int32(heartbeat_interval_ms);
DECLARE_int32(max_create_tablets_per_ts);
DECLARE_int32(raft_heartbeat_interval_ms);
DECLARE_int32(rpc_num_service_threads);
DECLARE_int32(rpc_service_queue_length);
DECLARE_int32(table_locations_ttl_ms);
DECLARE_string(location_mapping_cmd);
DECLARE_uint32(table_locations_cache_capacity_mb);

METRIC_DECLARE_entity(server);

METRIC_DECLARE_counter(rpcs_queue_overflow);
METRIC_DECLARE_counter(table_locations_cache_evictions);
METRIC_DECLARE_counter(table_locations_cache_hits);
METRIC_DECLARE_counter(table_locations_cache_inserts);
METRIC_DECLARE_counter(table_locations_cache_lookups);
METRIC_DECLARE_counter(table_locations_cache_misses);
METRIC_DECLARE_gauge_uint64(table_locations_cache_memory_usage);
METRIC_DECLARE_histogram(handler_latency_kudu_master_MasterService_GetTableLocations);

namespace kudu {
namespace master {

// Test the master GetTableLocations RPC. This can't be done in master-test,
// since it requires a running tablet server in order for the catalog manager to
// report tablet locations.
class TableLocationsTest : public KuduTest {
 public:

  void SetUp() override {
    KuduTest::SetUp();

    FLAGS_max_create_tablets_per_ts = 1000;
    SetUpConfig();

    InternalMiniClusterOptions opts;
    opts.num_tablet_servers = 3;

    cluster_.reset(new InternalMiniCluster(env_, opts));
    ASSERT_OK(cluster_->Start());

    // Create a client proxy to the master.
    MessengerBuilder bld("Client");
    ASSERT_OK(bld.Build(&client_messenger_));
    const auto& addr = cluster_->mini_master()->bound_rpc_addr();
    proxy_.reset(new MasterServiceProxy(client_messenger_, addr, addr.host()));
  }

  void TearDown() override {
    if (cluster_) {
      cluster_->Shutdown();
      cluster_.reset();
    }
    KuduTest::TearDown();
  }

  virtual void SetUpConfig() {}

  struct HashBucketSchema {
    vector<string> columns;
    int32_t num_buckets;
    uint32_t seed;
  };

  Status CreateTable(const string& table_name,
                     const Schema& schema,
                     const vector<KuduPartialRow>& split_rows,
                     const vector<pair<KuduPartialRow, KuduPartialRow>>& bounds,
                     const vector<vector<HashBucketSchema>>& range_hash_schema,
                     const vector<HashBucketSchema>& table_hash_schema);

  void CreateTable(const string& table_name, int num_splits);

  // Check that the master doesn't give back partial results while the table is being created.
  void CheckMasterTableCreation(const string &table_name, int tablet_locations_size);

  shared_ptr<Messenger> client_messenger_;
  unique_ptr<InternalMiniCluster> cluster_;
  unique_ptr<MasterServiceProxy> proxy_;
};

Status TableLocationsTest::CreateTable(const string& table_name,
                                       const Schema& schema,
                                       const vector<KuduPartialRow>& split_rows = {},
                                       const vector<pair<KuduPartialRow,
                                                         KuduPartialRow>>& bounds = {},
                                       const vector<vector<HashBucketSchema>>&
                                           range_hash_schema = {},
                                       const vector<HashBucketSchema>& table_hash_schema = {}) {

  CreateTableRequestPB req;
  CreateTableResponsePB resp;
  RpcController controller;

  req.set_name(table_name);
  RETURN_NOT_OK(SchemaToPB(schema, req.mutable_schema()));
  RowOperationsPBEncoder encoder(req.mutable_split_rows_range_bounds());
  for (const KuduPartialRow& row : split_rows) {
    encoder.Add(RowOperationsPB::SPLIT_ROW, row);
  }
  for (const auto& bound : bounds) {
    encoder.Add(RowOperationsPB::RANGE_LOWER_BOUND, bound.first);
    encoder.Add(RowOperationsPB::RANGE_UPPER_BOUND, bound.second);
  }

  for (const auto& hash_schemas : range_hash_schema) {
    auto* range_hash_schemas_pb = req.add_range_hash_schemas();
    for (const auto& hash_schema : hash_schemas) {
      auto* hash_schema_pb = range_hash_schemas_pb->add_hash_schemas();
      for (const string& col_name : hash_schema.columns) {
        hash_schema_pb->add_columns()->set_name(col_name);
      }
      hash_schema_pb->set_num_buckets(hash_schema.num_buckets);
      hash_schema_pb->set_seed(hash_schema.seed);
    }
  }

  if (!table_hash_schema.empty()) {
    auto* partition_schema_pb = req.mutable_partition_schema();
    for (const auto& hash_schema : table_hash_schema) {
      auto* hash_schema_pb = partition_schema_pb->add_hash_bucket_schemas();
      for (const string& col_name : hash_schema.columns) {
        hash_schema_pb->add_columns()->set_name(col_name);
      }
      hash_schema_pb->set_num_buckets(hash_schema.num_buckets);
      hash_schema_pb->set_seed(hash_schema.seed);
    }
  }

  return proxy_->CreateTable(req, &resp, &controller);
}

void TableLocationsTest::CheckMasterTableCreation(const string &table_name,
                                                  int tablet_locations_size) {
  GetTableLocationsRequestPB req;
  GetTableLocationsResponsePB resp;
  RpcController controller;
  req.set_max_returned_locations(1000);
  req.mutable_table()->set_table_name(table_name);

  for (int i = 1; ; i++) {
    if (i > 10) {
      FAIL() << "Create table timed out";
    }

    controller.Reset();
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    if (resp.has_error()) {
      ASSERT_EQ(MasterErrorPB::TABLET_NOT_RUNNING, resp.error().code());
      SleepFor(MonoDelta::FromMilliseconds(i * i * 100));
    } else {
      ASSERT_EQ(tablet_locations_size, resp.tablet_locations().size());
      break;
    }
  }
}

void TableLocationsTest::CreateTable(const string& table_name, int num_splits) {
  Schema schema({ ColumnSchema("key", INT32) }, 1);
  KuduPartialRow row(&schema);
  vector<KuduPartialRow> splits(num_splits, row);
  for (int i = 0; i < num_splits; i++) {
    ASSERT_OK(splits[i].SetInt32(0, i*1000));
  }

  ASSERT_OK(CreateTable(table_name, schema, splits));
  NO_FATALS(CheckMasterTableCreation(table_name, num_splits + 1));
}

// Test the tablet server location is properly set in the master GetTableLocations RPC.
class TableLocationsWithTSLocationTest : public TableLocationsTest {
 public:
  void SetUpConfig() override {
    const string location_cmd_path = JoinPathSegments(GetTestExecutableDirectory(),
                                                      "testdata/first_argument.sh");
    const string location = "/foo";
    FLAGS_location_mapping_cmd = Substitute("$0 $1", location_cmd_path, location);
  }
};

// Test that when the client requests table locations for a non-covered
// partition range, the master returns the first tablet previous to the begin
// partition key, as specified in the non-covering range partitions design
// document.
TEST_F(TableLocationsTest, TestGetTableLocations) {
  const string table_name = "test";
  Schema schema({ ColumnSchema("key", STRING) }, 1);
  KuduPartialRow row(&schema);

  vector<KuduPartialRow> splits(6, row);
  ASSERT_OK(splits[0].SetStringNoCopy(0, "aa"));
  ASSERT_OK(splits[1].SetStringNoCopy(0, "ab"));
  ASSERT_OK(splits[2].SetStringNoCopy(0, "ac"));
  ASSERT_OK(splits[3].SetStringNoCopy(0, "ca"));
  ASSERT_OK(splits[4].SetStringNoCopy(0, "cb"));
  ASSERT_OK(splits[5].SetStringNoCopy(0, "cc"));

  vector<pair<KuduPartialRow, KuduPartialRow>> bounds(2, { row, row });
  ASSERT_OK(bounds[0].first.SetStringNoCopy(0, "a"));
  ASSERT_OK(bounds[0].second.SetStringNoCopy(0, "b"));
  ASSERT_OK(bounds[1].first.SetStringNoCopy(0, "c"));
  ASSERT_OK(bounds[1].second.SetStringNoCopy(0, "d"));

  ASSERT_OK(CreateTable(table_name, schema, splits, bounds));

  NO_FATALS(CheckMasterTableCreation(table_name, 8));

  { // from "a"
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_partition_key_start("a");
    req.set_max_returned_locations(3);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(3, resp.tablet_locations().size());
    EXPECT_EQ("a", resp.tablet_locations(0).partition().partition_key_start());
    EXPECT_EQ("aa", resp.tablet_locations(1).partition().partition_key_start());
    EXPECT_EQ("ab", resp.tablet_locations(2).partition().partition_key_start());
  }

  { // from ""
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_partition_key_start("");
    req.set_max_returned_locations(3);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(3, resp.tablet_locations().size());
    EXPECT_EQ("a", resp.tablet_locations(0).partition().partition_key_start());
    EXPECT_EQ("aa", resp.tablet_locations(1).partition().partition_key_start());
    EXPECT_EQ("ab", resp.tablet_locations(2).partition().partition_key_start());

    auto get_tablet_location = [](MasterServiceProxy* proxy, const string& tablet_id) {
      rpc::RpcController rpc_old;
      GetTabletLocationsRequestPB req_old;
      GetTabletLocationsResponsePB resp_old;
      *req_old.add_tablet_ids() = tablet_id;
      ASSERT_OK(proxy->GetTabletLocations(req_old, &resp_old, &rpc_old));
      const auto& loc_old = resp_old.tablet_locations(0);
      ASSERT_GT(loc_old.deprecated_replicas_size(), 0);
      ASSERT_EQ(0, loc_old.interned_replicas_size());

      rpc::RpcController rpc_new;
      GetTabletLocationsRequestPB req_new;
      GetTabletLocationsResponsePB resp_new;
      *req_new.add_tablet_ids() = tablet_id;
      req_new.set_intern_ts_infos_in_response(true);
      ASSERT_OK(proxy->GetTabletLocations(req_new, &resp_new, &rpc_new));
      const auto& loc_new = resp_new.tablet_locations(0);
      ASSERT_GT(loc_new.interned_replicas_size(), 0);
      ASSERT_EQ(0, loc_new.deprecated_replicas_size());

      ASSERT_EQ(loc_old.tablet_id(), loc_new.tablet_id());
    };

    // Check that a UUID was returned for every replica
    for (const auto& loc : resp.tablet_locations()) {
      NO_FATALS(get_tablet_location(proxy_.get(), loc.tablet_id()));
      ASSERT_GT(loc.deprecated_replicas_size(), 0);
      ASSERT_EQ(0, loc.interned_replicas_size());
      for (const auto& replica : loc.deprecated_replicas()) {
        ASSERT_NE("", replica.ts_info().permanent_uuid());
      }
    }
  }

  { // from "", with TSInfo interning enabled.
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_partition_key_start("");
    req.set_max_returned_locations(3);
    req.set_intern_ts_infos_in_response(true);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(3, resp.tablet_locations().size());
    EXPECT_EQ("a", resp.tablet_locations(0).partition().partition_key_start());
    EXPECT_EQ("aa", resp.tablet_locations(1).partition().partition_key_start());
    EXPECT_EQ("ab", resp.tablet_locations(2).partition().partition_key_start());
    // Check that a UUID was returned for every replica
    for (const auto& loc : resp.tablet_locations()) {
      ASSERT_EQ(loc.deprecated_replicas_size(), 0);
      ASSERT_GT(loc.interned_replicas_size(), 0);
      for (const auto& replica : loc.interned_replicas()) {
        int idx = replica.ts_info_idx();
        ASSERT_GE(idx, 0);
        ASSERT_LE(idx, resp.ts_infos_size());
        ASSERT_NE("", resp.ts_infos(idx).permanent_uuid());
      }
    }
  }

  { // from "b"
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_partition_key_start("b");
    req.set_max_returned_locations(3);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(3, resp.tablet_locations().size());
    EXPECT_EQ("ac", resp.tablet_locations(0).partition().partition_key_start());
    EXPECT_EQ("c", resp.tablet_locations(1).partition().partition_key_start());
    EXPECT_EQ("ca", resp.tablet_locations(2).partition().partition_key_start());
  }

  { // from "z"
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_partition_key_start("z");
    req.set_max_returned_locations(3);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(1, resp.tablet_locations().size());
    EXPECT_EQ("cc", resp.tablet_locations(0).partition().partition_key_start());
  }

  { // Check that the tablet server location should be an empty string if not set.
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_max_returned_locations(1);
    req.set_intern_ts_infos_in_response(true);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(1, resp.tablet_locations().size());
    ASSERT_EQ(3, resp.tablet_locations(0).interned_replicas_size());
    const auto& loc = resp.tablet_locations(0);
    for (const auto& replica : loc.interned_replicas()) {
      EXPECT_EQ("", resp.ts_infos(replica.ts_info_idx()).location());
    }
  }
}

TEST_F(TableLocationsTest, TestRangeSpecificHashing) {
  const string table_name = "test";
  Schema schema({ ColumnSchema("key", STRING), ColumnSchema("val", STRING) }, 2);
  KuduPartialRow row(&schema);

  FLAGS_enable_per_range_hash_schemas = true; // enable for testing.

  vector<pair<KuduPartialRow, KuduPartialRow>> bounds(3, { row, row });
  ASSERT_OK(bounds[0].first.SetStringNoCopy(0, "a"));
  ASSERT_OK(bounds[0].second.SetStringNoCopy(0, "b"));
  ASSERT_OK(bounds[1].first.SetStringNoCopy(0, "c"));
  ASSERT_OK(bounds[1].second.SetStringNoCopy(0, "d"));
  ASSERT_OK(bounds[2].first.SetStringNoCopy(0, "e"));
  ASSERT_OK(bounds[2].second.SetStringNoCopy(0, "f"));

  vector<vector<HashBucketSchema>> range_hash_schema;
  vector<HashBucketSchema> hash_schema_4_by_2 = { { { "key" }, 4, 0 }, { { "val" }, 2, 0} };
  range_hash_schema.emplace_back(hash_schema_4_by_2);
  vector<HashBucketSchema> hash_schema_6 = { { { "key" }, 6, 2 } };
  range_hash_schema.emplace_back(hash_schema_6);

  // Table-wide hash schema, applied to range by default if no per-range schema is specified.
  vector<HashBucketSchema> table_hash_schema_5 = { { { "val" }, 5, 4 } };
  range_hash_schema.emplace_back(vector<HashBucketSchema>());

  ASSERT_OK(CreateTable(table_name, schema, {}, bounds, range_hash_schema, table_hash_schema_5));
  NO_FATALS(CheckMasterTableCreation(table_name, 19));

  GetTableLocationsRequestPB req;
  GetTableLocationsResponsePB resp;
  RpcController controller;
  req.mutable_table()->set_table_name(table_name);
  req.set_max_returned_locations(19);
  ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
  SCOPED_TRACE(SecureDebugString(resp));

  ASSERT_FALSE(resp.has_error());
  ASSERT_EQ(19, resp.tablet_locations().size());

  vector<string> partition_key_starts =  {
      string("\0\0\0\0" "\0\0\0\0" "a" "\0\0", 11),
      string("\0\0\0\0" "\0\0\0\1" "a" "\0\0", 11),
      string("\0\0\0\1" "\0\0\0\0" "a" "\0\0", 11),
      string("\0\0\0\1" "\0\0\0\1" "a" "\0\0", 11),
      string("\0\0\0\2" "\0\0\0\0" "a" "\0\0", 11),
      string("\0\0\0\2" "\0\0\0\1" "a" "\0\0", 11),
      string("\0\0\0\3" "\0\0\0\0" "a" "\0\0", 11),
      string("\0\0\0\3" "\0\0\0\1" "a" "\0\0", 11),
      string("\0\0\0\0" "c" "\0\0", 7),
      string("\0\0\0\1" "c" "\0\0", 7),
      string("\0\0\0\2" "c" "\0\0", 7),
      string("\0\0\0\3" "c" "\0\0", 7),
      string("\0\0\0\4" "c" "\0\0", 7),
      string("\0\0\0\5" "c" "\0\0", 7),
      string("\0\0\0\0" "e" "\0\0", 7),
      string("\0\0\0\1" "e" "\0\0", 7),
      string("\0\0\0\2" "e" "\0\0", 7),
      string("\0\0\0\3" "e" "\0\0", 7),
      string("\0\0\0\4" "e" "\0\0", 7)
  };
  // Sorting partition keys to match the tablets that are returned in sorted order.
  sort(partition_key_starts.begin(), partition_key_starts.end());
  ASSERT_EQ(partition_key_starts.size(), resp.tablet_locations_size());
  for (int i = 0; i < resp.tablet_locations_size(); i++) {
    EXPECT_EQ(partition_key_starts[i],
              resp.tablet_locations(i).partition().partition_key_start());
  }
}

TEST_F(TableLocationsWithTSLocationTest, TestGetTSLocation) {
  const string table_name = "test";
  Schema schema({ ColumnSchema("key", STRING) }, 1);
  ASSERT_OK(CreateTable(table_name, schema));

  NO_FATALS(CheckMasterTableCreation(table_name, 1));

  {
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    RpcController controller;
    req.mutable_table()->set_table_name(table_name);
    req.set_intern_ts_infos_in_response(true);
    ASSERT_OK(proxy_->GetTableLocations(req, &resp, &controller));
    SCOPED_TRACE(SecureDebugString(resp));

    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(1, resp.tablet_locations().size());
    ASSERT_EQ(3, resp.tablet_locations(0).interned_replicas_size());
    const auto& loc = resp.tablet_locations(0);
    for (const auto& replica : loc.interned_replicas()) {
      ASSERT_EQ("/foo", resp.ts_infos(replica.ts_info_idx()).location());
    }
  }
}

TEST_F(TableLocationsTest, GetTableLocationsBenchmark) {
  SKIP_IF_SLOW_NOT_ALLOWED();
  const int kNumSplits = FLAGS_benchmark_num_tablets - 1;
  const int kNumThreads = FLAGS_benchmark_num_threads;
  const auto kRuntime = MonoDelta::FromSeconds(FLAGS_benchmark_runtime_secs);

  const string table_name = "test";
  NO_FATALS(CreateTable(table_name, kNumSplits));

  // Make one proxy per thread, so each thread gets its own messenger and
  // reactor. If there were only one messenger, then only one reactor thread
  // would be used for the connection to the master, so this benchmark would
  // probably be testing the serialization and network code rather than the
  // master GTL code.
  vector<unique_ptr<MasterServiceProxy>> proxies;
  proxies.reserve(kNumThreads);
  for (int i = 0; i < kNumThreads; i++) {
    MessengerBuilder bld("Client");
    bld.set_num_reactors(1);
    shared_ptr<Messenger> msg;
    ASSERT_OK(bld.Build(&msg));
    const auto& addr = cluster_->mini_master()->bound_rpc_addr();
    proxies.emplace_back(new MasterServiceProxy(std::move(msg), addr, addr.host()));
  }

  std::atomic<bool> stop { false };
  vector<thread> threads;
  threads.reserve(kNumThreads);
  // If a thread encounters an error, the test is failed.
  vector<uint64_t> err_counters(kNumThreads, 0);
  for (auto i = 0; i < kNumThreads; ++i) {
    threads.emplace_back([&, i]() {
        GetTableLocationsRequestPB req;
        GetTableLocationsResponsePB resp;
        while (!stop) {
          RpcController controller;
          req.Clear();
          resp.Clear();
          req.mutable_table()->set_table_name(table_name);
          req.set_max_returned_locations(1000);
          req.set_intern_ts_infos_in_response(true);
          const auto s = proxies[i]->GetTableLocations(req, &resp, &controller);
          if (!s.ok()) {
            LOG(WARNING) << "GetTableLocations() failed: " << s.ToString();
            ++err_counters[i];
            break;
          }
          CHECK_EQ(resp.tablet_locations_size(), kNumSplits + 1);
        }
      });
  }

  SleepFor(kRuntime);
  stop = true;
  for (auto& t : threads) {
    t.join();
  }

  const auto& ent = cluster_->mini_master()->master()->metric_entity();
  auto counter = METRIC_rpcs_queue_overflow.Instantiate(ent);
  auto hist = METRIC_handler_latency_kudu_master_MasterService_GetTableLocations
      .Instantiate(ent);

  cluster_->Shutdown();

  LOG(INFO) << "RPC queue overflows: " << counter->value();
  const auto errors = accumulate(err_counters.begin(), err_counters.end(), 0UL);
  if (errors != 0) {
    FAIL() << Substitute("detected $0 errors", errors);
  }

  LOG(INFO) << Substitute(
      "GetTableLocations RPC: $0 req/sec",
      hist->histogram()->TotalCount() / kRuntime.ToSeconds());

  ostringstream ostr;
  ostr << "Stats on GetTableLocations RPC (times in microseconds): ";
  hist->histogram()->DumpHumanReadable(&ostr);
  LOG(INFO) << ostr.str();
}

// Similar to GetTableLocationsBenchmark, but calls the function directly within
// the mini-cluster process instead of issuing RPC.
TEST_F(TableLocationsTest, GetTableLocationsBenchmarkFunctionCall) {
  SKIP_IF_SLOW_NOT_ALLOWED();
  const string kUserName = "testuser";
  const size_t kNumSplits = FLAGS_benchmark_num_tablets - 1;
  const size_t kNumThreads = FLAGS_benchmark_num_threads;
  const auto kRuntime = MonoDelta::FromSeconds(FLAGS_benchmark_runtime_secs);

  const string table_name = "test";
  NO_FATALS(CreateTable(table_name, kNumSplits));

  CatalogManager* cm = cluster_->mini_master()->master()->catalog_manager();
  const boost::optional<const string&> username = kUserName;

  std::atomic<bool> stop(false);
  vector<thread> threads;
  threads.reserve(kNumThreads);
  vector<uint64_t> req_counters(kNumThreads, 0);
  vector<uint64_t> err_counters(kNumThreads, 0);
  for (size_t idx = 0; idx < kNumThreads; ++idx) {
    threads.emplace_back([&, idx]() {
      GetTableLocationsRequestPB req;
      google::protobuf::Arena arena;
      while (!stop) {
        req.Clear();
        arena.Reset();
        auto* resp = google::protobuf::Arena::CreateMessage<GetTableLocationsResponsePB>(&arena);
        req.mutable_table()->set_table_name(table_name);
        req.set_max_returned_locations(1000);
        req.set_intern_ts_infos_in_response(true);
        ++req_counters[idx];
        {
          CatalogManager::ScopedLeaderSharedLock l(cm);
          auto s = cm->GetTableLocations(&req, resp, username);
          if (!s.ok()) {
            ++err_counters[idx];
          }
        }
      }
    });
  }

  SleepFor(kRuntime);
  stop = true;
  for (auto& t : threads) {
    t.join();
  }
  cluster_->Shutdown();

  const auto errors = accumulate(err_counters.begin(), err_counters.end(), 0UL);
  if (errors != 0) {
    FAIL() << Substitute("detected $0 errors", errors);
  }

  const double total = accumulate(req_counters.begin(), req_counters.end(), 0UL);
  LOG(INFO) << Substitute(
      "GetTableLocations function call: $0 req/sec",
      total / kRuntime.ToSeconds());
}

// A small benchmark to see whether multiple clients refreshing table location
// might create the thundering herd pattern with their GetTableLocations()
// requests, overflowing the RPC service queue.
class RefreshTableLocationsTest : public KuduTest {
 public:
  void SetUp() override {
    KuduTest::SetUp();
    // Set the RPC queue size limit and the number of RPC service threads small
    // to register RPC queue overlows, if any. The idea is to make this test
    // use smaller numbers, while in real life there would be more clients
    // and regular length of the RPC service queue.
    FLAGS_rpc_num_service_threads = 2;
    FLAGS_rpc_service_queue_length = std::max(1, FLAGS_benchmark_num_threads / 3);
    // Make table location information expire faster so the clients would try
    // to refresh the information during the scenario, while not spending too
    // much time to run the scenario.
    FLAGS_table_locations_ttl_ms = 1500;
    {
      InternalMiniClusterOptions opts;
      opts.num_masters = 1;
      opts.num_tablet_servers = 3;
      cluster_.reset(new InternalMiniCluster(env_, std::move(opts)));
    }
    ASSERT_OK(cluster_->Start());
  }

  void TearDown() override {
    if (cluster_) {
      cluster_->Shutdown();
      cluster_.reset();
    }
    KuduTest::TearDown();
  }

 protected:
  unique_ptr<InternalMiniCluster> cluster_;
};

TEST_F(RefreshTableLocationsTest, ThunderingHerd) {
  SKIP_IF_SLOW_NOT_ALLOWED();
  const auto kNumThreads = FLAGS_benchmark_num_threads;
  const auto kRuntime =
      MonoDelta::FromMilliseconds(5 * FLAGS_table_locations_ttl_ms);

  // Create test tables and populate them with data. These tables are used by
  // the threads below, each running its own read test workload.
  for (auto idx = 0; idx < kNumThreads; ++idx) {
    TestWorkload w(cluster_.get());
    w.set_num_replicas(3);
    w.set_num_tablets(10);
    w.set_table_name(Substitute("test_table_$0", idx));
    w.set_num_write_threads(1);
    w.set_num_read_threads(0);
    w.Setup();
    w.Start();
    while (w.rows_inserted() < 100) {
      SleepFor(MonoDelta::FromMilliseconds(10));
    }
    w.StopAndJoin();
  }

  std::atomic<bool> stop(false);
  vector<thread> threads;
  threads.reserve(kNumThreads);
  for (auto idx = 0; idx < kNumThreads; ++idx) {
    // It's important to have multiple workloads to simulate multiple clients,
    // each with its own table locations cache.
    threads.emplace_back([&, idx]() {
      TestWorkload w(cluster_.get());
      w.set_table_name(Substitute("test_table_$0", idx));
      w.set_num_read_threads(1);
      w.set_num_write_threads(0);
      w.Setup();
      w.Start();
      while (!stop) {
        SleepFor(MonoDelta::FromMilliseconds(10));
      }
      w.StopAndJoin();
    });
  }

  // Make sure the reader threads refresh the expired information on table
  // locations few times during the scenario.
  SleepFor(kRuntime);
  stop = true;
  for (auto& t : threads) {
    t.join();
  }

  const auto& ent = cluster_->mini_master()->master()->metric_entity();
  auto counter = METRIC_rpcs_queue_overflow.Instantiate(ent);
  auto hist = METRIC_handler_latency_kudu_master_MasterService_GetTableLocations.
      Instantiate(ent);

  cluster_->Shutdown();

  LOG(INFO) << "RPC queue overflows: " << counter->value();
  LOG(INFO) << Substitute(
      "GetTableLocations RPC: $0 req/sec",
      hist->histogram()->TotalCount() / kRuntime.ToSeconds());

  ostringstream ostr;
  ostr << "Stats on GetTableLocations RPC (times in microseconds): ";
  hist->histogram()->DumpHumanReadable(&ostr);
  LOG(INFO) << ostr.str();
}

class TableLocationsCacheBaseTest : public KuduTest {
 public:
  TableLocationsCacheBaseTest(int num_tablet_servers,
                              int cache_capacity_mb)
      : num_tablet_servers_(num_tablet_servers),
        cache_capacity_mb_(cache_capacity_mb),
        schema_(itest::SimpleIntKeyKuduSchema()) {
  }

  void SetUp() override {
    KuduTest::SetUp();

    InternalMiniClusterOptions opts;
    opts.num_tablet_servers = num_tablet_servers_;

    // Setting the cache's capacity to anything > 0 enables caching of table
    // locations.
    FLAGS_table_locations_cache_capacity_mb = cache_capacity_mb_;

    cluster_.reset(new InternalMiniCluster(env_, opts));
    ASSERT_OK(cluster_->Start());
    ASSERT_OK(cluster_->CreateClient(nullptr, &client_));

    unique_ptr<KuduTableCreator> table_creator(client_->NewTableCreator());
    ASSERT_OK(table_creator->table_name(kTableName)
        .schema(&schema_)
        .set_range_partition_columns({ "key" })
        .num_replicas(std::min<int>(3, num_tablet_servers_))
        .add_hash_partitions({ "key" }, 5)
        .Create());
    metric_entity_ = cluster_->mini_master()->master()->metric_entity();
    ASSERT_NE(nullptr, metric_entity_.get());
  }

  void TearDown() override {
    if (cluster_) {
      cluster_->Shutdown();
      cluster_.reset();
    }
    KuduTest::TearDown();
  }

  int64_t GetCacheInserts() const {
    return METRIC_table_locations_cache_inserts.Instantiate(metric_entity_)->value();
  }

  int64_t GetCacheLookups() const {
    return METRIC_table_locations_cache_lookups.Instantiate(metric_entity_)->value();
  }

  int64_t GetCacheHits() const {
    return METRIC_table_locations_cache_hits.Instantiate(metric_entity_)->value();
  }

  int64_t GetCacheMisses() const {
    return METRIC_table_locations_cache_misses.Instantiate(metric_entity_)->value();
  }

  int64_t GetCacheEvictions() const {
    return METRIC_table_locations_cache_evictions.Instantiate(metric_entity_)->value();
  }

  uint64_t GetCacheMemoryUsage() const {
    return METRIC_table_locations_cache_memory_usage.Instantiate(metric_entity_, 0)->value();
  }

  // Insert about the 'rows_num' into the table named 'table_name', returning
  // the actual number of rows written.
  int64_t InsertRows(const string& table_name, int64_t rows_num = 10) {
    TestWorkload w(cluster_.get());
    w.set_table_name(table_name);
    w.set_schema(schema_);
    w.set_num_read_threads(0);
    w.set_num_write_threads(1);
    w.set_write_batch_size(1);
    w.set_write_pattern(TestWorkload::INSERT_SEQUENTIAL_ROWS);
    w.set_already_present_allowed(true);
    w.Setup();
    w.Start();
    while (w.rows_inserted() < rows_num) {
      SleepFor(MonoDelta::FromMilliseconds(10));
    }
    w.StopAndJoin();

    // Return the actual number of rows inserted.
    return w.rows_inserted();
  }

 protected:
  static constexpr const char* const kTableName = "table_locations_cache_test";
  const int num_tablet_servers_;
  const int cache_capacity_mb_;

  const KuduSchema schema_;

  unique_ptr<InternalMiniCluster> cluster_;
  client::sp::shared_ptr<KuduClient> client_;
  scoped_refptr<MetricEntity> metric_entity_;
};

class TableLocationsCacheDisabledTest : public TableLocationsCacheBaseTest {
 public:
  TableLocationsCacheDisabledTest()
      : TableLocationsCacheBaseTest(1 /*num_tablet_servers*/,
                                    0 /*cache_capacity_mb*/) {}
};

TEST_F(TableLocationsCacheDisabledTest, Basic) {
  // Verify that by default tablet location cache is disabled.
  google::CommandLineFlagInfo flag_info;
  ASSERT_TRUE(google::GetCommandLineFlagInfo(
      "table_locations_cache_capacity_mb", &flag_info));
  ASSERT_TRUE(flag_info.is_default);
  ASSERT_EQ("0", flag_info.default_value);

  InsertRows(kTableName);

  // When disabled, cache should not be used.
  EXPECT_EQ(0, GetCacheEvictions());
  EXPECT_EQ(0, GetCacheHits());
  EXPECT_EQ(0, GetCacheInserts());
  EXPECT_EQ(0, GetCacheLookups());
  EXPECT_EQ(0, GetCacheMisses());
  EXPECT_EQ(0, GetCacheMemoryUsage());
}

class TableLocationsCacheTest : public TableLocationsCacheBaseTest {
 public:
  TableLocationsCacheTest()
      : TableLocationsCacheBaseTest(1 /*num_tablet_servers*/,
                                    1 /*cache_capacity_mb*/) {}
};

// Verify that requests with different essential parameters produce
// different entries in the table locations cache.
TEST_F(TableLocationsCacheTest, DifferentRequestParameters) {
  MessengerBuilder bld("test_builder");
  shared_ptr<Messenger> messenger;
  ASSERT_OK(bld.Build(&messenger));
  const auto& addr = cluster_->mini_master()->bound_rpc_addr();
  MasterServiceProxy proxy(messenger, addr, addr.host());

  {
    // Issue one query many times -- there should be just one record inserted.
    const auto prev_cache_inserts = GetCacheInserts();
    for (auto i = 0; i < 5; ++i) {
      GetTableLocationsRequestPB req;
      GetTableLocationsResponsePB resp;
      req.mutable_table()->set_table_name(kTableName);
      req.set_max_returned_locations(1);
      RpcController ctl;
      ASSERT_OK(proxy.GetTableLocations(req, &resp, &ctl));
      ASSERT_TRUE(!resp.has_error());
    }
    ASSERT_EQ(prev_cache_inserts + 1, GetCacheInserts());
  }

  {
    // Issue a query with different value of 'max_returned_locations'. A new
    // entry should be added into the cache even if other parameters are the
    // same as in requests sent prior.
    const auto prev_cache_inserts = GetCacheInserts();
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    req.mutable_table()->set_table_name(kTableName);
    req.set_max_returned_locations(10);
    RpcController ctl;
    ASSERT_OK(proxy.GetTableLocations(req, &resp, &ctl));
    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(prev_cache_inserts + 1, GetCacheInserts());
  }

  {
    // Set 'replica_type_filter' parameter.
    const auto prev_cache_inserts = GetCacheInserts();
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    req.mutable_table()->set_table_name(kTableName);
    req.set_max_returned_locations(10);
    req.set_replica_type_filter(ReplicaTypeFilter::ANY_REPLICA);
    RpcController ctl;
    ASSERT_OK(proxy.GetTableLocations(req, &resp, &ctl));
    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(prev_cache_inserts + 1, GetCacheInserts());
  }

  {
    // Switch 'replica_type_filter' parameter to a different value.
    const auto prev_cache_inserts = GetCacheInserts();
    GetTableLocationsRequestPB req;
    GetTableLocationsResponsePB resp;
    req.mutable_table()->set_table_name(kTableName);
    req.set_max_returned_locations(10);
    req.set_replica_type_filter(ReplicaTypeFilter::VOTER_REPLICA);
    RpcController ctl;
    ASSERT_OK(proxy.GetTableLocations(req, &resp, &ctl));
    ASSERT_TRUE(!resp.has_error());
    ASSERT_EQ(prev_cache_inserts + 1, GetCacheInserts());
  }
}

// This scenario verifies basic functionality of the table locations cache.
TEST_F(TableLocationsCacheTest, Basic) {
  SKIP_IF_SLOW_NOT_ALLOWED();

  int64_t prev_cache_hits = 0;
  int64_t prev_cache_inserts = 0;
  int64_t prev_cache_lookups = 0;
  int64_t prev_cache_misses = 0;
  int64_t prev_cache_evictions = 0;

  {
    // Run a workload to prime the table locations cache.
    InsertRows(kTableName);
    const auto cache_lookups = GetCacheLookups();
    EXPECT_GT(cache_lookups, 0);
    const auto cache_misses = GetCacheMisses();
    const auto cache_inserts = GetCacheInserts();
    EXPECT_EQ(cache_misses, cache_inserts);
    const auto cache_hits = GetCacheHits();
    EXPECT_EQ(cache_lookups - cache_misses, cache_hits);

    // The cache usage should be non-zero since some entries are present.
    EXPECT_GT(GetCacheMemoryUsage(), 0);

    // Store the counters for the next round.
    prev_cache_hits = cache_hits;
    prev_cache_inserts = cache_inserts;
    prev_cache_lookups = cache_lookups;
    prev_cache_misses = cache_misses;
    prev_cache_evictions = GetCacheEvictions();
  }

  {
    // The cache has been primed -- now run another workload.
    InsertRows(kTableName);

    // No new cache inserts expected.
    EXPECT_EQ(prev_cache_inserts, GetCacheInserts());
    // No new cache misses expected.
    EXPECT_EQ(prev_cache_misses, GetCacheMisses());

    // There should be new cache lookups.
    const auto cache_lookups = GetCacheLookups();
    EXPECT_LT(prev_cache_lookups, cache_lookups);

    // There should be new cache hits.
    const auto cache_hits = GetCacheHits();
    EXPECT_LT(prev_cache_hits, cache_hits);

    // All new lookups should hit the cache.
    EXPECT_EQ(cache_lookups - prev_cache_lookups, cache_hits - prev_cache_hits);
  }

  {
    // Alter the test table and make sure the information isn't purged from the
    // table locations cache since the table locations haven't changed.
    const string new_table_name = string(kTableName) + "_renamed";

    unique_ptr<KuduTableAlterer> a0(client_->NewTableAlterer(kTableName));
    auto s = a0->RenameTo(new_table_name)->Alter();
    ASSERT_OK(s);

    // Wait for the master to get notified on the change.
    SleepFor(MonoDelta::FromMilliseconds(3 * FLAGS_heartbeat_interval_ms));

    // No new cache evictions are expected.
    EXPECT_EQ(prev_cache_evictions, GetCacheEvictions());

    InsertRows(new_table_name);

    EXPECT_EQ(prev_cache_evictions, GetCacheEvictions());
    // No new cache inserts nor misses are expected.
    EXPECT_EQ(prev_cache_inserts, GetCacheInserts());
    EXPECT_EQ(prev_cache_misses, GetCacheMisses());

    // Rename the table back.
    unique_ptr<KuduTableAlterer> a1(client_->NewTableAlterer(new_table_name));
    s = a1->RenameTo(kTableName)->Alter();
    ASSERT_OK(s);
  }

  {
    // Drop the test table and make sure the cached location information is
    // purged from the table locations cache.
    ASSERT_OK(client_->DeleteTable(kTableName));

    // Wait for the master to get notified on the change.
    SleepFor(MonoDelta::FromMilliseconds(3 * FLAGS_heartbeat_interval_ms));

    // All the previously inserted records should be purged.
    EXPECT_EQ(GetCacheEvictions(), GetCacheInserts());
    EXPECT_EQ(0, GetCacheMemoryUsage());
  }
}

class TableLocationsCacheTabletChangeTest :
    public TableLocationsCacheBaseTest {
 public:
  TableLocationsCacheTabletChangeTest()
      : TableLocationsCacheBaseTest(5 /*num_tablet_servers*/,
                                    1 /*cache_capacity_mb*/) {}
};

// Verify the behavior of the table locations cache when table's tablets change
// leadership roles or move around.
TEST_F(TableLocationsCacheTabletChangeTest, Basic) {
  SKIP_IF_SLOW_NOT_ALLOWED();

  {
    // Run a workload to prime the table locations cache.
    InsertRows(kTableName);
    const auto cache_lookups = GetCacheLookups();
    EXPECT_GT(cache_lookups, 0);
    const auto cache_misses = GetCacheMisses();
    EXPECT_GT(cache_misses, 0);
    EXPECT_EQ(cache_misses, GetCacheInserts());
    EXPECT_EQ(cache_lookups - cache_misses, GetCacheHits());
  }

  // Restart tablet servers to reshuffle leadership role of tablet replicas.
  for (auto idx = 0; idx < cluster_->num_tablet_servers(); ++idx) {
    cluster_->mini_tablet_server(idx)->Shutdown();
    SleepFor(MonoDelta::FromMilliseconds(static_cast<int64_t>(
        2 * FLAGS_leader_failure_max_missed_heartbeat_periods *
        FLAGS_raft_heartbeat_interval_ms) + 3 * FLAGS_heartbeat_interval_ms));
    ASSERT_OK(cluster_->mini_tablet_server(idx)->Restart());
  }

  {
    // All the previously inserted records should be purged.
    EXPECT_EQ(GetCacheEvictions(), GetCacheInserts());
    const auto prev_cache_inserts = GetCacheInserts();
    const auto prev_cache_misses = GetCacheMisses();

    InsertRows(kTableName);

    // New cache inserts and misses are expected since the previously inserted
    // records have been just purged.
    EXPECT_LT(prev_cache_inserts, GetCacheInserts());
    EXPECT_LT(prev_cache_misses, GetCacheMisses());
  }

  // Shutdown one tablet server and get the system some time to re-replicate
  // affected tablet replicas elsewhere.
  FLAGS_follower_unavailable_considered_failed_sec = 1;
  cluster_->mini_tablet_server(0)->Shutdown();
  SleepFor(MonoDelta::FromMilliseconds(static_cast<int64_t>(
      2 * FLAGS_leader_failure_max_missed_heartbeat_periods *
          FLAGS_raft_heartbeat_interval_ms +
      3 * 1000 * FLAGS_follower_unavailable_considered_failed_sec +
      3 * FLAGS_heartbeat_interval_ms)));

  // All the previously inserted records should be purged.
  EXPECT_EQ(GetCacheEvictions(), GetCacheInserts());
  EXPECT_EQ(0, GetCacheMemoryUsage());
}

class TableLocationsCacheMultiMasterTest : public KuduTest {
 public:
  TableLocationsCacheMultiMasterTest()
      : schema_(itest::SimpleIntKeyKuduSchema()) {
  }

  void SetUp() override {
    KuduTest::SetUp();

    ExternalMiniClusterOptions opts;
    opts.num_masters = 3;
    opts.num_tablet_servers = 3;
    opts.extra_master_flags = {
      "--table_locations_cache_capacity_mb=8",
      Substitute("--raft_heartbeat_interval_ms=$0",
          kRaftHeartbeatIntervalMs),
      Substitute("--leader_failure_max_missed_heartbeat_periods=$0",
          kMaxMissedHeartbeatPeriods),
    };
    cluster_.reset(new cluster::ExternalMiniCluster(std::move(opts)));
    ASSERT_OK(cluster_->Start());

    client::sp::shared_ptr<client::KuduClient> client;
    ASSERT_OK(cluster_->CreateClient(nullptr, &client));

    unique_ptr<KuduTableCreator> table_creator(client->NewTableCreator());
    ASSERT_OK(table_creator->table_name(kTableName)
        .schema(&schema_)
        .set_range_partition_columns({ "key" })
        .num_replicas(3)
        .add_hash_partitions({ "key" }, 5)
        .Create());
  }

  int64_t GetTableLocationCacheMetric(
      int master_idx,
      const MetricPrototype* metric_proto) const {
    CHECK_LT(master_idx, cluster_->num_masters());
    int64_t value;
    CHECK_OK(itest::GetInt64Metric(
        cluster_->master(master_idx)->bound_http_hostport(),
        &METRIC_ENTITY_server, "kudu.master", metric_proto, "value", &value));
    return value;
  }

  int64_t GetTableLocationCacheMetricAllMasters(
      const MetricPrototype* metric_proto) const {
    int64_t total = 0;
    for (auto idx = 0; idx < cluster_->num_masters(); ++idx) {
      const auto val = GetTableLocationCacheMetric(idx, metric_proto);
      total += val;
    }
    return total;
  }

  void CheckCacheMetricsReset(int master_idx) const {
    for (const auto* metric : {
         &METRIC_table_locations_cache_evictions,
         &METRIC_table_locations_cache_hits,
         &METRIC_table_locations_cache_inserts,
         &METRIC_table_locations_cache_lookups,
         &METRIC_table_locations_cache_misses, }) {
      SCOPED_TRACE(Substitute(
          "verifying value of '$0' metric ($1) for master at index $2",
          metric->name(), metric->label(), master_idx));
      const auto val = GetTableLocationCacheMetric(master_idx, metric);
      EXPECT_EQ(0, val);
    }
    const auto cache_memory_usage = GetTableLocationCacheMetric(
        master_idx, &METRIC_table_locations_cache_memory_usage);
    EXPECT_EQ(0, cache_memory_usage);
  }

 protected:
  static constexpr const char* const kTableName = "test_locations_cache_multi_master";
  static constexpr int32_t kRaftHeartbeatIntervalMs = 300;
  static constexpr int32_t kMaxMissedHeartbeatPeriods = 2;
  const KuduSchema schema_;
  std::unique_ptr<cluster::ExternalMiniCluster> cluster_;
};

// Verify that the table location cache is reset upon change once master
// starts its leadership role.
TEST_F(TableLocationsCacheMultiMasterTest, ResetCache) {
  SKIP_IF_SLOW_NOT_ALLOWED();

  // Make sure the cache's metrics are zeroed if no client activity has been
  // there yet.
  for (auto idx = 0; idx < cluster_->num_masters(); ++idx) {
    NO_FATALS(CheckCacheMetricsReset(idx));
  }

  TestWorkload w(cluster_.get());
  w.set_table_name(kTableName);
  w.set_schema(schema_);
  w.Setup();
  w.Start();
  while (w.rows_inserted() < 10) {
    SleepFor(MonoDelta::FromMilliseconds(10));
  }
  w.StopAndJoin();

  // Make sure some items are added into the cache.
  EXPECT_GT(GetTableLocationCacheMetricAllMasters(
      &METRIC_table_locations_cache_inserts), 0);

  int former_leader_master_idx;
  ASSERT_OK(cluster_->GetLeaderMasterIndex(&former_leader_master_idx));

  int leader_master_idx = -1;
  ASSERT_EVENTUALLY([&] {
    // Induce a change in master leadership (maybe, even few of them, up to the
    // number of masters in the cluster).
    for (auto idx = 0; idx < cluster_->num_masters(); ++idx) {
      ASSERT_OK(cluster_->master(idx)->Pause());
      // Make one master to stop sending heartbeats, and give the rest about
      // three heartbeat periods to elect a new leader in case if the stopped
      // master was a leader.
      SleepFor(MonoDelta::FromMilliseconds(
          2 * kRaftHeartbeatIntervalMs * kMaxMissedHeartbeatPeriods +
          3 * kRaftHeartbeatIntervalMs));
      ASSERT_OK(cluster_->master(idx)->Resume());
    }
    ASSERT_OK(cluster_->GetLeaderMasterIndex(&leader_master_idx));
    ASSERT_NE(former_leader_master_idx, leader_master_idx);
  });
  // An extra sanity check.
  ASSERT_NE(-1, leader_master_idx);

  // Make sure all the cache's metrics are reset once master just has become
  // a leader.
  NO_FATALS(CheckCacheMetricsReset(leader_master_idx));
}

} // namespace master
} // namespace kudu

max_stars_repo_name: attilabukor/kudu

max_issues_repo_name: attilabukor/kudu

max_forks_repo_name: attilabukor/kudu


text:
#ifndef OSRM_CUSTOMIZE_EDGE_BASED_GRAPH_HPP
#define OSRM_CUSTOMIZE_EDGE_BASED_GRAPH_HPP

#include "extractor/edge_based_edge.hpp"
#include "partitioner/edge_based_graph.hpp"
#include "partitioner/multi_level_graph.hpp"
#include "util/static_graph.hpp"
#include "util/typedefs.hpp"

#include "storage/shared_memory_ownership.hpp"

#include <boost/filesystem/path.hpp>

namespace osrm
{
namespace customizer
{

using EdgeBasedGraphEdgeData = partitioner::EdgeBasedGraphEdgeData;

struct MultiLevelEdgeBasedGraph
    : public partitioner::MultiLevelGraph<EdgeBasedGraphEdgeData, storage::Ownership::Container>
{
    using Base =
        partitioner::MultiLevelGraph<EdgeBasedGraphEdgeData, storage::Ownership::Container>;
    using Base::Base;
};

struct MultiLevelEdgeBasedGraphView
    : public partitioner::MultiLevelGraph<EdgeBasedGraphEdgeData, storage::Ownership::View>
{
    using Base = partitioner::MultiLevelGraph<EdgeBasedGraphEdgeData, storage::Ownership::View>;
    using Base::Base;
};

struct StaticEdgeBasedGraphEdge : MultiLevelEdgeBasedGraph::InputEdge
{
    using Base = MultiLevelEdgeBasedGraph::InputEdge;
    using Base::Base;
};
}
}

#endif

ext: hpp


text:
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <math.h>
#include <set>
#include <boost/multiprecision/gmp.hpp>
#include <boost/multiprecision/number.hpp>

using namespace std;
using namespace boost::multiprecision;

int target = 100;

int main(int argc, char** argv) {
  set<mpz_int> visited;

  for (int a = 2; a <= target; a++) {
    for (int b = 2; b <= target; b++) {
      mpz_int val = 1;
      for (int i = 0; i < b; i++) {
        val *= a;
      }
      visited.insert(val);
      cout << val << endl;
    }
  }

  cout << "Total unique elements: " << visited.size() << endl;

  return 0;
}

ext: cpp

max_stars_repo_name: DouglasSherk/project-euler

max_issues_repo_name: DouglasSherk/project-euler

max_forks_repo_name: DouglasSherk/project-euler


text:
/*
 Copyright (c) 2010, The Barbarian Group
 All rights reserved.

 Redistribution and use in source and binary forms, with or without modification, are permitted provided that
 the following conditions are met:

    * Redistributions of source code must retain the above copyright notice, this list of conditions and
	the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
	the following disclaimer in the documentation and/or other materials provided with the distribution.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
 WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.
*/

#include "cinder/ip/Flip.h"

#include <boost/preprocessor/seq.hpp>

namespace cinder { namespace ip {

template<typename T>
void flipVertical( SurfaceT<T> *surface )
{
	int32_t rowBytes = surface->getRowBytes();
	uint8_t *buffer = new uint8_t[rowBytes];

	int32_t lastRow = surface->getHeight() - 1;
	int32_t halfHeight = surface->getHeight() / 2;
	for( int32_t y = 0; y < halfHeight; ++y ) {
		memcpy( buffer, surface->getData( Vec2i( 0, y ) ), rowBytes );
		memcpy( surface->getData( Vec2i( 0, y ) ), surface->getData( Vec2i( 0, lastRow - y ) ), rowBytes );
		memcpy( surface->getData( Vec2i( 0, lastRow - y ) ), buffer, rowBytes );
	}

	delete [] buffer;
}

#define flip_PROTOTYPES(r,data,T)\
	template void flipVertical<T>( SurfaceT<T> *surface );

BOOST_PP_SEQ_FOR_EACH( flip_PROTOTYPES, ~, CHANNEL_TYPES )

} } // namespace cinder::ip

max_stars_repo_name: timmb/HarmonicMotion

max_issues_repo_name: timmb/HarmonicMotion

max_forks_repo_name: timmb/HarmonicMotion


text:
#include <Eigen/Dense>
#include <utility>
#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

#include "general_nonlinear_solver.h" // implements general_nonlinear_solver

//! \brief Implements a single step of the modified newton
//! \tparam argument type of argument to function f: such as double or vector etc...
//! \tparam function type for the function f, likely a lambda function
//! \tparam jacobian type for the jacobian df, likely a lambda function
//! \param[in] x previous value to use in step, also initial guess if needed
//! \param[out] x_next next step x_{k+1} of modified Newton
//! \param[in] f function handle for f(x) = 0
//! \param[in] df function handle for jacobian df of f
template <typename argument, class function, class jacobian>
void mod_newt_step_scalar(const argument & x, argument & x_next, function&& f, jacobian&& df) {
    argument y = x + f(x) / df(x);
    x_next = y - f(y) / df(x);
}

//! \brief Implements a single step of the modified newton
//! \tparam argument type of argument to function f: such as double or vector etc...
//! \tparam function type for the function f, likely a lambda function
//! \tparam jacobian type for the jacobian df, likely a lambda function
//! \param[in] x previous value to use in step, also initial guess if needed
//! \param[out] x_next next step x_{k+1} of modified Newton
//! \param[in] f function handle for f(x) = 0
//! \param[in] df function handle for jacobian df of f
template <typename argument, class function, class jacobian>
void mod_newt_step_system(const argument & x, argument & x_next, function&& f, jacobian&& df) {
    auto lu = df(x).lu();
    // reusing LU decomposition in the latter
    // the constructor performs LU decomposition in case argument = Vector
    argument y = x + lu.solve(f(x));
    x_next = y - lu.solve(f(y));
}

//! \brief Solve a scalar non-linear eq. with the modified Newton
void mod_newt_ord() {
    // Setting up values, functions and jacobian
    const double a = 0.123;
    auto f = [&a] (double x) { return atan(x) - a; }; // function, must see a
    auto df = [] (double x) { return 1. / (x*x+1.); };

    const double x_ex = tan(a); // exact solution

    // store solution and error at each time step
    std::vector<double> sol;
    std::vector<double> err;

    // Compute error and push back to err, used in general_nonlinear solver as breaking condition errf(x) < eps
    auto errf = [&err, x_ex] (double & x) { double e =  std::abs(x - x_ex); err.push_back(e); return e; };

    // Perform convergence study with Modified newton for scalar
    std::cout << std::endl << "*** Modified Newton method (scalar) ***" << std::endl << std::endl;
    std::cout << "Exact: " << x_ex << std::endl;
    // Initial guess and compute initial error
    double x_scalar = 5;
    sol.push_back(x_scalar);
    errf(x_scalar);
    // Lambda performing the next step, used to define a proper function handle to be passed to general_nonlinear_solver
    auto newt_scalar_step = [&sol, &f, &df] (double x, double & x_new) { mod_newt_step_scalar(x, x_new, f, df); sol.push_back(x_new); };
    // Actually perform the solution
    general_nonlinear_solver(newt_scalar_step, x_scalar, errf);
    // Print solution (final)
    std::cout << std::endl << "x^*_scalar = " << std::endl << x_scalar << std::endl;

    // Print table of solutions, errors and EOOC
    auto space = std::setw(15);
    std::cout << std::endl << space << "sol." << space << "err." << space << "order" << std::endl;
    for(unsigned i = 0; i < sol.size(); ++ i) {
        std::cout << space << sol.at(i) << space << err.at(i);
        if(i >= 3) std::cout << space << (log(err.at(i)) - log(err.at(i-1))) / (log(err.at(i-1)) - log(err.at(i-2)));
        std::cout << std::endl;
    }
}

//! \brief Solve a system non-linear eq. with the modified Newton
void mod_newt_sys() {
    //// Setup initial values
    using Vector = Eigen::VectorXd;
    using Matrix = Eigen::MatrixXd;

    Matrix A = Matrix::Random(4,4);
    Vector c = Vector::Random(4);

    A = (A*A.transpose()).eval(); // make sure A is symmetric
    c = c.cwiseAbs(); // make sure c is > 0 uniformly

    // Handler for function and jacobian in standard format. Must see matrix A and vector c
    auto F = [&A,&c] (const Vector & x) { Vector tmp =  A*x + c.cwiseProduct(x.array().exp().matrix()).eval(); return tmp; };
    std::function<Matrix(const Vector &)> dF = [&A, &c] (const Vector & x) { Matrix C = A; Vector temp = c.cwiseProduct(x.array().exp().matrix()); C += temp.asDiagonal(); return C; };

    // Container for errors
    std::vector<double> err;
    // Define lambda for breaking condition, which also stores the error of the previous step. Must see err
    auto rerr = [&err] (Vector & x) { if(err.size() > 0) return err.back(); else return 1.; };

    std::cout << std::endl << "*** Modified Newton method (system) ***" << std::endl << std::endl;
    Vector x_system = Vector::Zero(c.size()); // Initial guess

    // Refactor the step st. is compatible with general_nonlinear_solver,  must see F; dF and error vector
    auto newt_system_step = [&F, &dF, &err] (const Vector & x, Vector & x_new) { mod_newt_step_system(x, x_new, F, dF); double e = (x - x_new).norm()/ x_new.norm(); err.push_back(e); };
    // Actually performs computations
    general_nonlinear_solver(newt_system_step, x_system, rerr);

    // Sanity check
    std::cout << std::endl << "x^*_system = " << std::endl << x_system << std::endl;
    std::cout << std::endl << "F(x^*_system) = " << std::endl << F(x_system) << std::endl;

}

int main() {
    // Part 1: solve scalar
    mod_newt_ord();

    // Part 2: solve system
    mod_newt_sys();

}

max_stars_repo_path: Nummerical Methods for CSE/PS5/solutions_ps5/modnewt.cpp

max_stars_repo_name: valentinjacot/backupETHZ

max_issues_repo_name: valentinjacot/backupETHZ

max_forks_repo_name: valentinjacot/backupETHZ


text: #include <boost/pool/simple_segregated_storage.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/* Copyright (c) 2016 - 2020, the adamantine authors.
 *
 * This file is subject to the Modified BSD License and may not be distributed
 * without copyright and license information. Please refer to the file LICENSE
 * for the text and further information on this license.
 */

#define BOOST_TEST_MODULE ThermalOperator

#include <Geometry.hh>
#include <ThermalOperator.hh>

#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/numerics/matrix_tools.h>

#include <boost/property_tree/ptree.hpp>

#include "main.cc"

BOOST_AUTO_TEST_CASE(thermal_operator)
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Create the Geometry
  boost::property_tree::ptree geometry_database;
  geometry_database.put("length", 12);
  geometry_database.put("length_divisions", 4);
  geometry_database.put("height", 6);
  geometry_database.put("height_divisions", 5);
  adamantine::Geometry<2> geometry(communicator, geometry_database);
  // Create the DoFHandler
  dealii::FE_Q<2> fe(2);
  dealii::DoFHandler<2> dof_handler(geometry.get_triangulation());
  dof_handler.distribute_dofs(fe);
  dealii::AffineConstraints<double> affine_constraints;
  affine_constraints.close();
  dealii::QGauss<1> quad(3);

  // Create the MaterialProperty
  boost::property_tree::ptree mat_prop_database;
  mat_prop_database.put("property_format", "polynomial");
  mat_prop_database.put("n_materials", 1);
  mat_prop_database.put("material_0.solid.density", 1.);
  mat_prop_database.put("material_0.powder.density", 1.);
  mat_prop_database.put("material_0.liquid.density", 1.);
  mat_prop_database.put("material_0.solid.specific_heat", 1.);
  mat_prop_database.put("material_0.powder.specific_heat", 1.);
  mat_prop_database.put("material_0.liquid.specific_heat", 1.);
  mat_prop_database.put("material_0.solid.thermal_conductivity", 10.);
  mat_prop_database.put("material_0.powder.thermal_conductivity", 10.);
  mat_prop_database.put("material_0.liquid.thermal_conductivity", 10.);
  std::shared_ptr<adamantine::MaterialProperty<2>> mat_properties(
      new adamantine::MaterialProperty<2>(
          communicator, geometry.get_triangulation(), mat_prop_database));

  // Initialize the ThermalOperator
  adamantine::ThermalOperator<2, 2, dealii::MemorySpace::Host> thermal_operator(
      communicator, mat_properties);
  thermal_operator.reinit(dof_handler, affine_constraints, quad);
  thermal_operator.compute_inverse_mass_matrix(dof_handler, affine_constraints);
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dummy(
      thermal_operator.m());
  thermal_operator.evaluate_material_properties(dummy);
  BOOST_CHECK(thermal_operator.m() == 99);
  BOOST_CHECK(thermal_operator.m() == thermal_operator.n());

  // Check matrix-vector multiplications
  double const tolerance = 1e-15;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> src;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dst_1;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dst_2;

  dealii::MatrixFree<2, double> const &matrix_free =
      thermal_operator.get_matrix_free();
  matrix_free.initialize_dof_vector(src);
  matrix_free.initialize_dof_vector(dst_1);
  matrix_free.initialize_dof_vector(dst_2);

  src = 1.;
  thermal_operator.vmult(dst_1, src);
  BOOST_CHECK_CLOSE(dst_1.l1_norm(), 0., tolerance);

  thermal_operator.Tvmult(dst_2, src);
  BOOST_CHECK_CLOSE(dst_2.l1_norm(), dst_1.l1_norm(), tolerance);

  dst_2 = 1.;
  thermal_operator.vmult_add(dst_2, src);
  thermal_operator.vmult(dst_1, src);
  dst_1 += src;
  BOOST_CHECK_CLOSE(dst_1.l1_norm(), dst_2.l1_norm(), tolerance);

  dst_1 = 1.;
  thermal_operator.Tvmult_add(dst_1, src);
  BOOST_CHECK_CLOSE(dst_1.l1_norm(), dst_2.l1_norm(), tolerance);
}

BOOST_AUTO_TEST_CASE(spmv)
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Create the Geometry
  boost::property_tree::ptree geometry_database;
  geometry_database.put("length", 12);
  geometry_database.put("length_divisions", 4);
  geometry_database.put("height", 6);
  geometry_database.put("height_divisions", 5);
  adamantine::Geometry<2> geometry(communicator, geometry_database);
  // Create the DoFHandler
  dealii::FE_Q<2> fe(2);
  dealii::DoFHandler<2> dof_handler(geometry.get_triangulation());
  dof_handler.distribute_dofs(fe);
  dealii::AffineConstraints<double> affine_constraints;
  affine_constraints.close();
  dealii::QGauss<1> quad(3);

  // Create the MaterialProperty
  boost::property_tree::ptree mat_prop_database;
  mat_prop_database.put("property_format", "polynomial");
  mat_prop_database.put("n_materials", 1);
  mat_prop_database.put("material_0.solid.density", 1.);
  mat_prop_database.put("material_0.powder.density", 1.);
  mat_prop_database.put("material_0.liquid.density", 1.);
  mat_prop_database.put("material_0.solid.specific_heat", 1.);
  mat_prop_database.put("material_0.powder.specific_heat", 1.);
  mat_prop_database.put("material_0.liquid.specific_heat", 1.);
  mat_prop_database.put("material_0.solid.thermal_conductivity", 1.);
  mat_prop_database.put("material_0.powder.thermal_conductivity", 1.);
  mat_prop_database.put("material_0.liquid.thermal_conductivity", 1.);
  std::shared_ptr<adamantine::MaterialProperty<2>> mat_properties(
      new adamantine::MaterialProperty<2>(
          communicator, geometry.get_triangulation(), mat_prop_database));

  // Initialize the ThermalOperator
  adamantine::ThermalOperator<2, 2, dealii::MemorySpace::Host> thermal_operator(
      communicator, mat_properties);
  thermal_operator.reinit(dof_handler, affine_constraints, quad);
  thermal_operator.compute_inverse_mass_matrix(dof_handler, affine_constraints);
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dummy(
      thermal_operator.m());
  thermal_operator.evaluate_material_properties(dummy);
  BOOST_CHECK(thermal_operator.m() == 99);
  BOOST_CHECK(thermal_operator.m() == thermal_operator.n());

  // Build the matrix. This only works in serial.
  dealii::DynamicSparsityPattern dsp(dof_handler.n_dofs());
  dealii::DoFTools::make_sparsity_pattern(dof_handler, dsp, affine_constraints);
  dealii::SparsityPattern sparsity_pattern;
  sparsity_pattern.copy_from(dsp);
  dealii::SparseMatrix<double> sparse_matrix(sparsity_pattern);
  dealii::MatrixCreator::create_laplace_matrix(
      dof_handler, dealii::QGauss<2>(3), sparse_matrix);

  // Compare vmult using matrix free and building the matrix
  double const tolerance = 1e-12;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> src;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dst_1;
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> dst_2;

  dealii::MatrixFree<2, double> const &matrix_free =
      thermal_operator.get_matrix_free();
  matrix_free.initialize_dof_vector(src);
  matrix_free.initialize_dof_vector(dst_1);
  matrix_free.initialize_dof_vector(dst_2);

  for (unsigned int i = 0; i < thermal_operator.m(); ++i)
  {
    src = 0.;
    src[i] = 1;
    thermal_operator.vmult(dst_1, src);
    sparse_matrix.vmult(dst_2, src);
    for (unsigned int j = 0; j < thermal_operator.m(); ++j)
      BOOST_CHECK_CLOSE(dst_1[j], -dst_2[j], tolerance);
  }
}

max_stars_repo_name: stvdwtt/adamantine

max_issues_repo_name: stvdwtt/adamantine

max_forks_repo_name: stvdwtt/adamantine


text:
//
// The MIT License(MIT)
//
// Copyright(c) 2014 Demonsaw LLC
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#include <algorithm>
#include <iostream>
#include <boost/algorithm/string.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/lexical_cast.hpp>

#include "http_message.h"
#include "system/type.h"
#include "utility/boost.h"

namespace eja
{
	// Constructor
	http_message::http_message()
	{
		set_connection(http::keep_alive);
		set_type(http::text);
	}

	// Interface
	void http_message::clear()
	{
		m_body.clear();
		m_headers.clear();
		m_cookies.clear();

		set_connection(http::keep_alive);
		set_type(http::text);
	}

	// Utility
	std::string http_message::str() const
	{
		std::ostringstream sstream;

		for (const auto& pair : m_headers)
			sstream << pair.first << ": " << pair.second << http::eol;

		return sstream.str();
	}

	// Cookie
	bool http_message::is_cookie(const std::string& key) const
	{
		for (const auto& cookie : http::cookies)
		{
			if (boost::starts_with(key, cookie))
				return true;
		}

		return false;
	}

	void http_message::set_cookie(const std::string& data)
	{
		size_t pos = data.find_first_of(':');
		const std::string suffix = (pos != std::string::npos) ? data.substr(pos + 1) : data;

		boost::char_separator<char> separator(";");
		boost::string_tokenizer tokens(suffix, separator);

		for (const auto& line : tokens)
		{
			int pos = line.find_first_of('=');
			if (pos == std::string::npos)
				continue;

			key_type key = line.substr(0, pos);
			boost::trim(key);

			value_type value = line.substr(pos + 1);
			boost::trim(value);

			// TODO: Do we need to decode this?
			if (key == http::session)
				set_session(value);
			else
				m_cookies[key] = value;
		}
	}

	// Header
	void http_message::set_header(const std::string& data)
	{
		int pos = data.find_first_of(':');
		if (pos == std::string::npos)
			return;

		key_type key = data.substr(0, pos);
		boost::trim(key);

		value_type value = data.substr(pos + 1);
		boost::trim(value);

		m_headers[key] = value;
	}

	// Mutators
	void http_message::set_body(const char* body, const size_t body_size)
	{
		m_body.assign(body, body_size);;

		set_length(m_body.size());
	}
}

max_stars_repo_name: demonsaw/Code

max_issues_repo_name: demonsaw/Code

max_forks_repo_name: demonsaw/Code