nodes

This module provides functions for generating nodes used for solving PDEs with the RBF and RBF-FD method.

rbf.pde.nodes.disperse(nodes, domain, iterations=20, rho=None, fixed_nodes=None, neighbors=None, delta=0.1)

Disperses the nodes within the domain. The dispersion is analogous to electrostatic repulsion, where neighboring nodes exert a repulsive force on eachother. Each node steps in the direction of its net repulsive force with a step size proportional to the distance to its nearest neighbor. If a node is repelled into a boundary then it bounces back in.

Parameters:

  • nodes ((n, d) float array) – Initial node positions

  • domain ((p, d) float array and (q, d) int array) – Vertices of the domain and connectivity of the vertices.

  • iterations (int, optional) – Number of dispersion iterations.

  • rho (callable, optional) – Takes an (n, d) array as input and returns the repulsion force for a node at those position.

  • fixed_nodes ((k, d) float array, optional) – Nodes which do not move and only provide a repulsion force.

  • neighbors (int, optional) – The number of adjacent nodes used to determine repulsion forces for each node.

  • delta (float, optional) – The step size. Each node moves in the direction of the repulsion force by a distance delta times the distance to the nearest neighbor.

Returns:

(n, d) float array

rbf.pde.nodes.min_energy_nodes(n, domain, rho=None, build_rtree=False, start=0, **kwargs)

Generates nodes within a two or three dimensional. This first generates nodes with a rejection sampling algorithm, and then the nodes are dispersed to ensure a more even distribution.

Parameters:

  • n (int) – The number of nodes generated during rejection sampling. This is not necessarily equal to the number of nodes returned.

  • domain ((p, d) float array and (q, d) int array) – Vertices of the domain and connectivity of the vertices

  • rho (function, optional) – Node density function. Takes a (n, d) array of coordinates and returns an (n,) array of desired node densities at those coordinates. This function should be normalized to be between 0 and 1.

  • build_rtree (bool, optional) – If True, then an R-tree will be built to speed up computational geometry operations. This should be set to True if there are many (>10,000) simplices making up the domain.

  • start (int, optional) – The starting index for the Halton sequence, which is used to propose new points. Setting this value is akin to setting the seed for a random number generator.

  • **kwargs – Additional arguments passed to prepare_nodes

Returns:

  • (n, d) float array – Nodes positions

  • dict – The indices of nodes belonging to each group. There will always be a group called ‘interior’ containing the nodes that are not on the boundary. By default there is a group containing all the boundary nodes called ‘boundary:all’. If boundary_groups was specified, then those groups will be included in this dictionary and their names will be given a ‘boundary:’ prefix. If boundary_groups_with_ghosts was specified then those groups of ghost nodes will be included in this dictionary and their names will be given a ‘ghosts:’ prefix.

  • (n, d) float array – Outward normal vectors for each node. If a node is not on the boundary then its corresponding row will contain NaNs.

Notes

It is assumed that vert and smp define a closed domain. If this is not the case, then it is likely that an error message will be raised which says “ValueError: No intersection found for segment

Examples

make 9 nodes within the unit square

>>> vert = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
>>> smp = np.array([[0, 1], [1, 2], [2, 3], [3, 0]])
>>> out = min_energy_nodes(9, (vert, smp))

view the nodes

>>> out[0]
array([[ 0.50325675,  0.        ],
       [ 0.00605261,  1.        ],
       [ 1.        ,  0.51585247],
       [ 0.        ,  0.00956821],
       [ 1.        ,  0.99597894],
       [ 0.        ,  0.5026365 ],
       [ 1.        ,  0.00951112],
       [ 0.48867638,  1.        ],
       [ 0.54063894,  0.47960892]])

view the indices of nodes making each group

>>> out[1]
{'boundary:all': array([7, 6, 5, 4, 3, 2, 1, 0]),
 'interior': array([8])}

view the outward normal vectors for each node, note that the normal vector for the interior node is nan

