From 34c0632efe8dd2e489f2c5df5b0a3d2243a8ab88 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 11 Oct 2024 17:32:55 +0200 Subject: [PATCH 1/7] rename `LagrangeMesh` --> `Mesh` --- ext/IntiGmshExt.jl | 8 ++--- src/gmsh_api.jl | 2 +- src/mesh.jl | 75 ++++++++++++++++++++++------------------------ 3 files changed, 41 insertions(+), 44 deletions(-) diff --git a/ext/IntiGmshExt.jl b/ext/IntiGmshExt.jl index 388147eb..7dc13462 100644 --- a/ext/IntiGmshExt.jl +++ b/ext/IntiGmshExt.jl @@ -14,7 +14,7 @@ function Inti.import_mesh(filename = nothing; dim = 3) if isnothing(filename) gmsh.isInitialized() == 1 || error("gmsh is not initialized. Try gmsh.initialize() first.") - msh = Inti.LagrangeMesh{3,Float64}() + msh = Inti.Mesh{3,Float64}() _import_mesh!(msh) dim == 2 && (msh = Inti._convert_to_2d(msh)) # create iterators for the lagrange elements @@ -41,8 +41,8 @@ end """ _import_mesh!(msh) -Import the mesh from gmsh into `msh` as a [`LagrangeMesh`](@ref -Inti.LagrangeMesh). +Import the mesh from gmsh into `msh` as a [`Mesh`](@ref +Inti.Mesh). """ function _import_mesh!(msh) # NOTE: when importing the nodes, we will renumber them so that the first node has @@ -176,7 +176,7 @@ function _etype_to_type_tag(E::DataType) return gmsh.model.mesh.getElementType(family_name, order) end -function Inti.write_gmsh_model(msh::Inti.LagrangeMesh{N,Float64}; name = "") where {N} +function Inti.write_gmsh_model(msh::Inti.Mesh{N,Float64}; name = "") where {N} @assert N ∈ (2, 3) # lift the nodes to 3d if N == 2 nodes = N == 3 ? Inti.nodes(msh) : [SVector(x[1], x[2], 0) for x in Inti.nodes(msh)] diff --git a/src/gmsh_api.jl b/src/gmsh_api.jl index 7a61fe4b..6ddd0c29 100644 --- a/src/gmsh_api.jl +++ b/src/gmsh_api.jl @@ -1,7 +1,7 @@ """ import_mesh(filename = nothing; dim=3) -Open `filename` and create a [`LagrangeMesh`](@ref) from the `gmsh` model in it. +Open `filename` and create a [`Mesh`](@ref) from the `gmsh` model in it. If `filename` is `nothing`, the current `gmsh` model is used. Note that this assumes that the *Gmsh* API has been initialized through `gmsh.initialize`. diff --git a/src/mesh.jl b/src/mesh.jl index 04176dc5..df6b5bce 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -7,7 +7,7 @@ An abstract mesh structure in dimension `N` with primite data of type `T` (e.g. Concrete subtypes of `AbstractMesh` should implement [`ElementIterator`](@ref) for accessing the mesh elements. -See also: [`LagrangeMesh`](@ref) +See also: [`Mesh`](@ref) """ abstract type AbstractMesh{N,T} end @@ -58,7 +58,7 @@ Return the element types present in the `msh`. function element_types end """ - struct LagrangeMesh{N,T} <: AbstractMesh{N,T} + struct Mesh{N,T} <: AbstractMesh{N,T} Unstructured mesh is defined by a set of `nodes`` (of type `SVector{N,T}`), and a dictionary mapping element types to connectivity matrices. Each columns of a @@ -72,7 +72,7 @@ a given mesh using e.g. `view(msh,Γ)` or `msh[Γ]`, where `Γ` is a See [`elements`](@ref) for a way to iterate over the elements of a mesh. """ -struct LagrangeMesh{N,T} <: AbstractMesh{N,T} +struct Mesh{N,T} <: AbstractMesh{N,T} nodes::Vector{SVector{N,T}} # for each element type (key), return the connectivity matrix etype2mat::Dict{DataType,Matrix{Int}} @@ -83,8 +83,8 @@ struct LagrangeMesh{N,T} <: AbstractMesh{N,T} end # empty constructor -function LagrangeMesh{N,T}() where {N,T} - return LagrangeMesh{N,T}( +function Mesh{N,T}() where {N,T} + return Mesh{N,T}( SVector{N,T}[], Dict{DataType,Matrix{Int}}(), Dict{DataType,AbstractVector{<:ReferenceInterpolant}}(), @@ -92,16 +92,16 @@ function LagrangeMesh{N,T}() where {N,T} ) end -element_types(msh::LagrangeMesh) = keys(msh.etype2mat) +element_types(msh::Mesh) = keys(msh.etype2mat) -elements(msh::LagrangeMesh, E::DataType) = msh.etype2els[E] +elements(msh::Mesh, E::DataType) = msh.etype2els[E] # generic implementation function elements(msh::AbstractMesh) return Iterators.flatten(elements(msh, E) for E in element_types(msh)) end -nodes(msh::LagrangeMesh) = msh.nodes +nodes(msh::Mesh) = msh.nodes """ ent2etags(msh::AbstractMesh) @@ -109,9 +109,9 @@ nodes(msh::LagrangeMesh) = msh.nodes Return a dictionary mapping entities to a dictionary of element types to element tags. """ -ent2etags(msh::LagrangeMesh) = msh.ent2etags -etype2mat(msh::LagrangeMesh) = msh.etype2mat -etype2els(msh::LagrangeMesh) = msh.etype2els +ent2etags(msh::Mesh) = msh.ent2etags +etype2mat(msh::Mesh) = msh.etype2mat +etype2els(msh::Mesh) = msh.etype2els """ connectivity(msh::AbstractMesh,E::DataType) @@ -119,9 +119,9 @@ etype2els(msh::LagrangeMesh) = msh.etype2els Return the connectivity matrix for elements of type `E` in `msh`. The integer tags in the matrix refer to the points in `nodes(msh)` """ -connectivity(msh::LagrangeMesh, E::DataType) = msh.etype2mat[E] +connectivity(msh::Mesh, E::DataType) = msh.etype2mat[E] -function ent2nodetags(msh::LagrangeMesh, ent::EntityKey) +function ent2nodetags(msh::Mesh, ent::EntityKey) tags = Int[] for (E, t) in msh.ent2etags[ent] append!(tags, view(msh.etype2mat[E], :, t)) @@ -129,7 +129,7 @@ function ent2nodetags(msh::LagrangeMesh, ent::EntityKey) return tags end -entities(msh::LagrangeMesh) = keys(msh.ent2etags) +entities(msh::Mesh) = keys(msh.ent2etags) """ domain(msh::AbstractMesh) @@ -139,7 +139,7 @@ Return a [`Domain`] containing of all entities covered by the mesh. domain(msh::AbstractMesh) = Domain(entities(msh)) """ - dom2elt(m::LagrangeMesh,Ω,E)::Vector{Int} + dom2elt(m::Mesh,Ω,E)::Vector{Int} Compute the element indices `idxs` of the elements of type `E` composing `Ω`. """ @@ -152,8 +152,8 @@ function dom2elt(m::AbstractMesh, Ω::Domain, E::DataType) return idxs end -function Base.getindex(msh::LagrangeMesh{N,T}, Ω::Domain) where {N,T} - new_msh = LagrangeMesh{N,T}() +function Base.getindex(msh::Mesh{N,T}, Ω::Domain) where {N,T} + new_msh = Mesh{N,T}() (; nodes, etype2mat, etype2els, ent2etags) = new_msh foreach(k -> ent2etags[k] = Dict{DataType,Vector{Int}}(), keys(Ω)) glob2loc = Dict{Int,Int}() @@ -192,14 +192,14 @@ function Base.getindex(msh::LagrangeMesh{N,T}, Ω::Domain) where {N,T} end return new_msh end -Base.getindex(msh::LagrangeMesh, ent::EntityKey) = getindex(msh, Domain(ent)) +Base.getindex(msh::Mesh, ent::EntityKey) = getindex(msh, Domain(ent)) """ meshgen(Ω, n) meshgen(Ω, n_dict) meshgen(Ω; meshsize) -Generate a `LagrangeMesh` for the domain `Ω` where each curve is meshed using +Generate a `Mesh` for the domain `Ω` where each curve is meshed using `n` elements. Passing a dictionary allows for a finer control; in such cases, `n_dict[ent]` should return an integer for each entity `ent` in `Ω` of `geometric_dimension` one. @@ -222,7 +222,7 @@ function meshgen(Ω::Domain, args...; kwargs...) # 3d). Only makes sense if all entities have the same ambient dimension. N = ambient_dimension(first(Ω)) @assert all(p -> ambient_dimension(p) == N, entities(Ω)) "Entities must have the same ambient dimension" - mesh = LagrangeMesh{N,Float64}() + mesh = Mesh{N,Float64}() meshgen!(mesh, Ω, args...; kwargs...) return mesh end @@ -234,11 +234,11 @@ meshgen(e::EntityKey, args...; kwargs...) = meshgen(Domain(e), args...; kwargs.. Similar to [`meshgen`](@ref), but append entries to `mesh`. """ -function meshgen!(msh::LagrangeMesh, Ω::Domain, num_elements::Int) +function meshgen!(msh::Mesh, Ω::Domain, num_elements::Int) e1d = filter(k -> geometric_dimension(k) == 1, all_keys(Ω)) return meshgen!(msh, Ω, Dict(e => num_elements for e in e1d)) end -function meshgen!(msh::LagrangeMesh, Ω::Domain; meshsize::Real) +function meshgen!(msh::Mesh, Ω::Domain; meshsize::Real) # compute the length of each curve using an adaptive quadrature e1d = filter(k -> geometric_dimension(k) == 1, all_keys(Ω)) dict = Dict(k => ceil(Int, measure(k) / meshsize) for k in e1d) @@ -262,7 +262,7 @@ function meshgen!(msh::LagrangeMesh, Ω::Domain; meshsize::Real) end return meshgen!(msh, Ω, dict) end -function meshgen!(msh::LagrangeMesh, Ω::Domain, dict::Dict) +function meshgen!(msh::Mesh, Ω::Domain, dict::Dict) d = geometric_dimension(Ω) for k in keys(Ω) if d == 1 @@ -306,7 +306,7 @@ end # a simple mesh generation for GeometricEntity objects containing a # parametrization -function _meshgen!(mesh::LagrangeMesh, key::EntityKey, sz) +function _meshgen!(mesh::Mesh, key::EntityKey, sz) ent = global_get_entity(key) d = domain(ent) hasparametrization(ent) || error("$key has no parametrization") @@ -346,7 +346,7 @@ function _meshgen(f, d::HyperRectangle, sz::NTuple) end |> vec end -function _build_connectivity!(msh::LagrangeMesh{N,T}, tol = 1e-8) where {N,T} +function _build_connectivity!(msh::Mesh{N,T}, tol = 1e-8) where {N,T} nodes = msh.nodes connect_dict = msh.etype2mat # first build a naive connectivity matrix where duplicate points are present @@ -387,17 +387,14 @@ ElementIterator(m, E::DataType) = ElementIterator{E,typeof(m)}(m) # implement the interface for ElementIterator of lagrange elements on a generic # mesh. The elements are constructed on the flight based on the global nodes and # the connectivity list stored -function Base.length(iter::ElementIterator{E,<:LagrangeMesh}) where {E} +function Base.length(iter::ElementIterator{E,<:Mesh}) where {E} tags = iter.mesh.etype2mat[E]::Matrix{Int} _, Nel = size(tags) return Nel end -Base.size(iter::ElementIterator{E,<:LagrangeMesh}) where {E} = (length(iter),) +Base.size(iter::ElementIterator{E,<:Mesh}) where {E} = (length(iter),) -function Base.getindex( - iter::ElementIterator{E,<:LagrangeMesh}, - i::Int, -) where {E<:LagrangeElement} +function Base.getindex(iter::ElementIterator{E,<:Mesh}, i::Int) where {E<:LagrangeElement} tags = iter.mesh.etype2mat[E]::Matrix{Int} node_tags = view(tags, :, i) vtx = view(iter.mesh.nodes, node_tags) @@ -409,8 +406,8 @@ end # converting various element types to their 2d counterpart. These are needed # because some meshers like gmsh always create three-dimensional objects, so we # must convert after importing the mesh -function _convert_to_2d(mesh::LagrangeMesh{3,T}) where {T} - msh2d = LagrangeMesh{2,T}() +function _convert_to_2d(mesh::Mesh{3,T}) where {T} + msh2d = Mesh{2,T}() # create new dictionaries for elements and ent2etagsdict with 2d elements as keys new_etype2mat = msh2d.etype2mat new_ent2etags = msh2d.ent2etags @@ -438,7 +435,7 @@ function _convert_to_2d(::Type{LagrangeElement{R,N,SVector{3,T}}}) where {R,N,T} end _convert_to_2d(::Type{SVector{3,T}}) where {T} = SVector{2,T} -function remove_duplicate_nodes!(msh::LagrangeMesh, tol) +function remove_duplicate_nodes!(msh::Mesh, tol) nodes = copy(msh.nodes) new_nodes = empty!(msh.nodes) connect_dict = msh.etype2mat @@ -474,12 +471,12 @@ over elements of the submesh just like you would with a mesh. Construct `SubMesh`s using `view(parent,Ω::Domain)`. """ struct SubMesh{N,T} <: AbstractMesh{N,T} - parent::LagrangeMesh{N,T} + parent::Mesh{N,T} domain::Domain # etype2etags maps E => indices of elements in parent mesh contained in the # submesh etype2etags::Dict{DataType,Vector{Int}} - function SubMesh(mesh::LagrangeMesh{N,T}, Ω::Domain) where {N,T} + function SubMesh(mesh::Mesh{N,T}, Ω::Domain) where {N,T} etype2etags = Dict{DataType,Vector{Int}}() for E in element_types(mesh) # add the indices of the elements of type E in the submesh. Skip if @@ -491,8 +488,8 @@ struct SubMesh{N,T} <: AbstractMesh{N,T} end end -Base.view(m::LagrangeMesh, Ω::Domain) = SubMesh(m, Ω) -Base.view(m::LagrangeMesh, ent::EntityKey) = SubMesh(m, Domain(ent)) +Base.view(m::Mesh, Ω::Domain) = SubMesh(m, Ω) +Base.view(m::Mesh, ent::EntityKey) = SubMesh(m, Domain(ent)) Base.collect(msh::SubMesh) = msh.parent[msh.domain] @@ -546,7 +543,7 @@ function nodetags(msh::SubMesh) end return unique!(tags) end -nodetags(msh::LagrangeMesh) = collect(1:length(msh.nodes)) +nodetags(msh::Mesh) = collect(1:length(msh.nodes)) """ nodes(msh::SubMesh) From 48c56373d3f9171c0aadcc8b1b7d57ebfd8c9085 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 11 Oct 2024 17:38:26 +0200 Subject: [PATCH 2/7] rename `AbstractPDE` --> `AbstractDifferentialOperator` --- docs/src/pluto-examples/toy_example.md | 19 +++++++++ docs/src/tutorials/integral_operators.md | 4 +- src/bdim.jl | 2 +- src/kernels.jl | 54 ++++++++++++------------ src/vdim.jl | 8 +++- 5 files changed, 55 insertions(+), 32 deletions(-) create mode 100644 docs/src/pluto-examples/toy_example.md diff --git a/docs/src/pluto-examples/toy_example.md b/docs/src/pluto-examples/toy_example.md new file mode 100644 index 00000000..2ddca8fc --- /dev/null +++ b/docs/src/pluto-examples/toy_example.md @@ -0,0 +1,19 @@ + + +```julia +begin + import Pkg as _Pkg + haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"]) + using PlutoUI: with_terminal +end ; +``` + +# Toy example + +[![Pluto notebook](https://img.shields.io/badge/download-Pluto_notebook-blue)](../../build/pluto_examples/toy_example.jl) $\hspace{0.2cm}$ [![nbviewer](https://img.shields.io/badge/show-nbviewer-blue.svg)](../../build/pluto_examples/toy_example.html) + +All examples in Inti.jl are autogenerated by executing the `make.jl` script in the `docs` folder. The workflow uses [PlutoStaticHTML.jl](https://plutostatichtml.huijzer.xyz/stable/) and [PlutoSliderServer.jl](https://github.com/JuliaPluto/PlutoSliderServer.jl) via the intermediate of [ExampleJuggler.jl](https://github.com/j-fu/ExampleJuggler.jl) to generate (i) markdown files passed to Documenter.jl, and (ii) pluto notebook files downloadable from the example's page. + +```julia +using Inti +``` \ No newline at end of file diff --git a/docs/src/tutorials/integral_operators.md b/docs/src/tutorials/integral_operators.md index 217ecca9..44545679 100644 --- a/docs/src/tutorials/integral_operators.md +++ b/docs/src/tutorials/integral_operators.md @@ -20,7 +20,7 @@ handle custom kernels. ## Predefined kernels and integral operators To simplify the construction of integral operators for some commonly used PDEs, -Inti.jl defines a few [`AbstractPDE`](@ref)s types. For each of these PDEs, the +Inti.jl defines a few [`AbstractDifferentialOperator`](@ref)s types. For each of these PDEs, the package provides a [`SingleLayerKernel`](@ref), [`DoubleLayerKernel`](@ref), [`HyperSingularKernel`](@ref), and [`AdjointDoubleLayerKernel`](@ref) that can be used to construct the corresponding kernel functions, e.g.: @@ -173,7 +173,7 @@ up a custom kernel function, and how to build an integral operator from it. !!! note "Integral operators coming from PDEs" If the integral operator of interest arises from a PDE, it is recommended - to define a new [`AbstractPDE`](@ref) type, and implement the required + to define a new [`AbstractDifferentialOperator`](@ref) type, and implement the required methods for [`SingleLayerKernel`](@ref), [`DoubleLayerKernel`](@ref), [`AdjointDoubleLayerKernel`](@ref), and [`HyperSingularKernel`](@ref). This will enable the use of the high-level syntax for constructing boundary diff --git a/src/bdim.jl b/src/bdim.jl index f00b7f9d..edffa331 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -24,7 +24,7 @@ See [faria2021general](@cite) for more details on the method. ## Required: -- `pde` must be an [`AbstractPDE`](@ref) +- `pde` must be an [`AbstractDifferentialOperator`](@ref) - `Y` must be a [`Quadrature`](@ref) object of a closed surface - `X` is either inside, outside, or on `Y` - `S` and `D` are approximations to the single- and double-layer operators for diff --git a/src/kernels.jl b/src/kernels.jl index 747dc1c9..c80d421c 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -23,63 +23,63 @@ struct GenericKernel{T,F} <: AbstractKernel{T} end """ - abstract type AbstractPDE{N} + abstract type AbstractDifferentialOperator{N} -A partial differential equation in dimension `N`. `AbstractPDE` types are used -to define [`AbstractPDEKernel`s](@ref AbstractPDEKernel). +A partial differential equation in dimension `N`. `AbstractDifferentialOperator` types are used +to define [`AbstractDifferentialOperatorKernel`s](@ref AbstractDifferentialOperatorKernel). """ -abstract type AbstractPDE{N} end +abstract type AbstractDifferentialOperator{N} end -ambient_dimension(::AbstractPDE{N}) where {N} = N +ambient_dimension(::AbstractDifferentialOperator{N}) where {N} = N """ - abstract type AbstractPDEKernel{T,Op} <: AbstractKernel{T} + abstract type AbstractDifferentialOperatorKernel{T,Op} <: AbstractKernel{T} An [`AbstractKernel`](@ref) with an associated `pde::Op` field. """ -abstract type AbstractPDEKernel{T,Op} <: AbstractKernel{T} end +abstract type AbstractDifferentialOperatorKernel{T,Op} <: AbstractKernel{T} end """ - pde(K::AbstractPDEKernel) + pde(K::AbstractDifferentialOperatorKernel) -Return the underlying `AbstractPDE` associated with the kernel `K`. +Return the underlying `AbstractDifferentialOperator` associated with the kernel `K`. """ -pde(k::AbstractPDEKernel) = k.pde +pde(k::AbstractDifferentialOperatorKernel) = k.pde -parameters(k::AbstractPDEKernel) = parameters(pde(k)) +parameters(k::AbstractDifferentialOperatorKernel) = parameters(pde(k)) # convenient constructor for e.g. SingleLayerKernel(pde,Float64) or DoubleLayerKernel(pde,ComplexF64) function (::Type{K})( pde::Op, ::Type{T} = default_kernel_eltype(pde), -) where {T,Op,K<:AbstractPDEKernel} +) where {T,Op,K<:AbstractDifferentialOperatorKernel} return K{T,Op}(pde) end """ - struct SingleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} + struct SingleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} The free-space single-layer kernel (i.e. the fundamental solution) of an `Op <: -AbstractPDE`. +AbstractDifferentialOperator`. """ -struct SingleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} +struct SingleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} pde::Op end """ - struct DoubleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} + struct DoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} Given an operator `Op`, construct its free-space double-layer kernel. This corresponds to the `γ₁` trace of the [`SingleLayerKernel`](@ref). For operators such as [`Laplace`](@ref) or [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the source variable. """ -struct DoubleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} +struct DoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} pde::Op end """ - struct AdjointDoubleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} + struct AdjointDoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} Given an operator `Op`, construct its free-space adjoint double-layer kernel. This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the @@ -87,12 +87,12 @@ This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the target variable. """ -struct AdjointDoubleLayerKernel{T,Op} <: AbstractPDEKernel{T,Op} +struct AdjointDoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} pde::Op end """ - struct HyperSingularKernel{T,Op} <: AbstractPDEKernel{T,Op} + struct HyperSingularKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} Given an operator `Op`, construct its free-space hypersingular kernel. This corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the @@ -100,7 +100,7 @@ corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative respect to the target variable of the `DoubleLayerKernel`. """ -struct HyperSingularKernel{T,Op} <: AbstractPDEKernel{T,Op} +struct HyperSingularKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} pde::Op end @@ -108,7 +108,7 @@ end ################################# LAPLACE ###################################### ################################################################################ -struct Laplace{N} <: AbstractPDE{N} end +struct Laplace{N} <: AbstractDifferentialOperator{N} end """ Laplace(; dim) @@ -194,7 +194,7 @@ end ################################# Yukawa ####################################### ################################################################################ -struct Yukawa{N,K<:Real} <: AbstractPDE{N} +struct Yukawa{N,K<:Real} <: AbstractDifferentialOperator{N} λ::K end @@ -292,7 +292,7 @@ end ################################# Helmholtz #################################### ################################################################################ -struct Helmholtz{N,K} <: AbstractPDE{N} +struct Helmholtz{N,K} <: AbstractDifferentialOperator{N} k::K end @@ -406,7 +406,7 @@ function (HS::HyperSingularKernel{T,S})(target, source)::T where {T,S<:Helmholtz end ############################ STOKES ############################3 -struct Stokes{N,T} <: AbstractPDE{N} +struct Stokes{N,T} <: AbstractDifferentialOperator{N} μ::T end Stokes(; μ, dim = 3) = Stokes{dim}(μ) @@ -476,12 +476,12 @@ end ################################################################################ """ - struct Elastostatic{N,T} <: AbstractPDE{N} + struct Elastostatic{N,T} <: AbstractDifferentialOperator{N} Elastostatic equation in `N` dimensions: μΔu + (μ+λ)∇(∇⋅u) = 0. Note that the displacement u is a vector of length `N` since this is a vectorial problem. """ -struct Elastostatic{N,T} <: AbstractPDE{N} +struct Elastostatic{N,T} <: AbstractDifferentialOperator{N} μ::T λ::T end diff --git a/src/vdim.jl b/src/vdim.jl index 8e1c7354..142019c8 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -5,7 +5,7 @@ Compute a correction to the volume potential `V : Y → X` such that `V + δV` i more accurate approximation of the underlying volume potential operator. The correction is computed using the (volume) density interpolation method. -This function requires a `pde::AbstractPDE`, a target set `X`, a source +This function requires a `pde::AbstractDifferentialOperator`, a target set `X`, a source quadrature `Y`, a boundary quadrature `Y_boundary`, approximations `S : Y_boundary -> X` and `D : Y_boundary -> X` to the single- and double-layer potentials (correctly handling nearly-singular integrals), and a naive @@ -304,7 +304,11 @@ trace of `Pₙ`. Passing a point `center` will shift the monomials and solutions accordingly. """ -function polynomial_solutions_vdim(pde::AbstractPDE, order::Integer, center = nothing) +function polynomial_solutions_vdim( + pde::AbstractDifferentialOperator, + order::Integer, + center = nothing, +) N = ambient_dimension(pde) center = isnothing(center) ? zero(SVector{N,Float64}) : center # create empty arrays to store the monomials, solutions, and traces. For the From 56678c4e15ef7692896ed663c820df38f6077e2a Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 11 Oct 2024 18:12:03 +0200 Subject: [PATCH 3/7] fix typos in various docstrings - minus signs missing in def. of various opeartors - simplify naming --- src/kernels.jl | 146 ++++++++++++++++++++++++------------------------- 1 file changed, 71 insertions(+), 75 deletions(-) diff --git a/src/kernels.jl b/src/kernels.jl index c80d421c..d06fe82d 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -1,4 +1,4 @@ -const PREDEFINED_KERNELS = ["Laplace", "Helmholtz", "Stokes", "Yukawa"] +const PREDEFINED_OPERATORS = ["Laplace", "Helmholtz", "Stokes", "Yukawa"] """ abstract type AbstractKernel{T} @@ -13,73 +13,50 @@ abstract type AbstractKernel{T} end return_type(::AbstractKernel{T}) where {T} = T -""" - struct GenericKernel{T,F} <: AbstractKernel{T} - -An [`AbstractKernel`](@ref) with `kernel` of type `F`. -""" -struct GenericKernel{T,F} <: AbstractKernel{T} - kernel::F -end - """ abstract type AbstractDifferentialOperator{N} -A partial differential equation in dimension `N`. `AbstractDifferentialOperator` types are used -to define [`AbstractDifferentialOperatorKernel`s](@ref AbstractDifferentialOperatorKernel). +A partial differential operator in dimension `N`. + +`AbstractDifferentialOperator` types are used to define [`AbstractKernel`s](@ref +AbstractKernel) related to fundamental solutions of differential operators. """ abstract type AbstractDifferentialOperator{N} end ambient_dimension(::AbstractDifferentialOperator{N}) where {N} = N -""" - abstract type AbstractDifferentialOperatorKernel{T,Op} <: AbstractKernel{T} - -An [`AbstractKernel`](@ref) with an associated `pde::Op` field. -""" -abstract type AbstractDifferentialOperatorKernel{T,Op} <: AbstractKernel{T} end - -""" - pde(K::AbstractDifferentialOperatorKernel) - -Return the underlying `AbstractDifferentialOperator` associated with the kernel `K`. -""" -pde(k::AbstractDifferentialOperatorKernel) = k.pde - -parameters(k::AbstractDifferentialOperatorKernel) = parameters(pde(k)) - -# convenient constructor for e.g. SingleLayerKernel(pde,Float64) or DoubleLayerKernel(pde,ComplexF64) +# convenient constructor for e.g. SingleLayerKernel(op,Float64) or DoubleLayerKernel(op,ComplexF64) function (::Type{K})( - pde::Op, - ::Type{T} = default_kernel_eltype(pde), -) where {T,Op,K<:AbstractDifferentialOperatorKernel} - return K{T,Op}(pde) + op::Op, + ::Type{T} = default_kernel_eltype(op), +) where {T,Op,K<:AbstractKernel} + return K{T,Op}(op) end """ - struct SingleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} + struct SingleLayerKernel{T,Op} <: AbstractKernel{T,Op} The free-space single-layer kernel (i.e. the fundamental solution) of an `Op <: AbstractDifferentialOperator`. """ -struct SingleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} - pde::Op +struct SingleLayerKernel{T,Op} <: AbstractKernel{T,Op} + op::Op end """ - struct DoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} + struct DoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} Given an operator `Op`, construct its free-space double-layer kernel. This corresponds to the `γ₁` trace of the [`SingleLayerKernel`](@ref). For operators such as [`Laplace`](@ref) or [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the source variable. """ -struct DoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} - pde::Op +struct DoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} + op::Op end """ - struct AdjointDoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} + struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} Given an operator `Op`, construct its free-space adjoint double-layer kernel. This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the @@ -87,12 +64,12 @@ This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the target variable. """ -struct AdjointDoubleLayerKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} - pde::Op +struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} + op::Op end """ - struct HyperSingularKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} + struct HyperSingularKernel{T,Op} <: AbstractKernel{T,Op} Given an operator `Op`, construct its free-space hypersingular kernel. This corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the @@ -100,8 +77,8 @@ corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative respect to the target variable of the `DoubleLayerKernel`. """ -struct HyperSingularKernel{T,Op} <: AbstractDifferentialOperatorKernel{T,Op} - pde::Op +struct HyperSingularKernel{T,Op} <: AbstractKernel{T,Op} + op::Op end ################################################################################ @@ -113,12 +90,15 @@ struct Laplace{N} <: AbstractDifferentialOperator{N} end """ Laplace(; dim) -Laplace equation in `N` dimension: Δu = 0. +Laplace's differential operator in `dim` dimension: ``-Δu``. +``` + +Note the **negative sign** in the definition. """ Laplace(; dim) = Laplace{dim}() -function Base.show(io::IO, pde::Laplace) - return print(io, "Δu = 0") +function Base.show(io::IO, op::Laplace{N}) where {N} + return print(io, "Laplace operator in $N dimensions: -Δu") end default_kernel_eltype(::Laplace) = Float64 @@ -201,19 +181,28 @@ end """ Yukawa(; λ, dim) -Yukawa equation, also known as modified Helmholtz, in `N` dimensions: Δu - λ²u = -0. The parameter `λ` is a positive number. +Yukawa operator, also known as modified Helmholtz, in `dim` dimensions: ``-Δu + λ²u``. + +The parameter `λ` is a positive number. Note the **negative sign** in front of +the Laplacian. """ function Yukawa(; λ, dim) @assert λ > 0 "λ must be a positive number" return Yukawa{dim,typeof(λ)}(λ) end +""" + const ModifiedHelmholtz + +Type alias for the [`Yukawa`](@ref) operator. +""" const ModifiedHelmholtz = Yukawa -Base.show(io::IO, ::Yukawa) = print(io, "Δu + λ²u = 0") +function Base.show(io::IO, ::Yukawa{N}) where {N} + return print(io, "Yukawa operator in $N dimensions: -Δu + λ²u") +end -parameters(pde::Yukawa) = pde.λ +parameters(op::Yukawa) = op.λ default_kernel_eltype(::Yukawa) = Float64 default_density_eltype(::Yukawa) = Float64 @@ -262,7 +251,7 @@ end function (HS::HyperSingularKernel{T,<:Yukawa{N,K}})(target, source)::T where {N,T,K} x, y, nx, ny = coords(target), coords(source), normal(target), normal(source) - λ = parameters(pde(HS)) + λ = parameters(op(HS)) r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -299,17 +288,17 @@ end """ Helmholtz(; k, dim) -Helmholtz equation in `N` dimensions: Δu + k²u = 0. The parameter `k` can be a -real or complex number. +Helmholtz operator in `dim` dimensions: `-Δu - k²u`. -For purely imaginary wavenumbers, consider using the [`Yukawa`](@ref) kernel. +The parameter `k` can be a real or complex number. For purely imaginary +wavenumbers, consider using the [`Yukawa`](@ref) kernel. """ function Helmholtz(; k, dim) if k isa Complex @assert imag(k) ≥ 0 "k must have a non-negative imaginary part" if iszero(real(k)) - msg = """Purely imaginary wavenumber detected in Helmholtz equation. - Creating a modified Helmholtz (Yukawa) PDE instead.""" + msg = """Purely imaginary wavenumber detected in Helmholtz operator. + Creating a modified Helmholtz (Yukawa) op instead.""" @warn msg return Yukawa(; λ = imag(k), dim = dim) elseif iszero(imag(k)) @@ -319,12 +308,11 @@ function Helmholtz(; k, dim) return Helmholtz{dim,typeof(k)}(k) end -function Base.show(io::IO, ::Helmholtz) - # k = parameters(pde) - return print(io, "Δu + k² u = 0") +function Base.show(io::IO, ::Helmholtz{N}) where {N} + return print(io, "Helmholtz operator in $N dimensions: -Δu - k²u") end -parameters(pde::Helmholtz) = pde.k +parameters(op::Helmholtz) = op.k default_kernel_eltype(::Helmholtz) = ComplexF64 default_density_eltype(::Helmholtz) = ComplexF64 @@ -381,8 +369,8 @@ end # Hypersingular kernel function (HS::HyperSingularKernel{T,S})(target, source)::T where {T,S<:Helmholtz} x, y, nx, ny = coords(target), coords(source), normal(target), normal(source) - N = ambient_dimension(pde(HS)) - k = parameters(pde(HS)) + N = ambient_dimension(op(HS)) + k = parameters(op(HS)) r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -409,11 +397,17 @@ end struct Stokes{N,T} <: AbstractDifferentialOperator{N} μ::T end + +""" + Stokes(; μ, dim) + +Stokes operator in `dim` dimensions: ``[-μΔu + ∇p, ∇⋅u]``. +""" Stokes(; μ, dim = 3) = Stokes{dim}(μ) Stokes{N}(μ::T) where {N,T} = Stokes{N,T}(μ) -function Base.show(io::IO, pde::Stokes) - return println(io, "μΔu -∇p = 0, ∇⋅u = 0") +function Base.show(io::IO, op::Stokes{N}) where {N} + return println(io, "Stokes operator in $N dimensions: [-μΔu + ∇p, ∇⋅u]") end parameters(s::Stokes) = s.μ @@ -478,8 +472,10 @@ end """ struct Elastostatic{N,T} <: AbstractDifferentialOperator{N} -Elastostatic equation in `N` dimensions: μΔu + (μ+λ)∇(∇⋅u) = 0. Note that the -displacement u is a vector of length `N` since this is a vectorial problem. +Elastostatic operator in `N` dimensions: -μΔu - (μ+λ)∇(∇⋅u) + +Note that the displacement ``u`` is a vector of length `N` since this is a +vectorial problem. """ struct Elastostatic{N,T} <: AbstractDifferentialOperator{N} μ::T @@ -488,17 +484,17 @@ end Elastostatic(; μ, λ, dim) = Elastostatic{dim}(promote(μ, λ)...) Elastostatic{N}(μ::T, λ::T) where {N,T} = Elastostatic{N,T}(μ, λ) -function Base.show(io::IO, pde::Elastostatic) - return print(io, "μΔu + (μ+λ)∇(∇⋅u) = 0") +function Base.show(io::IO, op::Elastostatic) + return print(io, "Elastostatic operator in $N dimensions: -μΔu - (μ+λ)∇(∇⋅u)") end -parameters(pde::Elastostatic) = pde.μ, pde.λ +parameters(op::Elastostatic) = op.μ, op.λ default_kernel_eltype(::Elastostatic{N}) where {N} = SMatrix{N,N,Float64,N * N} default_density_eltype(::Elastostatic{N}) where {N} = SVector{N,Float64} function (SL::SingleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(pde(SL)) + μ, λ = parameters(op(SL)) ν = λ / (2 * (μ + λ)) x = coords(target) y = coords(source) @@ -514,7 +510,7 @@ function (SL::SingleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where { end function (DL::DoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(pde(DL)) + μ, λ = parameters(op(DL)) ν = λ / (2 * (μ + λ)) x = coords(target) y = coords(source) @@ -539,7 +535,7 @@ function (DL::DoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where { end function (ADL::AdjointDoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(pde(ADL)) + μ, λ = parameters(op(ADL)) ν = λ / (2 * (μ + λ)) x = coords(target) nx = normal(target) @@ -568,7 +564,7 @@ function (ADL::AdjointDoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T end function (HS::HyperSingularKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(pde(HS)) + μ, λ = parameters(op(HS)) ν = λ / (2 * (μ + λ)) x = coords(target) nx = normal(target) From cc5f2872e9878476c9b69f498efd5570b7401593 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Fri, 11 Oct 2024 19:04:33 +0200 Subject: [PATCH 4/7] bug fixes --- ext/IntiFMM3DExt.jl | 10 ++++---- src/kernels.jl | 59 ++++++++++++++++++--------------------------- src/mesh.jl | 4 +-- 3 files changed, 30 insertions(+), 43 deletions(-) diff --git a/ext/IntiFMM3DExt.jl b/ext/IntiFMM3DExt.jl index 93d7e666..461f996a 100644 --- a/ext/IntiFMM3DExt.jl +++ b/ext/IntiFMM3DExt.jl @@ -78,7 +78,7 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) # Helmholtz elseif K isa Inti.SingleLayerKernel{ComplexF64,<:Inti.Helmholtz{3}} charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant @. charges = 1 / (4 * π) * weights * x @@ -91,7 +91,7 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) normals[:, j] = Inti.normal(iop.source[j]) end dipvecs = similar(normals, ComplexF64) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n @@ -106,7 +106,7 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) xnormals[:, j] = Inti.normal(iop.target[j]) end charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant @. charges = 1 / (4 * π) * weights * x @@ -123,7 +123,7 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) ynormals[:, j] = Inti.normal(iop.source[j]) end dipvecs = similar(ynormals, ComplexF64) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n @@ -138,7 +138,7 @@ function Inti._assemble_fmm3d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) stoklet = Matrix{Float64}(undef, 3, n) return LinearMaps.LinearMap{SVector{3,Float64}}(m, n) do y, x # multiply by weights and constant - stoklet[:] = 1 / (4 * π * K.pde.μ) .* reinterpret(Float64, weights .* x) + stoklet[:] = 1 / (4 * π * K.op.μ) .* reinterpret(Float64, weights .* x) out = FMM3D.stfmm3d(rtol, sources; stoklet, targets, ppregt = 1) return copyto!(y, reinterpret(T, out.pottarg)) end diff --git a/src/kernels.jl b/src/kernels.jl index d06fe82d..74e432cb 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -34,29 +34,29 @@ function (::Type{K})( end """ - struct SingleLayerKernel{T,Op} <: AbstractKernel{T,Op} + struct SingleLayerKernel{T,Op} <: AbstractKernel{T} The free-space single-layer kernel (i.e. the fundamental solution) of an `Op <: AbstractDifferentialOperator`. """ -struct SingleLayerKernel{T,Op} <: AbstractKernel{T,Op} +struct SingleLayerKernel{T,Op} <: AbstractKernel{T} op::Op end """ - struct DoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} + struct DoubleLayerKernel{T,Op} <: AbstractKernel{T} Given an operator `Op`, construct its free-space double-layer kernel. This corresponds to the `γ₁` trace of the [`SingleLayerKernel`](@ref). For operators such as [`Laplace`](@ref) or [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the source variable. """ -struct DoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} +struct DoubleLayerKernel{T,Op} <: AbstractKernel{T} op::Op end """ - struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} + struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T} Given an operator `Op`, construct its free-space adjoint double-layer kernel. This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the @@ -64,12 +64,12 @@ This corresponds to the `transpose(γ₁,ₓ[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative of the fundamental solution respect to the target variable. """ -struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T,Op} +struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T} op::Op end """ - struct HyperSingularKernel{T,Op} <: AbstractKernel{T,Op} + struct HyperSingularKernel{T,Op} <: AbstractKernel{T} Given an operator `Op`, construct its free-space hypersingular kernel. This corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the @@ -77,7 +77,7 @@ corresponds to the `transpose(γ₁,ₓγ₁[G])`, where `G` is the [`Helmholtz`](@ref), this is simply the normal derivative respect to the target variable of the `DoubleLayerKernel`. """ -struct HyperSingularKernel{T,Op} <: AbstractKernel{T,Op} +struct HyperSingularKernel{T,Op} <: AbstractKernel{T} op::Op end @@ -104,8 +104,6 @@ end default_kernel_eltype(::Laplace) = Float64 default_density_eltype(::Laplace) = Float64 -parameters(::Laplace) = nothing - function (SL::SingleLayerKernel{T,Laplace{N}})( target, source, @@ -202,15 +200,13 @@ function Base.show(io::IO, ::Yukawa{N}) where {N} return print(io, "Yukawa operator in $N dimensions: -Δu + λ²u") end -parameters(op::Yukawa) = op.λ - default_kernel_eltype(::Yukawa) = Float64 default_density_eltype(::Yukawa) = Float64 function (SL::SingleLayerKernel{T,<:Yukawa{N,K}})(target, source)::T where {N,T,K} x = coords(target) y = coords(source) - λ = parameters(SL) + λ = SL.op.λ r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -223,7 +219,7 @@ end function (DL::DoubleLayerKernel{T,Yukawa{N,K}})(target, source)::T where {N,T,K} x, y, ny = coords(target), coords(source), normal(source) - λ = parameters(DL) + λ = DL.op.λ r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -237,7 +233,7 @@ end function (ADL::AdjointDoubleLayerKernel{T,<:Yukawa{N,K}})(target, source)::T where {N,T,K} x, y, nx = coords(target), coords(source), normal(target) - λ = parameters(ADL) + λ = ADL.op.λ r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -251,7 +247,7 @@ end function (HS::HyperSingularKernel{T,<:Yukawa{N,K}})(target, source)::T where {N,T,K} x, y, nx, ny = coords(target), coords(source), normal(target), normal(source) - λ = parameters(op(HS)) + λ = HS.op.λ r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -312,8 +308,6 @@ function Base.show(io::IO, ::Helmholtz{N}) where {N} return print(io, "Helmholtz operator in $N dimensions: -Δu - k²u") end -parameters(op::Helmholtz) = op.k - default_kernel_eltype(::Helmholtz) = ComplexF64 default_density_eltype(::Helmholtz) = ComplexF64 @@ -323,7 +317,7 @@ hankelh1(n, x::Complex) = SpecialFunctions.hankelh1(n, x) function (SL::SingleLayerKernel{T,<:Helmholtz{N}})(target, source)::T where {N,T} x = coords(target) y = coords(source) - k = parameters(SL) + k = SL.op.k r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -337,7 +331,7 @@ end # Double Layer Kernel function (DL::DoubleLayerKernel{T,<:Helmholtz{N}})(target, source)::T where {N,T} x, y, ny = coords(target), coords(source), normal(source) - k = parameters(DL) + k = DL.op.k r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -353,7 +347,7 @@ end # Adjoint double Layer Kernel function (ADL::AdjointDoubleLayerKernel{T,<:Helmholtz{N}})(target, source)::T where {N,T} x, y, nx = coords(target), coords(source), normal(target) - k = parameters(ADL) + k = ADL.op.k r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -367,10 +361,9 @@ function (ADL::AdjointDoubleLayerKernel{T,<:Helmholtz{N}})(target, source)::T wh end # Hypersingular kernel -function (HS::HyperSingularKernel{T,S})(target, source)::T where {T,S<:Helmholtz} +function (HS::HyperSingularKernel{T,<:Helmholtz{N}})(target, source)::T where {N,T} x, y, nx, ny = coords(target), coords(source), normal(target), normal(source) - N = ambient_dimension(op(HS)) - k = parameters(op(HS)) + k = HS.op.k r = x - y d = norm(r) d ≤ SAME_POINT_TOLERANCE && return zero(T) @@ -410,14 +403,12 @@ function Base.show(io::IO, op::Stokes{N}) where {N} return println(io, "Stokes operator in $N dimensions: [-μΔu + ∇p, ∇⋅u]") end -parameters(s::Stokes) = s.μ - default_kernel_eltype(::Stokes{N}) where {N} = SMatrix{N,N,Float64,N * N} default_density_eltype(::Stokes{N}) where {N} = SVector{N,Float64} # Single Layer function (SL::SingleLayerKernel{T,<:Stokes{N}})(target, source)::T where {N,T} - μ = parameters(SL) + μ = SL.op.μ x = coords(target) y = coords(source) r = x - y @@ -433,7 +424,6 @@ end # Double Layer Kernel function (DL::DoubleLayerKernel{T,<:Stokes{N}})(target, source)::T where {N,T} - μ = parameters(DL) x = coords(target) y = coords(source) ny = normal(source) @@ -449,7 +439,6 @@ end # Double Layer Kernel function (ADL::AdjointDoubleLayerKernel{T,<:Stokes{N}})(target, source)::T where {N,T} - μ = parameters(ADL) x = coords(target) nx = normal(target) y = coords(source) @@ -484,17 +473,15 @@ end Elastostatic(; μ, λ, dim) = Elastostatic{dim}(promote(μ, λ)...) Elastostatic{N}(μ::T, λ::T) where {N,T} = Elastostatic{N,T}(μ, λ) -function Base.show(io::IO, op::Elastostatic) +function Base.show(io::IO, op::Elastostatic{N}) where {N} return print(io, "Elastostatic operator in $N dimensions: -μΔu - (μ+λ)∇(∇⋅u)") end -parameters(op::Elastostatic) = op.μ, op.λ - default_kernel_eltype(::Elastostatic{N}) where {N} = SMatrix{N,N,Float64,N * N} default_density_eltype(::Elastostatic{N}) where {N} = SVector{N,Float64} function (SL::SingleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(op(SL)) + μ, λ = SL.op.μ, SL.op.λ ν = λ / (2 * (μ + λ)) x = coords(target) y = coords(source) @@ -510,7 +497,7 @@ function (SL::SingleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where { end function (DL::DoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(op(DL)) + μ, λ = DL.op.μ, DL.op.λ ν = λ / (2 * (μ + λ)) x = coords(target) y = coords(source) @@ -535,7 +522,7 @@ function (DL::DoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where { end function (ADL::AdjointDoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(op(ADL)) + μ, λ = ADL.op.μ, ADL.op.λ ν = λ / (2 * (μ + λ)) x = coords(target) nx = normal(target) @@ -564,7 +551,7 @@ function (ADL::AdjointDoubleLayerKernel{T,<:Elastostatic{N}})(target, source)::T end function (HS::HyperSingularKernel{T,<:Elastostatic{N}})(target, source)::T where {N,T} - μ, λ = parameters(op(HS)) + μ, λ = HS.op.μ, HS.op.λ ν = λ / (2 * (μ + λ)) x = coords(target) nx = normal(target) diff --git a/src/mesh.jl b/src/mesh.jl index df6b5bce..e79c07cd 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -60,8 +60,8 @@ function element_types end """ struct Mesh{N,T} <: AbstractMesh{N,T} -Unstructured mesh is defined by a set of `nodes`` (of type `SVector{N,T}`), and -a dictionary mapping element types to connectivity matrices. Each columns of a +Unstructured mesh defined by a set of `nodes`` (of type `SVector{N,T}`), and a +dictionary mapping element types to connectivity matrices. Each columns of a given connectivity matrix stores the integer tags of the nodes in the mesh comprising the element. From b5d38c72637c1562c7f33ffe87c12dbc6e03b619 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Sun, 13 Oct 2024 09:22:57 +0200 Subject: [PATCH 5/7] change `pde` to `op` in docs and `src` --- docs/src/examples/cavity2d_scattering.jl | 6 +- docs/src/examples/graded_mesh.jl | 10 +- docs/src/examples/helmholtz_scattering.jl | 18 +- docs/src/examples/lippmann_schwinger.jl | 6 +- docs/src/examples/poisson.jl | 10 +- docs/src/examples/poisson.md | 14 +- docs/src/examples/stokes_drag.jl | 10 +- docs/src/examples/transmission.jl | 16 +- docs/src/index.md | 6 +- .../pluto-examples/helmholtz_scattering.jl | 524 +++++++++--------- docs/src/pluto-examples/poisson.jl | 170 +++--- docs/src/tutorials/compression_methods.md | 4 +- docs/src/tutorials/getting_started.md | 10 +- docs/src/tutorials/integral_operators.md | 12 +- docs/src/tutorials/layer_potentials.md | 10 +- ext/IntiFMM2DExt.jl | 8 +- ext/IntiFMMLIB2DExt.jl | 8 +- src/api.jl | 42 +- src/bdim.jl | 14 +- src/nystrom.jl | 6 +- src/vdim.jl | 28 +- test/adaptive_test.jl | 16 +- test/convergence/dim_convergence.jl | 18 +- .../convergence/vdim_potential_convergence.jl | 6 +- test/convergence/vdim_potential_script.jl | 6 +- test/convergence/vdim_potential_script_3d.jl | 6 +- test/dim_test.jl | 28 +- test/fmm2d_test.jl | 12 +- test/fmm3d_test.jl | 18 +- test/hmatrices_test.jl | 6 +- test/integral_operator_test.jl | 4 +- test/quadrature_test.jl | 4 +- test/vdim_elastostatic_test.jl | 87 +++ 33 files changed, 620 insertions(+), 523 deletions(-) create mode 100644 test/vdim_elastostatic_test.jl diff --git a/docs/src/examples/cavity2d_scattering.jl b/docs/src/examples/cavity2d_scattering.jl index 7d3c234b..e748a530 100644 --- a/docs/src/examples/cavity2d_scattering.jl +++ b/docs/src/examples/cavity2d_scattering.jl @@ -54,9 +54,9 @@ Q = Inti.Quadrature(Γ_msh; qorder) println("Number of quadrature points: ", length(Q)) ## Setup the integral operators -pde = Inti.Helmholtz(; dim = 2, k) +op = Inti.Helmholtz(; dim = 2, k) S, D = Inti.single_double_layer(; - pde, + op target = Q, source = Q, correction = (method = :none,), @@ -76,7 +76,7 @@ L = I / 2 + LinearMap(D) - im * k * LinearMap(S) σ = gmres(L, g; restart = 1000, maxiter = 400, abstol = 1e-4, verbose = true) ## Plot a heatmap of the solution -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) u = (x) -> 𝒟[σ](x) - im * k * 𝒮[σ](x) xx = yy = range(-2, 2; step = meshsize) colorrange = (-2, 2) diff --git a/docs/src/examples/graded_mesh.jl b/docs/src/examples/graded_mesh.jl index 54598e05..2e88e095 100644 --- a/docs/src/examples/graded_mesh.jl +++ b/docs/src/examples/graded_mesh.jl @@ -37,9 +37,9 @@ Q = Inti.Quadrature(msh; qorder = 3) ## ## Create integral operators -pde = Inti.Helmholtz(; dim = 2, k) +op = Inti.Helmholtz(; dim = 2, k) S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :fmm, tol = 1e-6), @@ -47,11 +47,11 @@ S, D = Inti.single_double_layer(; ) L = I / 2 + D - im * k * S -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) ## Test with a source inside for validation if TEST - G = Inti.SingleLayerKernel(pde) + G = Inti.SingleLayerKernel(op) uₑ = x -> G(x, SVector(0.4, 0.4)) rhs = map(q -> uₑ(q.coords), Q) σ = gmres(L, rhs; restart = 100, maxiter = 100, abstol = 1e-8, verbose = true) @@ -76,7 +76,7 @@ if PLOT target = [SVector(x, y) for x in xx, y in yy] Spot, Dpot = Inti.single_double_layer(; - pde, + op, target = target, source = Q, compression = (method = :fmm, tol = 1e-4), diff --git a/docs/src/examples/helmholtz_scattering.jl b/docs/src/examples/helmholtz_scattering.jl index 25ee13b6..80827cd8 100644 --- a/docs/src/examples/helmholtz_scattering.jl +++ b/docs/src/examples/helmholtz_scattering.jl @@ -197,9 +197,9 @@ abs(Inti.integrate(x -> 1, Q) - 2π) # discrete approximation to the integral operators ``\mathrm{S}`` and # ``\mathrm{D}`` as follows: -pde = Inti.Helmholtz(; k, dim = 2) +op = Inti.Helmholtz(; k, dim = 2) S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), @@ -259,7 +259,7 @@ nothing #hide # Inti.single_double_layer_potential) to obtain the single- and double-layer # potentials, and then combine them as follows: -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) nothing #hide @@ -338,7 +338,7 @@ msh = Inti.meshgen(Γ; meshsize) Γ_msh = view(msh, Γ) Q = Inti.Quadrature(Γ_msh; qorder) S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), @@ -347,7 +347,7 @@ S, D = Inti.single_double_layer(; L = I / 2 + D - im * k * S rhs = map(q -> -uᵢ(q.coords), Q) σ = L \ rhs -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) vals = map(pt -> Inti.isinside(pt, Q) ? NaN : real(uₛ(pt) + uᵢ(pt)), Iterators.product(xx, yy)) @@ -445,9 +445,9 @@ nothing #hide # Next we assemble the integral operators, indicating that we wish to compress # them using hierarchical matrices: using HMatrices -pde = Inti.Helmholtz(; k, dim = 3) +op = Inti.Helmholtz(; k, dim = 3) S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :hmatrix, tol = 1e-4), @@ -497,7 +497,7 @@ end # As before, let us represent the solution using `IntegralPotential`s: -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) nothing #hide @@ -550,7 +550,7 @@ end target = Inti.nodes(Σ_msh) S, D = Inti.single_double_layer(; - pde, + op, target, source = Q, compression = (method = :hmatrix, tol = 1e-4), diff --git a/docs/src/examples/lippmann_schwinger.jl b/docs/src/examples/lippmann_schwinger.jl index 59860781..b3dec4dc 100644 --- a/docs/src/examples/lippmann_schwinger.jl +++ b/docs/src/examples/lippmann_schwinger.jl @@ -86,12 +86,12 @@ dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) # ## Volume Integral Operators and Volume Integral Equations using FMMLIB2D -pde = Inti.Helmholtz(; dim = 2, k = k₁) +op = Inti.Helmholtz(; dim = 2, k = k₁) # With quadratures constructed on the volume, we can define a discrete approximation # to the volume integral operator ``\mathcal{V}`` using VDIM. V_d2d = Inti.volume_potential(; - pde, + op, target = Ωₕ_quad, source = Ωₕ_quad, compression = (method = :fmm, tol = 1e-7), @@ -124,7 +124,7 @@ u, hist = gmres(L, rhs; log = true, abstol = 1e-7, verbose = false, restart = 200, maxiter = 200) @show hist -𝒱 = Inti.IntegralPotential(Inti.SingleLayerKernel(pde), Ωₕ_quad) +𝒱 = Inti.IntegralPotential(Inti.SingleLayerKernel(op), Ωₕ_quad) # The representation formula gives the solution in $\R^2 \setminus \Omega$: uˢ = (x) -> uⁱ(x) - k₁^2 * 𝒱[refr_map_d.*u](x) diff --git a/docs/src/examples/poisson.jl b/docs/src/examples/poisson.jl index dfaec803..3be084fe 100644 --- a/docs/src/examples/poisson.jl +++ b/docs/src/examples/poisson.jl @@ -107,18 +107,18 @@ fₑ = (x) -> -8 * uₑ(x) # ## Boundary and integral operators using FMM2D -pde = Inti.Laplace(; dim = 2) +op = Inti.Laplace(; dim = 2) ## Boundary operators S_b2b, D_b2b = Inti.single_double_layer(; - pde, + op, target = Γₕ_quad, source = Γₕ_quad, compression = (method = :fmm, tol = 1e-12), correction = (method = :dim,), ) S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ωₕ_quad, source = Γₕ_quad, compression = (method = :fmm, tol = 1e-12), @@ -127,14 +127,14 @@ S_b2d, D_b2d = Inti.single_double_layer(; ## Volume potentials V_d2d = Inti.volume_potential(; - pde, + op, target = Ωₕ_quad, source = Ωₕ_quad, compression = (method = :fmm, tol = 1e-12), correction = (method = :dim, interpolation_order), ) V_d2b = Inti.volume_potential(; - pde, + op, target = Γₕ_quad, source = Ωₕ_quad, compression = (method = :fmm, tol = 1e-12), diff --git a/docs/src/examples/poisson.md b/docs/src/examples/poisson.md index c65f5e3c..48555ac4 100644 --- a/docs/src/examples/poisson.md +++ b/docs/src/examples/poisson.md @@ -123,10 +123,10 @@ operator mapping to points on the boundary, i.e. operator: ```@example poisson using FMM2D #to accelerate the maps -pde = Inti.Laplace(; dim = 2) +op = Inti.Laplace(; dim = 2) # Newtonian potential mapping domain to boundary V_d2b = Inti.volume_potential(; - pde, + op, target = Γ_quad, source = Ω_quad, compression = (method = :fmm, tol = 1e-12), @@ -144,7 +144,7 @@ equation: ```@example poisson # Single and double layer operators on Γ S_b2b, D_b2b = Inti.single_double_layer(; - pde, + op, target = Γ_quad, source = Γ_quad, compression = (method = :fmm, tol = 1e-12), @@ -193,8 +193,8 @@ nothing # hide With the density function at hand, we can now reconstruct our approximate solution: ```@example poisson -G = Inti.SingleLayerKernel(pde) -dG = Inti.DoubleLayerKernel(pde) +G = Inti.SingleLayerKernel(op) +dG = Inti.DoubleLayerKernel(op) 𝒱 = Inti.IntegralPotential(G, Ω_quad) 𝒟 = Inti.IntegralPotential(dG, Γ_quad) u = (x) -> 𝒱[f](x) + 𝒟[σ](x) @@ -223,7 +223,7 @@ solution ``u`` at all the quadrature nodes of ``\Omega``: ```@example poisson V_d2d = Inti.volume_potential(; - pde, + op, target = Ω_quad, source = Ω_quad, compression = (method = :fmm, tol = 1e-8), @@ -236,7 +236,7 @@ our mesh nodes: ```@example poisson S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ω_quad, source = Γ_quad, compression = (method = :fmm, tol = 1e-8), diff --git a/docs/src/examples/stokes_drag.jl b/docs/src/examples/stokes_drag.jl index 34f2a0af..be2e16aa 100644 --- a/docs/src/examples/stokes_drag.jl +++ b/docs/src/examples/stokes_drag.jl @@ -63,10 +63,10 @@ Q = Inti.Quadrature(Γ_msh; qorder = 2) @show length(Q) @show abs(Inti.integrate(x -> 1, Q) - 4π * R^2) -## the pde and its integral kernels -pde = Inti.Stokes(; dim = 3, μ) -G = Inti.SingleLayerKernel(pde) -dG = Inti.DoubleLayerKernel(pde) +## the op and its integral kernels +op = Inti.Stokes(; dim = 3, μ) +G = Inti.SingleLayerKernel(op) +dG = Inti.DoubleLayerKernel(op) ## choice of a integral representation T = SVector{3,Float64} @@ -85,7 +85,7 @@ Smat = Inti.assemble_matrix(Sop) ## integral operators defined on the boundary S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), diff --git a/docs/src/examples/transmission.jl b/docs/src/examples/transmission.jl index 2f90a08d..20548403 100644 --- a/docs/src/examples/transmission.jl +++ b/docs/src/examples/transmission.jl @@ -44,12 +44,12 @@ msh = Inti.import_mesh(name; dim = 2) Q = Inti.Quadrature(Γ_msh; qorder) -pde₁ = Inti.Helmholtz(; k = k₁, dim = 2) -pde₂ = Inti.Helmholtz(; k = k₂, dim = 2) +op₁ = Inti.Helmholtz(; k = k₁, dim = 2) +op₂ = Inti.Helmholtz(; k = k₂, dim = 2) using FMMLIB2D S₁, D₁ = Inti.single_double_layer(; - pde = pde₁, + op = op₁, target = Q, source = Q, compression = (method = :fmm, tol = :1e-8), @@ -57,7 +57,7 @@ S₁, D₁ = Inti.single_double_layer(; ) K₁, N₁ = Inti.adj_double_layer_hypersingular(; - pde = pde₁, + op = op₁, target = Q, source = Q, compression = (method = :fmm, tol = :1e-8), @@ -65,7 +65,7 @@ K₁, N₁ = Inti.adj_double_layer_hypersingular(; ) S₂, D₂ = Inti.single_double_layer(; - pde = pde₂, + op = op₂, target = Q, source = Q, compression = (method = :fmm, tol = :1e-8), @@ -73,7 +73,7 @@ S₂, D₂ = Inti.single_double_layer(; ) K₂, N₂ = Inti.adj_double_layer_hypersingular(; - pde = pde₂, + op = op₂, target = Q, source = Q, compression = (method = :fmm, tol = :1e-8), @@ -120,8 +120,8 @@ nQ = size(Q, 1) sol = reshape(sol, nQ, 2) φ, ψ = sol[:, 1], sol[:, 2] -𝒮₁, 𝒟₁ = Inti.single_double_layer_potential(; pde = pde₁, source = Q) -𝒮₂, 𝒟₂ = Inti.single_double_layer_potential(; pde = pde₂, source = Q) +𝒮₁, 𝒟₁ = Inti.single_double_layer_potential(; op = op₁, source = Q) +𝒮₂, 𝒟₂ = Inti.single_double_layer_potential(; op = op₂, source = Q) v₁ = x -> 𝒟₁[φ](x) - 𝒮₁[ψ](x) v₂ = x -> -𝒟₂[φ](x) + 𝒮₂[ψ](x) diff --git a/docs/src/index.md b/docs/src/index.md index c14101fc..fdfa22bd 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -134,9 +134,9 @@ end msh = Inti.meshgen(Γ; meshsize = 0.1) Q = Inti.Quadrature(msh; qorder = 5) # create the integral operators -pde = Inti.Laplace(;dim=2) +op = Inti.Laplace(;dim=2) S, _ = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), @@ -148,7 +148,7 @@ g = map(q -> uₑ(q.coords), Q) # value at quad nodes # solve for σ σ = S \ g # use the single-layer potential to evaluate the solution -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) uₕ = x -> 𝒮[σ](x) ``` diff --git a/docs/src/pluto-examples/helmholtz_scattering.jl b/docs/src/pluto-examples/helmholtz_scattering.jl index c2c1b554..2207f88a 100644 --- a/docs/src/pluto-examples/helmholtz_scattering.jl +++ b/docs/src/pluto-examples/helmholtz_scattering.jl @@ -6,23 +6,23 @@ using InteractiveUtils # ╔═╡ 6b6d5e44-304d-4bed-912b-e4b89fce7c5a begin - import Pkg as _Pkg + import Pkg as _Pkg haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"]) - using PlutoUI: TableOfContents -end ; + using PlutoUI: TableOfContents +end; # ╔═╡ ef067551-aa44-4099-b32a-08debb81ee79 begin - using Inti - using LinearAlgebra - using StaticArrays - using Gmsh - using Meshes - using GLMakie - using SpecialFunctions - using GSL - using IterativeSolvers - using LinearMaps + using Inti + using LinearAlgebra + using StaticArrays + using Gmsh + using Meshes + using GLMakie + using SpecialFunctions + using GSL + using IterativeSolvers + using LinearMaps end # ╔═╡ a6c80f87-ff47-496f-8925-1275b58b02e1 @@ -122,11 +122,11 @@ and setup some of the (global) problem parameters: # ╔═╡ 4cabe108-93c6-4fbf-b8ec-37277abbbda1 begin - k = 4π - λ = 2π / k - qorder = 4 # quadrature order - gorder = 2 # order of geometrical approximation - nothing #hide + k = 4π + λ = 2π / k + qorder = 4 # quadrature order + gorder = 2 # order of geometrical approximation + nothing #hide end # ╔═╡ a33f066e-307c-42c4-836b-5b7d8236a893 @@ -138,22 +138,22 @@ We will use [Gmsh API](https://gmsh.info/doc/texinfo/gmsh.html#Gmsh-application- # ╔═╡ cb736b4d-bbd6-456c-9041-f69b2155a470 begin - function gmsh_circle(; name, meshsize, order = 1, radius = 1, center = (0, 0)) - try - gmsh.initialize() - gmsh.model.add("circle-mesh") - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - gmsh.model.occ.addDisk(center[1], center[2], 0, radius, radius) - gmsh.model.occ.synchronize() - gmsh.model.mesh.generate(1) - gmsh.model.mesh.setOrder(order) - gmsh.write(name) - finally - gmsh.finalize() - end - end - nothing #hide + function gmsh_circle(; name, meshsize, order = 1, radius = 1, center = (0, 0)) + try + gmsh.initialize() + gmsh.model.add("circle-mesh") + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + gmsh.model.occ.addDisk(center[1], center[2], 0, radius, radius) + gmsh.model.occ.synchronize() + gmsh.model.mesh.generate(1) + gmsh.model.mesh.setOrder(order) + gmsh.write(name) + finally + gmsh.finalize() + end + end + nothing #hide end # ╔═╡ 2ed4fd2d-b3c2-4211-bf3f-c551923e65cf @@ -163,9 +163,9 @@ Let us now use `gmsh_circle` to create a `circle.msh` file. As customary in wave # ╔═╡ 12e20620-042c-4040-bd64-af01bd0f2ad9 begin - name = joinpath(@__DIR__, "circle.msh") - meshsize = λ / 5 - gmsh_circle(; meshsize, order = gorder, name) + name = joinpath(@__DIR__, "circle.msh") + meshsize = λ / 5 + gmsh_circle(; meshsize, order = gorder, name) end # ╔═╡ 24454bd5-00cf-4282-926f-f1324141fe26 @@ -176,8 +176,8 @@ We can now import the file and parse the mesh and domain information into `Inti. # ╔═╡ ff73663c-7cf4-4b23-83b5-095dc66f711f # ╠═╡ show_logs = false begin - Inti.clear_entities!() # empty the entity cache - msh = Inti.import_mesh(name; dim = 2) + Inti.clear_entities!() # empty the entity cache + msh = Inti.import_mesh(name; dim = 2) end # ╔═╡ 176f1fef-f761-4ad0-8fea-982eab4fe24d @@ -196,10 +196,10 @@ for: # ╔═╡ 09a785ae-7286-4a47-b103-1e90efa20167 begin - Γ = Inti.boundary(Ω) - Γ_msh = view(msh, Γ) - Q = Inti.Quadrature(Γ_msh; qorder) - nothing #hide + Γ = Inti.boundary(Ω) + Γ_msh = view(msh, Γ) + Q = Inti.Quadrature(Γ_msh; qorder) + nothing #hide end # ╔═╡ 593216db-810e-47db-b6cb-aa6cf53e8590 @@ -231,15 +231,15 @@ With the [`Quadrature`](@ref Inti.Quadrature) constructed, we now can define dis # ╔═╡ 9b37d163-15ca-44af-9e00-7410956fc16c begin - pde = Inti.Helmholtz(; k, dim = 2) - S, D = Inti.single_double_layer(; - pde, - target = Q, - source = Q, - compression = (method = :none,), - correction = (method = :dim,), - ) - nothing #hide + op = Inti.Helmholtz(; k, dim = 2) + S, D = Inti.single_double_layer(; + op, + target = Q, + source = Q, + compression = (method = :none,), + correction = (method = :dim,), + ) + nothing #hide end # ╔═╡ 061b7503-bad0-411e-b126-eea4b89feb03 @@ -255,8 +255,8 @@ We can now combine `S` and `D` to form the combined-field operator: # ╔═╡ ebaf6ed5-8e01-4c5f-9d7f-8cdc32a68022 begin - L = I / 2 + D - im * k * S - nothing #hide + L = I / 2 + D - im * k * S + nothing #hide end # ╔═╡ 8ba9aadc-5b34-41cb-adba-1ebdc65a741d @@ -266,12 +266,12 @@ where `I` is the identity matrix. Assuming an incident field along the $x_1$ dir # ╔═╡ eaf2025e-1135-4ee6-9dea-25e2c0c09ec8 begin - uᵢ = x -> exp(im * k * x[1]) # plane-wave incident field - rhs = map(Q) do q - x = q.coords - return -uᵢ(x) - end - nothing #hide + uᵢ = x -> exp(im * k * x[1]) # plane-wave incident field + rhs = map(Q) do q + x = q.coords + return -uᵢ(x) + end + nothing #hide end # ╔═╡ 7fcd47c1-9d65-41e2-9363-141d53f0dbca @@ -287,8 +287,8 @@ We can now solve the integral equation using e.g. the backslash operator: # ╔═╡ 2925e1a3-1051-4b09-bbd9-461c6a76a1e3 begin - σ = L \ rhs - nothing #hide + σ = L \ rhs + nothing #hide end # ╔═╡ b26e0c64-b3f1-4b18-a6f5-8f789407a744 @@ -298,9 +298,9 @@ The variable `σ` contains the value of the approximate density at the quadratur # ╔═╡ fd0f2cf4-1b18-4c56-84a9-2339b3dfc10a begin - 𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) - uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) - nothing #hide + 𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) + uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) + nothing #hide end # ╔═╡ 12b3c26a-96b3-440b-9771-4a945fe14ce7 @@ -312,25 +312,25 @@ To assess the accuracy of the solution, we can compare it to the exact solution # ╔═╡ 185cfa1e-7c87-4243-8f5e-c5983932e135 begin - function circle_helmholtz_soundsoft(pt; radius = 1, k, θin) - x = pt[1] - y = pt[2] - r = sqrt(x^2 + y^2) - θ = atan(y, x) - u = 0.0 - r < radius && return u - c(n) = -exp(im * n * (π / 2 - θin)) * besselj(n, k * radius) / besselh(n, k * radius) - u = c(0) * besselh(0, k * r) - n = 1 - while (abs(c(n)) > 1e-12) - u += - c(n) * besselh(n, k * r) * exp(im * n * θ) + - c(-n) * besselh(-n, k * r) * exp(-im * n * θ) - n += 1 - end - return u - end - nothing #hide + function circle_helmholtz_soundsoft(pt; radius = 1, k, θin) + x = pt[1] + y = pt[2] + r = sqrt(x^2 + y^2) + θ = atan(y, x) + u = 0.0 + r < radius && return u + c(n) = -exp(im * n * (π / 2 - θin)) * besselj(n, k * radius) / besselh(n, k * radius) + u = c(0) * besselh(0, k * r) + n = 1 + while (abs(c(n)) > 1e-12) + u += + c(n) * besselh(n, k * r) * exp(im * n * θ) + + c(-n) * besselh(-n, k * r) * exp(-im * n * θ) + n += 1 + end + return u + end + nothing #hide end # ╔═╡ 4c37081c-4f68-48fe-be5b-bebf35f0c35e @@ -340,13 +340,13 @@ Here is the maximum error on some points located on a circle of radius `2`: # ╔═╡ 1daa0518-a98b-448c-bddd-1147a0f7ed4d begin - uₑ = x -> circle_helmholtz_soundsoft(x; k, radius = 1, θin = 0) # exact solution - er = maximum(0:0.01:2π) do θ - R = 2 - x = (R * cos(θ), R * sin(θ)) - return abs(uₛ(x) - uₑ(x)) - end - @info "maximum error = $er" + uₑ = x -> circle_helmholtz_soundsoft(x; k, radius = 1, θin = 0) # exact solution + er = maximum(0:0.01:2π) do θ + R = 2 + x = (R * cos(θ), R * sin(θ)) + return abs(uₛ(x) - uₑ(x)) + end + @info "maximum error = $er" end # ╔═╡ 3557cd23-43b2-4c79-a0e9-53c0d9650851 @@ -359,20 +359,22 @@ As we can see, the error is quite small! Let's use `Makie` to visualize the solu # ╔═╡ 56350a37-b2d8-463e-934f-e8a4fbf72c67 begin - xx = yy = range(-4; stop = 4, length = 200) - vals = - map(pt -> Inti.isinside(pt, Q) ? NaN : real(uₛ(pt) + uᵢ(pt)), Iterators.product(xx, yy)) - fig, ax, hm = heatmap( - xx, - yy, - vals; - colormap = :inferno, - interpolate = true, - axis = (aspect = DataAspect(), xgridvisible = false, ygridvisible = false), - ) - viz!(Γ_msh; color = :white, segmentsize = 5) - Colorbar(fig[1, 2], hm) - fig + xx = yy = range(-4; stop = 4, length = 200) + vals = map( + pt -> Inti.isinside(pt, Q) ? NaN : real(uₛ(pt) + uᵢ(pt)), + Iterators.product(xx, yy), + ) + fig, ax, hm = heatmap( + xx, + yy, + vals; + colormap = :inferno, + interpolate = true, + axis = (aspect = DataAspect(), xgridvisible = false, ygridvisible = false), + ) + viz!(Γ_msh; color = :white, segmentsize = 5) + Colorbar(fig[1, 2], hm) + fig end # ╔═╡ 471a5475-1ef5-4f40-bb45-a3e29ae5c0a8 @@ -384,43 +386,47 @@ Before moving on to the 3D example let us simply mention that, besides the fact # ╔═╡ 3f3f4794-f2bf-4171-a1b4-70fd38c01fc3 let - # vertices of an equilateral triangle centered at the origin with a vertex at (0,1) - a, b, c = SVector(0, 1), SVector(sqrt(3) / 2, -1 / 2), SVector(-sqrt(3) / 2, -1 / 2) - circle_f(center, radius) = s -> center + radius * SVector(cospi(2 * s[1]), sinpi(2 * s[1])) - disk1 = Inti.parametric_curve(circle_f(a, 1 / 2), 0, 1) - disk2 = Inti.parametric_curve(circle_f(b, 1 / 2), 0, 1) - disk3 = Inti.parametric_curve(circle_f(c, 1 / 2), 0, 1) - Γ = disk1 ∪ disk2 ∪ disk3 - msh = Inti.meshgen(Γ; meshsize) - Γ_msh = view(msh, Γ) - Q = Inti.Quadrature(Γ_msh; qorder) - S, D = Inti.single_double_layer(; - pde, - target = Q, - source = Q, - compression = (method = :none,), - correction = (method = :dim,), - ) - L = I / 2 + D - im * k * S - rhs = map(q -> -uᵢ(q.coords), Q) - σ = L \ rhs - 𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) - uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) - vals = - map(pt -> Inti.isinside(pt, Q) ? NaN : real(uₛ(pt) + uᵢ(pt)), Iterators.product(xx, yy)) - colorrange = (-2, 2) - fig, ax, hm = heatmap( - xx, - yy, - vals; - colormap = :inferno, - colorrange, - interpolate = true, - axis = (aspect = DataAspect(), xgridvisible = false, ygridvisible = false), - ) - viz!(Γ_msh; color = :black, segmentsize = 4) - Colorbar(fig[1, 2], hm) - fig + # vertices of an equilateral triangle centered at the origin with a vertex at (0,1) + a, b, c = SVector(0, 1), SVector(sqrt(3) / 2, -1 / 2), SVector(-sqrt(3) / 2, -1 / 2) + function circle_f(center, radius) + return s -> center + radius * SVector(cospi(2 * s[1]), sinpi(2 * s[1])) + end + disk1 = Inti.parametric_curve(circle_f(a, 1 / 2), 0, 1) + disk2 = Inti.parametric_curve(circle_f(b, 1 / 2), 0, 1) + disk3 = Inti.parametric_curve(circle_f(c, 1 / 2), 0, 1) + Γ = disk1 ∪ disk2 ∪ disk3 + msh = Inti.meshgen(Γ; meshsize) + Γ_msh = view(msh, Γ) + Q = Inti.Quadrature(Γ_msh; qorder) + S, D = Inti.single_double_layer(; + op, + target = Q, + source = Q, + compression = (method = :none,), + correction = (method = :dim,), + ) + L = I / 2 + D - im * k * S + rhs = map(q -> -uᵢ(q.coords), Q) + σ = L \ rhs + 𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) + uₛ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) + vals = map( + pt -> Inti.isinside(pt, Q) ? NaN : real(uₛ(pt) + uᵢ(pt)), + Iterators.product(xx, yy), + ) + colorrange = (-2, 2) + fig, ax, hm = heatmap( + xx, + yy, + vals; + colormap = :inferno, + colorrange, + interpolate = true, + axis = (aspect = DataAspect(), xgridvisible = false, ygridvisible = false), + ) + viz!(Γ_msh; color = :black, segmentsize = 4) + Colorbar(fig[1, 2], hm) + fig end # ╔═╡ b8fd284a-3865-4762-9602-6bc1eba464a4 @@ -440,28 +446,28 @@ The visualization is also more involved, and we will use the `Gmsh` API to creat # ╔═╡ b02fbb5b-ce97-48e7-ab2b-a8db6ca7c8d3 begin - function gmsh_sphere(; meshsize, order = gorder, radius = 1, visualize = false, name) - gmsh.initialize() - gmsh.model.add("sphere-scattering") - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) - sphere_tag = gmsh.model.occ.addSphere(0, 0, 0, radius) - xl, yl, zl = -2 * radius, -2 * radius, 0 - Δx, Δy = 4 * radius, 4 * radius - rectangle_tag = gmsh.model.occ.addRectangle(xl, yl, zl, Δx, Δy) - outDimTags, _ = - gmsh.model.occ.cut([(2, rectangle_tag)], [(3, sphere_tag)], -1, true, false) - gmsh.model.occ.synchronize() - gmsh.model.addPhysicalGroup(3, [sphere_tag], -1, "omega") - gmsh.model.addPhysicalGroup(2, [dt[2] for dt in outDimTags], -1, "sigma") - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(order) - visualize && gmsh.fltk.run() - gmsh.option.setNumber("Mesh.SaveAll", 1) # otherwise only the physical groups are saved - gmsh.write(name) - return gmsh.finalize() - end - nothing #hide + function gmsh_sphere(; meshsize, order = gorder, radius = 1, visualize = false, name) + gmsh.initialize() + gmsh.model.add("sphere-scattering") + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) + sphere_tag = gmsh.model.occ.addSphere(0, 0, 0, radius) + xl, yl, zl = -2 * radius, -2 * radius, 0 + Δx, Δy = 4 * radius, 4 * radius + rectangle_tag = gmsh.model.occ.addRectangle(xl, yl, zl, Δx, Δy) + outDimTags, _ = + gmsh.model.occ.cut([(2, rectangle_tag)], [(3, sphere_tag)], -1, true, false) + gmsh.model.occ.synchronize() + gmsh.model.addPhysicalGroup(3, [sphere_tag], -1, "omega") + gmsh.model.addPhysicalGroup(2, [dt[2] for dt in outDimTags], -1, "sigma") + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(order) + visualize && gmsh.fltk.run() + gmsh.option.setNumber("Mesh.SaveAll", 1) # otherwise only the physical groups are saved + gmsh.write(name) + return gmsh.finalize() + end + nothing #hide end # ╔═╡ 49f0edcd-9cb8-4722-ae56-b17ee48f8336 @@ -472,9 +478,9 @@ As before, lets write a file with our mesh, and import it into `Inti.jl`: # ╔═╡ a18b8726-fa97-4388-bb8b-53965a5d433b # ╠═╡ show_logs = false begin - name_sphere = joinpath(@__DIR__, "sphere.msh") - gmsh_sphere(; meshsize=(λ / 5), order = gorder, name=name_sphere, visualize = false) - msh_3d = Inti.import_mesh(name_sphere; dim = 3) + name_sphere = joinpath(@__DIR__, "sphere.msh") + gmsh_sphere(; meshsize = (λ / 5), order = gorder, name = name_sphere, visualize = false) + msh_3d = Inti.import_mesh(name_sphere; dim = 3) end # ╔═╡ 063e3185-8714-4955-b985-5413420a889b @@ -490,10 +496,10 @@ Since we created physical groups in `Gmsh`, we can use them to extract the relev # ╔═╡ 631328d5-c961-41e1-ac0e-181c3bb75112 begin - Ω_3d = Inti.Domain(e -> "omega" ∈ Inti.labels(e), Inti.entities(msh_3d)) - Σ_3d = Inti.Domain(e -> "sigma" ∈ Inti.labels(e), Inti.entities(msh_3d)) - Γ_3d = Inti.boundary(Ω_3d) - nothing #hide + Ω_3d = Inti.Domain(e -> "omega" ∈ Inti.labels(e), Inti.entities(msh_3d)) + Σ_3d = Inti.Domain(e -> "sigma" ∈ Inti.labels(e), Inti.entities(msh_3d)) + Γ_3d = Inti.boundary(Ω_3d) + nothing #hide end # ╔═╡ dbd6d298-0fbf-420f-9622-d1db1d550cc5 @@ -503,23 +509,23 @@ We can now create a quadrature as before # ╔═╡ 6096328b-e1af-4ec5-b72b-cd381c166cb7 begin - Γ_msh_3d = view(msh_3d, Γ_3d) - Q_3d = Inti.Quadrature(Γ_msh_3d; qorder) - nothing #hide + Γ_msh_3d = view(msh_3d, Γ_3d) + Q_3d = Inti.Quadrature(Γ_msh_3d; qorder) + nothing #hide end # ╔═╡ 498fa575-5cf7-4a7e-8d3b-c50d605aa57c begin - using HMatrices - pde_3d = Inti.Helmholtz(; k, dim = 3) - S_3d, D_3d = Inti.single_double_layer(; - pde = pde_3d, - target = Q_3d, - source = Q_3d, - compression = (method = :hmatrix, tol = 1e-4), - correction = (method = :dim,), - ) - nothing #hide + using HMatrices + op_3d = Inti.Helmholtz(; k, dim = 3) + S_3d, D_3d = Inti.single_double_layer(; + op = op_3d, + target = Q_3d, + source = Q_3d, + compression = (method = :hmatrix, tol = 1e-4), + correction = (method = :dim,), + ) + nothing #hide end # ╔═╡ b6cad085-eb74-4152-8850-226ed837febd @@ -540,8 +546,8 @@ Here is how much memory it would take to store the dense representation of these # ╔═╡ 54f07c21-339b-44c3-9ac2-22af163d55ae begin - mem = 2 * length(S_3d) * 16 / 1e9 # 16 bytes per complex number, 1e9 bytes per GB, two matrices - println("memory required to store S and D: $(mem) GB") + mem = 2 * length(S_3d) * 16 / 1e9 # 16 bytes per complex number, 1e9 bytes per GB, two matrices + println("memory required to store S and D: $(mem) GB") end # ╔═╡ 915d3636-70bd-4001-857f-3aae11eb2d2e @@ -562,8 +568,8 @@ We will use the generalized minimal residual (GMRES) iterative solver, for the l # ╔═╡ a2ddce1e-440e-49e2-8c3e-ec6960ff83f0 begin - L_3d = I / 2 + LinearMap(D_3d) - im * k * LinearMap(S_3d) - nothing #hide + L_3d = I / 2 + LinearMap(D_3d) - im * k * LinearMap(S_3d) + nothing #hide end # ╔═╡ f148c58a-50d3-4ac2-a1da-148563d533e7 @@ -574,13 +580,20 @@ We can now solve the linear system using GMRES solver: # ╔═╡ dcf006a3-5332-442a-9771-f427cd7189a5 # ╠═╡ show_logs = false begin - rhs_3d = map(Q_3d) do q - x = q.coords - return -uᵢ(x) - end - σ_3d, hist = - gmres(L_3d, rhs_3d; log = true, abstol = 1e-6, verbose = false, restart = 100, maxiter = 100) - @show hist + rhs_3d = map(Q_3d) do q + x = q.coords + return -uᵢ(x) + end + σ_3d, hist = gmres( + L_3d, + rhs_3d; + log = true, + abstol = 1e-6, + verbose = false, + restart = 100, + maxiter = 100, + ) + @show hist end # ╔═╡ 215677a6-af7b-4e5d-863f-aabdfb4a1400 @@ -590,9 +603,9 @@ As before, let us represent the solution using `IntegralPotential`s: # ╔═╡ 14a549ad-e23f-40a4-b7fa-ef31354805a9 begin - 𝒮_3d, 𝒟_3d = Inti.single_double_layer_potential(; pde=pde_3d, source = Q_3d) - uₛ_3d = x -> 𝒟_3d[σ_3d](x) - im * k * 𝒮_3d[σ_3d](x) - nothing #hide + 𝒮_3d, 𝒟_3d = Inti.single_double_layer_potential(; op = op_3d, source = Q_3d) + uₛ_3d = x -> 𝒟_3d[σ_3d](x) - im * k * 𝒮_3d[σ_3d](x) + nothing #hide end # ╔═╡ 832858dc-2475-4f6a-a820-a557b2883f0c @@ -602,31 +615,32 @@ To check the result, we compare against the exact solution obtained through a se # ╔═╡ e00c4d90-aa3b-41f9-b65e-7e9ace89a295 begin - sphbesselj(l, r) = sqrt(π / (2r)) * besselj(l + 1 / 2, r) - sphbesselh(l, r) = sqrt(π / (2r)) * besselh(l + 1 / 2, r) - sphharmonic(l, m, θ, ϕ) = GSL.sf_legendre_sphPlm(l, abs(m), cos(θ)) * exp(im * m * ϕ) - function sphere_helmholtz_soundsoft(xobs; radius = 1, k = 1, θin = 0, ϕin = 0) - x = xobs[1] - y = xobs[2] - z = xobs[3] - r = sqrt(x^2 + y^2 + z^2) - θ = acos(z / r) - ϕ = atan(y, x) - u = 0.0 - r < radius && return u - function c(l, m) - return -4π * im^l * sphharmonic(l, -m, θin, ϕin) * sphbesselj(l, k * radius) / sphbesselh(l, k * radius) - end - l = 0 - for l in 0:60 - for m in -l:l - u += c(l, m) * sphbesselh(l, k * r) * sphharmonic(l, m, θ, ϕ) - end - l += 1 - end - return u - end - nothing #hide + sphbesselj(l, r) = sqrt(π / (2r)) * besselj(l + 1 / 2, r) + sphbesselh(l, r) = sqrt(π / (2r)) * besselh(l + 1 / 2, r) + sphharmonic(l, m, θ, ϕ) = GSL.sf_legendre_sphPlm(l, abs(m), cos(θ)) * exp(im * m * ϕ) + function sphere_helmholtz_soundsoft(xobs; radius = 1, k = 1, θin = 0, ϕin = 0) + x = xobs[1] + y = xobs[2] + z = xobs[3] + r = sqrt(x^2 + y^2 + z^2) + θ = acos(z / r) + ϕ = atan(y, x) + u = 0.0 + r < radius && return u + function c(l, m) + return -4π * im^l * sphharmonic(l, -m, θin, ϕin) * sphbesselj(l, k * radius) / + sphbesselh(l, k * radius) + end + l = 0 + for l in 0:60 + for m in -l:l + u += c(l, m) * sphbesselh(l, k * r) * sphharmonic(l, m, θ, ϕ) + end + l += 1 + end + return u + end + nothing #hide end # ╔═╡ bb87f66d-5e68-49c5-9167-bdbab60109e5 @@ -636,13 +650,13 @@ We will compute the error on some point on the sphere of radius `2`: # ╔═╡ 3848717e-0c30-4c0b-82ec-b504c9210483 begin - uₑ_3d = (x) -> sphere_helmholtz_soundsoft(x; radius = 1, k = k, θin = π / 2, ϕin = 0) - er_3d = maximum(1:100) do _ - x̂ = rand(Inti.Point3D) |> normalize # an SVector of unit norm - x = 2 * x̂ - return abs(uₛ_3d(x) - uₑ_3d(x)) - end - @info "error with correction = $er_3d" + uₑ_3d = (x) -> sphere_helmholtz_soundsoft(x; radius = 1, k = k, θin = π / 2, ϕin = 0) + er_3d = maximum(1:100) do _ + x̂ = rand(Inti.Point3D) |> normalize # an SVector of unit norm + x = 2 * x̂ + return abs(uₛ_3d(x) - uₑ_3d(x)) + end + @info "error with correction = $er_3d" end # ╔═╡ 36ddb3be-7103-429f-9011-99c27c655de5 @@ -655,22 +669,22 @@ We see that, once again, the approximation is quite accurate. Let us now visuali # ╔═╡ b2b31b38-4c2f-4763-9273-971749c40d5c begin - Σ_msh = view(msh_3d, Σ_3d) - target = Inti.nodes(Σ_msh) - - S_viz, D_viz = Inti.single_double_layer(; - pde = pde_3d, - target, - source = Q_3d, - compression = (method = :hmatrix, tol = 1e-4), - # correction for the nearfield (for visual purposes, set to `:none` to disable) - correction = (method = :dim, maxdist = meshsize, target_location = :outside), - ) - - ui_eval_msh = uᵢ.(target) - us_eval_msh = D_viz * σ_3d - im * k * S_viz * σ_3d - u_eval_msh = ui_eval_msh + us_eval_msh - nothing #hide + Σ_msh = view(msh_3d, Σ_3d) + target = Inti.nodes(Σ_msh) + + S_viz, D_viz = Inti.single_double_layer(; + op = op_3d, + target, + source = Q_3d, + compression = (method = :hmatrix, tol = 1e-4), + # correction for the nearfield (for visual purposes, set to `:none` to disable) + correction = (method = :dim, maxdist = meshsize, target_location = :outside), + ) + + ui_eval_msh = uᵢ.(target) + us_eval_msh = D_viz * σ_3d - im * k * S_viz * σ_3d + u_eval_msh = ui_eval_msh + us_eval_msh + nothing #hide end # ╔═╡ 0cec09b1-8806-4aee-970d-cc0b1e4b841c @@ -680,15 +694,15 @@ Finalize, we use [`Meshes.viz`](@extref) to visualize the scattered field: # ╔═╡ dab9a19c-977e-4ab1-a3d9-f9970af69642 begin - nv = length(Inti.nodes(Γ_msh_3d)) - colorrange = extrema(real(u_eval_msh)) - colormap = :inferno - fig_3d = Figure(; size = (800, 500)) - ax_3d = Axis3(fig_3d[1, 1]; aspect = :data) - viz!(Γ_msh_3d; colorrange, colormap, color = zeros(nv), interpolate = true) - viz!(Σ_msh; colorrange, colormap, color = real(u_eval_msh)) - cb = Colorbar(fig_3d[1, 2]; label = "real(u)", colormap, colorrange) - fig_3d #hide + nv = length(Inti.nodes(Γ_msh_3d)) + colorrange = extrema(real(u_eval_msh)) + colormap = :inferno + fig_3d = Figure(; size = (800, 500)) + ax_3d = Axis3(fig_3d[1, 1]; aspect = :data) + viz!(Γ_msh_3d; colorrange, colormap, color = zeros(nv), interpolate = true) + viz!(Σ_msh; colorrange, colormap, color = real(u_eval_msh)) + cb = Colorbar(fig_3d[1, 2]; label = "real(u)", colormap, colorrange) + fig_3d #hide end # ╔═╡ 765e46fc-dac3-439a-ab07-2865298aed45 diff --git a/docs/src/pluto-examples/poisson.jl b/docs/src/pluto-examples/poisson.jl index ece579b9..e1d3222e 100644 --- a/docs/src/pluto-examples/poisson.jl +++ b/docs/src/pluto-examples/poisson.jl @@ -6,28 +6,28 @@ using InteractiveUtils # ╔═╡ 4a365b8b-fcb0-4507-883f-c49c252c4f3d begin - import Pkg as _Pkg + import Pkg as _Pkg haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"]) - using PlutoUI: TableOfContents -end ; + using PlutoUI: TableOfContents +end; # ╔═╡ 332bd3bb-c720-454e-80af-89ad65041773 begin - using Inti, Gmsh - meshsize = 0.1 - gmsh.initialize() - jellyfish = Inti.gmsh_curve(0, 2π; meshsize) do s - r = 1 + 0.3*cos(4*s + 2*sin(s)) - return r*Inti.Point2D(cos(s), sin(s)) - end - cl = gmsh.model.occ.addCurveLoop([jellyfish]) - surf = gmsh.model.occ.addPlaneSurface([cl]) - gmsh.model.occ.synchronize() - gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(2) - msh = Inti.import_mesh(; dim = 2) - gmsh.finalize() + using Inti, Gmsh + meshsize = 0.1 + gmsh.initialize() + jellyfish = Inti.gmsh_curve(0, 2π; meshsize) do s + r = 1 + 0.3 * cos(4 * s + 2 * sin(s)) + return r * Inti.Point2D(cos(s), sin(s)) + end + cl = gmsh.model.occ.addCurveLoop([jellyfish]) + surf = gmsh.model.occ.addPlaneSurface([cl]) + gmsh.model.occ.synchronize() + gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + gmsh.model.mesh.generate(2) + gmsh.model.mesh.setOrder(2) + msh = Inti.import_mesh(; dim = 2) + gmsh.finalize() end # ╔═╡ aca57c74-d084-40e7-9760-edf988a64915 @@ -117,19 +117,19 @@ We can now extract components of the mesh corresponding to the ``\Omega`` and # ╔═╡ 4e580f9f-6e22-4160-b262-ca941b6bfb8f begin - Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, msh) - Γ = Inti.boundary(Ω) - Ω_msh = view(msh, Ω) - Γ_msh = view(msh, Γ) - nothing #hide + Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, msh) + Γ = Inti.boundary(Ω) + Ω_msh = view(msh, Ω) + Γ_msh = view(msh, Γ) + nothing #hide end # ╔═╡ 2c4bace3-ef35-4500-829f-2f5ae6725249 begin - using Meshes, GLMakie - viz(Ω_msh; showsegments=true) - viz!(Γ_msh; color=:red) - Makie.current_figure() #hide + using Meshes, GLMakie + viz(Ω_msh; showsegments = true) + viz!(Γ_msh; color = :red) + Makie.current_figure() #hide end # ╔═╡ 901578c3-b7aa-4023-bb99-769cf5805f57 @@ -145,27 +145,23 @@ boundary: # ╔═╡ 85ce61fc-086d-4661-8f03-6a6ae4d55511 begin - Ω_quad = Inti.Quadrature(Ω_msh; qorder = 4) - Γ_quad = Inti.Quadrature(Γ_msh; qorder = 6) - nothing #hide + Ω_quad = Inti.Quadrature(Ω_msh; qorder = 4) + Γ_quad = Inti.Quadrature(Γ_msh; qorder = 6) + nothing #hide end # ╔═╡ 6e3ea607-2b70-485d-94f0-5626bab4f832 begin - using FMM2D #to accelerate the maps - pde = Inti.Laplace(; dim = 2) - # Newtonian potential mapping domain to boundary - V_d2b = Inti.volume_potential(; - pde, - target = Γ_quad, - source = Ω_quad, - compression = (method = :fmm, tol = 1e-12), - correction = ( - method = :dim, - maxdist = 5 * meshsize, - target_location = :on, - ), - ) + using FMM2D #to accelerate the maps + op = Inti.Laplace(; dim = 2) + # Newtonian potential mapping domain to boundary + V_d2b = Inti.volume_potential(; + op, + target = Γ_quad, + source = Ω_quad, + compression = (method = :fmm, tol = 1e-12), + correction = (method = :dim, maxdist = 5 * meshsize, target_location = :on), + ) end # ╔═╡ 731d9ae8-33a8-4f03-8cbc-5e40972996a2 @@ -186,7 +182,7 @@ equation: # ╔═╡ 878cfa05-066d-43df-b93a-2b16565f8b4e # Single and double layer operators on Γ S_b2b, D_b2b = Inti.single_double_layer(; - pde, + op, target = Γ_quad, source = Γ_quad, compression = (method = :fmm, tol = 1e-12), @@ -214,12 +210,12 @@ instead a manufactured solution ``u_e`` from which we will derive the functions # ╔═╡ 975b7c01-8147-44f8-a693-1185e7b8d63b begin - # Create a manufactured solution - uₑ = (x) -> cos(2 * x[1]) * sin(2 * x[2]) - fₑ = (x) -> 8 * cos(2 * x[1]) * sin(2 * x[2]) # -Δuₑ - g = map(q -> uₑ(q.coords), Γ_quad) - f = map(q -> fₑ(q.coords), Ω_quad) - nothing #hide + # Create a manufactured solution + uₑ = (x) -> cos(2 * x[1]) * sin(2 * x[2]) + fₑ = (x) -> 8 * cos(2 * x[1]) * sin(2 * x[2]) # -Δuₑ + g = map(q -> uₑ(q.coords), Γ_quad) + f = map(q -> fₑ(q.coords), Ω_quad) + nothing #hide end # ╔═╡ 700cfbbc-970a-466b-9f8c-748f7ff0bc6e @@ -230,15 +226,15 @@ homogeneous part of the solution: # ╔═╡ 6db682e4-e641-4272-ba90-1f10b2ff1150 begin - rhs = g - V_d2b*f - nothing #hide + rhs = g - V_d2b * f + nothing #hide end # ╔═╡ 6eb1d813-e792-4148-8d30-975c49e9dbc6 begin - using IterativeSolvers, LinearAlgebra - σ = gmres(-I/2 + D_b2b, rhs; abstol = 1e-8, verbose = true, restart = 1000) - nothing #hide + using IterativeSolvers, LinearAlgebra + σ = gmres(-I / 2 + D_b2b, rhs; abstol = 1e-8, verbose = true, restart = 1000) + nothing #hide end # ╔═╡ b21911c7-7276-44cb-a15f-64db3430a896 @@ -253,11 +249,11 @@ With the density function at hand, we can now reconstruct our approximate soluti # ╔═╡ e0e1fa9c-7a43-45b1-ad45-2510533e1aed begin - G = Inti.SingleLayerKernel(pde) - dG = Inti.DoubleLayerKernel(pde) - 𝒱 = Inti.IntegralPotential(G, Ω_quad) - 𝒟 = Inti.IntegralPotential(dG, Γ_quad) - u = (x) -> 𝒱[f](x) + 𝒟[σ](x) + G = Inti.SingleLayerKernel(op) + dG = Inti.DoubleLayerKernel(op) + 𝒱 = Inti.IntegralPotential(G, Ω_quad) + 𝒟 = Inti.IntegralPotential(dG, Γ_quad) + u = (x) -> 𝒱[f](x) + 𝒟[σ](x) end # ╔═╡ 2c2c943b-30ed-4244-a40c-f061050ca7b8 @@ -267,8 +263,8 @@ and evaluate it at any point in the domain: # ╔═╡ 2faa4311-59d9-4b85-b585-937391e94568 begin - x = Inti.Point2D(0.1,0.4) - println("error at $x: ", u(x)-uₑ(x)) + x = Inti.Point2D(0.1, 0.4) + println("error at $x: ", u(x) - uₑ(x)) end # ╔═╡ f3ee21f2-99e7-4be3-a2cc-930b6c4487f1 @@ -290,11 +286,11 @@ solution ``u`` at all the quadrature nodes of ``\Omega``: # ╔═╡ fccb6992-32c4-4bd2-a742-fb36906fb62a V_d2d = Inti.volume_potential(; - pde, + op, target = Ω_quad, source = Ω_quad, compression = (method = :fmm, tol = 1e-8), - correction = (method = :dim, ), + correction = (method = :dim,), ) # ╔═╡ 6489d40d-fdd3-488d-ae5d-4e4e489b669f @@ -305,11 +301,11 @@ our mesh nodes: # ╔═╡ 66a7946c-a601-4fc8-94a0-429d41b564ac S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ω_quad, source = Γ_quad, compression = (method = :fmm, tol = 1e-8), - correction = (method = :dim, maxdist = 2*meshsize, target_location = :inside), + correction = (method = :dim, maxdist = 2 * meshsize, target_location = :inside), ) # ╔═╡ a5c7fcd0-1966-49fe-b5bf-7d3c1d9f4aa7 @@ -320,10 +316,10 @@ manufactured: # ╔═╡ 5be59f68-cf21-4e7e-ac23-7bffd690dc03 begin - u_quad = V_d2d*f + D_b2d*σ - er_quad = u_quad - map(q -> uₑ(q.coords), Ω_quad) - println("maximum error at all quadrature nodes: ", norm(er_quad, Inf)) - nothing #hide + u_quad = V_d2d * f + D_b2d * σ + er_quad = u_quad - map(q -> uₑ(q.coords), Ω_quad) + println("maximum error at all quadrature nodes: ", norm(er_quad, Inf)) + nothing #hide end # ╔═╡ a61f336e-15ee-4bb3-a071-1115ac6f1be1 @@ -333,22 +329,22 @@ Lastly, let us visualize the solution and the error on the mesh nodes using [`qu # ╔═╡ ce03f9b7-c91c-4f53-bcda-49261a4bdcc2 begin - nodes = Inti.nodes(Ω_msh) - u_nodes = Inti.quadrature_to_node_vals(Ω_quad, u_quad) - er = u_nodes - map(uₑ, nodes) - colorrange = extrema(u_nodes) - fig = Figure(; size = (800, 300)) - ax = Axis(fig[1, 1]; aspect = DataAspect()) - viz!(Ω_msh; colorrange, color = u_nodes, interpolate = true) - cb = Colorbar(fig[1, 2]; label = "u", colorrange) - # plot error - log_er = log10.(abs.(er)) - colorrange = extrema(log_er) - colormap = :inferno - ax = Axis(fig[1, 3]; aspect = DataAspect()) - viz!(Ω_msh; colorrange, colormap, color = log_er, interpolate = true) - cb = Colorbar(fig[1, 4]; label = "log₁₀|u - uₑ|", colormap, colorrange) - fig #hide + nodes = Inti.nodes(Ω_msh) + u_nodes = Inti.quadrature_to_node_vals(Ω_quad, u_quad) + er = u_nodes - map(uₑ, nodes) + colorrange = extrema(u_nodes) + fig = Figure(; size = (800, 300)) + ax = Axis(fig[1, 1]; aspect = DataAspect()) + viz!(Ω_msh; colorrange, color = u_nodes, interpolate = true) + cb = Colorbar(fig[1, 2]; label = "u", colorrange) + # plot error + log_er = log10.(abs.(er)) + colorrange = extrema(log_er) + colormap = :inferno + ax = Axis(fig[1, 3]; aspect = DataAspect()) + viz!(Ω_msh; colorrange, colormap, color = log_er, interpolate = true) + cb = Colorbar(fig[1, 4]; label = "log₁₀|u - uₑ|", colormap, colorrange) + fig #hide end # ╔═╡ e63c776b-6b52-4e6e-aca8-3d67cd9a9c3f diff --git a/docs/src/tutorials/compression_methods.md b/docs/src/tutorials/compression_methods.md index e74c4d19..e1f802c7 100644 --- a/docs/src/tutorials/compression_methods.md +++ b/docs/src/tutorials/compression_methods.md @@ -46,8 +46,8 @@ geo = Inti.GeometricEntity("ellipsoid") Γ = Inti.boundary(Ω) Q = Inti.Quadrature(Γ; meshsize = 0.4, qorder = 5) # create the operator -pde = Inti.Helmholtz(; dim = 3, k = 2π) -K = Inti.SingleLayerKernel(pde) +op = Inti.Helmholtz(; dim = 3, k = 2π) +K = Inti.SingleLayerKernel(op) Sop = Inti.IntegralOperator(K, Q, Q) x = rand(eltype(Sop), length(Q)) rtol = 1e-8 diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md index 01da5b86..0c19a3b3 100644 --- a/docs/src/tutorials/getting_started.md +++ b/docs/src/tutorials/getting_started.md @@ -42,7 +42,7 @@ using Inti Inti.stack_weakdeps_env!() # add weak dependencies # PDE k = 2π -pde = Inti.Helmholtz(; dim = 2, k) +op = Inti.Helmholtz(; dim = 2, k) ``` Next, we generate the geometry of the problem. For this tutorial, we will @@ -138,7 +138,7 @@ To approximate ``S`` and ``D``, we can proceed as follows: ```@example getting_started S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), @@ -221,7 +221,7 @@ potentials defined as: to compute the solution ``u`` in the domain: ```@example getting_started -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) uₛ = x -> 𝒟[u](x) - 𝒮[g](x) ``` @@ -253,8 +253,8 @@ to the solution obtained numerically, as illustrated below: ```@example getting_started # build an exact solution -G = Inti.SingleLayerKernel(pde) -dG = Inti.DoubleLayerKernel(pde) +G = Inti.SingleLayerKernel(op) +dG = Inti.DoubleLayerKernel(op) xs = map(θ -> 0.5 * rand() * SVector(cos(θ), sin(θ)), 2π * rand(10)) cs = rand(ComplexF64, length(xs)) uₑ = q -> sum(c * G(x, q) for (x, c) in zip(xs, cs)) diff --git a/docs/src/tutorials/integral_operators.md b/docs/src/tutorials/integral_operators.md index 44545679..1dd6b63f 100644 --- a/docs/src/tutorials/integral_operators.md +++ b/docs/src/tutorials/integral_operators.md @@ -27,8 +27,8 @@ be used to construct the corresponding kernel functions, e.g.: ```@example integral_operators using Inti, StaticArrays, LinearAlgebra -pde = Inti.Helmholtz(; dim = 2, k = 2π) -G = Inti.SingleLayerKernel(pde) +op = Inti.Helmholtz(; dim = 2, k = 2π) +G = Inti.SingleLayerKernel(op) ``` Typically, we are not interested in the kernels themselves, but in the integral @@ -40,14 +40,14 @@ construct the four integral operators of Calderón calculus: Γ = Inti.parametric_curve(s -> SVector(cos(s), sin(s)), 0, 2π) |> Inti.Domain Q = Inti.Quadrature(Γ; meshsize = 0.1, qorder = 5) S, D = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), correction = (method = :dim,) ) K, N = Inti.adj_double_layer_hypersingular(; - pde, + op, target = Q, source = Q, compression = (method = :none,), @@ -85,14 +85,14 @@ something different: ```@example integral_operators using FMM2D # will load the extension Sfmm, Dfmm = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :fmm, tol = 1e-10), correction = (method = :dim, ) ) Kfmm, Nfmm = Inti.adj_double_layer_hypersingular(; - pde, + op, target = Q, source = Q, compression = (method = :fmm, tol = 1e-10), diff --git a/docs/src/tutorials/layer_potentials.md b/docs/src/tutorials/layer_potentials.md index fd6a4a64..077dec29 100644 --- a/docs/src/tutorials/layer_potentials.md +++ b/docs/src/tutorials/layer_potentials.md @@ -66,8 +66,8 @@ it is often more convenient to use the `single_layer_potential` when working with a supported PDE, e.g.: ```@example layer_potentials -pde = Inti.Laplace(; dim = 2) -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +op = Inti.Laplace(; dim = 2) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) ``` creates the single and double layer potentials for the Laplace equation in 2D. @@ -81,7 +81,7 @@ created through the Gmsh API. Do to so, let us first define the PDE:g using Inti, StaticArrays, LinearAlgebra, Meshes, GLMakie, Gmsh # define the PDE k = 4π -pde = Inti.Helmholtz(; dim = 2, k) +op = Inti.Helmholtz(; dim = 2, k) ``` We will now use the [`gmsh_curve`](@ref) function to create a smooth domain of a @@ -152,7 +152,7 @@ on a quadrature of ``\Gamma``: Γ = Inti.boundary(Ω) Q = Inti.Quadrature(view(msh,Γ); qorder = 5) # evaluate the layer potentials -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) +𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q) γ₀u = map(q -> u(q.coords), Q) γ₁u = map(q -> du(q.coords, q.normal), Q) uₕ = x -> 𝒮[γ₁u](x) - 𝒟[γ₀u](x) @@ -193,7 +193,7 @@ acceleration with a near-field correction to evaluate the layer potentials:: ```@example layer_potentials using FMM2D -S, D = Inti.single_double_layer(; pde, target, source = Q, +S, D = Inti.single_double_layer(; op, target, source = Q, compression = (method = :fmm, tol = 1e-12), correction = (method = :dim, target_location = :inside, maxdist = 0.2) ) diff --git a/ext/IntiFMM2DExt.jl b/ext/IntiFMM2DExt.jl index c0d74604..896a7d26 100644 --- a/ext/IntiFMM2DExt.jl +++ b/ext/IntiFMM2DExt.jl @@ -149,7 +149,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) # Helmholtz elseif K isa Inti.SingleLayerKernel{ComplexF64,<:Inti.Helmholtz{2}} charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant @. charges = weights * x @@ -181,7 +181,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) end dipvecs = similar(normals, Float64) dipstrs = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n @@ -219,7 +219,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) xnormals[:, j] = Inti.normal(iop.target[j]) end charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights @. charges = x * weights @@ -255,7 +255,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) end dipvecs = similar(ynormals, Float64) dipstrs = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n diff --git a/ext/IntiFMMLIB2DExt.jl b/ext/IntiFMMLIB2DExt.jl index 15b1860c..be83300d 100644 --- a/ext/IntiFMMLIB2DExt.jl +++ b/ext/IntiFMMLIB2DExt.jl @@ -153,7 +153,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) # Helmholtz elseif K isa Inti.SingleLayerKernel{ComplexF64,<:Inti.Helmholtz{2}} charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant @. charges = weights * x @@ -184,7 +184,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) end dipvecs = similar(normals, Float64) dipstrs = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n @@ -221,7 +221,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) xnormals[:, j] = Inti.normal(iop.target[j]) end charges = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights @. charges = x * weights @@ -258,7 +258,7 @@ function Inti._assemble_fmm2d(iop::Inti.IntegralOperator; rtol = sqrt(eps())) end dipvecs = similar(ynormals, Float64) dipstrs = Vector{ComplexF64}(undef, n) - zk = ComplexF64(K.pde.k) + zk = ComplexF64(K.op.k) return LinearMaps.LinearMap{ComplexF64}(m, n) do y, x # multiply by weights and constant for j in 1:n diff --git a/src/api.jl b/src/api.jl index d11b2453..e17c6921 100644 --- a/src/api.jl +++ b/src/api.jl @@ -14,11 +14,11 @@ Available correction methods for the singular and nearly-singular integrals in const CORRECTION_METHODS = [:none, :dim, :adaptive] """ - single_double_layer(; pde, target, source::Quadrature, compression, + single_double_layer(; op, target, source::Quadrature, compression, correction, derivative = false) Construct a discrete approximation to the single- and double-layer integral -operators for `pde`, mapping values defined on the quadrature nodes of `source` +operators for `op`, mapping values defined on the quadrature nodes of `source` to values defined on the nodes of `target`. If `derivative = true`, return instead the adjoint double-layer and hypersingular operators (which are the derivative of the single- and double-layer, respectively). @@ -57,7 +57,7 @@ integrals should be computed. The available options are: `target_location` is not needed. """ function single_double_layer(; - pde, + op, target, source, compression, @@ -66,8 +66,8 @@ function single_double_layer(; ) compression = _normalize_compression(compression, target, source) correction = _normalize_correction(correction, target, source) - G = derivative ? AdjointDoubleLayerKernel(pde) : SingleLayerKernel(pde) - dG = derivative ? HyperSingularKernel(pde) : DoubleLayerKernel(pde) + G = derivative ? AdjointDoubleLayerKernel(op) : SingleLayerKernel(op) + dG = derivative ? HyperSingularKernel(op) : DoubleLayerKernel(op) Sop = IntegralOperator(G, target, source) Dop = IntegralOperator(dG, target, source) # handle compression @@ -125,7 +125,7 @@ function single_double_layer(; glob_loc_near_trgs = glob_loc_near_trgs, ) δS, δD = bdim_correction( - pde, + op, target[glob_near_trgs], source, Sop_dim_mat, @@ -137,7 +137,7 @@ function single_double_layer(; ) else δS, δD = bdim_correction( - pde, + op, target, source, Smat, @@ -182,7 +182,7 @@ function single_double_layer(; end """ - adj_double_layer_hypersingular(; pde, target, source, compression, + adj_double_layer_hypersingular(; op, target, source, compression, correction) Similar to `single_double_layer`, but for the adjoint double-layer and @@ -190,14 +190,14 @@ hypersingular operators. See the documentation of [`single_double_layer`] for a description of the arguments. """ function adj_double_layer_hypersingular(; - pde, + op, target, source = target, compression, correction, ) return single_double_layer(; - pde, + op, target, source, compression, @@ -207,26 +207,26 @@ function adj_double_layer_hypersingular(; end """ - single_double_layer_potential(; pde, source) + single_double_layer_potential(; op, source) -Return the single- and double-layer potentials for `pde` as +Return the single- and double-layer potentials for `op` as [`IntegralPotential`](@ref)s. """ -function single_double_layer_potential(; pde, source) - G = SingleLayerKernel(pde) - dG = DoubleLayerKernel(pde) +function single_double_layer_potential(; op, source) + G = SingleLayerKernel(op) + dG = DoubleLayerKernel(op) 𝒮 = IntegralPotential(G, source) 𝒟 = IntegralPotential(dG, source) return 𝒮, 𝒟 end """ - volume_potential(; pde, target, source::Quadrature, compression, correction) + volume_potential(; op, target, source::Quadrature, compression, correction) Compute the volume potential operator for a given PDE. ## Arguments -- `pde`: The PDE (Partial Differential Equation) to solve. +- `op`: The PDE (Partial Differential Equation) to solve. - `target`: The target domain where the potential is computed. - `source`: The source domain where the potential is generated. - `compression`: The compression method to use for the potential operator. @@ -274,10 +274,10 @@ the specified compression method. If no compression is specified, the operator is returned as is. If a correction method is specified, the correction is computed and added to the compressed operator. """ -function volume_potential(; pde, target, source::Quadrature, compression, correction) +function volume_potential(; op, target, source::Quadrature, compression, correction) correction = _normalize_correction(correction, target, source) compression = _normalize_compression(compression, target, source) - G = SingleLayerKernel(pde) + G = SingleLayerKernel(op) V = IntegralOperator(G, target, source) # compress V if compression.method == :none @@ -312,7 +312,7 @@ function volume_potential(; pde, target, source::Quadrature, compression, correc # Advanced usage: Use previously constructed layer operators for VDIM if !haskey(correction, :S_b2d) || !haskey(correction, :D_b2d) S, D = single_double_layer(; - pde, + op, target, source = boundary, compression, @@ -324,7 +324,7 @@ function volume_potential(; pde, target, source::Quadrature, compression, correc end interpolation_order = correction.interpolation_order δV = vdim_correction( - pde, + op, target, source, boundary, diff --git a/src/bdim.jl b/src/bdim.jl index edffa331..fac21c40 100644 --- a/src/bdim.jl +++ b/src/bdim.jl @@ -10,9 +10,9 @@ Parameters associated with the density interpolation method used in end """ - bdim_correction(pde,X,Y,S,D; green_multiplier, kwargs...) + bdim_correction(op,X,Y,S,D; green_multiplier, kwargs...) -Given a `pde` and a (possibly innacurate) discretizations of its single and +Given a `op` and a (possibly innacurate) discretizations of its single and double-layer operators `S` and `D` (taking a vector of values on `Y` and returning a vector on of values on `X`), compute corrections `δS` and `δD` such that `S + δS` and `D + δD` are more accurate approximations of the underlying @@ -24,11 +24,11 @@ See [faria2021general](@cite) for more details on the method. ## Required: -- `pde` must be an [`AbstractDifferentialOperator`](@ref) +- `op` must be an [`AbstractDifferentialOperator`](@ref) - `Y` must be a [`Quadrature`](@ref) object of a closed surface - `X` is either inside, outside, or on `Y` - `S` and `D` are approximations to the single- and double-layer operators for - `pde` taking densities in `Y` and returning densities in `X`. + `op` taking densities in `Y` and returning densities in `X`. - `green_multiplier` (keyword argument) is a vector with the same length as `X` storing the value of `μ(x)` for `x ∈ X` in the Green identity `S\\[γ₁u\\](x) - D\\[γ₀u\\](x) + μ*u(x) = 0`. See [`_green_multiplier`](@ref). @@ -48,7 +48,7 @@ See [faria2021general](@cite) for more details on the method. """ function bdim_correction( - pde, + op, target, source::Quadrature, Sop, @@ -91,8 +91,8 @@ function bdim_correction( error("only 2D and 3D supported") end # compute traces of monopoles on the source mesh - G = SingleLayerKernel(pde, T) - γ₁G = AdjointDoubleLayerKernel(pde, T) + G = SingleLayerKernel(op, T) + γ₁G = AdjointDoubleLayerKernel(op, T) γ₀B = Dense{T}(undef, length(source), ns) γ₁B = Dense{T}(undef, length(source), ns) for k in 1:ns diff --git a/src/nystrom.jl b/src/nystrom.jl index 34ba9771..34accef6 100644 --- a/src/nystrom.jl +++ b/src/nystrom.jl @@ -109,7 +109,7 @@ end assemble_fmm(iop; atol) Set up a 2D or 3D FMM for evaluating the discretized integral operator `iop` -associated with the `pde`. In 2D the `FMM2D` or `FMMLIB2D` library is used +associated with the `op`. In 2D the `FMM2D` or `FMMLIB2D` library is used (whichever was most recently loaded) while in 3D `FMM3D` is used. !!! warning "FMMLIB2D" @@ -158,8 +158,8 @@ Helper function to help determine the constant σ in the Green identity S\\[γ point is inside a domain or not. """ function _green_multiplier(x::SVector, Q::Quadrature{N}) where {N} - pde = Laplace(; dim = N) - K = DoubleLayerKernel(pde) + op = Laplace(; dim = N) + K = DoubleLayerKernel(op) σ = sum(Q.qnodes) do q return K(x, q) * weight(q) end diff --git a/src/vdim.jl b/src/vdim.jl index 142019c8..467f8cca 100644 --- a/src/vdim.jl +++ b/src/vdim.jl @@ -1,11 +1,11 @@ """ - vdim_correction(pde,X,Y,Y_boundary,S,D,V; green_multiplier, kwargs...) + vdim_correction(op,X,Y,Y_boundary,S,D,V; green_multiplier, kwargs...) Compute a correction to the volume potential `V : Y → X` such that `V + δV` is a more accurate approximation of the underlying volume potential operator. The correction is computed using the (volume) density interpolation method. -This function requires a `pde::AbstractDifferentialOperator`, a target set `X`, a source +This function requires a `op::AbstractDifferentialOperator`, a target set `X`, a source quadrature `Y`, a boundary quadrature `Y_boundary`, approximations `S : Y_boundary -> X` and `D : Y_boundary -> X` to the single- and double-layer potentials (correctly handling nearly-singular integrals), and a naive @@ -28,7 +28,7 @@ See [anderson2024fast](@cite) for more details on the method. and rescaled to each element. """ function vdim_correction( - pde, + op, target, source::Quadrature, boundary::Quadrature, @@ -47,7 +47,7 @@ function vdim_correction( @assert eltype(Dop) == eltype(Sop) == T "eltype of Sop, Dop, and Vop must match" # figure out if we are dealing with a scalar or vector PDE m, n = length(target), length(source) - N = ambient_dimension(pde) + N = ambient_dimension(op) @assert ambient_dimension(source) == N "vdim only works for volume potentials" m, n = length(target), length(source) # a reasonable interpolation_order if not provided @@ -55,7 +55,7 @@ function vdim_correction( (interpolation_order = maximum(order, values(source.etype2qrule))) # by default basis centered at origin center = isnothing(center) ? zero(SVector{N,Float64}) : center - p, P, γ₁P, multiindices = polynomial_solutions_vdim(pde, interpolation_order, center) + p, P, γ₁P, multiindices = polynomial_solutions_vdim(op, interpolation_order, center) dict_near = etype_to_nearest_points(target, source; maxdist) R = _vdim_auxiliary_quantities( p, @@ -295,28 +295,28 @@ function vdim_mesh_center(msh::AbstractMesh) end """ - polynomial_solutions_vdim(pde, order[, center]) + polynomial_solutions_vdim(op, order[, center]) For every monomial term `pₙ` of degree `order`, compute a polynomial `Pₙ` such -that `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `pde`. +that `ℒ[Pₙ] = pₙ`, where `ℒ` is the differential operator associated with `op`. This function returns `{pₙ,Pₙ,γ₁Pₙ}`, where `γ₁Pₙ` is the generalized Neumann trace of `Pₙ`. Passing a point `center` will shift the monomials and solutions accordingly. """ function polynomial_solutions_vdim( - pde::AbstractDifferentialOperator, + op::AbstractDifferentialOperator, order::Integer, center = nothing, ) - N = ambient_dimension(pde) + N = ambient_dimension(op) center = isnothing(center) ? zero(SVector{N,Float64}) : center # create empty arrays to store the monomials, solutions, and traces. For the # neumann trace, we try to infer the concrete return type instead of simply # having a vector of `Function`. monomials = Vector{ElementaryPDESolutions.Polynomial{N,Float64}}() dirchlet_traces = Vector{ElementaryPDESolutions.Polynomial{N,Float64}}() - T = return_type(neumann_trace, typeof(pde), eltype(dirchlet_traces)) + T = return_type(neumann_trace, typeof(op), eltype(dirchlet_traces)) neumann_traces = Vector{T}() multiindices = Vector{MultiIndex{N}}() # iterate over N-tuples going from 0 to order @@ -325,8 +325,8 @@ function polynomial_solutions_vdim( # define the monomial basis functions, and the corresponding solutions. # TODO: adapt this to vectorial case p = ElementaryPDESolutions.Polynomial(I => 1 / factorial(MultiIndex(I))) - P = polynomial_solution(pde, p) - γ₁P = neumann_trace(pde, P) + P = polynomial_solution(op, p) + γ₁P = neumann_trace(op, P) push!(multiindices, MultiIndex(I)) push!(monomials, p) push!(dirchlet_traces, P) @@ -351,8 +351,8 @@ function polynomial_solution(::Laplace, p::ElementaryPDESolutions.Polynomial) return ElementaryPDESolutions.convert_coefs(P, Float64) end -function polynomial_solution(pde::Helmholtz, p::ElementaryPDESolutions.Polynomial) - k = pde.k +function polynomial_solution(op::Helmholtz, p::ElementaryPDESolutions.Polynomial) + k = op.k P = ElementaryPDESolutions.solve_helmholtz(p; k) return ElementaryPDESolutions.convert_coefs(P, Float64) end diff --git a/test/adaptive_test.jl b/test/adaptive_test.jl index 962b840f..9173447f 100644 --- a/test/adaptive_test.jl +++ b/test/adaptive_test.jl @@ -65,22 +65,22 @@ end Inti.Helmholtz(; k = 1.2, dim = N), # Inti.Stokes(; μ = 1.2, dim = N), ) - for pde in ops - @testset "Greens identity ($t) $(N)d $pde" begin + for op in ops + @testset "Greens identity ($t) $(N)d $op" begin xs = t == :interior ? ntuple(i -> 3, N) : ntuple(i -> 0.1, N) - T = Inti.default_density_eltype(pde) + T = Inti.default_density_eltype(op) c = rand(T) - u = (qnode) -> Inti.SingleLayerKernel(pde)(qnode, xs) * c - dudn = (qnode) -> Inti.AdjointDoubleLayerKernel(pde)(qnode, xs) * c + u = (qnode) -> Inti.SingleLayerKernel(op)(qnode, xs) * c + dudn = (qnode) -> Inti.AdjointDoubleLayerKernel(op)(qnode, xs) * c γ₀u = map(u, Γ_quad) γ₁u = map(dudn, Γ_quad) γ₀u_norm = norm(norm.(γ₀u, Inf), Inf) γ₁u_norm = norm(norm.(γ₁u, Inf), Inf) # single and double layer - G = Inti.SingleLayerKernel(pde) + G = Inti.SingleLayerKernel(op) S = Inti.IntegralOperator(G, Γ_quad) S0 = Inti.assemble_matrix(S) - dG = Inti.DoubleLayerKernel(pde) + dG = Inti.DoubleLayerKernel(op) D = Inti.IntegralOperator(dG, Γ_quad) D0 = Inti.assemble_matrix(D) e0 = norm(S0 * γ₁u - D0 * γ₀u - σ * γ₀u, Inf) / γ₀u_norm @@ -91,7 +91,7 @@ end δD = Inti.adaptive_correction(D; maxdist, tol = atol) Smat, Dmat = S0 + δS, D0 + δD e1 = norm(Smat * γ₁u - Dmat * γ₀u - σ * γ₀u, Inf) / γ₀u_norm - @testset "Single/double layer $(string(pde))" begin + @testset "Single/double layer $(string(op))" begin @test norm(e0, Inf) > 10 * norm(e1, Inf) @test norm(e1, Inf) < 10 * atol end diff --git a/test/convergence/dim_convergence.jl b/test/convergence/dim_convergence.jl index 25adf384..66a5f07b 100644 --- a/test/convergence/dim_convergence.jl +++ b/test/convergence/dim_convergence.jl @@ -13,10 +13,10 @@ t = :exterior σ = t == :interior ? 1 / 2 : -1 / 2 # For N = 3 one needs to use compression, e.g. `compression = :fmm`, for the operators below N = 2 -pde = Inti.Laplace(; dim = N) -# pde = Inti.Helmholtz(; dim = N, k = 2π) -# pde = Inti.Stokes(; dim = N, μ = 1.2) -@info "Greens identity ($t) $(N)d $pde" +op = Inti.Laplace(; dim = N) +# op = Inti.Helmholtz(; dim = N, k = 2π) +# op = Inti.Stokes(; dim = N, μ = 1.2) +@info "Greens identity ($t) $(N)d $op" Inti.clear_entities!() center = Inti.svector(i -> 0.1, N) radius = 1 @@ -53,10 +53,10 @@ for h in hh else center + Inti.svector(i -> 0.5 * radius, N) end - T = Inti.default_density_eltype(pde) + T = Inti.default_density_eltype(op) c = rand(T) - G = Inti.SingleLayerKernel(pde) - dG = Inti.DoubleLayerKernel(pde) + G = Inti.SingleLayerKernel(op) + dG = Inti.DoubleLayerKernel(op) u = (qnode) -> G(xs, qnode) * c dudn = (qnode) -> transpose(dG(xs, qnode)) * c γ₀u = map(u, Q) @@ -65,14 +65,14 @@ for h in hh γ₁u_norm = norm(norm.(γ₁u, Inf), Inf) # single and double layer S0, D0 = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), correction = (method = :none,), ) S1, D1 = Inti.single_double_layer(; - pde, + op, target = Q, source = Q, compression = (method = :none,), diff --git a/test/convergence/vdim_potential_convergence.jl b/test/convergence/vdim_potential_convergence.jl index 110d1aa7..b55e30d1 100644 --- a/test/convergence/vdim_potential_convergence.jl +++ b/test/convergence/vdim_potential_convergence.jl @@ -56,12 +56,12 @@ function test_volume_potential(; meshsize, bdry_qorder, interpolation_order) du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) f_d = map(q -> f(q.coords), Ωₕ_quad) - pde = k == 0 ? Inti.Laplace(; dim = 2) : Inti.Helmholtz(; dim = 2, k) + op = k == 0 ? Inti.Laplace(; dim = 2) : Inti.Helmholtz(; dim = 2, k) ## Boundary operators tbnd = @elapsed begin S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ωₕ_quad, source = Γₕ_quad, compression = (method = :hmatrix, tol = 1e-14), @@ -73,7 +73,7 @@ function test_volume_potential(; meshsize, bdry_qorder, interpolation_order) ## Volume potentials tvol = @elapsed begin V_d2d = Inti.volume_potential(; - pde, + op, target = Ωₕ_quad, source = Ωₕ_quad, compression = (method = :hmatrix, tol = 1e-14), diff --git a/test/convergence/vdim_potential_script.jl b/test/convergence/vdim_potential_script.jl index 496b9af1..e49d758d 100644 --- a/test/convergence/vdim_potential_script.jl +++ b/test/convergence/vdim_potential_script.jl @@ -60,12 +60,12 @@ u_b = map(q -> u(q.coords), Γₕ_quad) du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) f_d = map(q -> f(q.coords), Ωₕ_quad) -pde = k == 0 ? Inti.Laplace(; dim = 2) : Inti.Helmholtz(; dim = 2, k) +op = k == 0 ? Inti.Laplace(; dim = 2) : Inti.Helmholtz(; dim = 2, k) ## Boundary operators tbnd = @elapsed begin S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ωₕ_quad, source = Γₕ_quad, compression = (method = :fmm, tol = 1e-14), @@ -77,7 +77,7 @@ end ## Volume potentials tvol = @elapsed begin V_d2d = Inti.volume_potential(; - pde, + op, target = Ωₕ_quad, source = Ωₕ_quad, compression = (method = :fmm, tol = 1e-14), diff --git a/test/convergence/vdim_potential_script_3d.jl b/test/convergence/vdim_potential_script_3d.jl index 2a512ee7..5c59795f 100644 --- a/test/convergence/vdim_potential_script_3d.jl +++ b/test/convergence/vdim_potential_script_3d.jl @@ -60,12 +60,12 @@ u_b = map(q -> u(q.coords), Γₕ_quad) du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) f_d = map(q -> f(q.coords), Ωₕ_quad) -pde = k == 0 ? Inti.Laplace(; dim = 3) : Inti.Helmholtz(; dim = 3, k) +op = k == 0 ? Inti.Laplace(; dim = 3) : Inti.Helmholtz(; dim = 3, k) ## Boundary operators tbnd = @elapsed begin S_b2d, D_b2d = Inti.single_double_layer(; - pde, + op, target = Ωₕ_quad, source = Γₕ_quad, compression = (method = :fmm, tol = 1e-8), @@ -77,7 +77,7 @@ end ## Volume potentials tvol = @elapsed begin V_d2d = Inti.volume_potential(; - pde, + op, target = Ωₕ_quad, source = Ωₕ_quad, compression = (method = :fmm, tol = 1e-8), diff --git a/test/dim_test.jl b/test/dim_test.jl index 7c93ace1..6cd47dae 100644 --- a/test/dim_test.jl +++ b/test/dim_test.jl @@ -34,54 +34,54 @@ for N in dims Inti.Stokes(; μ = 1.2, dim = N), Inti.Elastostatic(; λ = 1, μ = 1, dim = N), ) - for pde in ops - @testset "Greens identity ($t) $(N)d $pde" begin + for op in ops + @testset "Greens identity ($t) $(N)d $op" begin xs = t == :interior ? ntuple(i -> 3, N) : ntuple(i -> 0.1, N) - T = Inti.default_density_eltype(pde) + T = Inti.default_density_eltype(op) c = rand(T) - u = (qnode) -> Inti.SingleLayerKernel(pde)(qnode, xs) * c - dudn = (qnode) -> Inti.AdjointDoubleLayerKernel(pde)(qnode, xs) * c + u = (qnode) -> Inti.SingleLayerKernel(op)(qnode, xs) * c + dudn = (qnode) -> Inti.AdjointDoubleLayerKernel(op)(qnode, xs) * c γ₀u = map(u, quad) γ₁u = map(dudn, quad) γ₀u_norm = norm(norm.(γ₀u, Inf), Inf) γ₁u_norm = norm(norm.(γ₁u, Inf), Inf) # single and double layer - G = Inti.SingleLayerKernel(pde) + G = Inti.SingleLayerKernel(op) S = Inti.IntegralOperator(G, quad) Smat = Inti.assemble_matrix(S) - dG = Inti.DoubleLayerKernel(pde) + dG = Inti.DoubleLayerKernel(op) D = Inti.IntegralOperator(dG, quad) Dmat = Inti.assemble_matrix(D) e0 = norm(Smat * γ₁u - Dmat * γ₀u - σ * γ₀u, Inf) / γ₀u_norm Sdim, Ddim = Inti.single_double_layer(; - pde, + op, target = quad, source = quad, compression = (method = :none,), correction = (method = :dim,), ) e1 = norm(Sdim * γ₁u - Ddim * γ₀u - σ * γ₀u, Inf) / γ₀u_norm - @testset "Single/double layer $(string(pde))" begin + @testset "Single/double layer $(string(op))" begin @test norm(e0, Inf) > norm(e1, Inf) @test norm(e1, Inf) < rtol1 end # adjoint double-layer and hypersingular. - pde isa Inti.Stokes && continue # TODO: implement hypersingular for Stokes? + op isa Inti.Stokes && continue # TODO: implement hypersingular for Stokes? - K = Inti.IntegralOperator(Inti.AdjointDoubleLayerKernel(pde), quad) + K = Inti.IntegralOperator(Inti.AdjointDoubleLayerKernel(op), quad) Kmat = Inti.assemble_matrix(K) - H = Inti.IntegralOperator(Inti.HyperSingularKernel(pde), quad) + H = Inti.IntegralOperator(Inti.HyperSingularKernel(op), quad) Hmat = Inti.assemble_matrix(H) e0 = norm(Kmat * γ₁u - Hmat * γ₀u - σ * γ₁u, Inf) / γ₁u_norm Kdim, Hdim = Inti.adj_double_layer_hypersingular(; - pde = pde, + op = op, target = quad, source = quad, compression = (method = :none,), correction = (method = :dim,), ) e1 = norm(Kdim * γ₁u - Hdim * γ₀u - σ * γ₁u, Inf) / γ₁u_norm - @testset "Adjoint double-layer/hypersingular $(string(pde))" begin + @testset "Adjoint double-layer/hypersingular $(string(op))" begin @test norm(e0, Inf) > norm(e1, Inf) @test norm(e1, Inf) < rtol2 end diff --git a/test/fmm2d_test.jl b/test/fmm2d_test.jl index 222a0508..8921f572 100644 --- a/test/fmm2d_test.jl +++ b/test/fmm2d_test.jl @@ -17,13 +17,13 @@ include("test_utils.jl") Γ₂ = Inti.external_boundary(Ω₂) Γ₂_quad = Inti.Quadrature(view(msh₂, Γ₂); qorder = 3) -for pde in (Inti.Laplace(; dim = 2), Inti.Helmholtz(; dim = 2, k = 1.2)) - @testset "PDE: $pde" begin +for op in (Inti.Laplace(; dim = 2), Inti.Helmholtz(; dim = 2, k = 1.2)) + @testset "PDE: $op" begin for K in ( - Inti.DoubleLayerKernel(pde), - Inti.SingleLayerKernel(pde), - Inti.AdjointDoubleLayerKernel(pde), - Inti.HyperSingularKernel(pde), + Inti.DoubleLayerKernel(op), + Inti.SingleLayerKernel(op), + Inti.AdjointDoubleLayerKernel(op), + Inti.HyperSingularKernel(op), ) for Γ_quad in (Γ₁_quad, Γ₂_quad) iop = Inti.IntegralOperator(K, Γ₁_quad, Γ_quad) diff --git a/test/fmm3d_test.jl b/test/fmm3d_test.jl index 4dab51ea..c5e2cabd 100644 --- a/test/fmm3d_test.jl +++ b/test/fmm3d_test.jl @@ -16,25 +16,25 @@ include("test_utils.jl") Γ₂ = Inti.external_boundary(Ω₂) Γ₂_quad = Inti.Quadrature(view(msh₂, Γ₂); qorder = 3) -for pde in ( +for op in ( Inti.Laplace(; dim = 3), Inti.Helmholtz(; dim = 3, k = 1.2), Inti.Stokes(; dim = 3, μ = 0.5), ) - @testset "PDE: $pde" begin + @testset "PDE: $op" begin for K in ( - Inti.DoubleLayerKernel(pde), - Inti.SingleLayerKernel(pde), - Inti.AdjointDoubleLayerKernel(pde), - Inti.HyperSingularKernel(pde), + Inti.DoubleLayerKernel(op), + Inti.SingleLayerKernel(op), + Inti.AdjointDoubleLayerKernel(op), + Inti.HyperSingularKernel(op), ) # TODO Stokes has only single and double layer implemented for now - (K isa Inti.AdjointDoubleLayerKernel && pde isa Inti.Stokes) && continue - (K isa Inti.HyperSingularKernel && pde isa Inti.Stokes) && continue + (K isa Inti.AdjointDoubleLayerKernel && op isa Inti.Stokes) && continue + (K isa Inti.HyperSingularKernel && op isa Inti.Stokes) && continue for Γ_quad in (Γ₁_quad, Γ₂_quad) iop = Inti.IntegralOperator(K, Γ₁_quad, Γ_quad) iop_fmm = Inti.assemble_fmm(iop; rtol = 1e-8) - x = rand(Inti.default_density_eltype(pde), size(iop, 2)) + x = rand(Inti.default_density_eltype(op), size(iop, 2)) yapprox = iop_fmm * x # test on a given index set idx_test = rand(1:size(iop, 1), 10) diff --git a/test/hmatrices_test.jl b/test/hmatrices_test.jl index 88f17894..8c255494 100644 --- a/test/hmatrices_test.jl +++ b/test/hmatrices_test.jl @@ -17,9 +17,9 @@ include("test_utils.jl") Γ_quad = Inti.Quadrature(view(msh, Γ); qorder = 3) W = [q.weight for q in Γ_quad] # test various PDEs and integral operators - for pde in (Inti.Laplace(; dim = 2), Inti.Helmholtz(; k = 1.2, dim = 2)) - @testset "PDE = $pde" begin - for K in (Inti.SingleLayerKernel(pde), Inti.DoubleLayerKernel(pde)) + for op in (Inti.Laplace(; dim = 2), Inti.Helmholtz(; k = 1.2, dim = 2)) + @testset "PDE = $op" begin + for K in (Inti.SingleLayerKernel(op), Inti.DoubleLayerKernel(op)) iop = Inti.IntegralOperator(K, Γ_quad) H = Inti.assemble_hmatrix(iop; atol = 1e-8) x = rand(eltype(iop), size(iop, 2)) diff --git a/test/integral_operator_test.jl b/test/integral_operator_test.jl index 55dc8788..a4335f5b 100644 --- a/test/integral_operator_test.jl +++ b/test/integral_operator_test.jl @@ -7,10 +7,10 @@ using Inti mesh = Inti.meshgen(Γ; meshsize = 0.4) quad = Inti.Quadrature(mesh; qorder = 1) - pde = Inti.Laplace(; dim = 3) + op = Inti.Laplace(; dim = 3) @test_throws ArgumentError Inti.single_double_layer(; - pde, + op, source = quad, target = quad, compression = (method = :none,), diff --git a/test/quadrature_test.jl b/test/quadrature_test.jl index 25915ae8..dded9f94 100644 --- a/test/quadrature_test.jl +++ b/test/quadrature_test.jl @@ -143,10 +143,10 @@ end mesh = Inti.meshgen(Γ; meshsize = 0.4) quad = Inti.Quadrature(mesh; qorder = 1) - pde = Inti.Laplace(; dim = 3) + op = Inti.Laplace(; dim = 3) @test_throws ArgumentError Inti.single_double_layer(; - pde, + op, source = quad, target = quad, compression = (method = :none,), diff --git a/test/vdim_elastostatic_test.jl b/test/vdim_elastostatic_test.jl new file mode 100644 index 00000000..b997357c --- /dev/null +++ b/test/vdim_elastostatic_test.jl @@ -0,0 +1,87 @@ +# # High-order convergence of vdim + +using Inti +using StaticArrays +using Gmsh +using LinearAlgebra +using HMatrices +using CairoMakie + +meshsize = 0.2 +meshorder = 1 +bdry_qorder = 8 +interpolation_order = 2 + +Inti.clear_entities!() +gmsh.initialize() +gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) +gmsh.option.setNumber("Mesh.MeshSizeMin", meshsize) +gmsh.model.occ.addDisk(0, 0, 0, 1, 1) +gmsh.model.occ.synchronize() +gmsh.model.mesh.generate(2) +gmsh.model.mesh.setOrder(meshorder) +msh = Inti.import_mesh(; dim = 2) +Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 2, Inti.entities(msh)) +gmsh.finalize() + +Γ = Inti.external_boundary(Ω) +Ωₕ = view(msh, Ω) +Γₕ = view(msh, Γ) +## + +VR_qorder = Inti.Triangle_VR_interpolation_order_to_quadrature_order(interpolation_order) + +# Use VDIM with the Vioreanu-Rokhlin quadrature rule for Ωₕ +Q = Inti.VioreanuRokhlin(; domain = :triangle, order = VR_qorder) +dict = Dict(E => Q for E in Inti.element_types(Ωₕ)) +Ωₕ_quad = Inti.Quadrature(Ωₕ, dict) +# Ωₕ_quad = Inti.Quadrature(Ωₕ; qorder = qorders[1]) +Γₕ_quad = Inti.Quadrature(Γₕ; qorder = bdry_qorder) + +## + +# build exact solution +μ = 1.0 +λ = 1.0 + +u = (x) -> SVector(1.0, 1.0) +du = (x, n) -> SVector(0.0, 0.0) +f = (x) -> SVector(0.0, 0.0) + +u_d = map(q -> u(q.coords), Ωₕ_quad) +u_b = map(q -> u(q.coords), Γₕ_quad) +du_b = map(q -> du(q.coords, q.normal), Γₕ_quad) +f_d = map(q -> f(q.coords), Ωₕ_quad) + +op = Inti.Elastostatic(; λ = λ, μ = μ, dim = 2) + +# m, d, n = Inti._polynomial_solutions_vec(op, 1) +# x = @SVector rand(2) +# m[1](x) +# d[1](x) +# n[1]((; coords = x, normal = x)) + +## Boundary operators +S_b2d, D_b2d = Inti.single_double_layer(; + op, + target = Ωₕ_quad, + source = Γₕ_quad, + compression = (method = :none, tol = 1e-14), + correction = (method = :dim, maxdist = 5 * meshsize, target_location = :inside), +) + +## Volume potentials +V_d2d = Inti.volume_potential(; + op, + target = Ωₕ_quad, + source = Ωₕ_quad, + compression = (method = :none, tol = 1e-14), + correction = (method = :dim, interpolation_order), +) + +vref = -u_d - D_b2d * u_b + S_b2d * du_b +vapprox = V_d2d * f_d +er = vref - vapprox + +ndofs = length(er) +@show ndofs, meshsize, norm(er, Inf), t From fe1bc809fb917752539f04e90f0f442f750a81b9 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Sun, 13 Oct 2024 13:30:28 +0200 Subject: [PATCH 6/7] run documentation on github worker --- .github/workflows/documentation.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index c781be80..70b33764 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -12,7 +12,7 @@ jobs: permissions: contents: write statuses: write - runs-on: self-hosted + runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 From 5f0b313fba5090fd87d3e70c6d1892ae3e2091cf Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Sun, 13 Oct 2024 13:30:45 +0200 Subject: [PATCH 7/7] add tolerance for codecov failure --- .github/codecov.yml | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 .github/codecov.yml diff --git a/.github/codecov.yml b/.github/codecov.yml new file mode 100644 index 00000000..b5edfb1f --- /dev/null +++ b/.github/codecov.yml @@ -0,0 +1,8 @@ +coverage: + status: + project: + default: # default is the status check's name, not default settings + target: auto + threshold: 5 + flags: + - unit \ No newline at end of file