Skip to content

Commit

Permalink
Build the RBF matrix from old points.
Browse files Browse the repository at this point in the history
  • Loading branch information
pvc1989 committed Oct 29, 2024
1 parent d80bc62 commit 0e115fe
Showing 1 changed file with 63 additions and 2 deletions.
65 changes: 63 additions & 2 deletions programming/mesh/cgns/shift_by_rbf.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from scipy import sparse
from scipy.spatial import KDTree
import sys
import argparse


X, Y, Z = 0, 1, 2
Expand Down Expand Up @@ -40,6 +41,8 @@ def build_dok_matrix(kdtree: KDTree, radius: float, verbose: bool) -> sparse.dok
if verbose and n_point <= 100:
print(f'i = {i_point:2d}, j = {j_point:2d}, d_ij = {distance_ij:.2e}')
dok_matrix[i_point, j_point] = distance_ij
if verbose:
print(f'point[{i_point}] has {len(j_neighbors)} neighbors')
if verbose and n_point <= 100:
print('points =\n', points)
print('matrix =\n', dok_matrix.todense())
Expand All @@ -49,7 +52,11 @@ def build_dok_matrix(kdtree: KDTree, radius: float, verbose: bool) -> sparse.dok
def test_on_random_points(n_point: int, radius: float, verbose: bool):
np.random.seed(123456789)
points = np.random.rand(n_point, 3)
if verbose:
print('building the KDTree ...')
kdtree = KDTree(points)
if verbose:
print('building the dok_matrix ...')
dok_matrix = build_dok_matrix(kdtree, radius, verbose)
sparsity = dok_matrix.nnz / (n_point * n_point)
print(f'nnz / (n * n) = {dok_matrix.nnz} / {n_point * n_point} = {sparsity:.2e}')
Expand All @@ -69,6 +76,60 @@ def test_on_random_points(n_point: int, radius: float, verbose: bool):
print(np.allclose(den_matrix.dot(u_coeffs), u_bc))


def build_rbf_matrix(old_points: np.ndarray, radius: float, verbose: bool) -> sparse.csc_matrix:
assert old_points.shape[1] == 3
print('shape =', old_points.shape)
n_point = len(old_points)
if verbose:
print('building the KDTree ...')
print(old_points)
kdtree = KDTree(old_points)
if verbose:
print('building the dok_matrix ...')
dok_matrix = build_dok_matrix(kdtree, radius, verbose)
sparsity = dok_matrix.nnz / (n_point * n_point)
print(f'nnz / (n * n) = {dok_matrix.nnz} / {n_point * n_point} = {sparsity:.2e}')
csc_matrix = dok_matrix.tocsc()
csc_file_name = f'csc_matrix_n={n_point}_r={radius:.2e}.npz'
sparse.save_npz(csc_file_name, csc_matrix)
return csc_matrix


def solve_rbf_system(rbf_matrix: sparse.csc_matrix, rhs_columns: np.ndarray, radius: float, verbose: bool) -> np.ndarray:
n_point = len(rhs_columns)
n_column = rhs_columns.shape[1]
sol_columns = np.ndarray((n_point + 4, n_column))
for i_column in range(n_column):
column_i = np.zeros((n_point + 4,))
column_i[0 : n_point] = rhs_columns[:, i_column]
if verbose:
print(f'solving column[{i_column}] ...')
sol_columns[:, i_column], exit_code = sparse.linalg.lgmres(
rbf_matrix, column_i, maxiter=20000, atol=0.0, rtol=1e-7)
print(f'column[{i_column}] solved with exit_code = {exit_code}')
sol_file_name = f'sol_columns_n={n_point}_r={radius:.2e}.npy'
np.save(sol_file_name, sol_columns)


if __name__ == '__main__':
n_point = int(sys.argv[1])
test_on_random_points(n_point, radius=0.1, verbose=True)
parser = argparse.ArgumentParser(
prog = f'python3 {sys.argv[0][:-3]}.py',
description='Shift interior points by RBF interpolation of boundary points.')
parser.add_argument('--mesh', type=str, help='the CGNS file of the mesh to be shifted')
parser.add_argument('--old_points', type=str, help='the NPY file of boundary points before shifting')
parser.add_argument('--new_points', type=str, help='the NPY file of boundary points after shifting')
parser.add_argument('--radius', type=float, help='the radius of the RBF basis')
parser.add_argument('--output', type=str, help='the output mesh file')
parser.add_argument('--verbose', default=True, action='store_true')
args = parser.parse_args()

old_points = np.load(args.old_points)
new_points = np.load(args.new_points)
assert old_points.shape == new_points.shape

# rbf_matrix = build_rbf_matrix(old_points, args.radius, args.verbose)
rbf_matrix = sparse.load_npz('./csc_matrix_n=270666_r=2.00e+00.npz')
rbf_coeffs = solve_rbf_system(rbf_matrix, new_points - old_points, args.radius, args.verbose)

# n_point = int(sys.argv[1])
# test_on_random_points(n_point, radius=0.1, verbose=True)

0 comments on commit 0e115fe

Please sign in to comment.