>>> out[2]
array([[  0.,  -1.],
       [  0.,   1.],
       [  1.,  -0.],
       [ -1.,  -0.],
       [  1.,  -0.],
       [ -1.,  -0.],
       [  1.,  -0.],
       [  0.,   1.],
       [ nan,  nan]])
rbf.pde.nodes.neighbor_argsort(nodes, m=None)

Returns a permutation array that sorts nodes so that each node and its m-1 nearest neighbors are close together in memory. This is done through the use of a KD Tree and the Reverse Cuthill-McKee algorithm.

Parameters:

  • nodes ((n, d) float array) –

  • m (int, optional) –

Returns:

(N,) int array

Examples

>>> nodes = np.array([[0.0, 1.0],
                      [2.0, 1.0],
                      [1.0, 1.0]])
>>> idx = neighbor_argsort(nodes, 2)
>>> nodes[idx]
array([[ 2., 1.],
       [ 1., 1.],
       [ 0., 1.]])
rbf.pde.nodes.poisson_disc_nodes(radius, domain, ntests=50, rmax_factor=1.5, build_rtree=False, **kwargs)

Generates nodes within a two or three dimensional domain. This first generate nodes with Poisson disc sampling, and then the nodes are dispersed to ensure a more even distribution. This function is considerably slower than min_energy_nodes but it has the advantage of directly specifying the node spacing.

Parameters:

  • radius (float or callable) – The radius for each disc. This is the minimum allowable distance between the nodes generated by Poisson disc sampling. This can be a float or a function that takes a (n, d) array of locations and returns an (n,) array of disc radii.

  • domain ((p, d) float array and (q, d) int array) – Vertices of the domain and connectivity of the vertices

  • build_rtree (bool, optional) – If True, then an R-tree will be built to speed up computational geometry operations. This should be set to True if there are many (>10,000) simplices making up the domain

  • **kwargs – Additional arguments passed to prepare_nodes

Returns:

  • (n, d) float array – Nodes positions

  • dict – The indices of nodes belonging to each group. There will always be a group called ‘interior’ containing the nodes that are not on the boundary. By default there is a group containing all the boundary nodes called ‘boundary:all’. If boundary_groups was specified, then those groups will be included in this dictionary and their names will be given a ‘boundary:’ prefix. If boundary_groups_with_ghosts was specified then those groups of ghost nodes will be included in this dictionary and their names will be given a ‘ghosts:’ prefix.

  • (n, d) float array – Outward normal vectors for each node. If a node is not on the boundary then its corresponding row will contain NaNs.

Notes

It is assumed that vert and smp define a closed domain. If this is not the case, then it is likely that an error message will be raised which says “ValueError: No intersection found for segment

rbf.pde.nodes.prepare_nodes(nodes, domain, rho=None, iterations=20, neighbors=None, dispersion_delta=0.1, pinned_nodes=None, snap_delta=0.5, boundary_groups=None, boundary_groups_with_ghosts=None, ghost_delta=0.5, include_vertices=False, orient_simplices=True)

Prepares a set of nodes for solving PDEs with the RBF and RBF-FD method. This includes: dispersing the nodes away from eachother to ensure a more even spacing, snapping nodes to the boundary, determining the normal vectors for each node, determining the group that each node belongs to, creating ghost nodes, sorting the nodes so that adjacent nodes are close in memory, and verifying that no two nodes are anomalously close to eachother.

The function returns a set of nodes, the normal vectors for each node, and a dictionary identifying which group each node belongs to.

Parameters:

  • nodes ((n, d) float arrary) – An initial sampling of nodes within the domain

  • domain ((p, d) float array and (q, d) int array) – Vertices of the domain and connectivity of the vertices

  • rho (function, optional) – Node density function. Takes a (n, d) array of coordinates and returns an (n,) array of desired node densities at those coordinates. This is used during the node dispersion step.

  • iterations (int, optional) – Number of dispersion iterations.

  • neighbors (int, optional) – Number of neighboring nodes to use when calculating the repulsion force. This defaults to 3 for 2D nodes and 4 for 3D nodes.

  • dispersion_delta (float, optional) – Scaling factor for the node step size in each iteration. The step size is equal to dispersion_delta times the distance to the nearest neighbor.

  • pinned_nodes ((k, d) array, optional) – Nodes which do not move and only provide a repulsion force. These nodes are included in the set of nodes returned by this function and they are in the group named “pinned”.

  • snap_delta (float, optional) – Controls the maximum snapping distance. The maximum snapping distance for each node is snap_delta times the distance to the nearest neighbor. This defaults to 0.5.

  • boundary_groups (dict, optional) – Dictionary defining the boundary groups. The keys are the names of the groups and the values are lists of simplex indices making up each group. This function will return a dictionary identifying which nodes belong to each boundary group. By default, there is a single group named ‘all’ for the entire boundary. Specifically, The default value is {‘all’:range(len(smp))}.

  • boundary_groups_with_ghosts (list of strs, optional) – List of boundary groups that will be given ghost nodes. By default, no boundary groups are given ghost nodes. The groups specified here must exist in boundary_groups.

  • ghost_delta (float, optional) – How far the ghost nodes should be from their corresponding boundary node. The distance is ghost_delta times the distance to the nearest neighbor.

  • include_vertices (bool, optional) – If True, then the vertices will be included in the output nodes. Each vertex will be assigned to the boundary group that its adjoining simplices are part of. If the simplices are in multiple groups, then the vertex will be assigned to the group containing the simplex that comes first in smp.

  • orient_simplices (bool, optional) – If False then it is assumed that the simplices are already oriented such that their normal vectors point outward.

Returns:

  • (m, d) float array – Nodes positions

  • dict – The indices of nodes belonging to each group. There will always be a group called ‘interior’ containing the nodes that are not on the boundary. By default there is a group containing all the boundary nodes called ‘boundary:all’. If boundary_groups was specified, then those groups will be included in this dictionary and their names will be given a ‘boundary:’ prefix. If boundary_groups_with_ghosts was specified then those groups of ghost nodes will be included in this dictionary and their names will be given a ‘ghosts:’ prefix.

  • (n, d) float array – Outward normal vectors for each node. If a node is not on the boundary then its corresponding row will contain NaNs.

Examples

''' 
This script demonstrates how to use `min_energy_nodes` to generate
nodes over a domain with spatially variable density. These nodes can
then be used to solve partial differential equation for this domain
using the spectral RBF method or the RBF-FD method.
'''
import numpy as np
import matplotlib.pyplot as plt
from rbf.pde.nodes import min_energy_nodes

# Define the problem domain with line segments.
vert = np.array([[2.0, 1.0], [1.0, 1.0], [1.0, 2.0], 
                 [0.0, 2.0], [0.0, 0.0], [2.0, 0.0]])
smp = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]])
# define which simplices belong to which boundary groups
boundary_groups = {'inside_corner': [0, 1],
                   'outside_corner': [2, 3, 4, 5]}
# define which boundary groups get ghost nodes
boundary_groups_with_ghosts = ['outside_corner']
# total number of nodes 
N = 500
# define a node density function. It takes an (N, D) array of positions
# and returns an (N,) array of normalized densities between 0 and 1
def rho(x):
  r = np.sqrt((x[:, 0] - 1.0)**2 + (x[:, 1] - 1.0)**2)
  return 0.2 + 0.8/((r/0.3)**4 + 1.0)

nodes, groups, normals = min_energy_nodes(
    N, 
    (vert, smp),
    rho=rho,
    boundary_groups=boundary_groups,
    boundary_groups_with_ghosts=boundary_groups_with_ghosts,
    include_vertices=True)

# plot the results
fig, ax = plt.subplots(figsize=(6, 6))
# plot the domain
for s in smp: 
  ax.plot(vert[s, 0], vert[s, 1], 'k-')

# plot the different node groups and their normal vectors
for i, (name, idx) in enumerate(groups.items()):
  ax.plot(nodes[idx, 0], nodes[idx, 1], 'C%s.' % i, label=name, ms=8)
  ax.quiver(nodes[idx, 0], nodes[idx, 1], 
            normals[idx, 0], normals[idx, 1], 
            color='C%s' % i)

ax.legend()
ax.set_aspect('equal')
fig.tight_layout()
plt.savefig('../figures/nodes.a.png')
plt.show()
_images/nodes.a.png 
'''
This script demonstrates generating nodes with `min_energy_nodes` and
it verifies that the density of the returned nodes accurately matches
the specified density.
'''
import logging
import numpy as np
import matplotlib.pyplot as plt

from rbf.basis import se
from rbf.pde.nodes import min_energy_nodes

logging.basicConfig(level=logging.DEBUG)


def desired_rho(x):
  '''
  desired node density
  '''
  r1 = np.sqrt((x[:, 0] - 0.25)**2 + (x[:, 1] - 0.25)**2)
  bump1 = 0.76 / ((r1/0.1)**2 + 1.0)

  r2 = np.sqrt((x[:, 0] - 0.25)**2 + (x[:, 1] - 0.75)**2)
  bump2 = 0.76 / ((r2/0.1)**4 + 1.0)

  r3 = np.sqrt((x[:, 0] - 0.75)**2 + (x[:, 1] - 0.75)**2)
  bump3 = 0.76 / ((r3/0.1)**8 + 1.0)

  r4 = np.sqrt((x[:, 0] - 0.75)**2 + (x[:, 1] - 0.25)**2)
  bump4 = 0.76 / ((r4/0.1)**16 + 1.0)

  out = 0.2 + bump1 + bump2 + bump3 + bump4
  return out


def actual_rho(x, nodes):
  '''
  compute the density of `nodes` and evaluate the density function at
  `x`. The output is normalize 1.0
  '''
  out = np.zeros(x.shape[0])
  for n in nodes:
    out += se(x, n[None, :], eps=0.01)[:, 0]
    
  out /= np.max(out)
  return out

vert = np.array([[0.0, 0.0],
                 [1.0, 0.0],
                 [1.0, 1.0],
                 [0.0, 1.0]])
smp = np.array([[0, 1], [1, 2], [2, 3], [3, 0]])                 

nodes = min_energy_nodes(10000, (vert, smp), rho=desired_rho)[0]

# plot the nodes
fig, ax = plt.subplots()
for s in smp:
  ax.plot(vert[s, 0], vert[s, 1], 'k-')

ax.plot(nodes[:, 0], nodes[:, 1], 'k.', ms=1)  
ax.set_aspect('equal')
ax.set_title('node positions')
fig.tight_layout()
plt.savefig('../figures/nodes.b.1.png')

fig, axs = plt.subplots(1, 2, figsize=(10, 4))

# evaluate and plot the node density
x = np.linspace(-0.1, 1.1, 250)
y = np.linspace(-0.1, 1.1, 250)
xg, yg = np.meshgrid(x, y)
x = xg.flatten()
y = yg.flatten()
points = np.array([x, y]).T

p = axs[0].tripcolor(points[:, 0], points[:, 1], desired_rho(points),
                     cmap='viridis', vmin=0.0, vmax=1.0)
for s in smp:
  axs[0].plot(vert[s, 0], vert[s, 1], 'w-')

axs[0].set_xlim(-0.1, 1.1)
axs[0].set_ylim(-0.1, 1.1)
axs[0].set_aspect('equal')
axs[0].set_title('desired node density')
fig.colorbar(p, ax=axs[0])

p = axs[1].tripcolor(points[:, 0], points[:, 1], actual_rho(points, nodes),
                     cmap='viridis', vmin=0.0, vmax=1.0)
for s in smp:
  axs[1].plot(vert[s, 0], vert[s, 1], 'w-')

axs[1].set_xlim(-0.1, 1.1)
axs[1].set_ylim(-0.1, 1.1)
axs[1].set_aspect('equal')
axs[1].set_title('actual node density')
fig.colorbar(p, ax=axs[1])
fig.tight_layout()
plt.savefig('../figures/nodes.b.2.png')
plt.show()
_images/nodes.b.1.png _images/nodes.b.2.png