FESpace, FEVector, FEMatrix improvements

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

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!

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

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))

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