diff --git a/CHANGELOG.md b/CHANGELOG.md index 343e833..9cd2288 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ - overhaul of interpolators: each finite element interpolation is composed of the three new structures NodalInterpolator, MomentInterpolator, FunctionalInterpolator + - improved show functions + - coffset function for FESpace + - first, last functions for FEMatrixBlocks and FEVectorBlocks + - improved submatrix function, submatrix function for FEMatrixBlock ## v1.0.0 April 7, 2025 - doc improvements diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index f15e2d4..df3983c 100644 --- a/src/ExtendableFEMBase.jl +++ b/src/ExtendableFEMBase.jl @@ -120,7 +120,7 @@ export interpolator_matrix export FEVectorBlock, FEVector export dot, norm, norms export FEMatrixBlock, FEMatrix, _addnz -export fill!, addblock!, addblock_matmul!, lrmatmul, mul!, add!, apply_penalties! +export fill!, addblock!, addblock_matmul!, lrmatmul, mul!, add!, apply_penalties!, coffset export submatrix export displace_mesh, displace_mesh! diff --git a/src/fematrix.jl b/src/fematrix.jl index 08d6755..6ce4d46 100644 --- a/src/fematrix.jl +++ b/src/fematrix.jl @@ -26,21 +26,14 @@ function Base.copy(FEMB::FEMatrixBlock{TvM, TiM, TvG, TiG, FETypeX, FETypeY, APT return FEMatrixBlock{TvM, TiM, TvG, TiG, FETypeX, FETypeY, APTX, APTY}(deepcopy(FEMB.name), copy(FEMB.FES), copy(FEMB.FESY), FEMB.offset, FEMB.offsetY, FEMB.last_index, FEMB.last_indexY, entries) end -function Base.show(io::IO, FEB::FEMatrixBlock; tol = 1.0e-14) - @printf(io, "\n") - for j in 1:(FEB.offset + 1):FEB.last_index - for k in 1:(FEB.offsetY + 1):FEB.last_indexY - if FEB.entries[j, k] > tol - @printf(io, " +%.1e", FEB.entries[j, k]) - elseif FEB.entries[j, k] < -tol - @printf(io, " %.1e", FEB.entries[j, k]) - else - @printf(io, " ********") - end - end - @printf(io, "\n") - end - return +""" +$(TYPEDSIGNATURES) + +Custom `show` function for `FEMatrixBlock` that prints its coordinates and the name. +""" +function Base.show(io::IO, FEB::FEMatrixBlock) + @printf(io, "[%d:%d,%d:%d]: %s", FEB.offset+1, FEB.last_index, FEB.offsetY+1, FEB.last_indexY, FEB.name) + return nothing end @@ -100,6 +93,9 @@ Base.getindex(FEF::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}, i::Int, Base.getindex(FEB::FEMatrixBlock, i::Int, j::Int) = FEB.entries[FEB.offset + i, FEB.offsetY + j] Base.getindex(FEB::FEMatrixBlock, i::Any, j::Any) = FEB.entries[FEB.offset .+ i, FEB.offsetY .+ j] Base.setindex!(FEB::FEMatrixBlock, v, i::Int, j::Int) = setindex!(FEB.entries, v, FEB.offset + i, FEB.offsetY + j) +Base.first(FEB::FEMatrixBlock) = (FEB.offset+1, FEB.offsetY+1) +Base.last(FEB::FEMatrixBlock) = (FEB.last_index, FEB.last_indexY) + """ @@ -147,14 +143,15 @@ Custom `show` function for `FEMatrix` that prints some information on its blocks function Base.show(io::IO, FEM::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}) where {TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal} println(io, "\nFEMatrix information") println(io, "====================") - println(io, " block | ndofsX | ndofsY | name (FETypeX, FETypeY) ") + println(io, " block | starts | ends | size | name ") for j in 1:length(FEM) n = mod(j - 1, nbrow) + 1 m = Int(ceil(j / nbrow)) - @printf(io, " [%d,%d] |", m, n) - @printf(io, " %6d |", FEM[j].FES.ndofs) - @printf(io, " %6d |", FEM[j].FESY.ndofs) - @printf(io, " %s (%s, %s)\n", FEM[j].name, FEM[j].FES.name, FEM[j].FESY.name) + @printf(io, " [%2d,%2d] |", m, n) + @printf(io, " (%6d,%6d) |", FEM[j].offset+1, FEM[j].offsetY+1) + @printf(io, " (%6d,%6d) |", FEM[j].last_index, FEM[j].last_indexY) + @printf(io, " (%6d,%6d) |", FEM[j].FES.ndofs, FEM[j].FESY.ndofs) + @printf(io, " %s \n", FEM[j].name) end return println(io, "\n nnzvals = $(length(FEM.entries.cscmatrix.nzval))") end @@ -166,19 +163,19 @@ FEMatrix{TvM,TiM}(name::String, FES::FESpace{TvG,TiG,FETypeX,APTX}) where {TvG,T Creates FEMatrix with one square block (FES,FES). """ -function FEMatrix(FES::FESpace; name = "auto") +function FEMatrix(FES::FESpace; name = :automatic) return FEMatrix{Float64, Int64}(FES; name = name) end -function FEMatrix{TvM}(FES::FESpace; name = "auto") where {TvM} +function FEMatrix{TvM}(FES::FESpace; name = :automatic) where {TvM} return FEMatrix{TvM, Int64}(FES; name = name) end -function FEMatrix{TvM, TiM}(FES::FESpace{TvG, TiG, FETypeX, APTX}; name = "auto") where {TvM, TiM, TvG, TiG, FETypeX, APTX} +function FEMatrix{TvM, TiM}(FES::FESpace{TvG, TiG, FETypeX, APTX}; name = :automatic) where {TvM, TiM, TvG, TiG, FETypeX, APTX} return FEMatrix{TvM, TiM}([FES], [FES]; name = name) end """ ```` -FEMatrix{TvM,TiM}(FESX, FESY; name = "auto") +FEMatrix{TvM,TiM}(FESX, FESY; name = :automatic) ```` Creates FEMatrix with one rectangular block (FESX,FESY) if FESX and FESY are single FESpaces, or @@ -209,12 +206,12 @@ end """ ```` -FEMatrix{TvM,TiM}(FESX, FESY; name = "auto") +FEMatrix{TvM,TiM}(FESX, FESY; name = :automatic) ```` Creates an FEMatrix with blocks coressponding to the ndofs of FESX (rows) and FESY (columns). """ -function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:FESpace{TvG, TiG}, 1}; entries = nothing, name = nothing, tags = nothing, tagsX = tags, tagsY = tagsX, npartitions = 1, kwargs...) where {TvM, TiM, TvG, TiG} +function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:FESpace{TvG, TiG}, 1}; entries = nothing, name = :automatic, tags = nothing, tagsX = tags, tagsY = tagsX, npartitions = 1, kwargs...) where {TvM, TiM, TvG, TiG} ndofsX, ndofsY = 0, 0 for j in 1:length(FESX) ndofsX += FESX[j].ndofs @@ -232,7 +229,7 @@ function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:F @assert size(entries) == (ndofsX, ndofsY) "size of given entries not matching number of dofs in given FE space(s)" end - if name === nothing + if name === :automatic || name === nothing name = "" end @@ -250,9 +247,9 @@ function FEMatrix{TvM, TiM}(FESX::Array{<:FESpace{TvG, TiG}, 1}, FESY::Array{<:F offsetY = 0 for k in 1:length(FESY) if (tagsX !== nothing) && (tagsY !== nothing) - blockname = name * " [$(tagsX[j]),$(tagsY[k])" + blockname = name * " $(tagsX[j])[$(FESX[j].name)] x $(tagsY[k])[$(FESY[k].name)]" else - blockname = name * " [$j,$k]" + blockname = name * " $(FESX[j].name) x $(FESY[k].name)" end Blocks[(j - 1) * length(FESY) + k] = FEMatrixBlock{TvM, TiM, TvG, TiG, eltype(FESX[j]), eltype(FESY[k]), assemblytype(FESX[j]), assemblytype(FESY[k])}(blockname, FESX[j], FESY[k], offset, offsetY, offset + FESX[j].ndofs, offsetY + FESY[k].ndofs, entries) @@ -623,21 +620,41 @@ $(TYPEDSIGNATURES) Generates an ExtendableSparseMatrix from the submatrix for the given row and col numbers """ -function submatrix(A::AbstractExtendableSparseMatrixCSC{Tv, Ti}, srows, scols) where {Tv, Ti} +function submatrix(A::AbstractExtendableSparseMatrixCSC{Tv, Ti}, srows, scols; factor = 1) where {Tv, Ti} cscmat::SparseMatrixCSC{Tv, Ti} = A.cscmatrix rows::Array{Ti, 1} = rowvals(cscmat) + valsA = cscmat.nzval + nrowsA = size(A,1) + ncolsA = size(A,2) S = ExtendableSparseMatrix{Tv, Ti}(length(srows), length(scols)) @assert maximum(srows) <= size(A, 1) "rows exceeds rowcount of A" @assert maximum(scols) <= size(A, 2) "cols exceeds colcount of A" - for col in 1:length(scols) - scol = scols[col] + ncols = length(scols) + nrows = length(srows) + newrows = zeros(Int, nrowsA) + newcols = zeros(Int, ncolsA) + newrows[srows] = 1:nrows + newcols[scols] = 1:ncols + minrow, maxrow = minimum(srows), maximum(srows) + + for scol in scols for r in nzrange(cscmat, scol) - j = findfirst(==(rows[r]), srows) - if j !== nothing - S[j, col] = A[rows[r], scol] + if newrows[rows[r]] > 0 + _addnz(S, newrows[rows[r]], newcols[scol], valsA[r], factor) + else + break end end end flush!(S) return S end + +""" +$(TYPEDSIGNATURES) + +Returns the FEMatrixBlock as an ExtendableSparseMatrix +""" +function submatrix(A::FEMatrixBlock{Tv, Ti}) where {Tv, Ti} + return submatrix(A.entries, A.offset+1:A.last_index, A.offsetY+1:A.last_indexY) +end diff --git a/src/fevector.jl b/src/fevector.jl index daefae0..cfcd314 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -22,20 +22,21 @@ function Base.copy(FEB::FEVectorBlock{T, Tv, Ti, FEType, APT}, entries) where {T return FEVectorBlock{T, Tv, Ti, FEType, APT}(deepcopy(FEB.name), copy(FEB.FES), FEB.offset, FEB.last_index, entries) end +""" +$(TYPEDSIGNATURES) + +Returns the number of components for the finite element in that block. +""" get_ncomponents(FB::FEVectorBlock) = get_ncomponents(get_FEType(FB.FES)) -function Base.show(io::IO, FEB::FEVectorBlock; tol = 1.0e-14) - return @printf(io, "\n") - #for j=1:FEB.offset+1:FEB.last_index - # if FEB.entries[j] > tol - # @printf(io, "[%d] = +%.3e", j, FEB.entries[j]); - # elseif FEB.entries[j] < -tol - # @printf(io, "[%d] %.3e", j, FEB.entries[j]); - # else - # @printf(io, "[%d] ********", j ); - # end - # @printf(io, "\n"); - #end +""" +$(TYPEDSIGNATURES) + +Custom `show` function for `FEVectorBlock` that prints some information and the view of that block. +""" +function Base.show(io::IO, ::MIME"text/plain", FEB::FEVectorBlock) + @printf(io, "block %s [%d:%d] = ", FEB.name, FEB.offset+1, FEB.last_index) + show(io, view(FEB)) end """ @@ -68,12 +69,14 @@ Base.setindex!(FEB::FEVectorBlock, v, i::AbstractArray) = (FEB.entries[FEB.offse Base.eltype(::FEVector{T}) where {T} = T Base.size(FEF::FEVector) = size(FEF.FEVectorBlocks) Base.size(FEB::FEVectorBlock) = FEB.last_index - FEB.offset +Base.first(FEB::FEVectorBlock) = FEB.offset+1 +Base.last(FEB::FEVectorBlock) = FEB.last_index """ $(TYPEDEF) -returns a view of the part of the full FEVector that coressponds to the block. +Returns a view of the part of the full FEVector that coressponds to the block. """ Base.view(FEB::FEVectorBlock) = view(FEB.entries, (FEB.offset + 1):FEB.last_index) @@ -88,7 +91,7 @@ end """ $(TYPEDEF) -returns a vector with the individual norms of all blocks +Returns a vector with the individual norms of all blocks. """ function norms(FEV::FEVector{T}, p::Real = 2) where {T} norms = zeros(T, length(FEV)) @@ -191,10 +194,12 @@ Custom `show` function for `FEVector` that prints some information on its blocks function Base.show(io::IO, FEF::FEVector) println(io, "\nFEVector information") println(io, "====================") - print(io, " block | ndofs \t| min / max \t| FEType \t\t (name/tag)") + print(io, " block | starts | ends | length | min / max \t| FEType \t\t (name/tag)") for j in 1:length(FEF) - @printf(io, "\n [%5d] | ", j) - @printf(io, " %6d\t|", FEF[j].FES.ndofs) + @printf(io, "\n [%5d] |", j) + @printf(io, " %6d |", FEF[j].offset+1) + @printf(io, " %6d |", FEF[j].last_index) + @printf(io, " %6d |", FEF[j].FES.ndofs) ext = extrema(view(FEF[j])) @printf(io, " %.2e/%.2e \t|", ext[1], ext[2]) if length(FEF.tags) >= j @@ -260,7 +265,7 @@ end """ $(TYPEDSIGNATURES) -Scalar product between two FEVEctorBlocks +Scalar product between two FEVEctorBlocks. """ function LinearAlgebra.dot(a::FEVectorBlock{T}, b::FEVectorBlock{T}) where {T} return dot(view(a), view(b)) diff --git a/src/finiteelements.jl b/src/finiteelements.jl index 167d21b..e51469b 100644 --- a/src/finiteelements.jl +++ b/src/finiteelements.jl @@ -95,6 +95,15 @@ ndofs(FES::FESpace) = FES.ndofs """ $(TYPEDSIGNATURES) +returns the offset between the degrees of freedom of each component +(i.e. the number of scalar degrees of freedom that influence a component, +vector-valued degrees of freedom are stored at the end). +""" +coffset(FES::FESpace) = FES.coffset + +""" +$(TYPEDSIGNATURES) + returns true if the finite element space is broken, false if not """ broken(FES::FESpace) = FES.broken @@ -215,12 +224,14 @@ function Base.show(io::IO, FES::FESpace{Tv, Ti, FEType, APT}) where {Tv, Ti, FET println(io, "\nFESpace information") println(io, "===================") println(io, " name = $(FES.name)") - println(io, " FEType = $FEType") + println(io, " FEType = $FEType ($(FES.broken ? "$APT, broken" : "$APT"))") println(io, " FEClass = $(supertype(FEType))") - println(io, " ndofs = $(FES.ndofs)\n") + println(io, " ndofs = $(FES.ndofs) $(FES.coffset !== FES.ndofs ? "(coffset = $(FES.coffset))" : "")\n") + println(io, " xgrid = $(FES.xgrid)") + println(io, " dofgrid = $(FES.dofgrid !== FES.xgrid ? FES.dofgrid : "xgrid")") println(io, "") println(io, "DofMaps") - println(io, "==========") + println(io, "=======") for tuple in FES.dofmaps println(io, "> $(tuple[1])") end