In [1]:
from HelmholtzSolverRAD import *
from ExampleBoundaries import SphereRAD
frequency = 40.0 # frequency [Hz]
k = frequencyToWavenumber(frequency)
interiorPoints = np.array([[0.0000, 0.0000],
[0.0000, 0.5000],
[0.0000, -0.5000],
[0.5000, 0.0000]], dtype=np.float32)
# Test Problem 1
# Dirichlet boundary condition with phi = sin(k*z)
#
solver = HelmholtzSolverRAD(*(SphereRAD()))
boundaryCondition = BoundaryCondition(solver.aElement.shape[0])
boundaryCondition.alpha.fill(1.0)
boundaryCondition.beta.fill(0.0)
boundaryCondition.f[:] = np.sin(k * solver.aCenters[:,1])
boundaryIncidence = BoundaryIncidence(solver.aElement.shape[0])
boundaryIncidence.phi.fill(0.0)
boundaryIncidence.v.fill(0.0)
interiorIncidentPhi = np.zeros(interiorPoints.shape[0], dtype=np.complex64)
boundarySolution = solver.solveInteriorBoundary(k, boundaryCondition, boundaryIncidence)
interiorPhi = solver.solveInterior(boundarySolution, interiorIncidentPhi, interiorPoints)
print "Test Problem 1"
print "==============\n"
print boundarySolution
printInteriorSolution(boundarySolution, interiorPhi)
# Test Problem 2
# von Neumann boundary condition such that phi = sin(k/sqrt(2) * x) * sin(k/sqrt(2) * y)
# Differentiate with respect to x and y to obtain outward normal
boundaryCondition.alpha.fill(0.0)
boundaryCondition.beta.fill(1.0)
for i in range(solver.aCenters.shape[0]):
z = solver.aCenters[i, 1]
n = Normal2D(solver.aVertex[solver.aElement[i,0],:],
solver.aVertex[solver.aElement[i,1],:])
boundaryCondition.f[i] = k * np.cos(k*z) * n[1]
boundarySolution = solver.solveInteriorBoundary(k, boundaryCondition, boundaryIncidence)
interiorPhi = solver.solveSamples(boundarySolution, interiorIncidentPhi, interiorPoints, 'interior')
print "\n\nTest Problem 2"
print "==============\n"
print boundarySolution
printInteriorSolution(boundarySolution, interiorPhi)
# Test Problem 3
# Dirichlet boundary condition, such that phi = sin(k/ sqrt(2) * x) * sin(k/sqrt(2) * y)
# Differentiate with respect to x and y to obtain outward normal
boundaryCondition.alpha.fill(1.0)
boundaryCondition.beta.fill(0.0)
zp = 0.25
for i in range(solver.aCenters.shape[0]):
r = solver.aCenters[i, 0]
z = solver.aCenters[i, 1]
# make input complex so proper sqrt is called
rpq = np.sqrt(0.0j + r**2 + (z - zp)**2)
boundaryCondition.f[i] = np.exp(1j * k * rpq) / (4.0 * np.pi * rpq)
boundaryIncidence.phi[i] = np.exp(1j * k * rpq) / (4.0 * np.pi * rpq)
n = Normal2D(solver.aVertex[solver.aElement[i,0],:],
solver.aVertex[solver.aElement[i,1],:])
drbdn = (r * n[0] + (z - zp) * n[1]) / rpq
boundaryIncidence.v[i] = drbdn * np.exp(1j * k * rpq) * (1j * k * rpq - 1.0) \
/ (4.0 * np.pi * rpq*rpq)
for i in range(interiorPoints.shape[0]):
r = interiorPoints[i, 0]
z = interiorPoints[i, 1]
# make input complex so proper sqrt is called
rpq = np.sqrt(0.0j + r**2 + (zp - z)**2)
interiorIncidentPhi[i] = np.exp(1j * k * rpq) / (4.0 * np.pi * rpq)
boundarySolution = solver.solveInteriorBoundary(k, boundaryCondition, boundaryIncidence)
interiorPhi = solver.solveInterior(boundarySolution, interiorIncidentPhi, interiorPoints)
print "\n\nTest Problem 3"
print "==============\n"
print boundarySolution
printInteriorSolution(boundarySolution, interiorPhi)
Copyright (C) 2017 Frank Jargstorff
This file is part of the AcousticBEM library.
AcousticBEM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
AcousticBEM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with AcousticBEM. If not, see http://www.gnu.org/licenses/.