Add test coverage for linear operators (#130)

* Add test coverage for linear operators * Fix ambiguity * More test cases
termi-official · Aug 8, 2024 · b087c9c · b087c9c
1 parent a4a0a8d
commit b087c9c
Show file tree

Hide file tree

Showing 5 changed files with 78 additions and 7 deletions.
diff --git a/src/discretization/operator.jl b/src/discretization/operator.jl
@@ -461,7 +461,7 @@ function update_operator!(op::LinearOperator, time)
 end
 
 function _update_linear_operator_on_subdomain!(b, sdh, element_cache, time)
-    ndofs = ndofs_per_cell(dh)
+    ndofs = ndofs_per_cell(sdh)
     bₑ = zeros(ndofs)
     @inbounds for cell in CellIterator(sdh)
         fill!(bₑ, 0)

diff --git a/src/discretization/rsafdq-operator.jl b/src/discretization/rsafdq-operator.jl
@@ -177,7 +177,7 @@ function setup_operator(f::RSAFDQ20223DFunction, solver::AbstractNonlinearSolver
 
     displacement_symbol = first(dh.field_names)
 
-    intorder = quadrature_order(structural_function, displacement_symbol)
+    intorder = default_quadrature_order(structural_function, displacement_symbol)
     qr = QuadratureRuleCollection(intorder)
     qr_face = FacetQuadratureRuleCollection(intorder)
 

diff --git a/src/solver/interface.jl b/src/solver/interface.jl
@@ -15,7 +15,7 @@ function setup_operator(f::AbstractQuasiStaticFunction, solver::AbstractNonlinea
 
     displacement_symbol = first(dh.field_names)
 
-    intorder = quadrature_order(f, displacement_symbol)::Int
+    intorder = default_quadrature_order(f, displacement_symbol)::Int
     qr = QuadratureRuleCollection(intorder)
     qr_face = FacetQuadratureRuleCollection(intorder)
 
@@ -57,7 +57,7 @@ end
 #     @assert length(dh.subdofhandlers) == 1 "Multiple subdomains not yet supported in the Newton solver."
 #     @assert length(dh.field_names) == 1 "Multiple fields not yet supported in the nonlinear solver."
 
-#     intorder = quadrature_order(problem, displacement_symbol)
+#     intorder = default_quadrature_order(problem, displacement_symbol)
 #     qr = QuadratureRuleCollection(intorder)
 #     qr_face = FacetQuadratureRuleCollection(intorder)
 
@@ -121,6 +121,6 @@ function create_system_vector(::Type{<:Vector{T}}, dh::DofHandler) where T
 end
 
 function create_quadrature_rule(f::AbstractSemidiscreteFunction, solver::AbstractSolver, field_name::Symbol)
-    intorder = quadrature_order(f, field_name)
+    intorder = default_quadrature_order(f, field_name)
     return QuadratureRuleCollection(intorder)
 end
diff --git a/src/utils.jl b/src/utils.jl
@@ -285,15 +285,15 @@ end
 @noinline check_subdomains(dh::Ferrite.AbstractDofHandler) = length(dh.subdofhandlers) == 1 || throw(ArgumentError("Using DofHandler with multiple subdomains is not currently supported"))
 @noinline check_subdomains(grid::Ferrite.AbstractGrid) = length(elementtypes(grid)) == 1 || throw(ArgumentError("Using mixed grid is not currently supported"))
 
-@inline function quadrature_order(f, fieldname)
+@inline function default_quadrature_order(f, fieldname)
     @unpack dh = f
     @assert fieldname ∈ dh.field_names "Field $fieldname not found in dof handler. Available fields are: $(dh.field_names)."
 
     for sdh in dh.subdofhandlers
         idx = findfirst(s->s==fieldname, sdh.field_names)
         idx === nothing && continue
         ip = sdh.field_interpolations[idx]
-        return 2*Ferrite.getorder(ip)
+        return max(2*Ferrite.getorder(ip)-1, 2)
     end
 end
 
@@ -319,6 +319,9 @@ struct DenseDataRange{DataVectorType, IndexVectorType}
     offsets::IndexVectorType
 end
 
+Base.size(v::DenseDataRange) = size(v.data)
+Base.getindex(v::DenseDataRange, i::Int) = getindex(v.data, i)
+
 Base.eltype(data::DenseDataRange) = eltype(data.data)
 
 @inline function get_data_for_index(r::DenseDataRange, i::Integer)
@@ -341,6 +344,13 @@ struct EAVector{T, EADataType <: AbstractVector{T}, IndexType <: AbstractVector{
     dof_to_element_map::DenseDataRange{DofMapType, IndexType}
 end
 
+Base.size(v::EAVector) = size(v.eadata)
+Base.getindex(v::EAVector, i::Int) = getindex(v.eadata, i)
+
+function Base.show(io::IO, mime::MIME"text/plain", data::EAVector{T, EADataType, IndexType}) where {T, EADataType, IndexType}
+    println(io, "EAVector{T=", T, ", EADataType=", EADataType, ", IndexType=", IndexType, "} with storate for ", size(data.eadata), " entries." )
+end
+
 @inline get_data_for_index(r::EAVector, i::Integer) = get_data_for_index(r.eadata, i)
 
 function EAVector(dh::DofHandler)

diff --git a/test/test_operators.jl b/test/test_operators.jl
@@ -71,6 +71,67 @@ using BlockArrays, SparseArrays
         @test Thunderbolt.getJ(bop_id) ≈ spdiagm(ones(6))
     end
 
+    @testset "Linear" begin
+        # Setup
+        grid = generate_grid(Quadrilateral, (2,2))
+        dh = DofHandler(grid)
+        add!(dh, :u, Lagrange{RefQuadrilateral,1}())
+        close!(dh)
+        qrc = QuadratureRuleCollection{2}()
+
+        @testset "Constant Cartesian" begin
+            cs = CartesianCoordinateSystem(grid)
+            protocol = AnalyticalTransmembraneStimulationProtocol(
+                AnalyticalCoefficient((x,t) -> 1.0, CoordinateSystemCoefficient(cs)),
+                [SVector((0.0, 1.0))]
+            )
+
+            linop = Thunderbolt.LinearOperator(
+                zeros(ndofs(dh)),
+                protocol,
+                qrc,
+                dh,
+            )
+            Thunderbolt.update_operator!(linop,0.0)
+            @test linop.b ≈ [0.25, 0.5, 1.0, 0.5, 0.25, 0.5, 0.5, 0.25, 0.25]
+
+            plinop = Thunderbolt.PEALinearOperator(
+                zeros(ndofs(dh)),
+                qrc,
+                protocol,
+                dh,
+            )
+            Thunderbolt.update_operator!(plinop,0.0)
+            @test linop.b ≈ plinop.b
+        end
+
+        @testset "Quadratic Cartesian" begin
+            cs = CartesianCoordinateSystem(grid)
+            protocol = AnalyticalTransmembraneStimulationProtocol(
+                AnalyticalCoefficient((x,t) -> norm(x)^2+1.0, CoordinateSystemCoefficient(cs)),
+                [SVector((0.0, 1.0))]
+            )
+
+            linop = Thunderbolt.LinearOperator(
+                zeros(ndofs(dh)),
+                protocol,
+                qrc,
+                dh,
+            )
+            Thunderbolt.update_operator!(linop,0.0)
+            @test linop.b ≈ [1.0/2, 5.0/6, 4.0/3, 5.0/6, 1.0/2, 5.0/6, 5.0/6, 1.0/2, 1.0/2]
+
+            plinop = Thunderbolt.PEALinearOperator(
+                zeros(ndofs(dh)),
+                qrc,
+                protocol,
+                dh,
+            )
+            Thunderbolt.update_operator!(plinop,0.0)
+            @test linop.b ≈ plinop.b
+        end
+    end
+
     @testset "Coupled" begin
         # TODO test with faulty blocks
         @testset "Block action" begin