FillBetweenPointVectorsTools Namespace Reference

Functions
void	fillBetweenPointVectorsGenerator (MeshBase &mesh, const std::vector< Point > &boundary_points_vec_1, const std::vector< Point > &boundary_points_vec_2, const unsigned int num_layers, const subdomain_id_type transition_layer_id, const boundary_id_type input_boundary_1_id, const boundary_id_type input_boundary_2_id, const boundary_id_type begin_side_boundary_id, const boundary_id_type end_side_boundary_id, const std::string type, const std::string name, const bool quad_elem=false, const Real bias_parameter=1.0, const Real sigma=3.0)
	Generates a 2D mesh with triangular elements for a region defined by two curves (sets of Points)
void	fillBetweenPointVectorsGenerator (MeshBase &mesh, const std::vector< Point > &boundary_points_vec_1, const std::vector< Point > &boundary_points_vec_2, const unsigned int num_layers, const subdomain_id_type transition_layer_id, const boundary_id_type external_boundary_id, const std::string type, const std::string name, const bool quad_elem=false)
	Generates a 2D mesh with triangular elements for a region defined by two curves (sets of Points)
void	elementsCreationFromNodesVectorsQuad (MeshBase &mesh, const std::vector< std::vector< Node * > > &nodes, const unsigned int num_layers, const std::vector< unsigned int > &node_number_vec, const subdomain_id_type transition_layer_id, const boundary_id_type input_boundary_1_id, const boundary_id_type input_boundary_2_id, const boundary_id_type begin_side_boundary_id, const boundary_id_type end_side_boundary_id)
	Generates a 2D mesh based on a 2D vector of Nodes using QUAD4 elements in the xy-plane.
void	elementsCreationFromNodesVectors (MeshBase &mesh, const std::vector< std::vector< Node * > > &nodes, const unsigned int num_layers, const std::vector< unsigned int > &node_number_vec, const subdomain_id_type transition_layer_id, const boundary_id_type input_boundary_1_id, const boundary_id_type input_boundary_2_id, const boundary_id_type begin_side_boundary_id, const boundary_id_type end_side_boundary_id)
	Generates a 2D mesh based on a 2D vector of Nodes using TRI3 elements in the xy-plane.
void	weightedInterpolator (const unsigned int vec_node_num, const std::vector< Point > &boundary_points_vec, std::vector< Real > &vec_index, std::vector< Real > &wt, std::vector< Real > &index, const Real sigma, std::unique_ptr< LinearInterpolation > &linear_vec_x, std::unique_ptr< LinearInterpolation > &linear_vec_y, std::unique_ptr< SplineInterpolation > &spline_vec_l)
	Generates weights, weighted indices and corresponding interpolation.
void	surrogateGenerator (std::vector< Real > &weighted_surrogate_index, std::vector< Real > &unweighted_surrogate_index, const std::vector< unsigned int > &node_number_vec, const std::vector< Real > &index, const std::vector< Real > &wt, const unsigned int boundary_node_num, const unsigned int i)
	Generates weighted surrogate index vectors on one side for points of a sublayer curve.
bool	needFlip (const std::vector< Point > &vec_pts_1, const std::vector< Point > &vec_pts_2)
	Decide whether one of the input vector of Points needs to be flipped to ensure correct transition layer shape.
bool	isBoundarySimpleClosedLoop (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &boundary_ordered_node_list, const Point origin_pt, const boundary_id_type bid)
	Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isBoundarySimpleClosedLoop (MeshBase &mesh, Real &max_node_radius, const Point origin_pt, const boundary_id_type bid)
	Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isBoundarySimpleClosedLoop (MeshBase &mesh, const Point origin_pt, const boundary_id_type bid)
	Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isBoundaryOpenSingleSegment (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &boundary_ordered_node_list, const Point origin_pt, const boundary_id_type bid)
	Decides whether a boundary of a given mesh is an open single-segment boundary.
bool	isCurveSimpleClosedLoop (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &ordered_node_list, const Point origin_pt)
	Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isCurveSimpleClosedLoop (MeshBase &mesh, Real &max_node_radius, const Point origin_pt)
	Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isCurveSimpleClosedLoop (MeshBase &mesh, const Point origin_pt)
	Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isCurveOpenSingleSegment (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &ordered_node_list, const Point origin_pt)
	Decides whether a curve contained in a given mesh is an open single-segment curve.
void	isClosedLoop (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &ordered_node_list, std::vector< std::pair< dof_id_type, dof_id_type > > &node_assm, std::vector< dof_id_type > &midpoint_node_list, const Point origin_pt, const std::string input_type, bool &is_closed_loop, const bool suppress_exception=false)
	Decides whether a series of nodes contained in a given mesh forms a closed loop with consecutive nodes's azimuthal angles change monotonically.
void	isClosedLoop (MeshBase &mesh, Real &max_node_radius, std::vector< dof_id_type > &ordered_node_list, std::vector< std::pair< dof_id_type, dof_id_type > > &node_assm, const Point origin_pt, const std::string input_type, bool &is_closed_loop, const bool suppress_exception=false)
	Decides whether a series of nodes contained in a given mesh forms a closed loop with consecutive nodes's azimuthal angles change monotonically.
bool	isExternalBoundary (MeshBase &mesh, const boundary_id_type bid)
	Decides whether a boundary of a given mesh is an external boundary.
bool	buildQuadElement (Elem *elem, Node *nd_0, Node *nd_1, Node *nd_2, Node *nd_3, const subdomain_id_type transition_layer_id)
	Creates an QUAD4 element that can be extruded in (0 0 1) direction.
bool	buildTriElement (Elem *elem, Node *nd_0, Node *nd_1, Node *nd_2, const subdomain_id_type transition_layer_id)
	Creates an TRI3 element that can be extruded in (0 0 1) direction.

Function Documentation

buildQuadElement()

bool FillBetweenPointVectorsTools::buildQuadElement	(	Elem *	elem,
		Node *	nd_0,
		Node *	nd_1,
		Node *	nd_2,
		Node *	nd_3,
		const subdomain_id_type	transition_layer_id
	)

Creates an QUAD4 element that can be extruded in (0 0 1) direction.

Parameters

elem	pointer of the element to be created
nd_0	pointer of element's first node
nd_1	pointer of element's second node
nd_2	pointer of element's third node
nd_3	pointer of element's fourth node
transition_layer_id	id of the subdomain that this element belongs to

Returns

whether the element is flipped to ensure (0 0 1) extrudability

Definition at line 748 of file FillBetweenPointVectorsTools.C.

{
  // Adjust the order of nodes in an element so that the mesh can be extruded in (0 0 1)
  // direction.
  elem->subdomain_id() = transition_layer_id;
  if (((*nd_1 - *nd_0).cross((*nd_2 - *nd_0)).unit())(2) > 0)
  {
    elem->set_node(0, nd_0);
    elem->set_node(1, nd_1);
    elem->set_node(2, nd_2);
    elem->set_node(3, nd_3);
    return false;
  }
  else
  {
    elem->set_node(0, nd_0);
    elem->set_node(3, nd_1);
    elem->set_node(2, nd_2);
    elem->set_node(1, nd_3);
    return true;
  }
}

Referenced by elementsCreationFromNodesVectorsQuad().

buildTriElement()

bool FillBetweenPointVectorsTools::buildTriElement	(	Elem *	elem,
		Node *	nd_0,
		Node *	nd_1,
		Node *	nd_2,
		const subdomain_id_type	transition_layer_id
	)

Creates an TRI3 element that can be extruded in (0 0 1) direction.

Parameters

elem	pointer of the element to be created
nd_0	pointer of element's first node
nd_1	pointer of element's second node
nd_2	pointer of element's third node
transition_layer_id	id of the subdomain that this element belongs to

Returns

whether the element is flipped to ensure (0 0 1) extrudability

Definition at line 777 of file FillBetweenPointVectorsTools.C.

{
  // Adjust the order of nodes in an element so that the mesh can be extruded in (0 0 1)
  // direction.
  elem->subdomain_id() = transition_layer_id;
  if (((*nd_1 - *nd_0).cross((*nd_2 - *nd_0)).unit())(2) > 0)
  {
    elem->set_node(0, nd_0);
    elem->set_node(1, nd_1);
    elem->set_node(2, nd_2);
    return false;
  }
  else
  {
    elem->set_node(0, nd_0);
    elem->set_node(2, nd_1);
    elem->set_node(1, nd_2);
    return true;
  }
}

Referenced by elementsCreationFromNodesVectors().

elementsCreationFromNodesVectors()

void FillBetweenPointVectorsTools::elementsCreationFromNodesVectors	(	MeshBase &	mesh,
		const std::vector< std::vector< Node * > > &	nodes,
		const unsigned int	num_layers,
		const std::vector< unsigned int > &	node_number_vec,
		const subdomain_id_type	transition_layer_id,
		const boundary_id_type	input_boundary_1_id,
		const boundary_id_type	input_boundary_2_id,
		const boundary_id_type	begin_side_boundary_id,
		const boundary_id_type	end_side_boundary_id
	)

Generates a 2D mesh based on a 2D vector of Nodes using TRI3 elements in the xy-plane.

Parameters

mesh	a reference MeshBase to contain the generated mesh
nodes	a 2D vector of nodes based on which the mesh is built
num_layers	number of sublayers of elements
node_number_vec	number of nodes on each sublayer
transition_layer_id	subdomain id of the generated mesh
input_boundary_1_id	id of the generated mesh's external boundary that originates from boundary_points_vec_1
input_boundary_2_id	id of the generated mesh's external boundary that originates from boundary_points_vec_2
begin_side_boundary_id	id of the generated mesh's external boundary that connects the first Points of the two input Point vectors
end_side_boundary_id	id of the generated mesh's external boundary that connects the last Points of the two input Point vectors

Definition at line 285 of file FillBetweenPointVectorsTools.C.

{
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
  for (const auto i : make_range(num_layers))
  {
    unsigned int nodes_up_it = 0;
    unsigned int nodes_down_it = 0;
    const unsigned int node_number_up = node_number_vec[i + 1];
    const unsigned int node_number_down = node_number_vec[i];
 
    while (nodes_up_it < node_number_up - 1 && nodes_down_it < node_number_down - 1 &&
           nodes_up_it + nodes_down_it < node_number_up + node_number_down - 3)
    {
      // Define the two possible options and chose the one with shorter distance
      Real dis1 = (*nodes[i + 1][nodes_up_it] - *nodes[i][nodes_down_it + 1]).norm();
      Real dis2 = (*nodes[i + 1][nodes_up_it + 1] - *nodes[i][nodes_down_it]).norm();
      if (MooseUtils::absoluteFuzzyGreaterThan(dis1, dis2))
      {
        Elem * elem = mesh.add_elem(new Tri3);
        bool is_elem_flip = buildTriElement(elem,
                                            nodes[i + 1][nodes_up_it],
                                            nodes[i][nodes_down_it],
                                            nodes[i + 1][nodes_up_it + 1],
                                            transition_layer_id);
        if (i == num_layers - 1)
          boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
        if (nodes_up_it == 0 && nodes_down_it == 0)
          boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
        nodes_up_it++;
      }
      else
      {
        Elem * elem = mesh.add_elem(new Tri3);
        bool is_elem_flip = buildTriElement(elem,
                                            nodes[i + 1][nodes_up_it],
                                            nodes[i][nodes_down_it],
                                            nodes[i][nodes_down_it + 1],
                                            transition_layer_id);
        if (i == 0)
          boundary_info.add_side(elem, 1, input_boundary_1_id);
        if (nodes_up_it == 0 && nodes_down_it == 0)
          boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
        nodes_down_it++;
      }
    }
    // Handle the end
    while (nodes_up_it < node_number_up - 1)
    {
      Elem * elem = mesh.add_elem(new Tri3);
      bool is_elem_flip = buildTriElement(elem,
                                          nodes[i + 1][nodes_up_it],
                                          nodes[i][nodes_down_it],
                                          nodes[i + 1][nodes_up_it + 1],
                                          transition_layer_id);
      nodes_up_it++;
      if (i == num_layers - 1)
        boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
      if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
        boundary_info.add_side(elem, 1, end_side_boundary_id);
    }
    while (nodes_down_it < node_number_down - 1)
    {
      Elem * elem = mesh.add_elem(new Tri3);
      bool is_elem_flip = buildTriElement(elem,
                                          nodes[i + 1][nodes_up_it],
                                          nodes[i][nodes_down_it],
                                          nodes[i][nodes_down_it + 1],
                                          transition_layer_id);
      nodes_down_it++;
      if (i == 0)
        boundary_info.add_side(elem, 1, input_boundary_1_id);
      if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
        boundary_info.add_side(elem, is_elem_flip ? 0 : 2, end_side_boundary_id);
    }
  }
}

