Skip to content

Commit

Permalink
implment symmetries and add symmetry to shapes and
Browse files Browse the repository at this point in the history
  • Loading branch information
arturgower committed Oct 27, 2021
1 parent c470e4b commit 0c02f1c
Show file tree
Hide file tree
Showing 21 changed files with 242 additions and 58 deletions.
4 changes: 2 additions & 2 deletions docs/src/library/base.md
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,9 @@ CapsuleParticle
iscongruent(::AbstractParticle,::AbstractParticle)
```

## [Source](@id source_base)
## [RegularSource](@id source_base)

Source types and functions for all physical mediums.
RegularSource types and functions for all physical mediums.

```@autodocs
Modules = [MultipleScattering]
Expand Down
4 changes: 2 additions & 2 deletions docs/src/manual/intro.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ Acoustic(1.0, 1.0 + 0.0im, 2)
```
At this step we have restricted the physics to acoustics, that is, solutions to the Helmholtz equation: $\nabla^2 u(x,y,\omega) + k^2 u(x,y,\omega) = 0$, where $k = \omega/c$, $\omega$ is the angular frequency and $c$ the sound speed of the medium.

### [Source wave](@id into_source)
### [RegularSource wave](@id into_source)

The host medium will determine the types of waves that can propagate. For example an incident plane wave $\mathrm e^{ \mathrm i k x}$ there is a convenient constructor
```jldoctest intro
Expand All @@ -29,7 +29,7 @@ julia> source = plane_source(host_medium; direction = [1.0,0.0]);
Often $\mathrm e^{ \mathrm i k x - \mathrm i \omega t}$ is considered to be a harmonic plane-wave travelling along the $x-$axis. We omit the part $ - \mathrm i \omega t$ as is common in frequency space.


We generally call the incident wave a source. See [Sources](@ref) for details, and see [Acoustic](@ref) for some functions for the `Acoustic` medium.
We generally call the incident wave a source. See [RegularSources](@ref) for details, and see [Acoustic](@ref) for some functions for the `Acoustic` medium.

### [Particles](@id into_particles)

Expand Down
2 changes: 1 addition & 1 deletion docs/src/manual/particles.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ end
```
[`Particle`](@ref) is a `struct` which define a scatterer, obstacle, or simply particle, which can scatter waves. See [Particle](@ref base_particle) for a list of relevant types and functions.

A `Particle` have two fields: `medium` and `shape`. The `medium` defines what kind of waves can propagate inside `Particle`, and what type of `Source` can be used to scatter waves from `Particle`. One example `medium` is [Acoustic](@ref). The `shape` completely defines the geometry and position of the particle, see [Shapes](@ref) for details.
A `Particle` have two fields: `medium` and `shape`. The `medium` defines what kind of waves can propagate inside `Particle`, and what type of `RegularSource` can be used to scatter waves from `Particle`. One example `medium` is [Acoustic](@ref). The `shape` completely defines the geometry and position of the particle, see [Shapes](@ref) for details.

For an example we can define a circular particle with acoustics medium:

Expand Down
10 changes: 5 additions & 5 deletions docs/src/manual/source.md
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
# Sources
# RegularSources

```@meta
DocTestSetup = quote
using MultipleScattering
end
```
[`Source`](@ref) is a `struct` which represents any source, also called an incident wave. See [Source](@ref source_base) for a list of relevant types and functions.
[`RegularSource`](@ref) is a `struct` which represents any source, also called an incident wave. See [RegularSource](@ref source_base) for a list of relevant types and functions.

## Acoustics

Expand Down Expand Up @@ -86,9 +86,9 @@ where `field_apply` is applied to the wave field at every point, the default is

To define a new source you will need to understand the internals below.

## Source internals
## RegularSource internals

The `struc` [`Source`](@ref) has three fields: `medium`, `field`, and `coef`, explained with examples below:
The `struc` [`RegularSource`](@ref) has three fields: `medium`, `field`, and `coef`, explained with examples below:
```jldoctest intro
julia> plane_wave = plane_source(Acoustic(1.0, 1.0, 2); direction = [1.0,0.0]);
Expand All @@ -101,7 +101,7 @@ julia> plane_wave.field(x,ω) # the value of the field
0.5403023058681398 + 0.8414709848078965im
```

To calculate the scattering from a particle due to a source, we need the coefficients $a_n(\mathbf x_0, \omega) = $ Source.coefficients(n,x0,ω). We use these coefficients to represent the source in a radial coordinate system. That is, for any origin $\mathbf x_0$, we need to represent the incident wave $u_{\text{in}}(\mathbf x)$ using a series or regular waves
To calculate the scattering from a particle due to a source, we need the coefficients $a_n(\mathbf x_0, \omega) = $ RegularSource.coefficients(n,x0,ω). We use these coefficients to represent the source in a radial coordinate system. That is, for any origin $\mathbf x_0$, we need to represent the incident wave $u_{\text{in}}(\mathbf x)$ using a series or regular waves

$u_{\text{in}}(\mathbf x) = \sum_n a_n(\mathbf x_0, \omega) \mathrm v_n(\mathbf x - \mathbf x_0),$

Expand Down
9 changes: 7 additions & 2 deletions src/MultipleScattering.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,10 @@ __precompile__()

module MultipleScattering

## Symmetries
export Symmetry, AbstractSymmetry, AbstractPlanarSymmetry, AbstractAzimuthalSymmetry
export WithoutSymmetry, PlanarSymmetry, PlanarAzimuthalSymmetry, AzimuthalSymmetry, RadialSymmetry

## Shapes
export Shape, Circle, Rectangle, Box, EmptyShape, Halfspace, Plate, TimeOfFlight, TimeOfFlightFromPoint

Expand All @@ -22,8 +26,8 @@ export Acoustic, AcousticCapsule, sound_hard, hard, rigid, zero_neumann, sound_s
## Particles
export AbstractParticle, Particle, CapsuleParticle, AbstractParticles

## Sources
export AbstractSource, Source, source_expand, regular_spherical_coefficients, self_test, constant_source, point_source, plane_source, regular_spherical_source, (*), (+)
## RegularSources
export AbstractSource, RegularSource, source_expand, regular_spherical_coefficients, self_test, constant_source, point_source, plane_source, regular_spherical_source, (*), (+)
export PlaneSource

## Main simulation and results
Expand Down Expand Up @@ -54,6 +58,7 @@ using ProgressMeter


# Generic machinery common to all physical models
include("types.jl")
include("shapes/shapes.jl")

include("physics/special_functions.jl") # Special functions missing from Base library
Expand Down
2 changes: 1 addition & 1 deletion src/impulse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ import Base.(*)
function *(a,impulse::ContinuousImpulse{T})::ContinuousImpulse{T} where {T}

if typeof(one(Complex{T})*a) != Complex{T}
throw(DomainError(a, "Multiplying impulse by $a would cause impulse type to change, please explicitly cast $a to same type as Source ($T)"))
throw(DomainError(a, "Multiplying impulse by $a would cause impulse type to change, please explicitly cast $a to same type as RegularSource ($T)"))
end

in_time(t) = a*impulse.in_time(t)
Expand Down
4 changes: 2 additions & 2 deletions src/physics/acoustics/acoustics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -256,11 +256,11 @@ See [`sound_soft`](@ref).
zero_dirichlet(host_medium::Acoustic{T,Dim}) where {T,Dim} = sound_soft(T, Dim)

"""
internal_field(x::AbstractVector, p::Particle{T,Dim,Acoustic{T,Dim}}, source::Source, ω::T, scattering_coefficients::AbstractVector{Complex{T}})
internal_field(x::AbstractVector, p::Particle{T,Dim,Acoustic{T,Dim}}, source::RegularSource, ω::T, scattering_coefficients::AbstractVector{Complex{T}})
The internal field for an acoustic particle in an acoustic medium. For a sphere and circular cylinder the result is exact, for everything else it is an approximation which assumes smooth fields.
"""
function internal_field(x::AbstractVector{T}, p::Particle{T,Dim,Acoustic{T,Dim}}, source::Source{T,Acoustic{T,Dim}}, ω::T, scattering_coefficients::AbstractVector{Complex{T}}) where {T,Dim}
function internal_field(x::AbstractVector{T}, p::Particle{T,Dim,Acoustic{T,Dim}}, source::RegularSource{T,Acoustic{T,Dim}}, ω::T, scattering_coefficients::AbstractVector{Complex{T}}) where {T,Dim}
if !(x p)
@error "Point $x is not inside the particle with shape $(p.shape)"
end
Expand Down
2 changes: 1 addition & 1 deletion src/physics/acoustics/concentric_capsule.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ function t_matrix(cap::CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, medium::A
return Diagonal(vcat(reverse(Tns), Tns[2:end]))
end

function internal_field(x::AbstractArray{T}, p::CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, source::Source{T,Acoustic{T,2}}, ω::T, scattering_coefficients::AbstractVector) where T
function internal_field(x::AbstractArray{T}, p::CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, source::RegularSource{T,Acoustic{T,2}}, ω::T, scattering_coefficients::AbstractVector) where T

Nh = Int((length(scattering_coefficients) - one(T))/T(2.0)) #shorthand

Expand Down
33 changes: 19 additions & 14 deletions src/physics/acoustics/source.jl
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
"""
point_source(medium::Acoustic, source_position, amplitude=1)::Source
point_source(medium::Acoustic, source_position, amplitude=1)::RegularSource
Create 2D [`Acoustic`](@ref) point [`Source`](@ref) (zeroth Hankel function of first type)
Create 2D [`Acoustic`](@ref) point [`RegularSource`](@ref) (zeroth Hankel function of first type)
"""
function point_source(medium::Acoustic{T,2}, source_position::AbstractVector, amplitude::Union{T,Complex{T},Function} = one(T))::Source{T,Acoustic{T,2}} where T <: AbstractFloat
function point_source(medium::Acoustic{T,2}, source_position::AbstractVector, amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,2}} where T <: AbstractFloat

# Convert to SVector for efficiency and consistency
source_position = SVector{2,T}(source_position)
Expand All @@ -23,11 +23,11 @@ function point_source(medium::Acoustic{T,2}, source_position::AbstractVector, am
return (amp(ω)*im)/4 * [hankelh1(-n,k*r) * exp(-im*n*θ) for n = -order:order]
end

return Source{T,Acoustic{T,2}}(medium, source_field, source_coef)
return RegularSource{T,Acoustic{T,2},WithoutSymmetry{2}}(medium, source_field, source_coef)
end

# If we replaced 3 with Dim below this could should work for all dimensions! Test carefully after changing.
function point_source(medium::Acoustic{T,3}, source_position, amplitude::Union{T,Complex{T},Function} = one(T))::Source{T,Acoustic{T,3}} where T <: AbstractFloat
function point_source(medium::Acoustic{T,3}, source_position, amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,3}} where T <: AbstractFloat

# Convert to SVector for efficiency and consistency
source_position = SVector{3,T}(source_position)
Expand Down Expand Up @@ -56,31 +56,33 @@ function point_source(medium::Acoustic{T,3}, source_position, amplitude::Union{T
return amp(ω) * U[1,:]
end

return Source{T,Acoustic{T,3}}(medium, source_field, source_coef)
return RegularSource{T,Acoustic{T,3},WithoutSymmetry{3}}(medium, source_field, source_coef)
end


function plane_source(medium::Acoustic{T,Dim}; position::AbstractArray{T} = SVector(zeros(T,Dim)...),
direction = SVector(one(T), zeros(T,Dim-1)...),
amplitude::Union{T,Complex{T},Function} = one(T))::Source{T,Acoustic{T,Dim}} where {T, Dim}
amplitude::Union{T,Complex{T},Function} = one(T))::RegularSource{T,Acoustic{T,Dim}} where {T, Dim}

plane_source(medium, position, direction, amplitude)
end

"""
plane_source(medium::Acoustic, source_position, source_direction=[1,0], amplitude=1)::Source
plane_source(medium::Acoustic, source_position, source_direction=[1,0], amplitude=1)::RegularSource
Create an [`Acoustic`](@ref) planar wave [`Source`](@ref)
Create an [`Acoustic`](@ref) planar wave [`RegularSource`](@ref)
"""
function plane_source(medium::Acoustic{T,2}, position::AbstractArray{T}, direction::AbstractArray{T} = SVector(one(T),zero(T)), amplitude::Union{T,Complex{T}} = one(T))::Source{T,Acoustic{T,2}} where {T}
function plane_source(medium::Acoustic{T,2}, position::AbstractArray{T}, direction::AbstractArray{T} = SVector(one(T),zero(T)), amplitude::Union{T,Complex{T}} = one(T))::RegularSource{T,Acoustic{T,2}} where {T}

# Convert to SVector for efficiency and consistency
position = SVector(position...)
direction = SVector((direction ./ norm(direction))...) # unit direction

S = (abs(dot(direction,azimuthalnormal(2))) == one(T)) ? PlanarAzimuthalSymmetry{2} : PlanarSymmetry{2}

# This pseudo-constructor is rarely called, so do some checks and conversions
if iszero(norm(direction))
throw(DomainError("Source direction must not have zero magnitude."))
throw(DomainError("RegularSource direction must not have zero magnitude."))
end

if typeof(amplitude) <: Number
Expand All @@ -97,7 +99,7 @@ function plane_source(medium::Acoustic{T,2}, position::AbstractArray{T}, directi
source_field(centre,ω) * [exp(im * n *(T(pi)/2 - θ)) for n = -order:order]
end

return Source{T,Acoustic{T,2}}(medium, source_field, source_coef)
return RegularSource{T,Acoustic{T,2},S}(medium, source_field, source_coef)
end

function plane_source(medium::Acoustic{T,3}, position::AbstractArray{T}, direction::AbstractArray{T} = SVector(zero(T),zero(T),one(T)), amplitude::Union{T,Complex{T}} = one(T)) where {T}
Expand All @@ -106,9 +108,11 @@ function plane_source(medium::Acoustic{T,3}, position::AbstractArray{T}, directi
position = SVector(position...)
direction = SVector( (direction ./ norm(direction))...) # unit direction

S = (abs(dot(direction,azimuthalnormal(3))) == one(T)) ? PlanarAzimuthalSymmetry{3} : PlanarSymmetry{3}

# This pseudo-constructor is rarely called, so do some checks and conversions
if iszero(norm(direction))
throw(DomainError("Source direction must not have zero magnitude."))
throw(DomainError("RegularSource direction must not have zero magnitude."))
end

if typeof(amplitude) <: Number
Expand All @@ -131,10 +135,11 @@ function plane_source(medium::Acoustic{T,3}, position::AbstractArray{T}, directi
for l = 0:order for m = -l:l]
end

return Source{T,Acoustic{T,3}}(medium, source_field, source_coef)
return RegularSource{T,Acoustic{T,3},S}(medium, source_field, source_coef)
end

function regular_spherical_coefficients(psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}) where {T,Dim}

source = plane_source(psource.medium;
amplitude = psource.amplitude[1],
position = psource.position,
Expand Down
2 changes: 1 addition & 1 deletion src/plot/plot.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,6 @@ end
end

"Plot the source field for a particular wavenumber"
@recipe function plot(s::Source, ω)
@recipe function plot(s::RegularSource, ω)
(FrequencySimulation(s), ω)
end
3 changes: 3 additions & 0 deletions src/shapes/halfspace.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,9 @@ Halfspace(normal::AbstractVector{T}, origin::AbstractVector{T} = zeros(T,length(
Halfspace(normal::NTuple{Dim,T}, origin::AbstractVector{T} = zeros(T,Dim)) where {T,Dim} = Halfspace{T,Dim}(normal ./ norm(normal), origin)

name(shape::Halfspace) = "Halfspace"
function Symmetry(shape::Halfspace{T,Dim}) where {T,Dim}
return (abs(dot(shape.normal,azimuthalnormal(Dim))) == one(T)) ? PlanarAzimuthalSymmetry{Dim}() : PlanarSymmetry{Dim}()
end

volume(shape::Halfspace) = Inf
outer_radius(hs::Halfspace{T}) where T = T(Inf)
Expand Down
3 changes: 3 additions & 0 deletions src/shapes/plate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,9 @@ Plate(normal::AbstractVector{T}, width::T, origin::AbstractVector{T} = zeros(T,l
Plate(normal::NTuple{Dim,T}, width::T, origin::AbstractVector{T} = zeros(T,Dim)) where {T,Dim} = Plate{T,Dim}(normal ./ norm(normal), width, origin)

name(shape::Plate) = "Plate"
function Symmetry(shape::Plate{T,Dim}) where {T,Dim}
return (abs(dot(shape.normal,azimuthalnormal(Dim))) == one(T)) ? PlanarAzimuthalSymmetry{Dim}() : PlanarSymmetry{Dim}()
end

volume(shape::Plate) = Inf
outer_radius(hs::Plate{T}) where T = T(Inf)
Expand Down
2 changes: 2 additions & 0 deletions src/shapes/shapes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@ Abstract idea which defines the external boundary of object.
"""
abstract type Shape{T<:AbstractFloat,Dim} end

Symmetry(s::Shape{T,Dim}) where {T,Dim} = WithoutSymmetry{Dim}()

"""
origin(shape::Shape)::SVector
Expand Down
1 change: 1 addition & 0 deletions src/shapes/sphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ Sphere(radius::T) where {T} = Sphere{T,3}(zeros(T,3), radius)

name(shape::Sphere) = "Sphere"
name(shape::Sphere{T,2}) where T = "Circle"
Symmetry(shape::Sphere{T,Dim}) where {T,Dim} = RadialSymmetry{Dim}()

outer_radius(sphere::Sphere) = sphere.radius
volume(shape::Sphere{T,3}) where T = 4//3 * π * shape.radius^3
Expand Down
2 changes: 1 addition & 1 deletion src/simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ After building, you can [`run`](@ref) the simulation to get a [`FrequencySimulat
mutable struct FrequencySimulation{T<:AbstractFloat,Dim,P<:PhysicalMedium} <: Simulation{T,Dim}
"Vector of particles, can be of different types."
particles::AbstractParticles
"Source wave, where source.medium is the background medium of the simulation."
"RegularSource wave, where source.medium is the background medium of the simulation."
source::AbstractSource{T,P}
end

Expand Down
Loading

0 comments on commit 0c02f1c

Please sign in to comment.