From 9fe6af40b526de245993b6e4be13ca144d768796 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 19 Sep 2024 15:04:49 +0200 Subject: [PATCH 01/17] Implemented the composed basis functionallity. This allows us to multiply, crossproduct or inner product a basis with a function or the normal vector. All functionallity is screened away from the user. Discussion needed about what can be made visible. This commit will support the later implemented composed operators. --- src/BEAST.jl | 3 + src/bases/composedbasis.jl | 91 +++++++++++++++++++++++++++ src/bases/local/localcomposedbasis.jl | 53 ++++++++++++++++ test/runtests.jl | 1 + test/test_composed_basis.jl | 34 ++++++++++ 5 files changed, 182 insertions(+) create mode 100644 src/bases/composedbasis.jl create mode 100644 src/bases/local/localcomposedbasis.jl create mode 100644 test/test_composed_basis.jl diff --git a/src/BEAST.jl b/src/BEAST.jl index 14daed8a..783ce83a 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -182,6 +182,9 @@ include("bases/stagedtimestep.jl") include("bases/timebasis.jl") include("bases/tensorbasis.jl") +include("bases/composedbasis.jl") +include("bases/local/localcomposedbasis.jl") + include("operator.jl") include("quadrature/quadstrats.jl") diff --git a/src/bases/composedbasis.jl b/src/bases/composedbasis.jl new file mode 100644 index 00000000..9b767d99 --- /dev/null +++ b/src/bases/composedbasis.jl @@ -0,0 +1,91 @@ +struct FunctionWrapper{T} + func::Function + function FunctionWrapper(f::Function;evalpoint = @SVector [0.0,0.0,0.0]) + new{typeof(f(evalpoint))}(f) + end +end + +function (f::FunctionWrapper{T})(x)::T where {T} + return f.func(x) +end + +function (f::FunctionWrapper{T})(x::CompScienceMeshes.MeshPointNM)::T where {T} + return f.func(cartesian(x)) +end +scalartype(F::FunctionWrapper{T}) where {T} = eltype(T) +function scalartype(::NormalVector) + @warn "The scallartype of the NormalVector is set at Float32, if used in combination with Float64 basis or operator this is no problem, in the case of Float16 it is." + return Float32 +end + +abstract type _BasisOperations{T} <: Space{T} end + +struct _BasisTimes{T} <: _BasisOperations{T} + el1 + el2 + function _BasisTimes(el1,el2) + new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) + end + function _BasisTimes{T}(el1,el2) where {T} + new{T}(el1,el2) + end +end + +struct _BasisCross{T} <: _BasisOperations{T} + el1 + el2 + function _BasisCross(el1,el2) + new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) + end + function _BasisCross{T}(el1,el2) where {T} + new{T}(el1,el2) + end +end + +struct _BasisDot{T} <: _BasisOperations{T} + el1 + el2 + function _BasisDot(el1,el2) + new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) + end + function _BasisDot{T}(el1,el2) where {T} + new{T}(el1,el2) + end +end + +#### wrapping of the functions +_BasisTimes(a::Function,b::Function) = _BasisTimes(FunctionWrapper(a),FunctionWrapper(b)) +_BasisTimes(a::Function,b) = _BasisTimes(FunctionWrapper(a),b) +_BasisTimes(a,b::Function) = _BasisTimes(a,FunctionWrapper(b)) + +_BasisCross(a::Function,b::Function) = _BasisCross(FunctionWrapper(a),FunctionWrapper(b)) +_BasisCross(a::Function,b) = _BasisCross(FunctionWrapper(a),b) +_BasisCross(a,b::Function) = _BasisCross(a,FunctionWrapper(b)) + +_BasisDot(a::Function,b::Function) = _BasisDot(FunctionWrapper(a),FunctionWrapper(b)) +_BasisDot(a::Function,b) = _BasisDot(FunctionWrapper(a),b) +_BasisDot(a,b::Function) = _BasisDot(a,FunctionWrapper(b)) + + +refspace(a::_BasisTimes{T}) where {T} = _LocalBasisTimes(T,refspace(a.el1),refspace(a.el2)) +refspace(a::_BasisCross{T}) where {T} = _LocalBasisCross(T,refspace(a.el1),refspace(a.el2)) +refspace(a::_BasisDot{T}) where {T} = _LocalBasisDot(T,refspace(a.el1),refspace(a.el2)) +refspace(a::Function) = FunctionWrapper(a) +refspace(a::FunctionWrapper) = a +refspace(a::NormalVector) = a + +numfunctions(a::Union{NormalVector,FunctionWrapper,Function}) = missing +numfunctions(a::_BasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) + +geometry(a::_BasisOperations) = coalesce(geometry(a.el1),geometry(a.el2)) +basisfunction(a::_BasisOperations,i) = coalesce(basisfunction(a.el1,i),basisfunction(a.el2,i)) +positions(a::_BasisOperations) = coalesce(positions(a.el1),positions(a.el2)) + +geometry(a::Union{NormalVector,FunctionWrapper,Function}) = missing +basisfunction(a::Union{NormalVector,FunctionWrapper,Function},i) = missing +positions(a::Union{NormalVector,FunctionWrapper,Function}) = missing + + +subset(a::T,I) where {T <: _BasisOperations} = T(subset(a.el1,I),subset(a.el2,I)) +subset(a::Union{NormalVector,FunctionWrapper,Function},I) = a + diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl new file mode 100644 index 00000000..7b169bec --- /dev/null +++ b/src/bases/local/localcomposedbasis.jl @@ -0,0 +1,53 @@ +abstract type _LocalBasisOperations{T} <: RefSpace{T,:None} end +numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) + +struct _LocalBasisTimes{T,U,V} <: _LocalBasisOperations{T} + el1::U + el2::V + function _LocalBasisTimes(::Type{T},el1::U,el2::V) where {T,U,V} + new{T,U,V}(el1,el2) + end +end + +struct _LocalBasisCross{T,U,V} <: _LocalBasisOperations{T} + el1::U + el2::V + function _LocalBasisCross(::Type{T},el1::U,el2::V) where {T,U,V} + new{T,U,V}(el1,el2) + end +end + +struct _LocalBasisDot{T,U,V} <: _LocalBasisOperations{T} + el1::U + el2::V + function _LocalBasisDot(::Type{T},el1::U,el2::V) where {T,U,V} + new{T,U,V}(el1,el2) + end +end + +function (op::U where {U <: _LocalBasisOperations})(p) + l = op.el1(p) + r = op.el2(p) + return _execute_operation(l,r,op) +end + +function (op::NormalVector)(x::CompScienceMeshes.MeshPointNM) + return normal(x) +end + +operation(a::_LocalBasisTimes) = * +operation(a::_LocalBasisCross) = × +operation(a::_LocalBasisDot) = ⋅ + +_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector,op::U) where {N,U} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) +_execute_operation(el1::SVector{N,<:NamedTuple},el2::U,op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) +_execute_operation(el1::SVector,el2::SVector{N,<:NamedTuple},op::U) where {N,U} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) +_execute_operation(el1::U,el2::SVector{N,<:NamedTuple},op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) +_execute_operation(el1::SVector,el2::SVector,op::U) where {U <: Union{_LocalBasisCross,_LocalBasisDot}} = operation(op)(el1,el2) +_execute_operation(el1::U,el2::V,op::_LocalBasisTimes) where {U <: Number,V <: Number}= el1*el2 +_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not (yet) supported" + +_execute_operation_named(a::NTuple{N},b,op) where {N} = ((value=op(a[1].value,b),), _execute_operation_named(Base.tail(a),b,op)...) +_execute_operation_named(a::NTuple{1},b,op) = ((value=op(a[1].value,b),),) +_execute_operation_named(b,a::NTuple{N},op) where {N} = ((value=op(b,a[1].value),), _execute_operation_named(b,Base.tail(a),op)...) +_execute_operation_named(b,a::NTuple{1},op) = ((value=op(b,a[1].value),),) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 956d5153..888b814e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,7 @@ include("test_variational.jl") include("test_handlers.jl") +include("test_composed_basis.jl") using TestItemRunner @run_package_tests diff --git a/test/test_composed_basis.jl b/test/test_composed_basis.jl new file mode 100644 index 00000000..078bf870 --- /dev/null +++ b/test/test_composed_basis.jl @@ -0,0 +1,34 @@ +######## test multiplied basis +using CompScienceMeshes +using LinearAlgebra +using BEAST +using Test + +Γ = meshcuboid(1.0,1.0,1.0,1.0) +X = raviartthomas(Γ) +f(x) = 2.0 +Y = BEAST._BasisTimes(X,f) +K = Maxwell3D.doublelayer(wavenumber=1.0) +m1 = assemble(K,Y,Y) +m2 = assemble(K,X,X) +@test m1 ≈ 4*m2 + +U = BEAST._BasisTimes(BEAST._BasisDot(x->(@SVector [1.0,1.0,1.0]),n),X) +m = assemble(K,U,U) + +L = lagrangecxd0(Γ) +Z = BEAST._BasisTimes(L,n) +m3 = assemble(K,Z,Z) +@test norm(m3) ≈ 0.08420116178577139 + + + +##### go in code +# using StaticArrays +# s = simplex((@SVector [1.0,0.0,0.0]),(@SVector [0.0,1.0,0.0]),(@SVector [0.0,0.0,0.0])) +# p = neighborhood(s,[0.5,0.2]) + +# rtref = BEAST.RTRefSpace{Float64}() +# f(x)= 2.0 + +# mref = BEAST._LocalBasisTimes(Float64,BEAST.FunctionWrapper(f),rtref) From d01cb68ab78723750636f4a19300a98aa8906765 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 19 Sep 2024 15:30:46 +0200 Subject: [PATCH 02/17] update_test_file --- test/test_composed_basis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_composed_basis.jl b/test/test_composed_basis.jl index 078bf870..e4ae5b24 100644 --- a/test/test_composed_basis.jl +++ b/test/test_composed_basis.jl @@ -6,8 +6,8 @@ using Test Γ = meshcuboid(1.0,1.0,1.0,1.0) X = raviartthomas(Γ) -f(x) = 2.0 -Y = BEAST._BasisTimes(X,f) +func(x) = 2.0 +Y = BEAST._BasisTimes(X,func) K = Maxwell3D.doublelayer(wavenumber=1.0) m1 = assemble(K,Y,Y) m2 = assemble(K,X,X) From 4f392862bebb9e99dc2f462aafebab00d553c424 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 31 Oct 2024 14:49:32 +0100 Subject: [PATCH 03/17] -) composed operator -) potentials -) displacement meshes -) trace operations --- examples/composedoperator.jl | 21 ++ src/BEAST.jl | 6 +- src/bases/composedbasis.jl | 100 +++--- src/bases/local/localcomposedbasis.jl | 6 +- src/composedoperators/composedoperator.jl | 413 ++++++++++++++++++++++ src/composedoperators/displacementmesh.jl | 29 ++ src/composedoperators/potentials.jl | 12 + src/composedoperators/smartmath.jl | 10 + src/composedoperators/trace.jl | 123 +++++++ 9 files changed, 672 insertions(+), 48 deletions(-) create mode 100644 examples/composedoperator.jl create mode 100644 src/composedoperators/composedoperator.jl create mode 100644 src/composedoperators/displacementmesh.jl create mode 100644 src/composedoperators/potentials.jl create mode 100644 src/composedoperators/smartmath.jl create mode 100644 src/composedoperators/trace.jl diff --git a/examples/composedoperator.jl b/examples/composedoperator.jl new file mode 100644 index 00000000..5cf9f764 --- /dev/null +++ b/examples/composedoperator.jl @@ -0,0 +1,21 @@ +##### user manual composed operator and potential + +using BEAST +using CompScienceMeshes + +Γ = meshcuboid(1.0,1.0,1.0,0.2); +X = lagrangecxd0(Γ) + +S = BEAST.PotentialIntegralOperator{2}(BEAST.HH3DGreen(1im),*,B->B) +St = BEAST.trace(S,1.0) +assemble(St,X,X) + +SB = Helmholtz3D.singlelayer(wavenumber=1.0) +assemble(SB,X,X) + +Y = raviartthomas(Γ) +K = BEAST.PotentialIntegralOperator{2}(BEAST.HH3DGradGreen(0*1im),×,B->B) +Kt = BEAST.ttrace(K,1.0;testfunction_tangential=true) +assemble(Kt,Y,Y) + +KB = Maxwell3D.doublelayer(wavenumber=0.0) \ No newline at end of file diff --git a/src/BEAST.jl b/src/BEAST.jl index ad32d754..a0013360 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -289,8 +289,10 @@ include("solvers/itsolver.jl") include("utils/plotlyglue.jl") - - +include("composedoperators/composedoperator.jl") +include("composedoperators/displacementmesh.jl") +include("composedoperators/potentials.jl") +include("composedoperators/trace.jl") const x̂ = point(1, 0, 0) const ŷ = point(0, 1, 0) diff --git a/src/bases/composedbasis.jl b/src/bases/composedbasis.jl index 9b767d99..c6817193 100644 --- a/src/bases/composedbasis.jl +++ b/src/bases/composedbasis.jl @@ -1,10 +1,9 @@ struct FunctionWrapper{T} func::Function - function FunctionWrapper(f::Function;evalpoint = @SVector [0.0,0.0,0.0]) - new{typeof(f(evalpoint))}(f) - end end - +# function FunctionWrapper(f::Function;evalpoint = @SVector [0.0,0.0,0.0]) +# FunctionWrapper{eltype(f(evalpoint))}(f) +# end function (f::FunctionWrapper{T})(x)::T where {T} return f.func(x) end @@ -13,79 +12,94 @@ function (f::FunctionWrapper{T})(x::CompScienceMeshes.MeshPointNM)::T where {T} return f.func(cartesian(x)) end scalartype(F::FunctionWrapper{T}) where {T} = eltype(T) -function scalartype(::NormalVector) - @warn "The scallartype of the NormalVector is set at Float32, if used in combination with Float64 basis or operator this is no problem, in the case of Float16 it is." - return Float32 -end +#scalartype(::NormalVector,s::Space) = scalartype(geometry(s)) +# function scalartype(::NormalVector) +# @warn "The scallartype of the NormalVector is set at Float32, if used in combination with Float64 basis or operator this is no problem, in the case of Float16 it is." +# return Float32 +# end abstract type _BasisOperations{T} <: Space{T} end struct _BasisTimes{T} <: _BasisOperations{T} el1 el2 - function _BasisTimes(el1,el2) - new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) - end - function _BasisTimes{T}(el1,el2) where {T} - new{T}(el1,el2) - end end +_BasisTimes(el1::NormalVector,el2) = _BasisTimes{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) +_BasisTimes(el2,el1::NormalVector) = _BasisTimes{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) +_BasisTimes(el1,el2) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) struct _BasisCross{T} <: _BasisOperations{T} el1 el2 - function _BasisCross(el1,el2) - new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) - end - function _BasisCross{T}(el1,el2) where {T} - new{T}(el1,el2) - end end - +_BasisCross(el1::NormalVector,el2) = _BasisCross{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) +_BasisCross(el2,el1::NormalVector) = _BasisCross{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) +_BasisCross(el1,el2) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) struct _BasisDot{T} <: _BasisOperations{T} el1 el2 - function _BasisDot(el1,el2) - new{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) - end - function _BasisDot{T}(el1,el2) where {T} - new{T}(el1,el2) - end end +_BasisDot(el1::NormalVector,el2) = _BasisDot{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) +_BasisDot(el2,el1::NormalVector) = _BasisDot{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) +_BasisDot(el1,el2) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) -#### wrapping of the functions -_BasisTimes(a::Function,b::Function) = _BasisTimes(FunctionWrapper(a),FunctionWrapper(b)) -_BasisTimes(a::Function,b) = _BasisTimes(FunctionWrapper(a),b) -_BasisTimes(a,b::Function) = _BasisTimes(a,FunctionWrapper(b)) +# #### wrapping of the functions +# _BasisTimes(a::Function,b::Function) = _BasisTimes(FunctionWrapper(a),FunctionWrapper(b)) +# _BasisTimes(a::Function,b) = _BasisTimes(FunctionWrapper(a),b) +# _BasisTimes(a,b::Function) = _BasisTimes(a,FunctionWrapper(b)) -_BasisCross(a::Function,b::Function) = _BasisCross(FunctionWrapper(a),FunctionWrapper(b)) -_BasisCross(a::Function,b) = _BasisCross(FunctionWrapper(a),b) -_BasisCross(a,b::Function) = _BasisCross(a,FunctionWrapper(b)) +# _BasisCross(a::Function,b::Function) = _BasisCross(FunctionWrapper(a),FunctionWrapper(b)) +# _BasisCross(a::Function,b) = _BasisCross(FunctionWrapper(a),b) +# _BasisCross(a,b::Function) = _BasisCross(a,FunctionWrapper(b)) -_BasisDot(a::Function,b::Function) = _BasisDot(FunctionWrapper(a),FunctionWrapper(b)) -_BasisDot(a::Function,b) = _BasisDot(FunctionWrapper(a),b) -_BasisDot(a,b::Function) = _BasisDot(a,FunctionWrapper(b)) +# _BasisDot(a::Function,b::Function) = _BasisDot(FunctionWrapper(a),FunctionWrapper(b)) +# _BasisDot(a::Function,b) = _BasisDot(FunctionWrapper(a),b) +# _BasisDot(a,b::Function) = _BasisDot(a,FunctionWrapper(b)) refspace(a::_BasisTimes{T}) where {T} = _LocalBasisTimes(T,refspace(a.el1),refspace(a.el2)) refspace(a::_BasisCross{T}) where {T} = _LocalBasisCross(T,refspace(a.el1),refspace(a.el2)) refspace(a::_BasisDot{T}) where {T} = _LocalBasisDot(T,refspace(a.el1),refspace(a.el2)) -refspace(a::Function) = FunctionWrapper(a) +# refspace(a::Function) = FunctionWrapper(a) refspace(a::FunctionWrapper) = a refspace(a::NormalVector) = a -numfunctions(a::Union{NormalVector,FunctionWrapper,Function}) = missing +numfunctions(a::Union{NormalVector,FunctionWrapper}) = missing numfunctions(a::_BasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) geometry(a::_BasisOperations) = coalesce(geometry(a.el1),geometry(a.el2)) basisfunction(a::_BasisOperations,i) = coalesce(basisfunction(a.el1,i),basisfunction(a.el2,i)) positions(a::_BasisOperations) = coalesce(positions(a.el1),positions(a.el2)) -geometry(a::Union{NormalVector,FunctionWrapper,Function}) = missing -basisfunction(a::Union{NormalVector,FunctionWrapper,Function},i) = missing -positions(a::Union{NormalVector,FunctionWrapper,Function}) = missing +geometry(a::Union{NormalVector,FunctionWrapper}) = missing +basisfunction(a::Union{NormalVector,FunctionWrapper},i) = missing +positions(a::Union{NormalVector,FunctionWrapper}) = missing subset(a::T,I) where {T <: _BasisOperations} = T(subset(a.el1,I),subset(a.el2,I)) -subset(a::Union{NormalVector,FunctionWrapper,Function},I) = a +subset(a::Union{NormalVector,FunctionWrapper},I) = a + +#### mathematical support will be called if there are no dedicated routines +×(n::NormalVector,s::Space) = _BasisCross(n,s) +×(s::Space,n::NormalVector) = _BasisCross(s,n) + +⋅(n::NormalVector,s::Space) = _BasisDot(n,s) +⋅(s::Space,n::NormalVector) = _BasisDot(s,n) + +*(n::NormalVector,s::Space) = _BasisTimes(n,s) +*(s::Space,n::NormalVector) = _BasisTimes(s,n) + + + +# function restrict(a::_LocalBasisOperations,cell1,cell2) +# s = sign(dot(normal(cell1),normal(cell2)))^number_of_normals(a) +# return s*restrict(inner_refspace(a),cell1,cell2) +# end + + +# number_of_normals(a::NormalVector) = 1 +# number_of_normals(a) = 0 +# number_of_normals(a::_LocalBasisOperations) = number_of_normals(a.el1) + number_of_normals(a.el2) +# #TODO save wrapped space in type, same with refspace. +# inner_refspace(a::) \ No newline at end of file diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl index 7b169bec..45c3718d 100644 --- a/src/bases/local/localcomposedbasis.jl +++ b/src/bases/local/localcomposedbasis.jl @@ -1,4 +1,4 @@ -abstract type _LocalBasisOperations{T} <: RefSpace{T,:None} end +abstract type _LocalBasisOperations{T} <: RefSpace{T} end numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) struct _LocalBasisTimes{T,U,V} <: _LocalBasisOperations{T} @@ -37,7 +37,7 @@ end operation(a::_LocalBasisTimes) = * operation(a::_LocalBasisCross) = × -operation(a::_LocalBasisDot) = ⋅ +operation(a::_LocalBasisDot) = (x,y) --> transpose(x)*y _execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector,op::U) where {N,U} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) _execute_operation(el1::SVector{N,<:NamedTuple},el2::U,op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) @@ -45,7 +45,7 @@ _execute_operation(el1::SVector,el2::SVector{N,<:NamedTuple},op::U) where {N,U} _execute_operation(el1::U,el2::SVector{N,<:NamedTuple},op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) _execute_operation(el1::SVector,el2::SVector,op::U) where {U <: Union{_LocalBasisCross,_LocalBasisDot}} = operation(op)(el1,el2) _execute_operation(el1::U,el2::V,op::_LocalBasisTimes) where {U <: Number,V <: Number}= el1*el2 -_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not (yet) supported" +_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not supported" _execute_operation_named(a::NTuple{N},b,op) where {N} = ((value=op(a[1].value,b),), _execute_operation_named(Base.tail(a),b,op)...) _execute_operation_named(a::NTuple{1},b,op) = ((value=op(a[1].value,b),),) diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl new file mode 100644 index 00000000..04b91e0e --- /dev/null +++ b/src/composedoperators/composedoperator.jl @@ -0,0 +1,413 @@ + +abstract type AbstractCombInt <: AbstractOperator end +abstract type Kernel{T} end +# abstract type _AbstractCombInt <: AbstractOperator end +#Make these structures, those are not yet pulled trough the smartmath function +struct CompDoubleInt{T} <: AbstractCombInt + tfunc + pairing1 + kernel::Kernel{T} + pairing2 + bfunc +end + +struct CompSingleInt{T} <: AbstractCombInt + tfunc + pairing1 + kernel::Kernel{T} + pairing2 + bfunc +end +scalartype(op::CompDoubleInt{T}) where {T} = T +scalartype(op::CompSingleInt{T}) where {T} = T + +#these structures are pulled trough the smartmath function +# struct _CompDoubleInt{T,K} <: _AbstractCombInt +# tfunc +# pairing1 +# kernel::K +# pairing2 +# bfunc +# end + +# struct _CompSingleInt{T,K} <: _AbstractCombInt +# tfunc +# pairing1 +# kernel::K +# pairing2 +# bfunc +# end + + +#the functionallity is pased to the basis at this point. +struct CompDoubleKern{T,U,V,W} <: IntegralOperator + pairing1::U + kernel::W + pairing2::V +end +function CompDoubleKern(p1,k::Kernel{T},p2) where {T} + CompDoubleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2) +end +struct CompSingleKern{T,U,V,W} <: LocalOperator + pairing1::U + kernel::W + pairing2::V +end +function CompSingleKern(p1,k::Kernel{T},p2) where {T} + CompSingleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2) +end +scalartype(op::CompDoubleKern{T}) where {T} = T +scalartype(op::CompSingleKern{T}) where {T} = T + +integralop(a::CompDoubleInt) = CompDoubleKern(a.pairing1,a.kernel,a.pairing2) +integralop(a::CompSingleInt) = CompSingleKern(a.pairing1,a.kernel,a.pairing2) + +### can be made more complicated later on +# function assemble!(op::AbstractCombInt, tfs::AbstractSpace, bfs::AbstractSpace, +# store, threading = Threading{:multi}; +# quadstrat=defaultquadstrat(op, tfs, bfs)) + +# return assemble!(smartmath(op), tfs,bfs, store, threading;quadstrat) +# end + +function assemble!(op::AbstractCombInt, tfs::AbstractSpace, bfs::AbstractSpace, + store, threading = Threading{:multi}; + quadstrat=defaultquadstrat(op, tfs, bfs)) + + assemble!(integralop(op),op.tfunc(tfs),op.bfunc(bfs),store,threading;quadstrat) + +end + +############################################################################################## mathematical support to build the operators +### not nesesarry, build potential structure yourself. + +# ######## support for symbolic basis +# #TODO add differentiation tools +# #Number of Rows: R +# #Number of columns : C +# abstract type SymbolicExpression{R,C} end +# struct SymbBasisFunction{R,1} <: SymbolicExpression{R,C} end +# SymbBasisFunction(i::Int) = SymbBasisFunction{i}() + +# struct SymbFunctionWrapper{R,C} <: SymbolicExpression{R,C} +# func +# end +# #TODO change this if allowed to the actual beast normalvector +# struct SymbolicNormal{3,1} <: SymbolicExpression{R,C} end + +# struct SymbMult{R,C} <: SymbolicExpression{R,C} +# el1 +# el2 +# function SymbMult(el1::SymbolicExpression{1,1},el2::SymbolicExpression{1,1}) +# return new{1,1}(el1,el2) +# end +# function SymbMult(el1::SymbolicExpression{1,1},el2::SymbolicExpression{R,C}) where {R,C} +# return new{R,C}(el1,el2) +# end +# function SymbMult(el1::SymbolicExpression{R,C},el2::SymbolicExpression{1,1}) where {R,C} +# return new{R,C}(el1,el2) +# end + +# end + +# struct SymbCross{3,1} <: SymbolicExpression{R,C} +# el1::SymbolicExpression{3,1} +# el2::SymbolicExpression{3,1} +# end +# struct SymbDot{R,C} <: SymbolicExpression{R,C} +# el1 +# el2 +# end + +# ### derivatives strongly restricted to the basis-function. This can be extended to SymbMult,... in the future. +# struct Symb3DCurl{3,1} <: SymbolicExpression{R,C} +# el::SymbBasisFunction +# end +# struct Symb2DCurl{R,C} <: SymbolicExpression{R,C} +# el::SymbBasisFunction +# function Symb2DCurl(el::SymbBasisFunction{3,1}) +# return new{1,1}(el) +# end +# function Symb2DCurl(el::SymbBasisFunction{1,1}) +# return new{3,1}(el) +# end +# end +# struct SymbDiv{1,1} <: SymbolicExpression{R,C} +# el::SymbBasisFunction +# end +# struct SymbGrad{3,1} <: SymbolicExpression{R,C} +# el::SymbBasisFunction +# end + +# # SymbMult(a::Function,b::Function) = SymbMult(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) +# # SymbMult(a::Function,b) = SymbMult(SymbFunctionWrapper(a),b) +# # SymbMult(a,b::Function) = SymbMult(a,SymbFunctionWrapper(b)) + +# # SymbCross(a::Function,b::Function) = SymbCross(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) +# # SymbCross(a::Function,b) = SymbCross(SymbFunctionWrapper(a),b) +# # SymbCross(a,b::Function) = SymbCross(a,SymbFunctionWrapper(b)) + +# # SymbDot(a::Function,b::Function) = SymbDot(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) +# # SymbDot(a::Function,b) = SymbDot(SymbFunctionWrapper(a),b) +# # SymbDot(a,b::Function) = SymbDot(a,SymbFunctionWrapper(b)) + + +# ##### casting symbolic representation to the actual basis +# ## for now the structures are used, when te user end of composed basis is finished it is better to use that, than a connection with the original basis functions is made +# function (a::SymbMult)(s::Space) +# BasisTimes(a.el1(s),a.el2(s)) +# end +# function (a::SymbCross)(s::Space) +# BasisCross(a.el1(s),a.el2(s)) +# end +# function (a::SymbDot)(s::Space) +# BasisDot(a.el1(s),a.el2(s)) +# end +# function (a::Symb3DCurl)(s::Space) +# return curl(a.el(s)) +# end +# function (a::Symb2DCurl)(s::Space) +# return curl(a.el(s)) +# end +# function (a::SymbDiv)(s::Space) +# return divergence(a.el(s)) +# end +# function (a::SymbGrad)(s::Space) +# return gradient(a.el(s)) +# end + +# function (a::SymbFunctionWrapper)(s::Space) +# return FunctionWrapper(a.func) +# end +# function (a::NormalVector)(s::Space) +# return a +# end +# function (a::SymbBasisFunction)(s::Space) +# return s +# end + +# ##### user build tools for Compbasis +# const B = SymbBasisFunction() +# #SymbolicElements = Union{SymbFunctionWrapper,NormalVector,SymbMult,SymbCross,SymbDot} +# *(a::SymbolicExpression,b::SymbolicExpression) = SymbMult(a,b) + + + +####### kernel definition (also posibilities to do symbolic actions) + + +#TODO fast computation of identity kernel +struct IdentityKernel <: Kernel{Int} end +struct NullKernel <: Kernel{Int} end +struct NxDotNy{T} <: Kernel{T} + mult::Kernel{T} +end +NxDotNy() = NxDotNy(IdentityKernel()) +NxDotNy(a::NxDotNy) = @error "dont do nx⋅ny*nx⋅ny (yet)" + +struct HH3DGreen{T} <: Kernel{T} + gamma::T +end +struct HH3DGradGreen{T} <: Kernel{T} + gamma::T +end + +function (op::HH3DGreen)(x,y) + gamma = op.gamma + + r = cartesian(x) - cartesian(y) + R = norm(r) + iR = 1/R + green = exp(-gamma*R)*(i4pi*iR) + green +end +function (op::HH3DGreen)(x::Union{SVector,Vector},y) + gamma = op.gamma + + r = x - cartesian(y) + R = norm(r) + iR = 1/R + green = exp(-gamma*R)*(i4pi*iR) + green +end + +function (op::HH3DGradGreen)(x,y) + gamma = op.gamma + + r = cartesian(x) - cartesian(y) + R = norm(r) + iR = 1/R + green = exp(-gamma*R)*(iR*i4pi) + gradgreen = -(gamma + iR) * green * (iR * r) + + return gradgreen +end +function (op::HH3DGradGreen)(x::Union{SVector,Vector},y) + gamma = op.gamma + + r = x - cartesian(y) + R = norm(r) + iR = 1/R + green = exp(-gamma*R)*(iR*i4pi) + gradgreen = -(gamma + iR) * green * (iR * r) + tt = gradgreen + + return tt +end + +##### integrand evaluation +function integrand(op::CompSingleKern,kerneldata,x,f,g,disp) + kernel = op.kernel(x,disp) + op1 = op.pairing1 + op2 = op.pairing2 + + op1(f.value,op2(kernel,g.value)) + +end + +function (op::Integrand{<:CompDoubleKern})(x,y,f,g) + kernel = op.operator.kernel(x,y) + op1 = op.operator.pairing1 + op2 = op.operator.pairing2 + _integrands(f,g) do fi, gi + op1(fi.value,op2(kernel,gi.value)) + end +end + +# _pairing(op1,op2,test::SVector{N,<:NamedTuple},trial::SVector{M,<:NamedTuple},kernel) where {N,M} = SMatrix{N,M}(_pairing(op1,op2,test.data,trial.data,kernel)) +# _pairing(op1,op2,test::SVector{N,<:NamedTuple},trial::SVector{M,<:NamedTuple},kernel::SVector) where {N,M} = SMatrix{N,M}(_pairing(op1,op2,test.data,trial.data,kernel)) + +# _pairing(op1,op2,test::NTuple{N},trial::NTuple{M},kernel) where {N,M} = (_pairing(op1,op2,test,(trial[1],),kernel)...,_pairing(op1,op2,test,Base.tail(trial),kernel)...) +# _pairing(op1,op2,test::NTuple{N},trial::NTuple{1},kernel) where {N} = (_pairing(op1,op2,(test[1],),trial,kernel)...,_pairing(op1,op2,Base.tail(test),trial,kernel)...) +# function _pairing(op1,op2,test::NTuple{1},trial::NTuple{1},kernel) +# # t = op2(kernel,trial[1].value) +# # u = op1(test[1].value,t) +# t = kernel*trial[1].value +# u = test[1].value*t +# return (u,) +# #(op1(test[1].value,op2(kernel,trial[1].value)),) +# end + + + +# function _pairing(op1,op2,test,trial,kern) +# display(op1) +# display(op2) +# display(test) +# display(trial) +# display(kern) +# @assert false +# end +defaultquadstrat(op::CompDoubleInt,tfs::Space,bfs::Space) = DoubleNumSauterQstrat(3,3,3,3,3,3) +defaultquadstrat(op::CompSingleInt,tfs::Space,bfs::Space) = SingleNumQStrat(3) +##### assembling the single integral +function assemble!(biop::CompSingleKern, tfs::Space, bfs::Space, store, + threading::Type{Threading{:multi}}; + quadstrat=defaultquadstrat(biop, tfs, bfs)) + + return assemble_local_mixed!(biop, tfs, bfs, store; quadstrat) +end + +function assemble_local_mixed!(biop::CompSingleKern, tfs::Space{T}, bfs::Space{T}, store; + quadstrat=defaultquadstrat(biop, tfs, bfs)) where {T} + + tol = sqrt(eps(T)) + + trefs = refspace(tfs) + brefs = refspace(bfs) + + tels, tad = assemblydata(tfs) + bels, bad = assemblydata(bfs) + + tgeo = geometry(tfs) + bgeo = geometry(bfs) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + num_trefs = numfunctions(trefs, tdom) + num_brefs = numfunctions(brefs, bdom) + + qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat) + + # store the bcells in an octree + tree = elementstree(bels) + + print("dots out of 10: ") + todo, done, pctg = length(tels), 0, 0 + for (p,tcell) in enumerate(tels) + + tc, ts = boundingbox(tcell.vertices) + pred = (c,s) -> boxesoverlap(c,s,tc,ts) + + for box in boxes(tree, pred) + for q in box + bcell = bels[q] + + if overlap(tcell, bcell) + + isct = intersection(tcell, bcell) + disp = displacement(displacementchart(geometry(tfs),p),displacementchart(geometry(bfs),q)) #number to multiply with test normal, if zero use the orientation of displacementvector + for cell in isct + volume(cell) < tol && continue + + P = restrict(brefs, bcell, cell) + Q = restrict(trefs, tcell, cell) + + qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat) + zlocal = cellinteractions(biop, trefs, brefs, cell, qr, disp) + zlocal = Q * zlocal * P' + + for i in 1 : num_trefs + for j in 1 : num_brefs + for (m,a) in tad[p,i] + for (n,b) in bad[q,j] + store(a * zlocal[i,j] * b, m, n) + end # next basis function this cell supports + end # next test function this cell supports + end # next refshape on basis side + end # next refshape on test side + + end # next cell in intersection + end # if overlap + end # next cell in the basis geometry + end # next box in the octree + + done += 1 + new_pctg = round(Int, done / todo * 100) + if new_pctg > pctg + 9 + print(".") + pctg = new_pctg + end + end # next cell in the test geometry + + println("") +end +kernelvals(::CompSingleKern,x) = nothing +function cellinteractions(biop::CompSingleKern, trefs::U, brefs::V, cell, qr, disp) where {U<:RefSpace{T},V<:RefSpace{T}} where {T} + + num_tshs = length(qr[1][3]) + num_bshs = length(qr[1][4]) + + zlocal = zeros(T, num_tshs, num_bshs) + for q in qr + + w, mp, tvals, bvals = q[1], q[2], q[3], q[4] + j = w * jacobian(mp) + kernel = kernelvals(biop, mp) + + for m in 1 : num_tshs + tval = tvals[m] + + for n in 1 : num_bshs + bval = bvals[n] + + igd = integrand(biop, kernel, mp, tval, bval,disp) + zlocal[m,n] += j * igd + + end + end + end + + return zlocal +end diff --git a/src/composedoperators/displacementmesh.jl b/src/composedoperators/displacementmesh.jl new file mode 100644 index 00000000..07174043 --- /dev/null +++ b/src/composedoperators/displacementmesh.jl @@ -0,0 +1,29 @@ +##### displacement meshes +abstract type DisplacementMesh{T} <: CompScienceMeshes.AbstractMesh{3,3,T} end + +""" +GlobalDisplacementMesh(mesh,epsilon) + every chart is displaced with epsilon allong the normal +""" +struct GlobalDisplacementMesh{T} <: DisplacementMesh{T} + mesh::Mesh{<:Any,<:Any,T} + epsilon +end +struct DisplacementChart + chart + epsilon +end + +displacementchart(m::GlobalDisplacementMesh,i) = DisplacementChart(chart(m.mesh,i),m.epsilon) +displacementchart(m::Mesh,i) = DisplacementChart(chart(m,i),-1.0) +CompScienceMeshes.chart(m::DisplacementMesh,i) = chart(m.mesh,i) +function displacement(a::DisplacementChart,b::DisplacementChart) + rel_orient = sign(dot(normal(a.chart),normal(b.chart))) + if a.epsilon*rel_orient ≈ b.epsilon + return 0 + end + testdisp = a.epsilon*normal(a.chart) + trialdisp = b.epsilon*normal(b.chart) + + return sign(dot(trialdisp-testdisp,normal(a.chart))) +end \ No newline at end of file diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl new file mode 100644 index 00000000..2ad45181 --- /dev/null +++ b/src/composedoperators/potentials.jl @@ -0,0 +1,12 @@ +#D: dimension of manifold over which integration is performed +struct PotentialIntegralOperator{D} + kernel + pairing2 + bfunc +end + +function potential(a::PotentialIntegralOperator,pts,u,s::Space) + +end + + diff --git a/src/composedoperators/smartmath.jl b/src/composedoperators/smartmath.jl new file mode 100644 index 00000000..e90828dd --- /dev/null +++ b/src/composedoperators/smartmath.jl @@ -0,0 +1,10 @@ +""" +support that can be called on double or single composed integral operator. The function can make connections with the rest of BEAST, simplify the results etc +For now it just returns the original structure. +""" +smartmath(x) = x + +smartmath(x::CompDoubleInt) = _CompDoubleInt(x.tfunc,t.pairing1,t.kernel,t.pairing2,t.bfunc) +smartmath(x::CompSingleInt) = _CompSingleInt(x.tfunc,t.pairing1,t.kernel,t.pairing2,t.bfunc) + +#TODO no smartmath, remove this section we assume the user is "smart" \ No newline at end of file diff --git a/src/composedoperators/trace.jl b/src/composedoperators/trace.jl new file mode 100644 index 00000000..4d937483 --- /dev/null +++ b/src/composedoperators/trace.jl @@ -0,0 +1,123 @@ + +#### trace vector + +# orientation is used when test and trial chart have same displacement, in this case the test-cart is displaced in the orientation*normal(test-chart) direction. +struct DisplacementVector{T} <: Kernel{T} + orientation +end + + +function (u::DisplacementVector)(x,disp) + if disp == 0 + return sign(u.orientation)*normal(x) + else + return sign(disp)*normal(x) + end +end + + +#### trace of a potential: +_trace(::Kernel,n;type=Float64,D) = NullKernel(), 0 +function _trace(::HH3DGradGreen,n,D;type=Float64) + if D==2 + return DisplacementVector{type}(n) , 0.5 + end + return NullKernel(), 0 +end + +""" +sign indicates if trace is taken allong the normal (+1) are opposite to the normal (-1) +""" +function ttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=false) where {D} + t,f = _trace(a.kernel,sign,D) + if testfunction_tangential + doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + else + + doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + end + if t==NullKernel() + return doubleint + else + #eventualy a scan to be non zero can be added here. + return doubleint + singleint + end +end + +function ntrace(a::PotentialIntegralOperator{D},sign) where {D} + t,f = _trace(a.kernel,sign,D) + doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + + if t==NullKernel() + return doubleint + else + #eventualy a scan to be non zero can be added here. + return doubleint + singleint + end +end + +function trace(a::PotentialIntegralOperator{D},sign) where {D} + t,f = _trace(a.kernel,sign,D) + doubleint = CompDoubleInt(B->B,*,a.kernel,a.pairing2,a.bfunc) + singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc) + + if t==NullKernel() + return doubleint + else + #eventualy a scan to be non zero can be added here. + return doubleint + singleint + end +end + +function strace(a::PotentialIntegralOperator{D},sign) where {D} + t,f = _trace(a.kernel,sign,D) + doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + + if t==NullKernel() + return doubleint + else + #eventualy a scan to be non zero can be added here. + return doubleint + singleint + end +end + +# the principal value of the trace. +function pvttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=false) where {D} + if testfunction_tangential + doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + else + doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + end + return doubleint +end + +function pvntrace(a::PotentialIntegralOperator{D},sign) where {D} + + doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + + return doubleint +end + +function pvtrace(a::PotentialIntegralOperator{D},sign) where {D} + doubleint = CompDoubleInt(B-> B,*,a.kernel,a.pairing2,a.bfunc) + + return doubleint +end + +function pvstrace(a::PotentialIntegralOperator{D},sign) where {D} + doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + return doubleint +end + +function volumetrace(a::PotentialIntegralOperator{D}) where {D} + CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) +end + + + + + From 4010c905169706a026f77e9fb53b97691813ed63 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Mon, 4 Nov 2024 17:21:08 +0100 Subject: [PATCH 04/17] - extended excitation posibility - cleanup --- src/BEAST.jl | 1 + src/composedoperators/analytic_excitation.jl | 23 +++ src/composedoperators/composedoperator.jl | 188 ++----------------- src/composedoperators/potentials.jl | 12 +- src/composedoperators/smartmath.jl | 10 - src/composedoperators/trace.jl | 22 +-- 6 files changed, 58 insertions(+), 198 deletions(-) create mode 100644 src/composedoperators/analytic_excitation.jl delete mode 100644 src/composedoperators/smartmath.jl diff --git a/src/BEAST.jl b/src/BEAST.jl index a0013360..5425ad7d 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -293,6 +293,7 @@ include("composedoperators/composedoperator.jl") include("composedoperators/displacementmesh.jl") include("composedoperators/potentials.jl") include("composedoperators/trace.jl") +include("composedoperators/analytic_excitation.jl") const x̂ = point(1, 0, 0) const ŷ = point(0, 1, 0) diff --git a/src/composedoperators/analytic_excitation.jl b/src/composedoperators/analytic_excitation.jl new file mode 100644 index 00000000..ada5e7b1 --- /dev/null +++ b/src/composedoperators/analytic_excitation.jl @@ -0,0 +1,23 @@ +import Base.* +struct NTimesTrace{T,F} <: Functional + field::F +end + +NTimesTrace(f::F) where {F} = NTimesTrace{scalartype(f), F}(f) +NTimesTrace{T}(f::F) where {T,F} = NTimesTrace{T,F}(f) +scalartype(s::NTimesTrace{T}) where {T} = T + +(ϕ::NTimesTrace)(p) = *(normal(p), ϕ.field(cartesian(p))) +integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) +*(::NormalVector, f) = NTimesTrace(f) + +struct Trace{T,F} <: Functional + field::F +end + +Trace(f::F) where {F} = Trace{scalartype(f), F}(f) +Trace{T}(f::F) where {T,F} = Trace{T,F}(f) +scalartype(s::Trace{T}) where {T} = T + +(ϕ::Trace)(p) = ϕ.field(cartesian(p)) +integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl index 04b91e0e..3c3fe5da 100644 --- a/src/composedoperators/composedoperator.jl +++ b/src/composedoperators/composedoperator.jl @@ -1,9 +1,9 @@ -abstract type AbstractCombInt <: AbstractOperator end +abstract type AbstractCompInt <: AbstractOperator end abstract type Kernel{T} end -# abstract type _AbstractCombInt <: AbstractOperator end +# abstract type _AbstractCompInt <: AbstractOperator end #Make these structures, those are not yet pulled trough the smartmath function -struct CompDoubleInt{T} <: AbstractCombInt +struct CompDoubleInt{T} <: AbstractCompInt tfunc pairing1 kernel::Kernel{T} @@ -11,7 +11,7 @@ struct CompDoubleInt{T} <: AbstractCombInt bfunc end -struct CompSingleInt{T} <: AbstractCombInt +struct CompSingleInt{T} <: AbstractCompInt tfunc pairing1 kernel::Kernel{T} @@ -21,23 +21,6 @@ end scalartype(op::CompDoubleInt{T}) where {T} = T scalartype(op::CompSingleInt{T}) where {T} = T -#these structures are pulled trough the smartmath function -# struct _CompDoubleInt{T,K} <: _AbstractCombInt -# tfunc -# pairing1 -# kernel::K -# pairing2 -# bfunc -# end - -# struct _CompSingleInt{T,K} <: _AbstractCombInt -# tfunc -# pairing1 -# kernel::K -# pairing2 -# bfunc -# end - #the functionallity is pased to the basis at this point. struct CompDoubleKern{T,U,V,W} <: IntegralOperator @@ -62,15 +45,8 @@ scalartype(op::CompSingleKern{T}) where {T} = T integralop(a::CompDoubleInt) = CompDoubleKern(a.pairing1,a.kernel,a.pairing2) integralop(a::CompSingleInt) = CompSingleKern(a.pairing1,a.kernel,a.pairing2) -### can be made more complicated later on -# function assemble!(op::AbstractCombInt, tfs::AbstractSpace, bfs::AbstractSpace, -# store, threading = Threading{:multi}; -# quadstrat=defaultquadstrat(op, tfs, bfs)) - -# return assemble!(smartmath(op), tfs,bfs, store, threading;quadstrat) -# end -function assemble!(op::AbstractCombInt, tfs::AbstractSpace, bfs::AbstractSpace, +function assemble!(op::AbstractCompInt, tfs::AbstractSpace, bfs::AbstractSpace, store, threading = Threading{:multi}; quadstrat=defaultquadstrat(op, tfs, bfs)) @@ -78,132 +54,16 @@ function assemble!(op::AbstractCombInt, tfs::AbstractSpace, bfs::AbstractSpace, end -############################################################################################## mathematical support to build the operators -### not nesesarry, build potential structure yourself. - -# ######## support for symbolic basis -# #TODO add differentiation tools -# #Number of Rows: R -# #Number of columns : C -# abstract type SymbolicExpression{R,C} end -# struct SymbBasisFunction{R,1} <: SymbolicExpression{R,C} end -# SymbBasisFunction(i::Int) = SymbBasisFunction{i}() - -# struct SymbFunctionWrapper{R,C} <: SymbolicExpression{R,C} -# func -# end -# #TODO change this if allowed to the actual beast normalvector -# struct SymbolicNormal{3,1} <: SymbolicExpression{R,C} end - -# struct SymbMult{R,C} <: SymbolicExpression{R,C} -# el1 -# el2 -# function SymbMult(el1::SymbolicExpression{1,1},el2::SymbolicExpression{1,1}) -# return new{1,1}(el1,el2) -# end -# function SymbMult(el1::SymbolicExpression{1,1},el2::SymbolicExpression{R,C}) where {R,C} -# return new{R,C}(el1,el2) -# end -# function SymbMult(el1::SymbolicExpression{R,C},el2::SymbolicExpression{1,1}) where {R,C} -# return new{R,C}(el1,el2) -# end - -# end - -# struct SymbCross{3,1} <: SymbolicExpression{R,C} -# el1::SymbolicExpression{3,1} -# el2::SymbolicExpression{3,1} -# end -# struct SymbDot{R,C} <: SymbolicExpression{R,C} -# el1 -# el2 -# end - -# ### derivatives strongly restricted to the basis-function. This can be extended to SymbMult,... in the future. -# struct Symb3DCurl{3,1} <: SymbolicExpression{R,C} -# el::SymbBasisFunction -# end -# struct Symb2DCurl{R,C} <: SymbolicExpression{R,C} -# el::SymbBasisFunction -# function Symb2DCurl(el::SymbBasisFunction{3,1}) -# return new{1,1}(el) -# end -# function Symb2DCurl(el::SymbBasisFunction{1,1}) -# return new{3,1}(el) -# end -# end -# struct SymbDiv{1,1} <: SymbolicExpression{R,C} -# el::SymbBasisFunction -# end -# struct SymbGrad{3,1} <: SymbolicExpression{R,C} -# el::SymbBasisFunction -# end - -# # SymbMult(a::Function,b::Function) = SymbMult(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) -# # SymbMult(a::Function,b) = SymbMult(SymbFunctionWrapper(a),b) -# # SymbMult(a,b::Function) = SymbMult(a,SymbFunctionWrapper(b)) - -# # SymbCross(a::Function,b::Function) = SymbCross(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) -# # SymbCross(a::Function,b) = SymbCross(SymbFunctionWrapper(a),b) -# # SymbCross(a,b::Function) = SymbCross(a,SymbFunctionWrapper(b)) - -# # SymbDot(a::Function,b::Function) = SymbDot(SymbFunctionWrapper(a),SymbFunctionWrapper(b)) -# # SymbDot(a::Function,b) = SymbDot(SymbFunctionWrapper(a),b) -# # SymbDot(a,b::Function) = SymbDot(a,SymbFunctionWrapper(b)) - - -# ##### casting symbolic representation to the actual basis -# ## for now the structures are used, when te user end of composed basis is finished it is better to use that, than a connection with the original basis functions is made -# function (a::SymbMult)(s::Space) -# BasisTimes(a.el1(s),a.el2(s)) -# end -# function (a::SymbCross)(s::Space) -# BasisCross(a.el1(s),a.el2(s)) -# end -# function (a::SymbDot)(s::Space) -# BasisDot(a.el1(s),a.el2(s)) -# end -# function (a::Symb3DCurl)(s::Space) -# return curl(a.el(s)) -# end -# function (a::Symb2DCurl)(s::Space) -# return curl(a.el(s)) -# end -# function (a::SymbDiv)(s::Space) -# return divergence(a.el(s)) -# end -# function (a::SymbGrad)(s::Space) -# return gradient(a.el(s)) -# end - -# function (a::SymbFunctionWrapper)(s::Space) -# return FunctionWrapper(a.func) -# end -# function (a::NormalVector)(s::Space) -# return a -# end -# function (a::SymbBasisFunction)(s::Space) -# return s -# end - -# ##### user build tools for Compbasis -# const B = SymbBasisFunction() -# #SymbolicElements = Union{SymbFunctionWrapper,NormalVector,SymbMult,SymbCross,SymbDot} -# *(a::SymbolicExpression,b::SymbolicExpression) = SymbMult(a,b) - ####### kernel definition (also posibilities to do symbolic actions) - -#TODO fast computation of identity kernel -struct IdentityKernel <: Kernel{Int} end -struct NullKernel <: Kernel{Int} end -struct NxDotNy{T} <: Kernel{T} - mult::Kernel{T} -end -NxDotNy() = NxDotNy(IdentityKernel()) -NxDotNy(a::NxDotNy) = @error "dont do nx⋅ny*nx⋅ny (yet)" +# struct NullKernel <: Kernel{Int} end +# struct NxDotNy{T} <: Kernel{T} +# mult::Kernel{T} +# end +# NxDotNy() = NxDotNy(IdentityKernel()) +# NxDotNy(a::NxDotNy) = @error "dont do nx⋅ny*nx⋅ny (yet)" struct HH3DGreen{T} <: Kernel{T} gamma::T @@ -274,30 +134,6 @@ function (op::Integrand{<:CompDoubleKern})(x,y,f,g) end end -# _pairing(op1,op2,test::SVector{N,<:NamedTuple},trial::SVector{M,<:NamedTuple},kernel) where {N,M} = SMatrix{N,M}(_pairing(op1,op2,test.data,trial.data,kernel)) -# _pairing(op1,op2,test::SVector{N,<:NamedTuple},trial::SVector{M,<:NamedTuple},kernel::SVector) where {N,M} = SMatrix{N,M}(_pairing(op1,op2,test.data,trial.data,kernel)) - -# _pairing(op1,op2,test::NTuple{N},trial::NTuple{M},kernel) where {N,M} = (_pairing(op1,op2,test,(trial[1],),kernel)...,_pairing(op1,op2,test,Base.tail(trial),kernel)...) -# _pairing(op1,op2,test::NTuple{N},trial::NTuple{1},kernel) where {N} = (_pairing(op1,op2,(test[1],),trial,kernel)...,_pairing(op1,op2,Base.tail(test),trial,kernel)...) -# function _pairing(op1,op2,test::NTuple{1},trial::NTuple{1},kernel) -# # t = op2(kernel,trial[1].value) -# # u = op1(test[1].value,t) -# t = kernel*trial[1].value -# u = test[1].value*t -# return (u,) -# #(op1(test[1].value,op2(kernel,trial[1].value)),) -# end - - - -# function _pairing(op1,op2,test,trial,kern) -# display(op1) -# display(op2) -# display(test) -# display(trial) -# display(kern) -# @assert false -# end defaultquadstrat(op::CompDoubleInt,tfs::Space,bfs::Space) = DoubleNumSauterQstrat(3,3,3,3,3,3) defaultquadstrat(op::CompSingleInt,tfs::Space,bfs::Space) = SingleNumQStrat(3) ##### assembling the single integral @@ -318,7 +154,7 @@ function assemble_local_mixed!(biop::CompSingleKern, tfs::Space{T}, bfs::Space{T tels, tad = assemblydata(tfs) bels, bad = assemblydata(bfs) - + tgeo = geometry(tfs) bgeo = geometry(bfs) diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl index 2ad45181..09b0934f 100644 --- a/src/composedoperators/potentials.jl +++ b/src/composedoperators/potentials.jl @@ -5,8 +5,18 @@ struct PotentialIntegralOperator{D} bfunc end -function potential(a::PotentialIntegralOperator,pts,u,s::Space) +struct PotentialIntegralOperatorKern{U,V} + kernel::U + pairing2::V +end +function integrand(a::PotentialIntegralOperatorKern,y,krn,f,x) + return a.pairing2(a.kernel(y,x),f) end +function potential(op::PotentialIntegralOperator, points, coeffs, basis; + type=SVector{3,ComplexF64}, + quadstrat=defaultquadstrat(op, basis)) + return potential(PotentialIntegralOperatorKern(op.kernel,op.pairing2),points,coeffs,op.bfunc(basis);type,quadstrat) +end diff --git a/src/composedoperators/smartmath.jl b/src/composedoperators/smartmath.jl deleted file mode 100644 index e90828dd..00000000 --- a/src/composedoperators/smartmath.jl +++ /dev/null @@ -1,10 +0,0 @@ -""" -support that can be called on double or single composed integral operator. The function can make connections with the rest of BEAST, simplify the results etc -For now it just returns the original structure. -""" -smartmath(x) = x - -smartmath(x::CompDoubleInt) = _CompDoubleInt(x.tfunc,t.pairing1,t.kernel,t.pairing2,t.bfunc) -smartmath(x::CompSingleInt) = _CompSingleInt(x.tfunc,t.pairing1,t.kernel,t.pairing2,t.bfunc) - -#TODO no smartmath, remove this section we assume the user is "smart" \ No newline at end of file diff --git a/src/composedoperators/trace.jl b/src/composedoperators/trace.jl index 4d937483..b25fe668 100644 --- a/src/composedoperators/trace.jl +++ b/src/composedoperators/trace.jl @@ -17,12 +17,12 @@ end #### trace of a potential: -_trace(::Kernel,n;type=Float64,D) = NullKernel(), 0 +_trace(::Kernel,n;type=Float64,D) = 0, 0 function _trace(::HH3DGradGreen,n,D;type=Float64) if D==2 return DisplacementVector{type}(n) , 0.5 end - return NullKernel(), 0 + return 0, 0 end """ @@ -32,13 +32,13 @@ function ttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=fal t,f = _trace(a.kernel,sign,D) if testfunction_tangential doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) else doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) end - if t==NullKernel() + if t==0 return doubleint else #eventualy a scan to be non zero can be added here. @@ -49,9 +49,9 @@ end function ntrace(a::PotentialIntegralOperator{D},sign) where {D} t,f = _trace(a.kernel,sign,D) doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) - if t==NullKernel() + if t==0 return doubleint else #eventualy a scan to be non zero can be added here. @@ -62,9 +62,9 @@ end function trace(a::PotentialIntegralOperator{D},sign) where {D} t,f = _trace(a.kernel,sign,D) doubleint = CompDoubleInt(B->B,*,a.kernel,a.pairing2,a.bfunc) - singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc) + t!=0 && singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc) - if t==NullKernel() + if t==0 return doubleint else #eventualy a scan to be non zero can be added here. @@ -75,9 +75,9 @@ end function strace(a::PotentialIntegralOperator{D},sign) where {D} t,f = _trace(a.kernel,sign,D) doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) - if t==NullKernel() + if t==0 return doubleint else #eventualy a scan to be non zero can be added here. From 73f50297cdcbcf85f5de9e731c123a50fcc92be6 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 13 Nov 2024 14:10:04 +0100 Subject: [PATCH 05/17] update composed operator + vp multi -trace example --- examples/vp_global_multi_trace.jl | 194 +++++++++++++++++++ src/bases/composedbasis.jl | 25 ++- src/bases/local/localcomposedbasis.jl | 75 +++++-- src/bases/ndspace.jl | 16 +- src/bases/rtspace.jl | 2 +- src/composedoperators/analytic_excitation.jl | 22 +-- src/composedoperators/composedoperator.jl | 4 +- src/composedoperators/displacementmesh.jl | 15 +- src/composedoperators/trace.jl | 29 ++- src/identityop.jl | 26 +++ test/runtests.jl | 3 + test/test_analytic_excitation.jl | 41 ++++ test/test_composed_basis.jl | 18 +- test/test_composed_operator.jl | 53 +++++ 14 files changed, 452 insertions(+), 71 deletions(-) create mode 100644 examples/vp_global_multi_trace.jl create mode 100644 test/test_analytic_excitation.jl create mode 100644 test/test_composed_operator.jl diff --git a/examples/vp_global_multi_trace.jl b/examples/vp_global_multi_trace.jl new file mode 100644 index 00000000..6ee8e33b --- /dev/null +++ b/examples/vp_global_multi_trace.jl @@ -0,0 +1,194 @@ +using BEAST +using CompScienceMeshes +using LinearAlgebra +using SparseArrays +using PlotlyJS +using StaticArrays + +ω = 10.0^8*2*pi +ϵ0 = 8.854e-12 +μ0 = 4π*1e-7 +κ0 = ω*√(μ0*ϵ0) + +ϵr = [1,2,2] +μr = [2,1,2] + +h = 0.5 +fn = joinpath(@__DIR__, "assets/rectangular_torus.geo") +Γ11 = CompScienceMeshes.meshgeo(fn; dim=2, h) +Γ11 = Mesh([point(-z,y,x) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) +function onlys(a) + c = sort.(a) + b = [i for i in a if length(findall(==(sort(i)),c))==1] + return b +end + +#from paul to test something +# Γ12 = -Mesh([point(-x-4,y,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) +# Γ1 = weld(Γ11,Γ12;glueop=onlys) +# Γ2 = Mesh([point(-x-4,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) +# Γ3 = -Mesh([point(x,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) + +Γ12 = -Mesh([point(x,y,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) +Γ1 = weld(Γ11,Γ12;glueop=onlys) +Γ2 = Mesh([point(x,-y-4,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) +Γ3 = -Mesh([point(x,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) + + +Γ = [Γ1,Γ2,Γ3] # the meshes + +κ = [ω*√(μ0*ϵ0*ϵri*μri) for (ϵri,μri) in zip(ϵr,μr)] + +########RHS definition (lorentz gauge:) +x0 = @SVector [0.0,0.0,0.0] # used to check the type of the incidnet fields +A(x) = x[1]*sqrt(ϵ0*μ0)exp(-1im*κ0*x[3])* (@SVector [0.0,0.0,1.0]) +curlA(x) = -(@SVector [0.0,1.0,0.0])*sqrt(ϵ0*μ0)*exp(-1.0im*κ0 *x[3]) +phi(x) = x[1]*exp(-1.0im*κ0 *x[3]) +gradphi(x) = (@SVector [1.0,0.0,-1im*κ0*x[1]])*exp(-1.0im*κ0 *x[3]) + +#################################################### +# +# From here on do nothing +# +#################################################### +@hilbertspace a1 a2 +@hilbertspace b1 b2 + +#### potentials + operators +G0 = BEAST.HH3DGreen(1im*κ0) +G = [BEAST.HH3DGreen(1im*k) for k in κ] +dG0 = BEAST.HH3DGradGreen(1im*κ0) +dG = [BEAST.HH3DGradGreen(1im*k) for k in κ] + + U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) + U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0) + U22 = BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->n*b),-1.0) + +DDL0 = U11[a1,b1] + U12[a1,b2] + U22[a2,b2] + + U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] + U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for g in G] + U22 = [BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->n*b),1.0) for g in dG] + +DDL = [U11i[a1,b1] + er*mr*U12i[a1,b2] + U22i[a2,b2] for (U11i,U12i,U22i,er,mr) in zip(U11,U12,U22,ϵr,μr)] + + U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) + U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) + U22 = BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0) + +NDL0 = U11[a1,b1] + U21[a2,b1] + U22[a2,b2] + + U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] + U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G] + U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG] + +NDL = [U11i[a1,b1] + er*mr*U21i[a2,b1] + U22i[a2,b2] for (U11i,U21i,U22i,er,mr) in zip(U11,U21,U22,ϵr,μr)] + + U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) + U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0) + U21 = -BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->b),-1.0) + U22 = -κ0^2* BEAST.trace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) + +NDSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2] + + U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G] + U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG] + U21 = [-BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->b),1.0) for g in dG] + U22 = [-k^2*BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for (g,k) in zip(G,κ)] + +NDSL = [mr*U11i[a1,b1] + 1/er*U12i[a1,b2] + 1/er*U21i[a2,b1] + 1/(er^2*mr)*U22i[a2,b2] for (U11i,U12i,U21i,U22i,er,mr) in zip(U11,U12,U21,U22,ϵr,μr)] + + U11 = -κ0^2 * BEAST.pvstrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0)-BEAST.CompDoubleInt(B->divergence(n×B),*,G0,*,B->divergence(B)) + U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->n*b),-1.0) + U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) + U22 = BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0) + +DNSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2] + + U11 = [-k^2 * BEAST.pvstrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0)-BEAST.CompDoubleInt(B->divergence(n×B),*,g,*,B->divergence(B)) for (g,k) in zip(G,κ)] + U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->n*b),1.0) for g in dG] + U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] + U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for (g,k) in zip(G,κ)] + +DNSL = [1/mr*U11i[a1,b1] + er*U12i[a1,b2] + er*U21i[a2,b1] + er^2*mr*U22i[a2,b2] for (U11i,U12i,U21i,U22i,er,mr) in zip(U11,U12,U21,U22,ϵr,μr)] + +### definition per domain bilinear forms + +A0 = DNSL0[a1,b1] + NDSL0[a2,b2] + DDL0[a2,b1] + NDL0[a1,b2] +Ad = [DNSLi[a1,b1] + NDSLi[a2,b2] + DDLi[a2,b1] + NDLi[a1,b2] for (DNSLi,NDSLi,DDLi,NDLi) in zip(DNSL,NDSL,DDL,NDL)] + +p = BEAST.hilbertspace(:p, length(κ)) +q = BEAST.hilbertspace(:q, length(κ)) + +LHSA = A0[p,q] + BEAST.Variational.DirectProductKernel(Ad)[p,q] + +#################################################### +# +# RHS definition +# +#################################################### + +### wrapping incident field + +w_A = BEAST.FunctionWrapper{typeof(A(x0))}(A) +w_curlA = BEAST.FunctionWrapper{typeof(curlA(x0))}(curlA) +w_phi = BEAST.FunctionWrapper{typeof(phi(x0))}(phi) +w_gradphi = BEAST.FunctionWrapper{typeof(gradphi(x0))}(gradphi) + +### traces RHS + +t_nxA = BEAST.CrossTraceMW(w_A) +t_nxcurlA = BEAST.CrossTraceMW(w_curlA) +t_ndA = BEAST.NDotTrace(w_A) +t_div_rescA = BEAST.Trace(w_phi) + +t_d = t_nxA[a1] + -1im*κ0^2/ω*t_div_rescA[a2] +t_n = t_nxcurlA[a1] + t_ndA[a2] +RHSAi = t_d[a2] + t_n[a1] + +t_phi = BEAST.Trace(w_phi) +t_ngradphi = BEAST.NDotTrace(w_gradphi) + +RHSA = RHSAi[p] + +### Spaces +RT = raviartthomas.(Γ) +ND = Ref(n) .× RT +MND = Ref(n) .× (Ref(n) .× (Ref(n) .× RT)) +L0 = lagrangecxd0.(Γ) +L1 = lagrangec0d1.(Γ) + +Xd = [RTi×L1i for (RTi,L1i) in zip(RT,L1)] +Xn = [RTi×L0i for (RTi,L0i) in zip(RT,L0)] + +Yn = [NDi×L1i for (NDi,L1i) in zip(MND,L1)] +Yd = [NDi×L0i for (NDi,L0i) in zip(ND,L0)] + +Qₕ = [BEAST.DirectProductSpace([Xdi,Xni]) for (Xdi,Xni) in zip(Xd,Xn)] +Pₕ = [BEAST.DirectProductSpace([Yni,Ydi]) for (Ydi,Yni) in zip(Yd,Yn)] + +deq = BEAST.discretise(LHSA==RHSA, (p .∈ Pₕ)..., (q .∈ Qₕ)...) +# Z = assemble(LHSA,BEAST.DirectProductSpace(Pₕ),BEAST.DirectProductSpace(Qₕ)) +# B = assemble(RHSA,BEAST.DirectProductSpace(Pₕ)) +u = solve(deq) + +# ### to check aginst previous solution +# using JLD2 +# d = load("examples/0.5Z_#1.jld2")["dict"] +# #comparing against previous solution +# X_server = load("examples/RT.jld2")["dict"] + + +# resc = zeros(size(Z)) + LinearAlgebra.I + +# for i in 1:length(X_server) +# if X_server.fns[i] != RT[1].fns[i] +# resc[i,i] *= -1 +# end +# end + +# for i in 1:length(X_server) +# if X_server.fns[i] != RT[1].fns[i] +# resc[i+length(Xd[1]),i+length(Xd[1])] *= -1 +# end +# end diff --git a/src/bases/composedbasis.jl b/src/bases/composedbasis.jl index c6817193..3ba24c2b 100644 --- a/src/bases/composedbasis.jl +++ b/src/bases/composedbasis.jl @@ -12,7 +12,7 @@ function (f::FunctionWrapper{T})(x::CompScienceMeshes.MeshPointNM)::T where {T} return f.func(cartesian(x)) end scalartype(F::FunctionWrapper{T}) where {T} = eltype(T) -#scalartype(::NormalVector,s::Space) = scalartype(geometry(s)) +#scalartype(::NormalVector,s::Space) = vertextype(geometry(s)) # function scalartype(::NormalVector) # @warn "The scallartype of the NormalVector is set at Float32, if used in combination with Float64 basis or operator this is no problem, in the case of Float16 it is." # return Float32 @@ -24,24 +24,28 @@ struct _BasisTimes{T} <: _BasisOperations{T} el1 el2 end -_BasisTimes(el1::NormalVector,el2) = _BasisTimes{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) -_BasisTimes(el2,el1::NormalVector) = _BasisTimes{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) -_BasisTimes(el1,el2) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisTimes(el1::NormalVector,el2::Space) = _BasisTimes{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2) +_BasisTimes(el2::Space,el1::NormalVector) = _BasisTimes{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1) +_BasisTimes(el1::Space,el2::FunctionWrapper) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisTimes(el1::FunctionWrapper,el2::Space) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) struct _BasisCross{T} <: _BasisOperations{T} el1 el2 end -_BasisCross(el1::NormalVector,el2) = _BasisCross{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) -_BasisCross(el2,el1::NormalVector) = _BasisCross{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) -_BasisCross(el1,el2) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisCross(el1::NormalVector,el2::Space) = _BasisCross{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2) +_BasisCross(el2::Space,el1::NormalVector) = _BasisCross{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1) +_BasisCross(el1::FunctionWrapper,el2::Space) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisCross(el1::Space,el2::FunctionWrapper) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) + struct _BasisDot{T} <: _BasisOperations{T} el1 el2 end -_BasisDot(el1::NormalVector,el2) = _BasisDot{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el1,el2) -_BasisDot(el2,el1::NormalVector) = _BasisDot{promote_type(scalartype(geometry(el2)),scalartype(el2))}(el2,el1) -_BasisDot(el1,el2) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisDot(el1::NormalVector,el2::Space) = _BasisDot{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2) +_BasisDot(el2::Space,el1::NormalVector) = _BasisDot{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1) +_BasisDot(el1::Space,el2::FunctionWrapper) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) +_BasisDot(el1::FunctionWrapper,el2::Space) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) # #### wrapping of the functions # _BasisTimes(a::Function,b::Function) = _BasisTimes(FunctionWrapper(a),FunctionWrapper(b)) @@ -65,6 +69,7 @@ refspace(a::FunctionWrapper) = a refspace(a::NormalVector) = a numfunctions(a::Union{NormalVector,FunctionWrapper}) = missing +numfunctions(a::Union{NormalVector,FunctionWrapper},simp) = missing numfunctions(a::_BasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) geometry(a::_BasisOperations) = coalesce(geometry(a.el1),geometry(a.el2)) diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl index 45c3718d..772cd33f 100644 --- a/src/bases/local/localcomposedbasis.jl +++ b/src/bases/local/localcomposedbasis.jl @@ -1,6 +1,21 @@ abstract type _LocalBasisOperations{T} <: RefSpace{T} end numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2)) +numfunctions(a::_LocalBasisOperations,simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp)) +struct TriangleSupport end +struct TetraSupport end + +support(a::_LocalBasisOperations) = coalesce(support(a.el1),support(a.el2)) +support(a::NormalVector) = missing +support(a::FunctionWrapper) = missing +support(a::LagrangeRefSpace{T,N,3}) where {T,N} = TriangleSupport() +support(a::LagrangeRefSpace{T,N,4}) where {T,N} = TetraSupport() +restrict(a::_LocalBasisOperations,bcell,cell) = restrict(_refspace(a),bcell,cell) + +_refspace(a::_LocalBasisOperations) = coalesce(_refspace(a.el1) , _refspace(a.el2)) +_refspace(a::NormalVector) = missing +_refspace(a::FunctionWrapper) = missing +_refspace(a::RefSpace) = a struct _LocalBasisTimes{T,U,V} <: _LocalBasisOperations{T} el1::U el2::V @@ -28,26 +43,60 @@ end function (op::U where {U <: _LocalBasisOperations})(p) l = op.el1(p) r = op.el2(p) - return _execute_operation(l,r,op) + #return _execute_operation(l,r,op) + _modbasis(l,r) do a,b + operation(op)(a,b) + end end function (op::NormalVector)(x::CompScienceMeshes.MeshPointNM) return normal(x) end +function _modbasis_genleft(::Type{U}, ::Type{V}) where {U<:SVector{N}, V} where {N} + ex = :(SVector{N}(())) + for n in 1:N + # push!(ex.args[2].args, :(dot(a[$n], b[$m]))) + push!(ex.args[2].args, :(value = f(a[$n].value, b),)) + end + + return ex +end +function _modbasis_genright(::Type{U}, ::Type{V}) where {V<:SVector{N}, U} where {N} + ex = :(SVector{N}(())) + + for n in 1:N + # push!(ex.args[2].args, :(dot(a[$n], b[$m]))) + push!(ex.args[2].args, :(value = f(a, b[$n].value),)) + end + + return ex +end +@generated function _modbasis(f, a::SVector{N,T}, b::U) where {N, T <:NamedTuple,U} + ex = _modbasis_genleft(a,b) + + return ex +end +@generated function _modbasis(f, a::U, b::SVector{N,T}) where {N,T<:NamedTuple,U} + ex = _modbasis_genright(a,b) + + return ex +end + operation(a::_LocalBasisTimes) = * operation(a::_LocalBasisCross) = × operation(a::_LocalBasisDot) = (x,y) --> transpose(x)*y -_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector,op::U) where {N,U} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) -_execute_operation(el1::SVector{N,<:NamedTuple},el2::U,op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) -_execute_operation(el1::SVector,el2::SVector{N,<:NamedTuple},op::U) where {N,U} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) -_execute_operation(el1::U,el2::SVector{N,<:NamedTuple},op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) -_execute_operation(el1::SVector,el2::SVector,op::U) where {U <: Union{_LocalBasisCross,_LocalBasisDot}} = operation(op)(el1,el2) -_execute_operation(el1::U,el2::V,op::_LocalBasisTimes) where {U <: Number,V <: Number}= el1*el2 -_execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not supported" - -_execute_operation_named(a::NTuple{N},b,op) where {N} = ((value=op(a[1].value,b),), _execute_operation_named(Base.tail(a),b,op)...) -_execute_operation_named(a::NTuple{1},b,op) = ((value=op(a[1].value,b),),) -_execute_operation_named(b,a::NTuple{N},op) where {N} = ((value=op(b,a[1].value),), _execute_operation_named(b,Base.tail(a),op)...) -_execute_operation_named(b,a::NTuple{1},op) = ((value=op(b,a[1].value),),) \ No newline at end of file +# _execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector,op::U) where {N,U} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) +# _execute_operation(el1::SVector{N,<:NamedTuple},el2::U,op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) +# _execute_operation(el1::SVector,el2::SVector{N,<:NamedTuple},op::U) where {N,U} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) +# _execute_operation(el1::U,el2::SVector{N,<:NamedTuple},op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) +# _execute_operation(el1::SVector,el2::SVector,op::U) where {U <: Union{_LocalBasisCross,_LocalBasisDot}} = operation(op)(el1,el2) +# _execute_operation(el1::U,el2::V,op::_LocalBasisTimes) where {U <: Number,V <: Number}= el1*el2 +# _execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not supported" + +# _execute_operation_named(a::NTuple{N},b,op) where {N} = ((value=op(a[1].value,b),), _execute_operation_named(Base.tail(a),b,op)...) +# _execute_operation_named(a::NTuple{1},b,op) = ((value=op(a[1].value,b),),) +# _execute_operation_named(b,a::NTuple{N},op) where {N} = ((value=op(b,a[1].value),), _execute_operation_named(b,Base.tail(a),op)...) +# _execute_operation_named(b,a::NTuple{1},op) = ((value=op(b,a[1].value),),) + diff --git a/src/bases/ndspace.jl b/src/bases/ndspace.jl index 0393e5bb..8ff9f9f0 100644 --- a/src/bases/ndspace.jl +++ b/src/bases/ndspace.jl @@ -40,12 +40,20 @@ function nedelec(surface, edges=skeleton(surface,1)) NDBasis(surface, fns, pos) end +# function LinearAlgebra.cross(::NormalVector, s::NDBasis) +# # @assert CompScienceMeshes.isoriented(s.geo) +# RTBasis(s.geo, s.fns, s.pos) +# end + function LinearAlgebra.cross(::NormalVector, s::NDBasis) - # @assert CompScienceMeshes.isoriented(s.geo) - RTBasis(s.geo, s.fns, s.pos) + @assert CompScienceMeshes.isoriented(s.geo) + fns = similar(s.fns) + for (i,fn) in pairs(s.fns) + fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn] + end + RTBasis(s.geo, fns, s.pos) end - function curl(space::NDBasis) - divergence(n × space) + -divergence(n × space) end diff --git a/src/bases/rtspace.jl b/src/bases/rtspace.jl index b6a280a6..6f21271d 100644 --- a/src/bases/rtspace.jl +++ b/src/bases/rtspace.jl @@ -264,7 +264,7 @@ function LinearAlgebra.cross(::NormalVector, s::RTBasis) @assert CompScienceMeshes.isoriented(s.geo) fns = similar(s.fns) for (i,fn) in pairs(s.fns) - fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn] + fns[i] = [Shape(sh.cellid, sh.refid, sh.coeff) for sh in fn] end NDBasis(s.geo, fns, s.pos) end diff --git a/src/composedoperators/analytic_excitation.jl b/src/composedoperators/analytic_excitation.jl index ada5e7b1..31263c3f 100644 --- a/src/composedoperators/analytic_excitation.jl +++ b/src/composedoperators/analytic_excitation.jl @@ -1,15 +1,15 @@ -import Base.* -struct NTimesTrace{T,F} <: Functional - field::F -end +# import Base.* +# struct NTimesTrace{T,F} <: Functional +# field::F +# end -NTimesTrace(f::F) where {F} = NTimesTrace{scalartype(f), F}(f) -NTimesTrace{T}(f::F) where {T,F} = NTimesTrace{T,F}(f) -scalartype(s::NTimesTrace{T}) where {T} = T +# NTimesTrace(f::F) where {F} = NTimesTrace{scalartype(f), F}(f) +# NTimesTrace{T}(f::F) where {T,F} = NTimesTrace{T,F}(f) +# scalartype(s::NTimesTrace{T}) where {T} = T -(ϕ::NTimesTrace)(p) = *(normal(p), ϕ.field(cartesian(p))) -integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) -*(::NormalVector, f) = NTimesTrace(f) +# (ϕ::NTimesTrace)(p) = *(normal(p), ϕ.field(cartesian(p))) +# integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) +# *(::NormalVector, f) = NTimesTrace(f) struct Trace{T,F} <: Functional field::F @@ -20,4 +20,4 @@ Trace{T}(f::F) where {T,F} = Trace{T,F}(f) scalartype(s::Trace{T}) where {T} = T (ϕ::Trace)(p) = ϕ.field(cartesian(p)) -integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) +integrand(::Trace, g, ϕ) = *(g.value, ϕ) diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl index 3c3fe5da..0e0f0255 100644 --- a/src/composedoperators/composedoperator.jl +++ b/src/composedoperators/composedoperator.jl @@ -134,8 +134,8 @@ function (op::Integrand{<:CompDoubleKern})(x,y,f,g) end end -defaultquadstrat(op::CompDoubleInt,tfs::Space,bfs::Space) = DoubleNumSauterQstrat(3,3,3,3,3,3) -defaultquadstrat(op::CompSingleInt,tfs::Space,bfs::Space) = SingleNumQStrat(3) +defaultquadstrat(op::CompDoubleInt,tfs::Space,bfs::Space) = DoubleNumSauterQstrat(5,5,5,5,5,5) +defaultquadstrat(op::CompSingleInt,tfs::Space,bfs::Space) = SingleNumQStrat(6) ##### assembling the single integral function assemble!(biop::CompSingleKern, tfs::Space, bfs::Space, store, threading::Type{Threading{:multi}}; diff --git a/src/composedoperators/displacementmesh.jl b/src/composedoperators/displacementmesh.jl index 07174043..e6079ce4 100644 --- a/src/composedoperators/displacementmesh.jl +++ b/src/composedoperators/displacementmesh.jl @@ -13,10 +13,21 @@ struct DisplacementChart chart epsilon end - +CompScienceMeshes.mesh(m::GlobalDisplacementMesh) = m.mesh displacementchart(m::GlobalDisplacementMesh,i) = DisplacementChart(chart(m.mesh,i),m.epsilon) displacementchart(m::Mesh,i) = DisplacementChart(chart(m,i),-1.0) -CompScienceMeshes.chart(m::DisplacementMesh,i) = chart(m.mesh,i) +CompScienceMeshes.chart(m::DisplacementMesh,i) = chart(mesh(m),i) +CompScienceMeshes.cells(m::DisplacementMesh) = cells(mesh(m)) +CompScienceMeshes.celltype(m::DisplacementMesh) = CompScienceMeshes.celltype(mesh(m)) +CompScienceMeshes.indextype(m::DisplacementMesh,a) = CompScienceMeshes.indextype(mesh(m),a) +CompScienceMeshes.celltype(m::DisplacementMesh,a) = CompScienceMeshes.celltype(mesh(m),a) +CompScienceMeshes.indices(m::DisplacementMesh,n) = CompScienceMeshes.indices(mesh(m),n) +CompScienceMeshes.vertices(m::DisplacementMesh) = CompScienceMeshes.vertices(mesh(m)) +CompScienceMeshes.numvertices(m::DisplacementMesh) = CompScienceMeshes.numvertices(mesh(m)) +CompScienceMeshes.numcells(m::DisplacementMesh) = CompScienceMeshes.numcells(mesh(m)) +CompScienceMeshes.vertextype(m::DisplacementMesh) = CompScienceMeshes.vertextype(mesh(m)) + + function displacement(a::DisplacementChart,b::DisplacementChart) rel_orient = sign(dot(normal(a.chart),normal(b.chart))) if a.epsilon*rel_orient ≈ b.epsilon diff --git a/src/composedoperators/trace.jl b/src/composedoperators/trace.jl index b25fe668..c2afe6ad 100644 --- a/src/composedoperators/trace.jl +++ b/src/composedoperators/trace.jl @@ -17,39 +17,38 @@ end #### trace of a potential: -_trace(::Kernel,n;type=Float64,D) = 0, 0 -function _trace(::HH3DGradGreen,n,D;type=Float64) - if D==2 - return DisplacementVector{type}(n) , 0.5 - end - return 0, 0 +_trace(::Kernel,n,o;type=Float64) = 0, 0 +function _trace(::HH3DGradGreen,n,::Val{2};type=Float64) + return DisplacementVector{type}(n) , 0.5 end """ sign indicates if trace is taken allong the normal (+1) are opposite to the normal (-1) """ function ttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=false) where {D} - t,f = _trace(a.kernel,sign,D) + t,f = _trace(a.kernel,sign,Val{D}()) if testfunction_tangential doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) else doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && (singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) end if t==0 + return doubleint else + #eventualy a scan to be non zero can be added here. return doubleint + singleint end end function ntrace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,D) + t,f = _trace(a.kernel,sign,Val{D}()) doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) if t==0 return doubleint @@ -60,9 +59,9 @@ function ntrace(a::PotentialIntegralOperator{D},sign) where {D} end function trace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,D) + t,f = _trace(a.kernel,sign,Val{D}()) doubleint = CompDoubleInt(B->B,*,a.kernel,a.pairing2,a.bfunc) - t!=0 && singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc)) if t==0 return doubleint @@ -73,9 +72,9 @@ function trace(a::PotentialIntegralOperator{D},sign) where {D} end function strace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,D) + t,f = _trace(a.kernel,sign,Val{D}()) doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) if t==0 return doubleint diff --git a/src/identityop.jl b/src/identityop.jl index c55d1dce..93ced4be 100644 --- a/src/identityop.jl +++ b/src/identityop.jl @@ -86,6 +86,8 @@ function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg}, return qd, A end +support(::LinearRefSpaceTriangle) = TriangleSupport() +support(::LinearRefSpaceTetr) = TetraSupport() defaultquadstrat(::LocalOperator, ::GWPDivRefSpace{<:Real,D1}, ::GWPDivRefSpace{<:Real,D2}) where {D1,D2} = SingleNumQStrat(7) @@ -110,3 +112,27 @@ function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs: end return A end + + +defaultquadstrat(::LocalOperator, ::TriangleSupport, ::TriangleSupport) = SingleNumQStrat(6) +defaultquadstrat(::LocalOperator, ::TetraSupport, ::TetraSupport) = SingleNumQStrat(3) +defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b) = defaultquadstrat(op,support(a),support(b)) +defaultquadstrat(op::LocalOperator, a, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) +defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) + +function quaddata(op::LocalOperator, g, f, tels, bels, + qs::SingleNumQStrat) + gs = support(g) + fs = support(f) + quaddata(op,g,f,gs,fs,tels,bels,qs) +end + +function quaddata(op::LocalOperator, g, f,supg::TriangleSupport,supf::TriangleSupport, tels, bels, + qs::SingleNumQStrat) + + u, w = trgauss(qs.quad_rule) + qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] + A = _alloc_workspace(qd, g, f, tels, bels) + + return qd, A +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index de253fbd..9d544756 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -80,6 +80,9 @@ include("test_handlers.jl") include("test_gridfunction.jl") include("test_composed_basis.jl") +include("test_composed_operator.jl") +include("test_analytic_excitation.jl") + using TestItemRunner @run_package_tests diff --git a/test/test_analytic_excitation.jl b/test/test_analytic_excitation.jl new file mode 100644 index 00000000..08ca71e5 --- /dev/null +++ b/test/test_analytic_excitation.jl @@ -0,0 +1,41 @@ +using BEAST +using CompScienceMeshes +using LinearAlgebra +using Test +using StaticArrays + +# define the plane wave excitation analytically +Γ = meshcuboid(1.0,1.0,1.0,0.2) +Y = unitfunctionc0d1(Γ) +X = raviartthomas(Γ) + +### test crosstracemw +a = (@SVector [1.0,0.0,0.0]) +e(x) = a*exp(-1.0im*dot((@SVector [1.0,0.0,0.0]),x))# not physical but to test the traces. +ex = BEAST.CrossTraceMW(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,ComplexF64}}(e))) +# probeer testen te verzinnen die integraal nul geven. + +rhs = assemble(ex, n× X) + +# loop projector +∂Γ = boundary(Γ) +edges = setminus(skeleton(Γ,1), ∂Γ) +verts = setminus(skeleton(Γ,0), skeleton(∂Γ,0)) +Λ = Matrix(connectivity(verts, edges, sign)) + +PΛ = Λ * pinv(Λ'*Λ) * Λ' + +z = PΛ * rhs +@test isapprox(norm(z), 0.0;atol = sqrt(eps())) + +### test NDotTraceMW + +a = (@SVector [1.0,0.0,0.0]) +e(x) = a +ex = BEAST.NDotTrace(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,Float64}}(e))) +# probeer testen te verzinnen die integraal nul geven. + +rhs = assemble(ex, Y) +@test isapprox(norm(rhs), 0.0;atol = sqrt(eps())) + + diff --git a/test/test_composed_basis.jl b/test/test_composed_basis.jl index e4ae5b24..527e2f3b 100644 --- a/test/test_composed_basis.jl +++ b/test/test_composed_basis.jl @@ -3,32 +3,24 @@ using CompScienceMeshes using LinearAlgebra using BEAST using Test +using StaticArrays Γ = meshcuboid(1.0,1.0,1.0,1.0) X = raviartthomas(Γ) -func(x) = 2.0 +func = BEAST.FunctionWrapper{Float64}(x->2.0) Y = BEAST._BasisTimes(X,func) K = Maxwell3D.doublelayer(wavenumber=1.0) m1 = assemble(K,Y,Y) m2 = assemble(K,X,X) @test m1 ≈ 4*m2 -U = BEAST._BasisTimes(BEAST._BasisDot(x->(@SVector [1.0,1.0,1.0]),n),X) -m = assemble(K,U,U) L = lagrangecxd0(Γ) Z = BEAST._BasisTimes(L,n) m3 = assemble(K,Z,Z) @test norm(m3) ≈ 0.08420116178577139 +m1 = assemble(NCross(),X,X) +m2 = assemble(Identity(),X×n,X) - -##### go in code -# using StaticArrays -# s = simplex((@SVector [1.0,0.0,0.0]),(@SVector [0.0,1.0,0.0]),(@SVector [0.0,0.0,0.0])) -# p = neighborhood(s,[0.5,0.2]) - -# rtref = BEAST.RTRefSpace{Float64}() -# f(x)= 2.0 - -# mref = BEAST._LocalBasisTimes(Float64,BEAST.FunctionWrapper(f),rtref) +@test m1 ≈ m2 diff --git a/test/test_composed_operator.jl b/test/test_composed_operator.jl new file mode 100644 index 00000000..f5a88a1b --- /dev/null +++ b/test/test_composed_operator.jl @@ -0,0 +1,53 @@ +######## test multiplied basis +using CompScienceMeshes +using LinearAlgebra +using BEAST +using Test +using StaticArrays + + +# different options to compute selfpatch trace +Γ = meshcuboid(1.0,1.0,1.0,0.4) +Γd = BEAST.GlobalDisplacementMesh(Γ,1.0) +X = raviartthomas(Γ) +Xd = raviartthomas(Γd) + +∇G = BEAST.HH3DGradGreen(0.0) +Z = BEAST.PotentialIntegralOperator{2}(∇G,×,b->b) +K = BEAST.ttrace(Z,-1.0;testfunction_tangential=false) #-1.0 trace taken allong the test normal (from outside to inside) +M = assemble(K,X,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) + + + +K_compare = Maxwell3D.doublelayer(wavenumber = 0.0) - 0.5*NCross() #trace from outside to inside +M_compare = assemble(K_compare,X,X;quadstrat = [BEAST.SingleNumQStrat(6),BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)]) + +M_displaced = assemble(K,Xd,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) + +@test M_displaced ≈ M +@test M_displaced ≈ M_compare + +#################################################### +# +# Material interactions (global multi-trace context) +# +#################################################### + +Γ_mirror = -mirrormesh(Γ,(@SVector [1.0,0.0,0.0]),(@SVector [0.0,0.0,0.0])) +Γd_mirror = BEAST.GlobalDisplacementMesh(Γ_mirror,-1.0) +Xm = raviartthomas(Γ_mirror) +Xdm = raviartthomas(Γd_mirror) + +K = ttrace(Z,0.0) #number does not matter here, trace part should just be spawned, this number only matters fr selfpatch or if user is not carefull with displacement meshes +M = assemble(K,Xm,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) + +K_compare = Maxwell3D.doublelayer(wavenumber = 0.0) + 0.5*NCross() #trace from outside to inside +M_compare = assemble(K_compare,Xm,X;quadstrat = [BEAST.SingleNumQStrat(6),BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)]) + +M_displaced = assemble(K,Xdm,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) + +@test M_displaced ≈ M +@test M_displaced ≈ M_compare + + + From f0b03ee14687b79a8e48bd17054f763b7596eb41 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 13 Nov 2024 14:38:24 +0100 Subject: [PATCH 06/17] debug tests --- src/bases/ndspace.jl | 2 +- test/test_analytic_excitation.jl | 8 ++++---- test/test_basis.jl | 4 ++-- test/test_dipole.jl | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/bases/ndspace.jl b/src/bases/ndspace.jl index 8ff9f9f0..ad0f0dab 100644 --- a/src/bases/ndspace.jl +++ b/src/bases/ndspace.jl @@ -55,5 +55,5 @@ function LinearAlgebra.cross(::NormalVector, s::NDBasis) end function curl(space::NDBasis) - -divergence(n × space) + divergence(n × (n × (n × space))) end diff --git a/test/test_analytic_excitation.jl b/test/test_analytic_excitation.jl index 08ca71e5..b3d5e3de 100644 --- a/test/test_analytic_excitation.jl +++ b/test/test_analytic_excitation.jl @@ -11,8 +11,8 @@ X = raviartthomas(Γ) ### test crosstracemw a = (@SVector [1.0,0.0,0.0]) -e(x) = a*exp(-1.0im*dot((@SVector [1.0,0.0,0.0]),x))# not physical but to test the traces. -ex = BEAST.CrossTraceMW(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,ComplexF64}}(e))) +eloc(x) = a*exp(-1.0im*dot((@SVector [1.0,0.0,0.0]),x))# not physical but to test the traces. +ex = BEAST.CrossTraceMW(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,ComplexF64}}(eloc))) # probeer testen te verzinnen die integraal nul geven. rhs = assemble(ex, n× X) @@ -31,8 +31,8 @@ z = PΛ * rhs ### test NDotTraceMW a = (@SVector [1.0,0.0,0.0]) -e(x) = a -ex = BEAST.NDotTrace(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,Float64}}(e))) +eloc(x) = a +ex = BEAST.NDotTrace(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,Float64}}(eloc))) # probeer testen te verzinnen die integraal nul geven. rhs = assemble(ex, Y) diff --git a/test/test_basis.jl b/test/test_basis.jl index 64a7d0ef..8e83f4b7 100644 --- a/test/test_basis.jl +++ b/test/test_basis.jl @@ -251,8 +251,8 @@ for i in eachindex(crl.fns) crli = sort(crl.fns[i], by=sh->(sh.cellid, sh.refid)) nxgradi = sort(nxgrad.fns[i], by=sh->(sh.cellid, sh.refid)) for j in eachindex(crl.fns[i]) - # @test crl.fns[i][j] == nxgrad.fns[i][j] - @test crli[j].coeff == -nxgradi[j].coeff + @test crl.fns[i][j] == nxgrad.fns[i][j] + #@test crli[j].coeff == -nxgradi[j].coeff end end diff --git a/test/test_dipole.jl b/test/test_dipole.jl index f11e277f..f96f287b 100644 --- a/test/test_dipole.jl +++ b/test/test_dipole.jl @@ -57,7 +57,7 @@ for U in [Float32,Float64] nf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(𝓚), pts, j_EFIE, X) ff_E_EFIE = potential(MWFarField3D(𝓣), pts, j_EFIE, X) ff_H_EFIE = potential(BEAST.MWDoubleLayerFarField3D(𝓚), pts, j_EFIE, X) - ff_H_EFIE_rotated = potential(n × BEAST.MWDoubleLayerFarField3D(𝓚), pts, -j_EFIE, n × X) + ff_H_EFIE_rotated = -potential(n × BEAST.MWDoubleLayerFarField3D(𝓚), pts, -j_EFIE, n × X) ff_H_EFIE_doublerotated = potential(n × BEAST.MWDoubleLayerRotatedFarField3D(n × 𝓚), pts, -j_EFIE, X) From 55d041cdeeb2f46b5a6154cd5144ccebd1010e52 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 13 Nov 2024 14:42:42 +0100 Subject: [PATCH 07/17] cleanup --- examples/vp_global_multi_trace.jl | 26 ------------- src/bases/composedbasis.jl | 40 ++------------------ src/bases/local/localcomposedbasis.jl | 19 ++-------- src/composedoperators/analytic_excitation.jl | 13 ------- src/composedoperators/composedoperator.jl | 2 +- 5 files changed, 7 insertions(+), 93 deletions(-) diff --git a/examples/vp_global_multi_trace.jl b/examples/vp_global_multi_trace.jl index 6ee8e33b..206e240e 100644 --- a/examples/vp_global_multi_trace.jl +++ b/examples/vp_global_multi_trace.jl @@ -23,12 +23,6 @@ function onlys(a) return b end -#from paul to test something -# Γ12 = -Mesh([point(-x-4,y,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) -# Γ1 = weld(Γ11,Γ12;glueop=onlys) -# Γ2 = Mesh([point(-x-4,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) -# Γ3 = -Mesh([point(x,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) - Γ12 = -Mesh([point(x,y,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) Γ1 = weld(Γ11,Γ12;glueop=onlys) Γ2 = Mesh([point(x,-y-4,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11))) @@ -172,23 +166,3 @@ deq = BEAST.discretise(LHSA==RHSA, (p .∈ Pₕ)..., (q .∈ Qₕ)...) # B = assemble(RHSA,BEAST.DirectProductSpace(Pₕ)) u = solve(deq) -# ### to check aginst previous solution -# using JLD2 -# d = load("examples/0.5Z_#1.jld2")["dict"] -# #comparing against previous solution -# X_server = load("examples/RT.jld2")["dict"] - - -# resc = zeros(size(Z)) + LinearAlgebra.I - -# for i in 1:length(X_server) -# if X_server.fns[i] != RT[1].fns[i] -# resc[i,i] *= -1 -# end -# end - -# for i in 1:length(X_server) -# if X_server.fns[i] != RT[1].fns[i] -# resc[i+length(Xd[1]),i+length(Xd[1])] *= -1 -# end -# end diff --git a/src/bases/composedbasis.jl b/src/bases/composedbasis.jl index 3ba24c2b..de2e135c 100644 --- a/src/bases/composedbasis.jl +++ b/src/bases/composedbasis.jl @@ -1,9 +1,7 @@ struct FunctionWrapper{T} func::Function end -# function FunctionWrapper(f::Function;evalpoint = @SVector [0.0,0.0,0.0]) -# FunctionWrapper{eltype(f(evalpoint))}(f) -# end + function (f::FunctionWrapper{T})(x)::T where {T} return f.func(x) end @@ -12,11 +10,7 @@ function (f::FunctionWrapper{T})(x::CompScienceMeshes.MeshPointNM)::T where {T} return f.func(cartesian(x)) end scalartype(F::FunctionWrapper{T}) where {T} = eltype(T) -#scalartype(::NormalVector,s::Space) = vertextype(geometry(s)) -# function scalartype(::NormalVector) -# @warn "The scallartype of the NormalVector is set at Float32, if used in combination with Float64 basis or operator this is no problem, in the case of Float16 it is." -# return Float32 -# end + abstract type _BasisOperations{T} <: Space{T} end @@ -47,24 +41,10 @@ _BasisDot(el2::Space,el1::NormalVector) = _BasisDot{promote_type(eltype(vertexty _BasisDot(el1::Space,el2::FunctionWrapper) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) _BasisDot(el1::FunctionWrapper,el2::Space) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2) -# #### wrapping of the functions -# _BasisTimes(a::Function,b::Function) = _BasisTimes(FunctionWrapper(a),FunctionWrapper(b)) -# _BasisTimes(a::Function,b) = _BasisTimes(FunctionWrapper(a),b) -# _BasisTimes(a,b::Function) = _BasisTimes(a,FunctionWrapper(b)) - -# _BasisCross(a::Function,b::Function) = _BasisCross(FunctionWrapper(a),FunctionWrapper(b)) -# _BasisCross(a::Function,b) = _BasisCross(FunctionWrapper(a),b) -# _BasisCross(a,b::Function) = _BasisCross(a,FunctionWrapper(b)) - -# _BasisDot(a::Function,b::Function) = _BasisDot(FunctionWrapper(a),FunctionWrapper(b)) -# _BasisDot(a::Function,b) = _BasisDot(FunctionWrapper(a),b) -# _BasisDot(a,b::Function) = _BasisDot(a,FunctionWrapper(b)) - - refspace(a::_BasisTimes{T}) where {T} = _LocalBasisTimes(T,refspace(a.el1),refspace(a.el2)) refspace(a::_BasisCross{T}) where {T} = _LocalBasisCross(T,refspace(a.el1),refspace(a.el2)) refspace(a::_BasisDot{T}) where {T} = _LocalBasisDot(T,refspace(a.el1),refspace(a.el2)) -# refspace(a::Function) = FunctionWrapper(a) + refspace(a::FunctionWrapper) = a refspace(a::NormalVector) = a @@ -94,17 +74,3 @@ subset(a::Union{NormalVector,FunctionWrapper},I) = a *(n::NormalVector,s::Space) = _BasisTimes(n,s) *(s::Space,n::NormalVector) = _BasisTimes(s,n) - - -# function restrict(a::_LocalBasisOperations,cell1,cell2) -# s = sign(dot(normal(cell1),normal(cell2)))^number_of_normals(a) -# return s*restrict(inner_refspace(a),cell1,cell2) -# end - - -# number_of_normals(a::NormalVector) = 1 -# number_of_normals(a) = 0 -# number_of_normals(a::_LocalBasisOperations) = number_of_normals(a.el1) + number_of_normals(a.el2) - -# #TODO save wrapped space in type, same with refspace. -# inner_refspace(a::) \ No newline at end of file diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl index 772cd33f..447539b7 100644 --- a/src/bases/local/localcomposedbasis.jl +++ b/src/bases/local/localcomposedbasis.jl @@ -43,7 +43,7 @@ end function (op::U where {U <: _LocalBasisOperations})(p) l = op.el1(p) r = op.el2(p) - #return _execute_operation(l,r,op) + _modbasis(l,r) do a,b operation(op)(a,b) end @@ -56,7 +56,7 @@ end function _modbasis_genleft(::Type{U}, ::Type{V}) where {U<:SVector{N}, V} where {N} ex = :(SVector{N}(())) for n in 1:N - # push!(ex.args[2].args, :(dot(a[$n], b[$m]))) + push!(ex.args[2].args, :(value = f(a[$n].value, b),)) end @@ -66,7 +66,7 @@ function _modbasis_genright(::Type{U}, ::Type{V}) where {V<:SVector{N}, U} where ex = :(SVector{N}(())) for n in 1:N - # push!(ex.args[2].args, :(dot(a[$n], b[$m]))) + push!(ex.args[2].args, :(value = f(a, b[$n].value),)) end @@ -87,16 +87,3 @@ operation(a::_LocalBasisTimes) = * operation(a::_LocalBasisCross) = × operation(a::_LocalBasisDot) = (x,y) --> transpose(x)*y -# _execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector,op::U) where {N,U} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) -# _execute_operation(el1::SVector{N,<:NamedTuple},el2::U,op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1.data,el2,operation(op))) -# _execute_operation(el1::SVector,el2::SVector{N,<:NamedTuple},op::U) where {N,U} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) -# _execute_operation(el1::U,el2::SVector{N,<:NamedTuple},op::_LocalBasisTimes) where {N,U <: Number} = SVector{N}(_execute_operation_named(el1,el2.data,operation(op))) -# _execute_operation(el1::SVector,el2::SVector,op::U) where {U <: Union{_LocalBasisCross,_LocalBasisDot}} = operation(op)(el1,el2) -# _execute_operation(el1::U,el2::V,op::_LocalBasisTimes) where {U <: Number,V <: Number}= el1*el2 -# _execute_operation(el1::SVector{N,<:NamedTuple},el2::SVector{M,<:NamedTuple},op) where {N,M} = @error "multiplication of basisses not supported" - -# _execute_operation_named(a::NTuple{N},b,op) where {N} = ((value=op(a[1].value,b),), _execute_operation_named(Base.tail(a),b,op)...) -# _execute_operation_named(a::NTuple{1},b,op) = ((value=op(a[1].value,b),),) -# _execute_operation_named(b,a::NTuple{N},op) where {N} = ((value=op(b,a[1].value),), _execute_operation_named(b,Base.tail(a),op)...) -# _execute_operation_named(b,a::NTuple{1},op) = ((value=op(b,a[1].value),),) - diff --git a/src/composedoperators/analytic_excitation.jl b/src/composedoperators/analytic_excitation.jl index 31263c3f..48163f17 100644 --- a/src/composedoperators/analytic_excitation.jl +++ b/src/composedoperators/analytic_excitation.jl @@ -1,16 +1,3 @@ -# import Base.* -# struct NTimesTrace{T,F} <: Functional -# field::F -# end - -# NTimesTrace(f::F) where {F} = NTimesTrace{scalartype(f), F}(f) -# NTimesTrace{T}(f::F) where {T,F} = NTimesTrace{T,F}(f) -# scalartype(s::NTimesTrace{T}) where {T} = T - -# (ϕ::NTimesTrace)(p) = *(normal(p), ϕ.field(cartesian(p))) -# integrand(::NTimesTrace, g, ϕ) = *(g.value, ϕ) -# *(::NormalVector, f) = NTimesTrace(f) - struct Trace{T,F} <: Functional field::F end diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl index 0e0f0255..a0c4453e 100644 --- a/src/composedoperators/composedoperator.jl +++ b/src/composedoperators/composedoperator.jl @@ -1,7 +1,7 @@ abstract type AbstractCompInt <: AbstractOperator end abstract type Kernel{T} end -# abstract type _AbstractCompInt <: AbstractOperator end + #Make these structures, those are not yet pulled trough the smartmath function struct CompDoubleInt{T} <: AbstractCompInt tfunc From a42bac5f3c2951388c716dbf87b15ed5ba6c33e2 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Thu, 14 Nov 2024 14:33:13 +0100 Subject: [PATCH 08/17] make it posible to assemble an operator withous making it a linear combination of operators --- src/solvers/solver.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/solver.jl b/src/solvers/solver.jl index 92542222..6bbfb93a 100644 --- a/src/solvers/solver.jl +++ b/src/solvers/solver.jl @@ -257,8 +257,8 @@ function assemble(bf::BilForm, X::DirectProductSpace, Y::DirectProductSpace; y = op[end](op[1:end-1]..., y) end - a = term.coeff * term.kernel - z = materialize(a, x, y) + a = term.kernel + z = term.coeff * materialize(a, x, y) Smap = lift(z, Block(term.test_id), Block(term.trial_id), U, V) T = promote_type(T, eltype(Smap)) From e6ca0af1efd78adb16118334422e530db960d60f Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Fri, 15 Nov 2024 12:04:56 +0100 Subject: [PATCH 09/17] nx on other side to for RT and ND space --- src/bases/ndspace.jl | 9 ++++++++- src/bases/rtspace.jl | 8 ++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/bases/ndspace.jl b/src/bases/ndspace.jl index ad0f0dab..10042d8d 100644 --- a/src/bases/ndspace.jl +++ b/src/bases/ndspace.jl @@ -53,7 +53,14 @@ function LinearAlgebra.cross(::NormalVector, s::NDBasis) end RTBasis(s.geo, fns, s.pos) end - +function LinearAlgebra.cross(s::NDBasis, ::NormalVector) + @assert CompScienceMeshes.isoriented(s.geo) + fns = similar(s.fns) + for (i,fn) in pairs(s.fns) + fns[i] = [Shape(sh.cellid, sh.refid, sh.coeff) for sh in fn] + end + RTBasis(s.geo, fns, s.pos) +end function curl(space::NDBasis) divergence(n × (n × (n × space))) end diff --git a/src/bases/rtspace.jl b/src/bases/rtspace.jl index 6f21271d..b4e4d42f 100644 --- a/src/bases/rtspace.jl +++ b/src/bases/rtspace.jl @@ -268,3 +268,11 @@ function LinearAlgebra.cross(::NormalVector, s::RTBasis) end NDBasis(s.geo, fns, s.pos) end +function LinearAlgebra.cross(s::RTBasis,::NormalVector) + @assert CompScienceMeshes.isoriented(s.geo) + fns = similar(s.fns) + for (i,fn) in pairs(s.fns) + fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn] + end + NDBasis(s.geo, fns, s.pos) +end \ No newline at end of file From 187f2613345e49d8eb9ccb150434fe4a2fd78379 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Fri, 22 Nov 2024 11:44:42 +0100 Subject: [PATCH 10/17] correcting momintegrals for nonconforming overlap for small triangles --- src/quadrature/nonconformingoverlapqrule.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/quadrature/nonconformingoverlapqrule.jl b/src/quadrature/nonconformingoverlapqrule.jl index 12c563b0..651be0ad 100644 --- a/src/quadrature/nonconformingoverlapqrule.jl +++ b/src/quadrature/nonconformingoverlapqrule.jl @@ -1,7 +1,10 @@ struct NonConformingOverlapQRule{S} conforming_qstrat::S end - +function tangent_rank(p::Simplex{U,D}) where {U,D} + G = [dot(p.tangents[i], p.tangents[j]) for i in 1:D, j in 1:D] + return rank(G) == D +end function momintegrals!(op, test_local_space, basis_local_space, @@ -23,6 +26,9 @@ function momintegrals!(op, test_charts = [ch for ch in test_charts if volume(ch) .> 1e6 * eps(T)] bsis_charts = [ch for ch in bsis_charts if volume(ch) .> 1e6 * eps(T)] + test_charts = [ch for ch in test_charts if tangent_rank(ch)] + bsis_charts = [ch for ch in bsis_charts if tangent_rank(ch)] + isempty(test_charts) && return isempty(bsis_charts) && return From 938b70d1b4b9d2d741199d1e03442ac2dd5c50af Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Fri, 22 Nov 2024 11:47:13 +0100 Subject: [PATCH 11/17] debug tangent rank --- src/quadrature/nonconformingoverlapqrule.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quadrature/nonconformingoverlapqrule.jl b/src/quadrature/nonconformingoverlapqrule.jl index 651be0ad..ec8fc8ff 100644 --- a/src/quadrature/nonconformingoverlapqrule.jl +++ b/src/quadrature/nonconformingoverlapqrule.jl @@ -1,7 +1,7 @@ struct NonConformingOverlapQRule{S} conforming_qstrat::S end -function tangent_rank(p::Simplex{U,D}) where {U,D} +function tangent_rank(p::CompScienceMeshes.Simplex{U,D}) where {U,D} G = [dot(p.tangents[i], p.tangents[j]) for i in 1:D, j in 1:D] return rank(G) == D end From f8ba6d53893829614c050e585d3517e2ceea0c9d Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Fri, 29 Nov 2024 11:39:56 +0100 Subject: [PATCH 12/17] apply comments --- src/composedoperators/composedoperator.jl | 28 ++--- src/composedoperators/potentials.jl | 8 +- src/composedoperators/trace.jl | 135 +++++++++++++++------- src/identityop.jl | 95 ++++++++------- 4 files changed, 167 insertions(+), 99 deletions(-) diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl index a0c4453e..13a5e9f3 100644 --- a/src/composedoperators/composedoperator.jl +++ b/src/composedoperators/composedoperator.jl @@ -5,17 +5,17 @@ abstract type Kernel{T} end #Make these structures, those are not yet pulled trough the smartmath function struct CompDoubleInt{T} <: AbstractCompInt tfunc - pairing1 + op1 kernel::Kernel{T} - pairing2 + op2 bfunc end struct CompSingleInt{T} <: AbstractCompInt tfunc - pairing1 + op1 kernel::Kernel{T} - pairing2 + op2 bfunc end scalartype(op::CompDoubleInt{T}) where {T} = T @@ -24,17 +24,17 @@ scalartype(op::CompSingleInt{T}) where {T} = T #the functionallity is pased to the basis at this point. struct CompDoubleKern{T,U,V,W} <: IntegralOperator - pairing1::U + op1::U kernel::W - pairing2::V + op2::V end function CompDoubleKern(p1,k::Kernel{T},p2) where {T} CompDoubleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2) end struct CompSingleKern{T,U,V,W} <: LocalOperator - pairing1::U + op1::U kernel::W - pairing2::V + op2::V end function CompSingleKern(p1,k::Kernel{T},p2) where {T} CompSingleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2) @@ -42,8 +42,8 @@ end scalartype(op::CompDoubleKern{T}) where {T} = T scalartype(op::CompSingleKern{T}) where {T} = T -integralop(a::CompDoubleInt) = CompDoubleKern(a.pairing1,a.kernel,a.pairing2) -integralop(a::CompSingleInt) = CompSingleKern(a.pairing1,a.kernel,a.pairing2) +integralop(a::CompDoubleInt) = CompDoubleKern(a.op1,a.kernel,a.op2) +integralop(a::CompSingleInt) = CompSingleKern(a.op1,a.kernel,a.op2) function assemble!(op::AbstractCompInt, tfs::AbstractSpace, bfs::AbstractSpace, @@ -118,8 +118,8 @@ end ##### integrand evaluation function integrand(op::CompSingleKern,kerneldata,x,f,g,disp) kernel = op.kernel(x,disp) - op1 = op.pairing1 - op2 = op.pairing2 + op1 = op.op1 + op2 = op.op2 op1(f.value,op2(kernel,g.value)) @@ -127,8 +127,8 @@ end function (op::Integrand{<:CompDoubleKern})(x,y,f,g) kernel = op.operator.kernel(x,y) - op1 = op.operator.pairing1 - op2 = op.operator.pairing2 + op1 = op.operator.op1 + op2 = op.operator.op2 _integrands(f,g) do fi, gi op1(fi.value,op2(kernel,gi.value)) end diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl index 09b0934f..d3764b90 100644 --- a/src/composedoperators/potentials.jl +++ b/src/composedoperators/potentials.jl @@ -1,22 +1,22 @@ #D: dimension of manifold over which integration is performed struct PotentialIntegralOperator{D} kernel - pairing2 + op2 bfunc end struct PotentialIntegralOperatorKern{U,V} kernel::U - pairing2::V + op2::V end function integrand(a::PotentialIntegralOperatorKern,y,krn,f,x) - return a.pairing2(a.kernel(y,x),f) + return a.op2(a.kernel(y,x),f) end function potential(op::PotentialIntegralOperator, points, coeffs, basis; type=SVector{3,ComplexF64}, quadstrat=defaultquadstrat(op, basis)) - return potential(PotentialIntegralOperatorKern(op.kernel,op.pairing2),points,coeffs,op.bfunc(basis);type,quadstrat) + return potential(PotentialIntegralOperatorKern(op.kernel,op.op2),points,coeffs,op.bfunc(basis);type,quadstrat) end diff --git a/src/composedoperators/trace.jl b/src/composedoperators/trace.jl index c2afe6ad..ad5b1736 100644 --- a/src/composedoperators/trace.jl +++ b/src/composedoperators/trace.jl @@ -1,39 +1,45 @@ - -#### trace vector - -# orientation is used when test and trial chart have same displacement, in this case the test-cart is displaced in the orientation*normal(test-chart) direction. struct DisplacementVector{T} <: Kernel{T} - orientation + interior::Bool end function (u::DisplacementVector)(x,disp) if disp == 0 - return sign(u.orientation)*normal(x) + return (u.interior - !u.interior) * normal(x) else return sign(disp)*normal(x) end end -#### trace of a potential: -_trace(::Kernel,n,o;type=Float64) = 0, 0 -function _trace(::HH3DGradGreen,n,::Val{2};type=Float64) - return DisplacementVector{type}(n) , 0.5 +""" + _trace(kernel,interior,Val{N}) + + computes the jump contribution of a kernel, where N is the dimension of the potential integral. +""" +_trace(::Kernel,interior,_) = 0, 0 +function _trace(::HH3DGradGreen{T},interior::Bool,::Val{2}) where {T} + return DisplacementVector{real(T)}(interior) , real(T)(0.5) end """ -sign indicates if trace is taken allong the normal (+1) are opposite to the normal (-1) + ttrace(Potential,interior::Bool;testfunction_tangential=false) + +Compute the tangential trace, -n×n× Potential, of a potential operator mapping it to a boundary operator. + +This function assumes the normalvector on the mesh to point outwards. +Global multi-trace structures are supported when this function is called. +A small acceleration can be gained by setting testfunction_tangential to true. """ -function ttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=false) where {D} - t,f = _trace(a.kernel,sign,Val{D}()) +function ttrace(a::PotentialIntegralOperator{D},interior::Bool;testfunction_tangential=false) where {D} + t,f = _trace(a.kernel,interior,Val{D}()) if testfunction_tangential - doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && (singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) + doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.op2,a.bfunc)) else - doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && (singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) + doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) + t!=0 && (singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.op2,a.bfunc)) end if t==0 @@ -44,11 +50,18 @@ function ttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=fal return doubleint + singleint end end +""" + ntrace(Potential,interior::Bool) + +Compute the normal trace, n⋅ Potential, of a potential operator mapping it to a boundary operator. -function ntrace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,Val{D}()) - doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && (singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) +This function assumes the normalvector on the mesh to point outwards. +Global multi-trace structures are supported when this function is called. +""" +function ntrace(a::PotentialIntegralOperator{D},interior::Bool) where {D} + t,f = _trace(a.kernel,interior,Val{D}()) + doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.op2,a.bfunc)) if t==0 return doubleint @@ -57,11 +70,18 @@ function ntrace(a::PotentialIntegralOperator{D},sign) where {D} return doubleint + singleint end end +""" + trace(Potential,interior::Bool) -function trace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,Val{D}()) - doubleint = CompDoubleInt(B->B,*,a.kernel,a.pairing2,a.bfunc) - t!=0 && (singleint = f*CompSingleInt(B->B,*,t,a.pairing2,a.bfunc)) +Compute the trace of a potential operator mapping it to a boundary operator. + +This function assumes the normalvector on the mesh to point outwards. +Global multi-trace structures are supported when this function is called. +""" +function trace(a::PotentialIntegralOperator{D},interior::Bool) where {D} + t,f = _trace(a.kernel,interior,Val{D}()) + doubleint = CompDoubleInt(B->B,*,a.kernel,a.op2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B->B,*,t,a.op2,a.bfunc)) if t==0 return doubleint @@ -70,11 +90,18 @@ function trace(a::PotentialIntegralOperator{D},sign) where {D} return doubleint + singleint end end +""" + rtrace(Potential,interior::Bool) + +Compute the rotated trace, n × Potential, of a potential operator mapping it to a boundary operator. -function strace(a::PotentialIntegralOperator{D},sign) where {D} - t,f = _trace(a.kernel,sign,Val{D}()) - doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) - t!=0 && (singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.pairing2,a.bfunc)) +This function assumes the normalvector on the mesh to point outwards. +Global multi-trace structures are supported when this function is called. +""" +function rtrace(a::PotentialIntegralOperator{D},interior::Bool) where {D} + t,f = _trace(a.kernel,interior,Val{D}()) + doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) + t!=0 && (singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.op2,a.bfunc)) if t==0 return doubleint @@ -83,37 +110,67 @@ function strace(a::PotentialIntegralOperator{D},sign) where {D} return doubleint + singleint end end +""" + pvttrace(Potential,interior::Bool;testfunction_tangential=false) +Compute the principal value of the tangential trace, -n×n× Potential. + +A small acceleration can be gained by setting testfunction_tangential to true. +""" # the principal value of the trace. -function pvttrace(a::PotentialIntegralOperator{D},sign;testfunction_tangential=false) where {D} +function pvttrace(a::PotentialIntegralOperator{D};testfunction_tangential=false) where {D} if testfunction_tangential - doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) else - doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) end return doubleint end -function pvntrace(a::PotentialIntegralOperator{D},sign) where {D} +""" + pvntrace(Potential,interior::Bool) + +Compute the principal value of the normal trace, n⋅ Potential. - doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) +This function assumes the normalvector on the mesh to point outwards. +""" +function pvntrace(a::PotentialIntegralOperator{D}) where {D} + + doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) return doubleint end -function pvtrace(a::PotentialIntegralOperator{D},sign) where {D} - doubleint = CompDoubleInt(B-> B,*,a.kernel,a.pairing2,a.bfunc) +""" + pvtrace(Potential,interior::Bool) + +Compute the principal value of the trace. +""" +function pvtrace(a::PotentialIntegralOperator{D}) where {D} + doubleint = CompDoubleInt(B-> B,*,a.kernel,a.op2,a.bfunc) return doubleint end -function pvstrace(a::PotentialIntegralOperator{D},sign) where {D} - doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) +""" + pvrtrace(Potential,interior::Bool) + +Compute the principal value of the rotated trace, n× Potential. + +This function assumes the normalvector on the mesh to point outwards. +""" +function pvrtrace(a::PotentialIntegralOperator{D}) where {D} + doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) return doubleint end +""" + volumetrace(Potential) + +Mapping the potential to a volume operator. +""" function volumetrace(a::PotentialIntegralOperator{D}) where {D} - CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.pairing2,a.bfunc) + CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc) end diff --git a/src/identityop.jl b/src/identityop.jl index 93ced4be..c38bde36 100644 --- a/src/identityop.jl +++ b/src/identityop.jl @@ -22,25 +22,25 @@ end const LinearRefSpaceTriangle = Union{RTRefSpace, RT2RefSpace, NDRefSpace, ND2RefSpace, BDMRefSpace, NCrossBDMRefSpace} defaultquadstrat(::LocalOperator, ::LinearRefSpaceTriangle, ::LinearRefSpaceTriangle) = SingleNumQStrat(6) -function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels, - qs::SingleNumQStrat) +# function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels, +# qs::SingleNumQStrat) - u, w = trgauss(qs.quad_rule) - qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] - A = _alloc_workspace(qd, g, f, tels, bels) +# u, w = trgauss(qs.quad_rule) +# qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] +# A = _alloc_workspace(qd, g, f, tels, bels) - return qd, A -end +# return qd, A +# end defaultquadstrat(::LocalOperator, ::subReferenceSpace, ::subReferenceSpace) = SingleNumQStrat(6) -function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels, - qs::SingleNumQStrat) +# function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels, +# qs::SingleNumQStrat) - u, w = trgauss(qs.quad_rule) - qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] - A = _alloc_workspace(qd, g, f, tels, bels) - return qd, A -end +# u, w = trgauss(qs.quad_rule) +# qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] +# A = _alloc_workspace(qd, g, f, tels, bels) +# return qd, A +# end const LinearRefSpaceTetr = Union{NDLCCRefSpace, NDLCDRefSpace, BDM3DRefSpace} defaultquadstrat(::LocalOperator, ::LinearRefSpaceTetr, ::LinearRefSpaceTetr) = SingleNumQStrat(3) @@ -65,14 +65,14 @@ function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg}, end defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,3}, ::LagrangeRefSpace{T,D2,3}) where {T,D1,D2} = SingleNumQStrat(6) -function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg}, - f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat) +# function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg}, +# f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat) - u, w = trgauss(qs.quad_rule) - qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)] - A = _alloc_workspace(qd, g, f, tels, bels) - return qd, A -end +# u, w = trgauss(qs.quad_rule) +# qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)] +# A = _alloc_workspace(qd, g, f, tels, bels) +# return qd, A +# end defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,4}, ::LagrangeRefSpace{T,D2,4}) where {T,D1,D2} = SingleNumQStrat(6) function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg}, @@ -93,15 +93,15 @@ defaultquadstrat(::LocalOperator, ::GWPDivRefSpace{<:Real,D1}, ::GWPDivRefSpace{<:Real,D2}) where {D1,D2} = SingleNumQStrat(7) -function quaddata(op::LocalOperator, g::GWPDivRefSpace, f::GWPDivRefSpace, - tels::Vector, bels::Vector, - qs::SingleNumQStrat) +# function quaddata(op::LocalOperator, g::GWPDivRefSpace, f::GWPDivRefSpace, +# tels::Vector, bels::Vector, +# qs::SingleNumQStrat) - u, w = trgauss(qs.quad_rule) - qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] - A = _alloc_workspace(qd, g, f, tels, bels) - return qd, A -end +# u, w = trgauss(qs.quad_rule) +# qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] +# A = _alloc_workspace(qd, g, f, tels, bels) +# return qd, A +# end function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs::SingleNumQStrat) @@ -114,21 +114,32 @@ function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs: end -defaultquadstrat(::LocalOperator, ::TriangleSupport, ::TriangleSupport) = SingleNumQStrat(6) -defaultquadstrat(::LocalOperator, ::TetraSupport, ::TetraSupport) = SingleNumQStrat(3) -defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b) = defaultquadstrat(op,support(a),support(b)) -defaultquadstrat(op::LocalOperator, a, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) -defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) +defaultquadstrat(::LocalOperator, _, _) = SingleNumQStrat(6) +# defaultquadstrat(::LocalOperator, ::TetraSupport, ::TetraSupport) = SingleNumQStrat(3) +# defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b) = defaultquadstrat(op,support(a),support(b)) +# defaultquadstrat(op::LocalOperator, a, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) +# defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b)) -function quaddata(op::LocalOperator, g, f, tels, bels, - qs::SingleNumQStrat) - gs = support(g) - fs = support(f) - quaddata(op,g,f,gs,fs,tels,bels,qs) -end +# function quaddata(op::LocalOperator, g, f, tels, bels, +# qs::SingleNumQStrat) +# gs = support(g) +# fs = support(f) +# quaddata(op,g,f,gs,fs,tels,bels,qs) +# end + +# function quaddata(op::LocalOperator, g, f,supg::TriangleSupport,supf::TriangleSupport, tels, bels, +# qs::SingleNumQStrat) + +# u, w = trgauss(qs.quad_rule) +# qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] +# A = _alloc_workspace(qd, g, f, tels, bels) + +# return qd, A +# end -function quaddata(op::LocalOperator, g, f,supg::TriangleSupport,supf::TriangleSupport, tels, bels, - qs::SingleNumQStrat) +function quaddata(op::LocalOperator, g, f, + tels::Vector{<:Simplex{U,2}}, bels::Vector{<:Simplex{U,2}}, + qs::SingleNumQStrat) where {U} u, w = trgauss(qs.quad_rule) qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] From f4263cbb98206bbe75eefaf5834ff43958534b7f Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Fri, 29 Nov 2024 14:22:51 +0100 Subject: [PATCH 13/17] update + test pass --- examples/vp_global_multi_trace.jl | 28 ++++++++++++++-------------- src/identityop.jl | 16 ++++++++-------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/examples/vp_global_multi_trace.jl b/examples/vp_global_multi_trace.jl index 206e240e..553f7a2d 100644 --- a/examples/vp_global_multi_trace.jl +++ b/examples/vp_global_multi_trace.jl @@ -54,53 +54,53 @@ G = [BEAST.HH3DGreen(1im*k) for k in κ] dG0 = BEAST.HH3DGradGreen(1im*κ0) dG = [BEAST.HH3DGradGreen(1im*k) for k in κ] - U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) - U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0) + U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) + U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0) U22 = BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->n*b),-1.0) DDL0 = U11[a1,b1] + U12[a1,b2] + U22[a2,b2] - U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] - U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for g in G] + U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] + U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for g in G] U22 = [BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->n*b),1.0) for g in dG] DDL = [U11i[a1,b1] + er*mr*U12i[a1,b2] + U22i[a2,b2] for (U11i,U12i,U22i,er,mr) in zip(U11,U12,U22,ϵr,μr)] - U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) + U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) U22 = BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0) NDL0 = U11[a1,b1] + U21[a2,b1] + U22[a2,b2] - U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] + U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G] U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG] NDL = [U11i[a1,b1] + er*mr*U21i[a2,b1] + U22i[a2,b2] for (U11i,U21i,U22i,er,mr) in zip(U11,U21,U22,ϵr,μr)] - U11 = -BEAST.strace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) - U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0) + U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) + U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0) U21 = -BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->b),-1.0) U22 = -κ0^2* BEAST.trace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0) NDSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2] - U11 = [-BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G] - U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG] + U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G] + U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG] U21 = [-BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->b),1.0) for g in dG] U22 = [-k^2*BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for (g,k) in zip(G,κ)] NDSL = [mr*U11i[a1,b1] + 1/er*U12i[a1,b2] + 1/er*U21i[a2,b1] + 1/(er^2*mr)*U22i[a2,b2] for (U11i,U12i,U21i,U22i,er,mr) in zip(U11,U12,U21,U22,ϵr,μr)] - U11 = -κ0^2 * BEAST.pvstrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0)-BEAST.CompDoubleInt(B->divergence(n×B),*,G0,*,B->divergence(B)) - U12 = BEAST.strace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->n*b),-1.0) + U11 = -κ0^2 * BEAST.pvrtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b))-BEAST.CompDoubleInt(B->divergence(n×B),*,G0,*,B->divergence(B)) + U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->n*b),-1.0) U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0) U22 = BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0) DNSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2] - U11 = [-k^2 * BEAST.pvstrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0)-BEAST.CompDoubleInt(B->divergence(n×B),*,g,*,B->divergence(B)) for (g,k) in zip(G,κ)] - U12 = [BEAST.strace(BEAST.PotentialIntegralOperator{2}(g,×,b->n*b),1.0) for g in dG] + U11 = [-k^2 * BEAST.pvrtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b))-BEAST.CompDoubleInt(B->divergence(n×B),*,g,*,B->divergence(B)) for (g,k) in zip(G,κ)] + U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->n*b),1.0) for g in dG] U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG] U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for (g,k) in zip(G,κ)] diff --git a/src/identityop.jl b/src/identityop.jl index c38bde36..0209bae7 100644 --- a/src/identityop.jl +++ b/src/identityop.jl @@ -33,14 +33,14 @@ defaultquadstrat(::LocalOperator, ::LinearRefSpaceTriangle, ::LinearRefSpaceTria # end defaultquadstrat(::LocalOperator, ::subReferenceSpace, ::subReferenceSpace) = SingleNumQStrat(6) -# function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels, -# qs::SingleNumQStrat) +function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels, + qs::SingleNumQStrat) -# u, w = trgauss(qs.quad_rule) -# qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] -# A = _alloc_workspace(qd, g, f, tels, bels) -# return qd, A -# end + u, w = trgauss(qs.quad_rule) + qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)] + A = _alloc_workspace(qd, g, f, tels, bels) + return qd, A +end const LinearRefSpaceTetr = Union{NDLCCRefSpace, NDLCDRefSpace, BDM3DRefSpace} defaultquadstrat(::LocalOperator, ::LinearRefSpaceTetr, ::LinearRefSpaceTetr) = SingleNumQStrat(3) @@ -138,7 +138,7 @@ defaultquadstrat(::LocalOperator, _, _) = SingleNumQStrat(6) # end function quaddata(op::LocalOperator, g, f, - tels::Vector{<:Simplex{U,2}}, bels::Vector{<:Simplex{U,2}}, + tels::Vector{<:CompScienceMeshes.Simplex{U,2}}, bels::Vector{<:CompScienceMeshes.Simplex{U,2}}, qs::SingleNumQStrat) where {U} u, w = trgauss(qs.quad_rule) From a7103d83c9ef394e9bd7a1589fcd3baa28070032 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Tue, 3 Dec 2024 14:01:50 +0100 Subject: [PATCH 14/17] update potentials + test --- src/composedoperators/potentials.jl | 4 ++-- test/test_composed_operator.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl index d3764b90..026d7c12 100644 --- a/src/composedoperators/potentials.jl +++ b/src/composedoperators/potentials.jl @@ -10,8 +10,8 @@ struct PotentialIntegralOperatorKern{U,V} op2::V end -function integrand(a::PotentialIntegralOperatorKern,y,krn,f,x) - return a.op2(a.kernel(y,x),f) +function integrand(a::PotentialIntegralOperatorKern,krn,y,f,x) + return a.op2(a.kernel(y,x),f.value) end function potential(op::PotentialIntegralOperator, points, coeffs, basis; type=SVector{3,ComplexF64}, diff --git a/test/test_composed_operator.jl b/test/test_composed_operator.jl index f5a88a1b..71f2cec2 100644 --- a/test/test_composed_operator.jl +++ b/test/test_composed_operator.jl @@ -14,7 +14,7 @@ Xd = raviartthomas(Γd) ∇G = BEAST.HH3DGradGreen(0.0) Z = BEAST.PotentialIntegralOperator{2}(∇G,×,b->b) -K = BEAST.ttrace(Z,-1.0;testfunction_tangential=false) #-1.0 trace taken allong the test normal (from outside to inside) +K = BEAST.ttrace(Z,false;testfunction_tangential=false) #-1.0 trace taken allong the test normal (from outside to inside) M = assemble(K,X,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) @@ -38,7 +38,7 @@ M_displaced = assemble(K,Xd,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6 Xm = raviartthomas(Γ_mirror) Xdm = raviartthomas(Γd_mirror) -K = ttrace(Z,0.0) #number does not matter here, trace part should just be spawned, this number only matters fr selfpatch or if user is not carefull with displacement meshes +K = ttrace(Z,true) #number does not matter here, trace part should just be spawned, this number only matters fr selfpatch or if user is not carefull with displacement meshes M = assemble(K,Xm,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)]) K_compare = Maxwell3D.doublelayer(wavenumber = 0.0) + 0.5*NCross() #trace from outside to inside From 2c91be1ea3760d422c041295f13d39f31b0ae1b4 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 4 Dec 2024 14:34:24 +0100 Subject: [PATCH 15/17] correct duallagrangecxd0 and duallagrangec0d1 --- src/bases/lagrange.jl | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index cbc3a617..552f21cb 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -296,10 +296,15 @@ function singleduallagd0(fine, F, v; interpolatory=false) T = coordtype(fine) fn = Shape{T}[] + vol = T(0.0) + for cellid in F + ptch = chart(fine, cellid) + vol += volume(ptch) + end for cellid in F # cell = cells(fine)[cellid] ptch = chart(fine, cellid) - coeff = interpolatory ? T(1.0) : 1 / volume(ptch) / length(F) + coeff = interpolatory ? T(1.0) : 1 / vol refid = 1 push!(fn, Shape(cellid, refid, coeff)) end @@ -498,7 +503,6 @@ end function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) - T = coordtype(mesh) num_faces = dimension(mesh)+1 @@ -559,7 +563,9 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) fine_idcs = cells(refined)[c] local_id = something(findfirst(isequal(v), fine_idcs),0) @assert local_id != 0 - shape = Shape(c, local_id, 1/n/2) + #shape = Shape(c, local_id, 1/n/2) + shape = Shape(c, local_id, 1/2) + push!(fns[i], shape) end end @@ -572,7 +578,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) fine_idcs = cells(refined)[c] local_id = something(findfirst(isequal(v), fine_idcs),0) @assert local_id != 0 - shape = Shape(c, local_id, 1/n/2) + shape = Shape(c, local_id, 1/(n/2)) push!(fns[i], shape) end end @@ -1045,4 +1051,22 @@ function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) for f in mesh pos[nV + nE + (f-1)*nf + 1: nV + nE + f*nf] .= Ref(cartesian(center(chart(mesh, f)))) end return LagrangeBasis{order,0,localdim}(mesh, fns, pos) +end + +function _surface(X) + C = unitfunctioncxd0(X.geo) + G = assemble(Identity(),X,X) + u = Vector(assemble(Identity(),X,C)[1:end,1]) + dot(u,G\u) +end +@testitem "duallagrangec0d1" begin + using CompScienceMeshes + Γ = meshcuboid(1.0,1.0,1.0,0.5) + @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 6.0 +end + +@testitem "duallagrangecxd0" begin + using CompScienceMeshes + Γ = meshcuboid(1.0,1.0,1.0,0.5) + @test BEAST._surface(duallagrangecxd0(Γ)) ≈ 6.0 end \ No newline at end of file From 55bfc664626f8c06758690482a03c9bc05d2c324 Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 11 Dec 2024 15:48:22 +0100 Subject: [PATCH 16/17] update postprocessing potential for composed operators --- src/composedoperators/potentials.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl index 026d7c12..0dbcbcef 100644 --- a/src/composedoperators/potentials.jl +++ b/src/composedoperators/potentials.jl @@ -20,3 +20,6 @@ function potential(op::PotentialIntegralOperator, points, coeffs, basis; return potential(PotentialIntegralOperatorKern(op.kernel,op.op2),points,coeffs,op.bfunc(basis);type,quadstrat) end +defaultquadstrat(op::PotentialIntegralOperator, basis) = SingleNumQStrat(6) +quaddata(op::PotentialIntegralOperatorKern,rs,els,qs::SingleNumQStrat) = quadpoints(rs,els,(qs.quad_rule,)) +quadrule(op::PotentialIntegralOperatorKern,refspace,p,y,q,el,qdata,qs::SingleNumQStrat) = qdata[1,q] \ No newline at end of file From 898495034fa617975e9766655539148bab35b8bb Mon Sep 17 00:00:00 2001 From: PaulOlyslager Date: Wed, 11 Dec 2024 15:59:09 +0100 Subject: [PATCH 17/17] update potentials: add kernelvals --- src/composedoperators/potentials.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl index 0dbcbcef..e7752b7c 100644 --- a/src/composedoperators/potentials.jl +++ b/src/composedoperators/potentials.jl @@ -22,4 +22,5 @@ end defaultquadstrat(op::PotentialIntegralOperator, basis) = SingleNumQStrat(6) quaddata(op::PotentialIntegralOperatorKern,rs,els,qs::SingleNumQStrat) = quadpoints(rs,els,(qs.quad_rule,)) -quadrule(op::PotentialIntegralOperatorKern,refspace,p,y,q,el,qdata,qs::SingleNumQStrat) = qdata[1,q] \ No newline at end of file +quadrule(op::PotentialIntegralOperatorKern,refspace,p,y,q,el,qdata,qs::SingleNumQStrat) = qdata[1,q] +kernelvals(op::PotentialIntegralOperatorKern,y,p) = nothing \ No newline at end of file