Referenced by fillBetweenPointVectorsGenerator().

elementsCreationFromNodesVectorsQuad()

void FillBetweenPointVectorsTools::elementsCreationFromNodesVectorsQuad	(	MeshBase &	mesh,
		const std::vector< std::vector< Node * > > &	nodes,
		const unsigned int	num_layers,
		const std::vector< unsigned int > &	node_number_vec,
		const subdomain_id_type	transition_layer_id,
		const boundary_id_type	input_boundary_1_id,
		const boundary_id_type	input_boundary_2_id,
		const boundary_id_type	begin_side_boundary_id,
		const boundary_id_type	end_side_boundary_id
	)

Generates a 2D mesh based on a 2D vector of Nodes using QUAD4 elements in the xy-plane.

Parameters

mesh	a reference MeshBase to contain the generated mesh
nodes	a 2D vector of nodes based on which the mesh is built
num_layers	number of sublayers of elements
node_number_vec	number of nodes on each sublayer
transition_layer_id	subdomain id of the generated mesh
input_boundary_1_id	id of the generated mesh's external boundary that originates from boundary_points_vec_1
input_boundary_2_id	id of the generated mesh's external boundary that originates from boundary_points_vec_2
begin_side_boundary_id	id of the generated mesh's external boundary that connects the first Points of the two input Point vectors
end_side_boundary_id	id of the generated mesh's external boundary that connects the last Points of the two input Point vectors

Definition at line 250 of file FillBetweenPointVectorsTools.C.

{
  const unsigned int node_number = node_number_vec.front();
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
  for (const auto i : make_range(num_layers))
    for (unsigned int j = 1; j < node_number; j++)
    {
      Elem * elem = mesh.add_elem(new Quad4);
      bool is_elem_flip = buildQuadElement(elem,
                                           nodes[i][j - 1],
                                           nodes[i + 1][j - 1],
                                           nodes[i + 1][j],
                                           nodes[i][j],
                                           transition_layer_id);
      if (i == 0)
        boundary_info.add_side(elem, is_elem_flip ? 0 : 3, input_boundary_1_id);
      if (i == num_layers - 1)
        boundary_info.add_side(elem, is_elem_flip ? 2 : 1, input_boundary_2_id);
      if (j == 1)
        boundary_info.add_side(elem, is_elem_flip ? 3 : 0, begin_side_boundary_id);
      if (j == node_number - 1)
        boundary_info.add_side(elem, is_elem_flip ? 1 : 2, end_side_boundary_id);
    }
}

Referenced by fillBetweenPointVectorsGenerator().

fillBetweenPointVectorsGenerator() [1/2]

void FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator	(	MeshBase &	mesh,
		const std::vector< Point > &	boundary_points_vec_1,
		const std::vector< Point > &	boundary_points_vec_2,
		const unsigned int	num_layers,
		const subdomain_id_type	transition_layer_id,
		const boundary_id_type	external_boundary_id,
		const std::string	type,
		const std::string	name,
		const bool	quad_elem = false
	)

Generates a 2D mesh with triangular elements for a region defined by two curves (sets of Points)

Parameters

mesh	a reference MeshBase to contain the generated mesh
boundary_points_vec_1	vector of Points of Side #1
boundary_points_vec_2	vector of Points of Side #2
num_layers	number of sublayers of elements
transition_layer_id	subdomain id of the generated mesh
external_boundary_id	id of the generated mesh's external boundary
type	type of the MOOSE object that calls this method
name	name of the MOOSE object that calls this method
quad_elem	whether the QUAD4 elements are used to construct the mesh

Definition at line 225 of file FillBetweenPointVectorsTools.C.

{
  fillBetweenPointVectorsGenerator(mesh,
                                   boundary_points_vec_1,
                                   boundary_points_vec_2,
                                   num_layers,
                                   transition_layer_id,
                                   external_boundary_id,
                                   external_boundary_id,
                                   external_boundary_id,
                                   external_boundary_id,
                                   type,
                                   name,
                                   quad_elem);
}

fillBetweenPointVectorsGenerator() [2/2]

void FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator	(	MeshBase &	mesh,
		const std::vector< Point > &	boundary_points_vec_1,
		const std::vector< Point > &	boundary_points_vec_2,
		const unsigned int	num_layers,
		const subdomain_id_type	transition_layer_id,
		const boundary_id_type	input_boundary_1_id,
		const boundary_id_type	input_boundary_2_id,
		const boundary_id_type	begin_side_boundary_id,
		const boundary_id_type	end_side_boundary_id,
		const std::string	type,
		const std::string	name,
		const bool	quad_elem = false,
		const Real	bias_parameter = 1.0,
		const Real	sigma = 3.0
	)

Generates a 2D mesh with triangular elements for a region defined by two curves (sets of Points)

Parameters

