Skip to content
Snippets Groups Projects
Commit a418ba0e authored by Tizian Wenzel's avatar Tizian Wenzel
Browse files

Added example using Neumann boundary conditions.

parent 383d5343
Branches main
No related tags found
No related merge requests found
File moved
# Example file to demonstrate the use of the VKOGA-PDE code for
# a PDE that uses both Dirichlet and Neumann boundary conditions.
# We use \Omega = [0, 2]^2, L = -LaPlace and a Wendland kernel.
# Some imports
import numpy as np
from vkoga_pde.kernels_PDE import wendland33_laplace
from vkoga_pde.vkoga_PDE import VKOGA_PDE
from matplotlib import pyplot as plt
import math
np.random.seed(1)
## Set number of points that will be sampled
N1 = int(5e3)
N2 = int(1e2)
N3 = int(1e2)
## Create some data to describe domain (X1), Dirichlet (X2) and Neumann (X3) boundary
dim = 2
X1 = np.random.rand(N1, dim)
X2 = np.random.rand(2*N2, dim)
X2[0:N2, 1] = 0
X2[N2:, 1] = 1
X3 = np.random.rand(2*N3, dim)
X3[:N3, 0] = 1
X3[N3:, 0] = 0
n_X3 = np.zeros((2*N3, dim))
n_X3[:N3, 0] = 1
n_X3[:N3, 0] = 0
n_X3[N3:, 0] = -1
n_X3[N3:, 1] = 0
# Rescale the domain a bit: Kernel width parameters are not implemented correctly, thus just scale the domain
scale_factor = 2
X1 = scale_factor * X1
X2 = scale_factor * X2
X3 = scale_factor * X3
# Visualize the domain
plt.figure(1)
plt.clf()
plt.plot(X1[:, 0], X1[:, 1], 'm.')
plt.plot(X2[:, 0], X2[:, 1], 'r.')
plt.plot(X3[:, 0], X3[:, 1], 'b.')
## Define functions u (solution of PDE), f (-Laplace u = f) and g (for Neumann)
u = lambda x: 20 * x[:, [1]] * x[:, [0]] + 1
f = lambda x: np.zeros((x.shape[0], 1))
g = lambda x: 20 * x[:, [1]] * np.sign(x[:, [0]] - .5)
# Define right hand side values
y1 = f(X1) # Application of -Laplace
y2 = u(X2) # Dirichlet boundary, given by solution
y3 = g(X3) # Neumann boundary
## Initialize and run models
kernel = wendland33_laplace(dim=2) # or also try this une
maxIter = 100
model_pde = VKOGA_PDE(kernel=kernel, greedy_type='f_greedy')
_ = model_pde.fit(X1, y1, X2, y2, X3=X3, y3=y3, n_X3=n_X3, maxIter=maxIter)#, y1_sol=u(X1))
## Plot the solution
N_plots = 1000 # don't plot too many points, otherwise plot is slow
fig = plt.figure(105)
ax = fig.add_subplot(projection='3d')
ax.scatter(X1[:N_plots, 0], X1[:N_plots, 1], model_pde.predict_s(X1[:N_plots, :]))
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
#ax.set_zlim(-2,2)
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment