From 0c02f1c1d51a56f10a11251aa543b24447f94623 Mon Sep 17 00:00:00 2001 From: arturgower Date: Wed, 27 Oct 2021 14:35:36 +0100 Subject: [PATCH] implment symmetries and add symmetry to shapes and --- docs/src/library/base.md | 4 +- docs/src/manual/intro.md | 4 +- docs/src/manual/particles.md | 2 +- docs/src/manual/source.md | 10 +-- src/MultipleScattering.jl | 9 +- src/impulse.jl | 2 +- src/physics/acoustics/acoustics.jl | 4 +- src/physics/acoustics/concentric_capsule.jl | 2 +- src/physics/acoustics/source.jl | 33 ++++--- src/plot/plot.jl | 2 +- src/shapes/halfspace.jl | 3 + src/shapes/plate.jl | 3 + src/shapes/shapes.jl | 2 + src/shapes/sphere.jl | 1 + src/simulation.jl | 2 +- src/source.jl | 58 +++++++------ src/types.jl | 95 +++++++++++++++++++++ test/print_tests.jl | 2 +- test/runtests.jl | 3 + test/source_test.jl | 2 +- test/symmetrytests.jl | 57 +++++++++++++ 21 files changed, 242 insertions(+), 58 deletions(-) create mode 100644 src/types.jl create mode 100644 test/symmetrytests.jl diff --git a/docs/src/library/base.md b/docs/src/library/base.md index e7905e35..3a7a96ae 100644 --- a/docs/src/library/base.md +++ b/docs/src/library/base.md @@ -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] diff --git a/docs/src/manual/intro.md b/docs/src/manual/intro.md index 78a18146..32ab5d6d 100644 --- a/docs/src/manual/intro.md +++ b/docs/src/manual/intro.md @@ -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 @@ -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) diff --git a/docs/src/manual/particles.md b/docs/src/manual/particles.md index a7297640..6f1a92bf 100644 --- a/docs/src/manual/particles.md +++ b/docs/src/manual/particles.md @@ -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: diff --git a/docs/src/manual/source.md b/docs/src/manual/source.md index a1af136f..7cc74649 100644 --- a/docs/src/manual/source.md +++ b/docs/src/manual/source.md @@ -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 @@ -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]); @@ -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),$ diff --git a/src/MultipleScattering.jl b/src/MultipleScattering.jl index d2d3f515..2fd65097 100644 --- a/src/MultipleScattering.jl +++ b/src/MultipleScattering.jl @@ -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 @@ -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 @@ -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 diff --git a/src/impulse.jl b/src/impulse.jl index 31e120e6..45da4ab4 100644 --- a/src/impulse.jl +++ b/src/impulse.jl @@ -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) diff --git a/src/physics/acoustics/acoustics.jl b/src/physics/acoustics/acoustics.jl index 35bd537a..b79020f2 100644 --- a/src/physics/acoustics/acoustics.jl +++ b/src/physics/acoustics/acoustics.jl @@ -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 diff --git a/src/physics/acoustics/concentric_capsule.jl b/src/physics/acoustics/concentric_capsule.jl index e18ed3fa..0f3c5e46 100644 --- a/src/physics/acoustics/concentric_capsule.jl +++ b/src/physics/acoustics/concentric_capsule.jl @@ -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 diff --git a/src/physics/acoustics/source.jl b/src/physics/acoustics/source.jl index ce4c259c..83265faf 100644 --- a/src/physics/acoustics/source.jl +++ b/src/physics/acoustics/source.jl @@ -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) @@ -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) @@ -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 @@ -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} @@ -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 @@ -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, diff --git a/src/plot/plot.jl b/src/plot/plot.jl index 25f0132c..a22806fd 100644 --- a/src/plot/plot.jl +++ b/src/plot/plot.jl @@ -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 diff --git a/src/shapes/halfspace.jl b/src/shapes/halfspace.jl index bb8e96cd..2f50144f 100644 --- a/src/shapes/halfspace.jl +++ b/src/shapes/halfspace.jl @@ -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) diff --git a/src/shapes/plate.jl b/src/shapes/plate.jl index 15cd5cb6..7230fc87 100644 --- a/src/shapes/plate.jl +++ b/src/shapes/plate.jl @@ -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) diff --git a/src/shapes/shapes.jl b/src/shapes/shapes.jl index 56594a27..4a1639ba 100644 --- a/src/shapes/shapes.jl +++ b/src/shapes/shapes.jl @@ -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 diff --git a/src/shapes/sphere.jl b/src/shapes/sphere.jl index e5219546..8cc37fc8 100644 --- a/src/shapes/sphere.jl +++ b/src/shapes/sphere.jl @@ -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 diff --git a/src/simulation.jl b/src/simulation.jl index 9e869abf..d199f92b 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -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 diff --git a/src/source.jl b/src/source.jl index 1f4a45fd..51b9e42d 100644 --- a/src/source.jl +++ b/src/source.jl @@ -67,53 +67,56 @@ function PlaneSource(medium::PhysicalMedium{T,Dim,FieldDim}; PlaneSource(medium,direction,position,amplitude) end +Symmetry(s::PlaneSource{T,Dim}) where {T,Dim} = (abs(dot(s.direction,azimuthalnormal(Dim))) == norm(s.direction)) ? PlanarAzimuthalSymmetry{Dim}() : PlanarSymmetry{Dim}() """ - Source(medium::P, field::Function, coef::Function) + RegularSource(medium::P, field::Function, coef::Function) -Is a struct type which describes the source field that drives/forces the whole system. It is also described as an incident wave. It has three fields `Source.medium`, `Source.field`, and `Source.coef`. +Is a struct type which describes the source field that drives/forces the whole system. It is also known as an incident wave. It has three fields `RegularSource.medium`, `RegularSource.field`, and `RegularSource.coef`. The source field at the position 'x' and angular frequency 'ω' is given by ```julia-repl x = [1.0,0.0] ω = 1.0 -Source.field(x,ω) +RegularSource.field(x,ω) ``` -The field `Source.coef` +The field `RegularSource.coef` regular_basis_function(medium::Acoustic{T,2}, ω::T) """ -struct Source{T<:AbstractFloat,P<:PhysicalMedium} <: AbstractSource{T,P} +struct RegularSource{T<:AbstractFloat,P<:PhysicalMedium,S<:AbstractSymmetry} <: AbstractSource{T,P} medium::P "Use: field(x,ω)" field::Function "Use: coefficients(n,x,ω)" coefficients::Function # Enforce that the Types are the same - function Source{T,P}(medium::P,field::Function,coefficients::Function) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}} - s = new{T,P}(medium,field,coefficients) + function RegularSource{T,P,S}(medium::P,field::Function,coefficients::Function) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim},S<:AbstractSymmetry{Dim}} + s = new{T,P,S}(medium,field,coefficients) self_test(s) return s end end -function Source(medium::P,field::Function,coefficients::Function) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}} - return Source{T,P}(medium,field,coefficients) +function RegularSource(medium::P,field::Function,coefficients::Function; symmetry::AbstractSymmetry{Dim} = WithoutSymmetry{Dim}()) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}} + return RegularSource{T,P,typeof(symmetry)}(medium,field,coefficients) end +Symmetry(reg::RegularSource{T,P,S}) where {T,P,S} = S() + """ - regular_spherical_coefficients(source::Source) + regular_spherical_coefficients(source::RegularSource) return a function which can calculate the coefficients of a regular spherical wave basis. """ -function regular_spherical_coefficients(source::Source) +function regular_spherical_coefficients(source::RegularSource) source.coefficients end """ Check that the source functions return the correct types """ -function self_test(source::Source{T,P}) where {P,T} +function self_test(source::RegularSource{T,P}) where {P,T} ω = one(T) # choose rand postion, hopefully not the source position/origin @@ -136,11 +139,11 @@ function self_test(source::Source{T,P}) where {P,T} return true end -field(s::Source) = s.field -field(s::Source{T},x::AbstractArray{T}, ω::T) where T = field(s)(x,ω) +field(s::RegularSource) = s.field +field(s::RegularSource{T},x::AbstractArray{T}, ω::T) where T = field(s)(x,ω) function constant_source(medium::P, num::Complex{T} = zero(Float64) * im) where {P,T} - return Source{T,P}(medium, (x,ω) -> num, (order,x,ω) -> [num]) + return RegularSource{T,P}(medium, (x,ω) -> num, (order,x,ω) -> [num]) end """ @@ -148,7 +151,10 @@ end ``c_n = ```regular_coefficients[n]`, ``x_o=```position`, and let ``v_n(kx)`` represent the regular spherical basis with wavenumber ``k`` at position ``x``. The above returns a `source` where `source.field(x,ω) =```\\sum_n c_n v_n(kx -k x_0)`` """ -function regular_spherical_source(medium::PhysicalMedium{T,Dim},regular_coefficients::AbstractVector{CT}; amplitude::Union{T,Complex{T}} = one(T), position::AbstractArray{T} = SVector(zeros(T,Dim)...)) where {T,Dim,CT<:Union{T,Complex{T}}} +function regular_spherical_source(medium::PhysicalMedium{T,Dim},regular_coefficients::AbstractVector{CT}; + symmetry::AbstractSymmetry{Dim} = WithoutSymmetry{Dim}(), + amplitude::Union{T,Complex{T}} = one(T), + position::AbstractArray{T} = SVector(zeros(T,Dim)...)) where {T,Dim,CT<:Union{T,Complex{T}}} coeff_order = basislength_to_basisorder(typeof(medium),length(regular_coefficients)) @@ -169,11 +175,11 @@ function regular_spherical_source(medium::PhysicalMedium{T,Dim},regular_coeffici return amplitude .* sum(regular_coefficients[n] .* V[n,1:len2] for n in 1:len1) end - return Source(medium, source_field, source_coef) + return RegularSource(medium, source_field, source_coef; symmetry = symmetry) end """ - source_expand(Source, centre; basis_order = 4) + source_expand(RegularSource, centre; basis_order = 4) Returns a function of `(x,ω)` which approximates the value of the source at `(x,ω)`. That is, the source is written in terms of a regular basis expansion centred at `centre`. """ @@ -191,26 +197,30 @@ function source_expand(source::AbstractSource{T}, centre::AbstractVector{T}; bas end import Base.(+) -function +(s1::Source{T,P},s2::Source{T,P})::Source{T,P} where {P,T} +function +(s1::RegularSource{T,P},s2::RegularSource{T,P})::RegularSource{T,P} where {P,T} if s1.medium != s2.medium error("Can not add sources from different physical mediums.") end field(x,ω) = s1.field(x,ω) + s2.field(x,ω) coef(n,centre,ω) = s1.coefficients(n,centre,ω) + s2.coefficients(n,centre,ω) - Source{T,P}(s1.medium,field,coef) + sym1 = Symmetry(s1) + sym2 = Symmetry(s2) + S = typeof(Symmetry(sym1,sym2)) + + RegularSource{T,P,S}(s1.medium,field,coef) end import Base.(*) -function *(a,s::Source{T,P})::Source{T,P} where {P,T} +function *(a,s::RegularSource{T,P,S})::RegularSource{T,P,S} where {P,T,S} if typeof(one(Complex{T})*a) != Complex{T} - error("Multiplying source by $a would cause source type to change, please explicitly cast $a to same type as Source ($T)") + error("Multiplying source by $a would cause source type to change, please explicitly cast $a to same type as RegularSource ($T)") end field(x,ω) = a * s.field(x,ω) coef(n,centre,ω) = a .* s.coefficients(n,centre,ω) - Source{T,P}(s.medium,field,coef) + RegularSource{T,P,S}(s.medium,field,coef) end -*(s::Source,a) = *(a,s) +*(s::RegularSource,a) = *(a,s) diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 00000000..a36cb528 --- /dev/null +++ b/src/types.jl @@ -0,0 +1,95 @@ +""" + AbstractSymmetry + +An abstract types which dictates the symmetry of the setup. That is, the symmetry shared between the incident wave and the shape of the material. +""" +abstract type AbstractSymmetry{Dim} end + +""" + WithoutSymmetry + +A type used to describe materials and incident waves which share no common symmetry. This will lead to the most general regular wave expansion for the eignvectors. +""" +struct WithoutSymmetry{Dim} <: AbstractSymmetry{Dim} end + +""" +An incident plane-wave and halfspace material will result in all fields being plane-waves. +""" +abstract type AbstractPlanarSymmetry{Dim} <: AbstractSymmetry{Dim} end + +""" + PlanarSymmetry + +A type used to describe materials and incident waves that both share a planar symmetry. +""" +struct PlanarSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end + +""" +For spatial dimension > 2, we can consider problems that have azimuthal symmetry. For example, a plane-wave incident on a sphere. +""" +abstract type AbstractAzimuthalSymmetry{Dim} <: AbstractSymmetry{Dim} end + +""" + AzimuthalSymmetry + +A type used to describe materials and incident waves in 3 spatial dimensions that share a symmetry about the azimuth. +""" +struct AzimuthalSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end +AzimuthalSymmetry() = AzimuthalSymmetry{3}() + +""" + RadialSymmetry + +A type used to describe materials and incident waves in that are both radially symmetric. That is, the material is a sphere, and the incident wave is radially symmetric around the center of this spherical material. +""" +struct RadialSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end + +""" + PlanarAzimuthalSymmetry + +A type used to describe a materail and incident wave which have both [`PlanarSymmetry`](@ref) and [`AzimuthalSymmetry`](@ref). +""" +struct PlanarAzimuthalSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end +PlanarAzimuthalSymmetry() = PlanarAzimuthalSymmetry{3}() + +function azimuthalnormal(dim::Int) + if dim == 2 + return [1,0] + elseif dim == 3 + return [0,0,1] + else + return [zeros(dim-1);1] + end +end + +# If T1 is a shape, source, or material, first calculate its symmetry +""" + Symmetry(T) + +returns the symmetry of the object T. Some possible symmetries include [`PlanarSymmetry`](@ref), [`RadialSymmetry`](@ref), and [`WithoutSymmetry`](@ref). +""" + + +""" + Symmetry(T1,T2) + +returns the combined symmetry of the two objects T1 and T2. +""" +Symmetry(T1,T2) = Symmetry(Symmetry(T1),Symmetry(T2)) + +Symmetry(sym1::AbstractSymmetry{Dim},sym2::AbstractSymmetry{Dim}) where Dim = WithoutSymmetry{Dim}() +Symmetry(sym1::S,sym2::S) where S<:AbstractSymmetry = sym1 + +Symmetry(sym1::AbstractPlanarSymmetry{Dim},sym2::AbstractPlanarSymmetry{Dim}) where Dim = PlanarSymmetry{Dim}() +Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = PlanarAzimuthalSymmetry{Dim}() + +Symmetry(sym1::AzimuthalSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() +Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::AzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() + +Symmetry(sym1::AbstractAzimuthalSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() +Symmetry(sym1::RadialSymmetry{Dim},sym2::AbstractAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() + +Symmetry(sym1::RadialSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = RadialSymmetry{Dim}() + +Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() +Symmetry(sym1::RadialSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}() diff --git a/test/print_tests.jl b/test/print_tests.jl index 1da0f167..0d36ef35 100644 --- a/test/print_tests.jl +++ b/test/print_tests.jl @@ -17,7 +17,7 @@ @test circle == show_and_eval(circle) @test particle == show_and_eval(particle) - # Currently can't do this with FrequencySimulations because of Source type + # Currently can't do this with FrequencySimulations because of RegularSource type @test show(devnull, simulation) == nothing end diff --git a/test/runtests.jl b/test/runtests.jl index 42474ad7..914a110f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,9 @@ include("fourier_tests.jl") include("acoustic_physics_tests.jl") include("source_test.jl") + +include("symmetrytests.jl") + include("boundary_conditions_tests.jl") include("end_to_end_tests.jl") diff --git a/test/source_test.jl b/test/source_test.jl index 92d306c0..6ff309a4 100644 --- a/test/source_test.jl +++ b/test/source_test.jl @@ -39,7 +39,7 @@ @test_throws(DimensionMismatch,PlaneSource(a3_host; direction = direction)) @test_throws(DimensionMismatch,PlaneSource(a3_host; amplitude = [1.0,2.0])) - # Test equivalence between plane-Sources + # Test equivalence between plane-sources psource = PlaneSource(a2_host; direction = direction, position = position, amplitude = amplitude) source = plane_source(a2_host; direction = direction, position = position, amplitude = amplitude) diff --git a/test/symmetrytests.jl b/test/symmetrytests.jl new file mode 100644 index 00000000..b59d815b --- /dev/null +++ b/test/symmetrytests.jl @@ -0,0 +1,57 @@ +@testset "Symmetries" begin + + rectangle = Box([[0.0,0.0], [2.0, 3.0]]) + circle = Sphere(rand(2), 2.0) + h2Dazi = Halfspace([1.0,0.0], [0.0,0.0]) + + shapes2D = [rectangle,circle,h2Dazi] + + a2 = Acoustic(0.1,0.1 + 0.0im,2) + + source_position = [0.0,1.0] + s1 = point_source(a2, source_position) + s2 = plane_source(a2; position = source_position, direction = [1.0,0.0]) + + # NOTE: should use the function regular_spherical_source, but for 2D one of the translation matrices has not been implemented yet + s3 = RegularSource{Float64,Acoustic{Float64,2},RadialSymmetry{2}}(a2, (x,ω) -> norm(x)+1.0im, (order,centre,ω) -> [1.0+1.0im for m = -order:order]) + + sources2D = [s1,s2,s3] + + syms = [Symmetry(sh,s) for sh in shapes2D, s in sources2D] + @test syms == AbstractSymmetry{2}[ + WithoutSymmetry{2}() WithoutSymmetry{2}() WithoutSymmetry{2}(); + WithoutSymmetry{2}() AzimuthalSymmetry{2}() RadialSymmetry{2}(); + WithoutSymmetry{2}() PlanarAzimuthalSymmetry{2}() AzimuthalSymmetry{2}() + ] + + check_commute = [Symmetry(sh,s) == Symmetry(s,sh) for s in sources2D, sh in shapes2D] + @test count(check_commute) == length(sources2D) * length(shapes2D) + + + sphere = Sphere(rand(3), 2.0) + p = Plate([1.0,1.0,1.0], 1.0) + hazi = Halfspace([0.0,0.0,1.0]) + + shapes3D = [p,sphere,hazi] + + a3 = Acoustic(3, ρ = 1.3, c = 1.5) + + psource = plane_source(a3; direction = [0.0,0.0,1.0]) + pointsource = point_source(a3, [0.0,0.0,0.0]) + radsource = regular_spherical_source(a3, [1.0+0.0im]; + symmetry = RadialSymmetry{3}() + ); + + sources3D = [pointsource,psource,radsource]; + + syms = [Symmetry(sh,s) for sh in shapes3D, s in sources3D]; + @test syms == AbstractSymmetry{3}[ + WithoutSymmetry{3}() PlanarSymmetry{3}() WithoutSymmetry{3}(); + WithoutSymmetry{3}() AzimuthalSymmetry{3}() RadialSymmetry{3}(); + WithoutSymmetry{3}() PlanarAzimuthalSymmetry{3}() AzimuthalSymmetry{3}() + ] + + check_commute = [Symmetry(sh,s) == Symmetry(s,sh) for s in sources3D, sh in shapes3D] + @test count(check_commute) == length(sources3D) * length(shapes3D) + +end