mesh	a reference MeshBase to contain the generated mesh
boundary_points_vec_1	vector of Points of Side #1
boundary_points_vec_2	vector of Points of Side #2
num_layers	number of sublayers of elements
transition_layer_id	subdomain id of the generated mesh
input_boundary_1_id	id of the generated mesh's external boundary that originates from boundary_points_vec_1
input_boundary_2_id	id of the generated mesh's external boundary that originates from boundary_points_vec_2
begin_side_boundary_id	id of the generated mesh's external boundary that connects the first Points of the two input Point vectors
end_side_boundary_id	id of the generated mesh's external boundary that connects the last Points of the two input Point vectors
type	type of the MOOSE object that calls this method for error messages
name	name of the MOOSE object that calls this method for error messages
quad_elem	whether the QUAD4 elements are used to construct the mesh
bias_parameter	a parameter to control bias options (1) >0, biasing factor (2) <=0 automatic biasing
sigma	Gaussian parameter used for Gaussian blurring of local node density; only used if bias_parameter <= 0

Definition at line 30 of file FillBetweenPointVectorsTools.C.

{
  if (!MooseMeshUtils::isCoPlanar(
          boundary_points_vec_1, Point(0.0, 0.0, 1.0), Point(0.0, 0.0, 0.0)) ||
      !MooseMeshUtils::isCoPlanar(
          boundary_points_vec_2, Point(0.0, 0.0, 1.0), Point(0.0, 0.0, 0.0)))
    mooseError("In ",
               type,
               " ",
               name,
               ", the input vectors of points for "
               "FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator "
               "must be in XY plane.");
 
  const unsigned int vec_1_node_num = boundary_points_vec_1.size();
  const unsigned int vec_2_node_num = boundary_points_vec_2.size();
 
  if (vec_1_node_num < 2 || vec_2_node_num < 2)
    mooseError("In ",
               type,
               " ",
               name,
               ", the two input vectors of points for "
               "FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator "
               "must respectively contain at least two elements.");
 
  if (quad_elem && boundary_points_vec_1.size() != boundary_points_vec_2.size())
    mooseError("In ",
               type,
               " ",
               name,
               ", QUAD4 elements option can only be selected when the two input vectors of Points "
               "have the same length. In the current instance, the first vector has ",
               boundary_points_vec_1.size(),
               " points and the second ",
               boundary_points_vec_2.size(),
               " points");
 
  std::vector<Point> possibly_reoriented_boundary_points_vec_2;
  const std::vector<Point> * oriented_boundary_points_vec_2 = &boundary_points_vec_2;
 
  if (needFlip(boundary_points_vec_1, boundary_points_vec_2))
  {
    possibly_reoriented_boundary_points_vec_2.assign(boundary_points_vec_2.rbegin(),
                                                     boundary_points_vec_2.rend());
    oriented_boundary_points_vec_2 = &possibly_reoriented_boundary_points_vec_2;
 
    // This isn't worth warning about.  The way
    // MooseMeshUtils::makeOrderedNodeList works, we can end up
    // finding a flip necessary on one element numbering and
    // unnecessary on another.
    //
    // mooseWarning(
    //     "In FillBetweenPointVectorsTools, one of the vector of Points must be flipped to ensure "
    //     "correct transition layer shape.");
  }
 
  std::vector<Real> vec_1_index; // Unweighted index
  std::vector<Real> vec_2_index; // Unweighted index
 
  std::vector<Real> wt_1;
  std::vector<Real> index_1; // Weighted index
  std::vector<Real> wt_2;
  std::vector<Real> index_2; // Weighted index
 
  // create interpolations
  std::unique_ptr<LinearInterpolation> linear_vec_1_x;
  std::unique_ptr<LinearInterpolation> linear_vec_1_y;
  std::unique_ptr<SplineInterpolation> spline_vec_1_l;
  std::unique_ptr<LinearInterpolation> linear_vec_2_x;
  std::unique_ptr<LinearInterpolation> linear_vec_2_y;
  std::unique_ptr<SplineInterpolation> spline_vec_2_l;
 
  weightedInterpolator(vec_1_node_num,
                       boundary_points_vec_1,
                       vec_1_index,
                       wt_1,
                       index_1,
                       sigma,
                       linear_vec_1_x,
                       linear_vec_1_y,
                       spline_vec_1_l);
  weightedInterpolator(vec_2_node_num,
                       *oriented_boundary_points_vec_2,
                       vec_2_index,
                       wt_2,
                       index_2,
                       sigma,
                       linear_vec_2_x,
                       linear_vec_2_y,
                       spline_vec_2_l);
 
  // If the two input vectors have different sizes
  // The node numbers of the intermediate layers change linearly
  const Real increment = ((Real)vec_2_node_num - (Real)vec_1_node_num) / (Real)(num_layers);
  // Number of nodes in each sublayer
  std::vector<unsigned int> node_number_vec;
  // 2D vector of nodes
  std::vector<std::vector<Node *>> nodes(num_layers + 1);
  // Node counter
  unsigned int node_counter = 0;
 
  for (const auto i : make_range(num_layers + 1))
  {
    // calculate number of nodes in each sublayer
    node_number_vec.push_back(
        (unsigned int)(vec_1_node_num + (long)(increment * i + 0.5 - (increment < 0))));
    // Reserve memory for new nodes
    nodes[i] = std::vector<Node *>(node_number_vec[i]);
 
    // Calculate vectors of weighted surrogated index for side #1
    std::vector<Real> weighted_surrogate_index_1;
    std::vector<Real> unweighted_surrogate_index_1;
 
    surrogateGenerator(weighted_surrogate_index_1,
                       unweighted_surrogate_index_1,
                       node_number_vec,
                       wt_1,
                       vec_1_index,
                       vec_1_node_num,
                       i);
 
    // Calculate vectors of weighted surrogated index for side #2
    std::vector<Real> weighted_surrogate_index_2;
    std::vector<Real> unweighted_surrogate_index_2;
 
    surrogateGenerator(weighted_surrogate_index_2,
                       unweighted_surrogate_index_2,
                       node_number_vec,
                       wt_2,
                       vec_2_index,
                       vec_2_node_num,
                       i);
 
    for (const auto j : make_range(node_number_vec[i]))
    {
      // Create surrogate Points on side #1 for Point #j on the sublayer
      Point surrogate_pos_1 = Point(linear_vec_1_x->sample(weighted_surrogate_index_1[j]),
                                    linear_vec_1_y->sample(weighted_surrogate_index_1[j]),
                                    0.0);
      // Create surrogate Points on side #2 for Point #j on the sublayer
      Point surrogate_pos_2 = Point(linear_vec_2_x->sample(weighted_surrogate_index_2[j]),
                                    linear_vec_2_y->sample(weighted_surrogate_index_2[j]),
                                    0.0);
      const Real l_ratio = bias_parameter <= 0.0
                               ? std::pow(spline_vec_2_l->sample(weighted_surrogate_index_2[j]) /
                                              spline_vec_1_l->sample(weighted_surrogate_index_1[j]),
                                          1.0 / ((Real)num_layers - 1.0))
                               : bias_parameter;
      const Real index_factor =
          MooseUtils::absoluteFuzzyEqual(l_ratio, 1.0)
              ? (Real)i / (Real)num_layers
              : (1.0 - std::pow(l_ratio, (Real)i)) / (1.0 - std::pow(l_ratio, (Real)num_layers));
      Point tmp_point = surrogate_pos_2 * index_factor + surrogate_pos_1 * (1.0 - index_factor);
      nodes[i][j] = mesh.add_point(tmp_point, j + node_counter);
    }
    node_counter += node_number_vec[i];
  }
  // Create triangular elements based on the 2D Node vector
  if (quad_elem)
    elementsCreationFromNodesVectorsQuad(mesh,
                                         nodes,
                                         num_layers,
                                         node_number_vec,
                                         transition_layer_id,
                                         input_boundary_1_id,
                                         input_boundary_2_id,
                                         begin_side_boundary_id,
                                         end_side_boundary_id);
  else
    elementsCreationFromNodesVectors(mesh,
                                     nodes,
                                     num_layers,
                                     node_number_vec,
                                     transition_layer_id,
                                     input_boundary_1_id,
                                     input_boundary_2_id,
                                     begin_side_boundary_id,
                                     end_side_boundary_id);
}

Referenced by fillBetweenPointVectorsGenerator(), FillBetweenCurvesGenerator::generate(), FillBetweenPointVectorsGenerator::generate(), and FillBetweenSidesetsGenerator::generate().

isBoundaryOpenSingleSegment()

bool FillBetweenPointVectorsTools::isBoundaryOpenSingleSegment	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	boundary_ordered_node_list,
		const Point	origin_pt,
		const boundary_id_type	bid
	)

Decides whether a boundary of a given mesh is an open single-segment boundary.

Parameters

mesh	input mesh that contains the boundary to be examined
max_node_radius	the maximum radius of the nodes on the boundary
boundary_ordered_node_list	the ordered node ids on the given boundary
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
bid	ID of the boundary to be examined

Returns

whether the boundary is an open single-segment boundary

Definition at line 562 of file FillBetweenPointVectorsTools.C.

{
  try
  {
    isBoundarySimpleClosedLoop(mesh, max_node_radius, boundary_ordered_node_list, origin_pt, bid);
  }
  catch (MooseException & e)
  {
    if (((std::string)e.what())
            .compare("This mesh generator does not work for the provided external boundary as it "
                     "is not a closed loop.") != 0)
      throw MooseException("The provided boundary is not an open single-segment boundary.");
    else
      return true;
  }
 
  throw MooseException("The provided boundary is closed loop, which is not supported.");
}

Referenced by FillBetweenSidesetsGenerator::generate().

isBoundarySimpleClosedLoop() [1/3]

bool FillBetweenPointVectorsTools::isBoundarySimpleClosedLoop	(	MeshBase &	mesh,
		const Point	origin_pt,
		const boundary_id_type	bid
	)

Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the boundary to be examined
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
bid	ID of the boundary to be examined

Returns

whether the boundary is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 555 of file FillBetweenPointVectorsTools.C.

{
  Real dummy_max_node_radius;
  return isBoundarySimpleClosedLoop(mesh, dummy_max_node_radius, origin_pt, bid);
}

isBoundarySimpleClosedLoop() [2/3]

bool FillBetweenPointVectorsTools::isBoundarySimpleClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		const Point	origin_pt,
		const boundary_id_type	bid
	)

Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the boundary to be examined
max_node_radius	the maximum radius of the nodes on the boundary
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
bid	ID of the boundary to be examined

Returns

whether the boundary is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 544 of file FillBetweenPointVectorsTools.C.

{
  std::vector<dof_id_type> dummy_boundary_ordered_node_list;
  return isBoundarySimpleClosedLoop(
      mesh, max_node_radius, dummy_boundary_ordered_node_list, origin_pt, bid);
}

isBoundarySimpleClosedLoop() [3/3]

bool FillBetweenPointVectorsTools::isBoundarySimpleClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	boundary_ordered_node_list,
		const Point	origin_pt,
		const boundary_id_type	bid
	)

Decides whether a boundary of a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the boundary to be examined
max_node_radius	the maximum radius of the nodes on the boundary
boundary_ordered_node_list	the ordered node ids on the given boundary
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
bid	ID of the boundary to be examined

Returns

whether the boundary is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 499 of file FillBetweenPointVectorsTools.C.

{
  // This has no communication and expects elem_ptr to find any
  // element, so it only works on serialized meshes
  libMesh::MeshSerializer serial(mesh);
 
  max_node_radius = 0.0;
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto side_list_tmp = boundary_info.build_side_list();
  std::vector<std::pair<dof_id_type, dof_id_type>> boundary_node_assm;
  std::vector<dof_id_type> boundary_midpoint_node_list;
  for (const auto i : index_range(side_list_tmp))
  {
    if (std::get<2>(side_list_tmp[i]) == bid)
    {
      // store two nodes of each side
      const auto elem = mesh.elem_ptr(std::get<0>(side_list_tmp[i]));
      const auto side = elem->side_ptr(std::get<1>(side_list_tmp[i]));
      boundary_node_assm.push_back(std::make_pair(side->node_id(0), side->node_id(1)));
      // see if there is a midpoint
      const auto & side_type = elem->side_type(std::get<1>(side_list_tmp[i]));
      if (side_type == EDGE3)
        boundary_midpoint_node_list.push_back(
            elem->node_id(elem->n_vertices() + std::get<1>(side_list_tmp[i])));
      else
        boundary_midpoint_node_list.push_back(DofObject::invalid_id);
    }
  }
  bool is_closed_loop;
  isClosedLoop(mesh,
               max_node_radius,
               boundary_ordered_node_list,
               boundary_node_assm,
               boundary_midpoint_node_list,
               origin_pt,
               "external boundary",
               is_closed_loop);
  return is_closed_loop;
}

Referenced by isBoundaryOpenSingleSegment(), isBoundarySimpleClosedLoop(), and isBoundarySimpleClosedLoop().

isClosedLoop() [1/2]

void FillBetweenPointVectorsTools::isClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	ordered_node_list,
		std::vector< std::pair< dof_id_type, dof_id_type > > &	node_assm,
		const Point	origin_pt,
		const std::string	input_type,
		bool &	is_closed_loop,
		const bool	suppress_exception = false
	)

Decides whether a series of nodes contained in a given mesh forms a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the node series to be examined
max_node_radius	the maximum radius of the nodes in the series
ordered_node_list	the ordered node ids on the given curve
node_assm	the assembly of the node id pairs that correspond vertices of a series of sides
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
is_closed_loop	whether the series of nodes form a closed loop with consecutive nodes's azimuthal angles change monotonically
suppress_exception	whether to suppress the exceptions thrown by this function if the boundary is not closed

Definition at line 726 of file FillBetweenPointVectorsTools.C.

{
  std::vector<dof_id_type> dummy_midpoint_node_list(node_assm.size(), DofObject::invalid_id);
  isClosedLoop(mesh,
               max_node_radius,
               ordered_node_list,
               node_assm,
               dummy_midpoint_node_list,
               origin_pt,
               input_type,
               is_closed_loop,
               suppress_exception);
}

isClosedLoop() [2/2]

void FillBetweenPointVectorsTools::isClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	ordered_node_list,
		std::vector< std::pair< dof_id_type, dof_id_type > > &	node_assm,
		std::vector< dof_id_type > &	midpoint_node_list,
		const Point	origin_pt,
		const std::string	input_type,
		bool &	is_closed_loop,
		const bool	suppress_exception = false
	)

Decides whether a series of nodes contained in a given mesh forms a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the node series to be examined
max_node_radius	the maximum radius of the nodes in the series
ordered_node_list	the ordered node ids on the given curve
node_assm	the assembly of the node id pairs that correspond to vertices of the sides on the closed loop
midpoint_node_list	the node ids of the midpoints of the sides for quadratic sides
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)
is_closed_loop	whether the series of nodes form a closed loop with consecutive nodes's azimuthal angles changing monotonically
suppress_exception	whether to suppress the exceptions thrown by this function if the boundary is not closed

Definition at line 664 of file FillBetweenPointVectorsTools.C.

{
  // This has no communication and expects node_ptr to find any
  // node, so it only works on serialized meshes
  libMesh::MeshSerializer serial(mesh);
 
  std::vector<dof_id_type> dummy_elem_list = std::vector<dof_id_type>(node_assm.size(), 0);
  std::vector<dof_id_type> ordered_dummy_elem_list;
  is_closed_loop = false;
  MooseMeshUtils::makeOrderedNodeList(
      node_assm, dummy_elem_list, midpoint_node_list, ordered_node_list, ordered_dummy_elem_list);
  // If the code ever gets here, node_assm is empty.
  // If the ordered_node_list front and back are not the same, the boundary is not a loop.
  // This is not done inside the loop just for some potential applications in the future.
  if (ordered_node_list.front() != ordered_node_list.back())
  {
    // This is invalid type #2
    if (!suppress_exception)
      throw MooseException("This mesh generator does not work for the provided ",
                           input_type,
                           " as it is not a closed loop.");
  }
  // It the curve is a loop, check if azimuthal angles change monotonically
  else
  {
    // Utilize cross product here.
    // If azimuthal angles change monotonically,
    // the z components of the cross products are always negative or positive.
    std::vector<Real> ordered_node_azi_list;
    for (const auto i : make_range(ordered_node_list.size() - 1))
    {
      ordered_node_azi_list.push_back(
          (*mesh.node_ptr(ordered_node_list[i]) - origin_pt)
              .cross(*mesh.node_ptr(ordered_node_list[i + 1]) - origin_pt)(2));
      // Use this opportunity to calculate maximum radius
      max_node_radius =
          std::max((*mesh.node_ptr(ordered_node_list[i]) - origin_pt).norm(), max_node_radius);
    }
    std::sort(ordered_node_azi_list.begin(), ordered_node_azi_list.end());
    if (ordered_node_azi_list.front() * ordered_node_azi_list.back() < 0.0)
    {
      // This is invalid type #3
      if (!suppress_exception)
        throw MooseException(
            "This mesh generator does not work for the provided ",
            input_type,
            " as azimuthal angles of consecutive nodes do not change monotonically.");
    }
    else
      is_closed_loop = true;
  }
}

Referenced by CircularBoundaryCorrectionGenerator::generate(), isBoundarySimpleClosedLoop(), isClosedLoop(), and isCurveSimpleClosedLoop().

isCurveOpenSingleSegment()

bool FillBetweenPointVectorsTools::isCurveOpenSingleSegment	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	ordered_node_list,
		const Point	origin_pt
	)

Decides whether a curve contained in a given mesh is an open single-segment curve.

Parameters

mesh	input mesh that contains the curve to be examined
max_node_radius	the maximum radius of the nodes on the curve
ordered_node_list	the ordered node ids on the given curve
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)

Returns

whether the boundary is an open single-segment boundary

Definition at line 641 of file FillBetweenPointVectorsTools.C.

{
  try
  {
    isCurveSimpleClosedLoop(mesh, max_node_radius, ordered_node_list, origin_pt);
  }
  catch (MooseException & e)
  {
    if (((std::string)e.what())
            .compare("This mesh generator does not work for the provided curve as it is not a "
                     "closed loop.") != 0)
      throw MooseException("The provided curve is not an open single-segment boundary.");
    else
      return true;
  }
  throw MooseException("The provided curve is closed loop, which is not supported.");
  return false;
}

Referenced by FillBetweenCurvesGenerator::generate().

isCurveSimpleClosedLoop() [1/3]

bool FillBetweenPointVectorsTools::isCurveSimpleClosedLoop	(	MeshBase &	mesh,
		const Point	origin_pt
	)

Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the curve to be examined
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)

Returns

whether the curve is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 634 of file FillBetweenPointVectorsTools.C.

{
  Real dummy_max_node_radius;
  return isCurveSimpleClosedLoop(mesh, dummy_max_node_radius, origin_pt);
}

isCurveSimpleClosedLoop() [2/3]

bool FillBetweenPointVectorsTools::isCurveSimpleClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		const Point	origin_pt
	)

Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the curve to be examined
max_node_radius	the maximum radius of the nodes on the curve
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)

Returns

whether the curve is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 627 of file FillBetweenPointVectorsTools.C.

{
  std::vector<dof_id_type> dummy_ordered_node_list;
  return isCurveSimpleClosedLoop(mesh, max_node_radius, dummy_ordered_node_list, origin_pt);
}

isCurveSimpleClosedLoop() [3/3]

bool FillBetweenPointVectorsTools::isCurveSimpleClosedLoop	(	MeshBase &	mesh,
		Real &	max_node_radius,
		std::vector< dof_id_type > &	ordered_node_list,
		const Point	origin_pt
	)

Decides whether a curve contained in a given mesh is a closed loop with consecutive nodes's azimuthal angles change monotonically.

Parameters

mesh	input mesh that contains the curve to be examined
max_node_radius	the maximum radius of the nodes on the curve
ordered_node_list	the ordered node ids on the given curve
origin_pt	origin position of the given mesh (used for azimuthal angle calculation)

Returns

whether the curve is a closed loop with consecutive nodes's azimuthal angles change monotonically

Definition at line 607 of file FillBetweenPointVectorsTools.C.

{
  // This has no communication and expects to loop over all elements
  // on every processor, so it only works on serialized meshes
  libMesh::MeshSerializer serial(mesh);
 
  max_node_radius = 0.0;
  std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
  for (auto it = mesh.active_elements_begin(); it != mesh.active_elements_end(); it++)
    node_assm.push_back(std::make_pair((*it)->node_id(0), (*it)->node_id(1)));
  bool is_closed_loop;
  isClosedLoop(
      mesh, max_node_radius, ordered_node_list, node_assm, origin_pt, "curve", is_closed_loop);
  return is_closed_loop;
}

Referenced by isCurveOpenSingleSegment(), isCurveSimpleClosedLoop(), and isCurveSimpleClosedLoop().

isExternalBoundary()

bool FillBetweenPointVectorsTools::isExternalBoundary	(	MeshBase &	mesh,
		const boundary_id_type	bid
	)

Decides whether a boundary of a given mesh is an external boundary.

Parameters

mesh	input mesh that contains the boundary to be examined
bid	ID of the boundary to be examined

Returns

whether the boundary is the external boundary of the given mesh

Definition at line 586 of file FillBetweenPointVectorsTools.C.

{
  // This has no communication and expects elem_ptr to find any
  // element, so it only works on serialized meshes
  libMesh::MeshSerializer serial(mesh);
 
  if (!mesh.is_prepared())
    mesh.find_neighbors();
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  auto side_list = boundary_info.build_side_list();
  for (const auto i : index_range(side_list))
  {
    if (std::get<2>(side_list[i]) == bid)
      if (mesh.elem_ptr(std::get<0>(side_list[i]))->neighbor_ptr(std::get<1>(side_list[i])) !=
          nullptr)
        return false;
  }
  return true;
}

Referenced by FillBetweenSidesetsGenerator::generate().

needFlip()

bool FillBetweenPointVectorsTools::needFlip	(	const std::vector< Point > &	vec_pts_1,
		const std::vector< Point > &	vec_pts_2
	)

Decide whether one of the input vector of Points needs to be flipped to ensure correct transition layer shape.

Parameters

vec_pts_1	first vector of points to be examined
vec_pts_2	second vector of points to be examined

Returns

whether one of the vectors needs to be flipped to ensure correct transition layer shape

Definition at line 477 of file FillBetweenPointVectorsTools.C.

{
  const Real th1 =
      acos((vec_pts_1.back() - vec_pts_1.front()) * (vec_pts_2.front() - vec_pts_1.front()) /
           (vec_pts_1.back() - vec_pts_1.front()).norm() /
           (vec_pts_2.front() - vec_pts_1.front()).norm());
  const Real th2 = acos(
      (vec_pts_2.back() - vec_pts_1.back()) * (vec_pts_1.front() - vec_pts_1.back()) /
      (vec_pts_2.back() - vec_pts_1.back()).norm() / (vec_pts_1.front() - vec_pts_1.back()).norm());
  const Real th3 = acos(
      (vec_pts_2.front() - vec_pts_2.back()) * (vec_pts_1.back() - vec_pts_2.back()) /
      (vec_pts_2.front() - vec_pts_2.back()).norm() / (vec_pts_1.back() - vec_pts_2.back()).norm());
  const Real th4 =
      acos((vec_pts_1.front() - vec_pts_2.front()) * (vec_pts_2.back() - vec_pts_2.front()) /
           (vec_pts_1.front() - vec_pts_2.front()).norm() /
           (vec_pts_2.back() - vec_pts_2.front()).norm());
  if (MooseUtils::absoluteFuzzyEqual(th1 + th2 + th3 + th4, 2 * M_PI))
    return false;
  return true;
}

Referenced by fillBetweenPointVectorsGenerator().

surrogateGenerator()

void FillBetweenPointVectorsTools::surrogateGenerator	(	std::vector< Real > &	weighted_surrogate_index,
		std::vector< Real > &	unweighted_surrogate_index,
		const std::vector< unsigned int > &	node_number_vec,
		const std::vector< Real > &	index,
		const std::vector< Real > &	wt,
		const unsigned int	boundary_node_num,
		const unsigned int	i
	)

Generates weighted surrogate index vectors on one side for points of a sublayer curve.

Parameters

weighted_surrogate_index	weighted surrogate index vector
unweighted_surrogate_index	unweighted surrogate index vector
node_number_vec	number of points/nodes on the sublayer curve
index	weighted index vector of the points on the side
wt	weight vector of the points on the side
boundary_node_num	number of points/nodes on the side
i	index of the sublayer

Definition at line 436 of file FillBetweenPointVectorsTools.C.

{
  // First element is trivial
  weighted_surrogate_index.push_back(0.0);
  unweighted_surrogate_index.push_back(0.0);
  for (unsigned int j = 1; j < node_number_vec[i]; j++)
  {
    // uniform interval for unweighted index
    unweighted_surrogate_index.push_back((Real)j / ((Real)node_number_vec[i] - 1.0));
    // >
    const auto it_0 =
        std::upper_bound(index.begin(), index.end(), unweighted_surrogate_index[j - 1]);
    // >=
    const auto it_1 = std::lower_bound(index.begin(), index.end(), unweighted_surrogate_index[j]);
    //
    const auto it_dist = std::distance(it_0, it_1);
    //
    const auto it_start = std::distance(index.begin(), it_0);
 
    if (it_0 == it_1)
      weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1] +
                                         wt[it_start - 1] / ((Real)node_number_vec[i] - 1.0));
    else
    {
      weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1]);
      weighted_surrogate_index[j] += (*it_0 - unweighted_surrogate_index[j - 1]) * wt[it_start - 1];
      weighted_surrogate_index[j] +=
          (unweighted_surrogate_index[j] - *(it_1 - 1)) * wt[it_start + it_dist - 1];
      for (unsigned int k = 1; k < it_dist; k++)
        weighted_surrogate_index[j] += wt[it_start + k - 1] / ((Real)boundary_node_num - 1.0);
    }
  }
}

Referenced by fillBetweenPointVectorsGenerator().

weightedInterpolator()

void FillBetweenPointVectorsTools::weightedInterpolator	(	const unsigned int	vec_node_num,
		const std::vector< Point > &	boundary_points_vec,
		std::vector< Real > &	vec_index,
		std::vector< Real > &	wt,
		std::vector< Real > &	index,
		const Real	sigma,
		std::unique_ptr< LinearInterpolation > &	linear_vec_x,
		std::unique_ptr< LinearInterpolation > &	linear_vec_y,
		std::unique_ptr< SplineInterpolation > &	spline_vec_l
	)

Generates weights, weighted indices and corresponding interpolation.

Parameters

vec_node_num	number of points/nodes on the given curve
boundary_points_vec	vector of points on the given curve
vec_index	unweighted index vector
wt	weights based on point-to-point distances
index	weighted index vector
sigma	Gaussian parameter used for Gaussian blurring of local node density
linear_vec_x	linear interpolation of x coordinates based on weighted index
linear_vec_y	linear interpolation of y coordinates based on weighted index
spline_vec_l	spline interpolation of inter-node length based on weighted index

Definition at line 372 of file FillBetweenPointVectorsTools.C.

{
  std::vector<Real> pos_x;
  std::vector<Real> pos_y;
  std::vector<Real> dist_vec;
  std::vector<Real> pos_l;
 
  for (const auto i : make_range(vec_node_num))
  {
    // Unweighted, the index interval is just uniform
    // Normalized range 0~1
    vec_index.push_back((Real)i / ((Real)vec_node_num - 1.0));
    // X and Y coordinates cooresponding to the index
    pos_x.push_back(boundary_points_vec[i](0));
    pos_y.push_back(boundary_points_vec[i](1));
    // Use Point-to-Point distance as unnormalized weight
    if (i > 0)
    {
      wt.push_back((boundary_points_vec[i] - boundary_points_vec[i - 1]).norm());
      dist_vec.push_back(wt.back());
    }
    // Accumulated unnormalized weights to get unnormalized weighted index
    index.push_back(std::accumulate(wt.begin(), wt.end(), 0.0));
  }
  const Real dist_vec_total = index.back(); // Total accumulated distances
  const Real wt_norm_factor = dist_vec_total / ((Real)vec_node_num - 1.0); // Normalization factor
  // Normalization for both weights and weighted indices
  std::transform(
      wt.begin(), wt.end(), wt.begin(), [wt_norm_factor](Real & c) { return c / wt_norm_factor; });
  std::transform(index.begin(),
                 index.end(),
                 index.begin(),
                 [dist_vec_total](Real & c) { return c / dist_vec_total; });
  // Use Gaussian blurring to smoothen local density
  for (const auto i : make_range(vec_node_num))
  {
    Real gaussian_factor(0.0);
    Real sum_tmp(0.0);
    // Use interval as parameter now, consider distance in the future
    for (const auto j : make_range(vec_node_num - 1))
    {
      // dis_vec and index are off by 0.5
      const Real tmp_factor =
          exp(-((Real)(i - j) - 0.5) * ((Real)(i - j) - 0.5) / 2.0 / sigma / sigma);
      gaussian_factor += tmp_factor;
      sum_tmp += tmp_factor * dist_vec[j];
    }
    pos_l.push_back(sum_tmp / gaussian_factor);
  }
  // Interpolate positions based on weighted indices
  linear_vec_x = std::make_unique<LinearInterpolation>(index, pos_x);
  linear_vec_y = std::make_unique<LinearInterpolation>(index, pos_y);
  spline_vec_l = std::make_unique<SplineInterpolation>(index, pos_l);
}

Referenced by fillBetweenPointVectorsGenerator().