version bump for ExtendableGrids, FESpace now has a dof grid that is …

…used for dof numbering which allows

to handle elements that live only on subgrids; copy functions for FESpace, FEVector, FEMatrix that do not copy the grid

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 authors = ["Christian Merdon <[email protected]>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
@@ -20,7 +20,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 [compat]
 DiffResults = "1"
 DocStringExtensions = "0.8,0.9"
-ExtendableGrids = "1"
+ExtendableGrids = "1.3"
 ExtendableSparse = "^1.2"
 ForwardDiff = "^0.10.35"
 Term = "2"

src/ExtendableFEMBase.jl

@@ -42,8 +42,9 @@ export OperatorPair, OperatorTriple
 include("finiteelements.jl")
 export DofMap
 export CellDofs, FaceDofs, EdgeDofs, BFaceDofs, BEdgeDofs
+export CellDofsParent, FaceDofsParent, EdgeDofsParent, BFaceDofsParent, BEdgeDofsParent
 export DofMapTypes
-export Dofmap4AssemblyType, ItemGeometries4DofMap, EffAT4AssemblyType
+export Dofmap4AssemblyType, ItemGeometries4DofMap, EffAT4AssemblyType, ParentDofmap4Dofmap
 export AbstractFiniteElement
 export FESpace, FESpaces, get_AT, get_FEType
 

src/dofmaps.jl

@@ -11,6 +11,11 @@ abstract type FaceDofs <: DofMap end
 abstract type EdgeDofs <: DofMap end
 abstract type BFaceDofs <: DofMap end
 abstract type BEdgeDofs <: DofMap end
+abstract type CellDofsParent <: DofMap end
+abstract type FaceDofsParent <: DofMap end
+abstract type EdgeDofsParent <: DofMap end
+abstract type BFaceDofsParent <: DofMap end
+abstract type BEdgeDofsParent <: DofMap end
 
 const DofMapTypes{Ti} = Union{VariableTargetAdjacency{Ti}, SerialVariableTargetAdjacency{Ti}, Array{Ti, 2}}
 
@@ -46,6 +51,12 @@ Dofmap4AssemblyType(::Type{ON_BFACES}) = BFaceDofs
 Dofmap4AssemblyType(::Type{<:ON_EDGES}) = EdgeDofs
 Dofmap4AssemblyType(::Type{ON_BEDGES}) = BEdgeDofs
 
+ParentDofmap4Dofmap(::Type{CellDofs}) = CellDofsParent
+ParentDofmap4Dofmap(::Type{FaceDofs}) = FaceDofsParent
+ParentDofmap4Dofmap(::Type{EdgeDofs}) = EdgeDofsParent
+ParentDofmap4Dofmap(::Type{BFaceDofs}) = BFaceDofsParent
+ParentDofmap4Dofmap(::Type{BEdgeDofs}) = BEdgeDofsParent
+
 EffAT4AssemblyType(::Type{ON_CELLS}, ::Type{ON_CELLS}) = ON_CELLS
 EffAT4AssemblyType(::Type{ON_CELLS}, ::Type{<:ON_FACES}) = ON_FACES
 EffAT4AssemblyType(::Type{ON_CELLS}, ::Type{ON_BFACES}) = ON_BFACES
@@ -197,13 +208,12 @@ function Dofmap4AssemblyType(FES::FESpace, AT::Type{<:AssemblyType})
 	return FES[Dofmap4AssemblyType(EffAT4AssemblyType(assemblytype(FES), AT))]
 end
 
-
 function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<:DofMap}) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
 	## Beware: Automatic broken DofMap generation currently only reliable for CellDofs
 
 	@debug "Generating $DM for $(FES.name)"
 	## prepare dofmap patterns
-	xgrid = FES.xgrid
+	xgrid = FES.dofgrid
 	EG = xgrid[UCG4DofMap(DM)]
 	ncomponents::Int = get_ncomponents(FEType)
 	need_faces::Bool = false
@@ -368,15 +378,50 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{<
 			end
 		end
 	end
-	# save dofmap
+
 	FES[DM] = xItemDofs
+
+	if FES.dofgrid !== FES.xgrid
+		## assume parent relation between xgrid and dofgrid
+		@assert FES.dofgrid[ParentGrid] == FES.xgrid "xgrid is not the parent grid of dofgrid !!!"
+		@assert FES.dofgrid[ParentGridRelation] == SubGrid "dofgrid is not a subgrid of xgrid !!!"
+		## lift subgrid dofmap to parent xgrid to allow assembly on common parent grid
+		## by constructing a variable target adjacency with empty dofs for parent items in xgrid
+		## that are not in the dofgrid --> this allows to use the dofmap in assembly over full xgrid
+		xItemDofsLifted = VariableTargetAdjacency(Ti)
+		if DM == CellDofs
+			parentitems = FES.dofgrid[CellParents]
+		elseif DM == FaceDofs
+			parentitems = FES.dofgrid[FaceParents]
+		elseif DM == EdgeDofs
+			parentitems = FES.dofgrid[EdgeParents]
+		elseif DM == BFaceDofs
+			parentitems = FES.dofgrid[BFaceParents]
+		elseif DM == BEdgeDofs
+			parentitems = FES.dofgrid[BEdgeParents]
+		end
+		DMP = ParentDofmap4Dofmap(DM)
+		subitem = 1
+		for i = 1 : num_sources(xItemDofs)
+			t = parentitems[i]
+			for j = subitem : t-1
+				append!(xItemDofsLifted, [])
+				subitem += 1
+			end
+			append!(xItemDofsLifted, view(xItemDofs, :, i))
+			subitem += 1
+		end
+
+		FES[DMP] = xItemDofsLifted
+	end
+	return FES[DM]
 end
 
 
 function init_broken_dofmap!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Union{Type{BFaceDofs}, Type{FaceDofs}}) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
 
 	## prepare dofmap patterns
-	xgrid = FES.xgrid
+	xgrid = FES.dofgrid
 	cellEG = xgrid[UniqueCellGeometries]
 	EG = xgrid[UCG4DofMap(DM)]
 	ncomponents::Int = get_ncomponents(FEType)
@@ -405,15 +450,15 @@ function init_broken_dofmap!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Union{Type{B
 		end
 	end
 
-	xFaceCells = FES.xgrid[FaceCells]
-	xCellNodes = FES.xgrid[CellNodes]
+	xFaceCells = xgrid[FaceCells]
+	xCellNodes = xgrid[CellNodes]
 	xCellDofs = FES[CellDofs]
-	xFaceNodes = FES.xgrid[SuperItemNodes4DofMap(DM)]
+	xFaceNodes = xgrid[SuperItemNodes4DofMap(DM)]
 	xCellGeometries = xgrid[CellGeometries]
 	xFaceGeometries = xgrid[ItemGeometries4DofMap(DM)]
 	xFaceDofs = VariableTargetAdjacency(Ti)
 	if DM == BFaceDofs
-		xRealFace = FES.xgrid[BFaceFaces]
+		xRealFace = xgrid[BFaceFaces]
 		nfaces = length(xRealFace)
 	elseif DM == FaceDofs
 		nfaces = num_sources(xFaceNodes)

src/feevaluator_h1.jl

@@ -265,7 +265,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Lapl
 end
 
 
-# CURLSCALAR H1
+# CURLSCALAR H1 : R1 -> R2
 function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:CurlScalar, <:AbstractH1FiniteElement})
 	L2GAinv = _update_trafo!(FEBE)
 	subset = _update_subset!(FEBE)
@@ -285,7 +285,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl
 	return nothing
 end
 
-# CURL2D H1
+# CURL2D H1 : R2 -> R1
 function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl2D, <:AbstractH1FiniteElement})
 	L2GAinv = _update_trafo!(FEBE)
 	subset = _update_subset!(FEBE)

src/fematrix.jl

@@ -22,6 +22,9 @@ struct FEMatrixBlock{TvM, TiM, TvG, TiG, FETypeX, FETypeY, APTX, APTY} <: Abstra
 	entries::AbstractSparseMatrix{TvM, TiM} # shares with parent object
 end
 
+function Base.copy(FEMB::FEMatrixBlock{TvM, TiM, TvG, TiG, FETypeX, FETypeY, APTX, APTY}, entries) where {TvM, TiM, TvG, TiG, FETypeX, FETypeY, APTX, APTY}
+    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 = 1e-14)
 	@printf(io, "\n")
@@ -50,6 +53,11 @@ struct FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal} <: AbstractSparseMatr
 	entries::AbstractSparseMatrix{TvM, TiM}
 end
 
+function Base.copy(FEV::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}) where {TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}
+    entries = deepcopy(FEV.entries)
+    return FEVector{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}([copy(B, entries) for B in FEV.FEMatrixBlocks], entries)
+end
+
 #Add value to matrix if it is nonzero
 @inline function _addnz(ESM::ExtendableSparseMatrix, i, j, v::Tv, fac) where Tv
 	if v != zero(Tv)

src/fevector.jl

@@ -18,6 +18,10 @@ struct FEVectorBlock{T, Tv, Ti, FEType, APT} <: AbstractArray{T, 1}
 	entries::Array{T, 1} # shares with parent object
 end
 
+function Base.copy(FEB::FEVectorBlock{T,Tv,Ti,FEType,APT}, entries) where {T,Tv,Ti,FEType,APT}
+    return FEVectorBlock{T,Tv,Ti,FEType,APT}(deepcopy(FEB.name), copy(FEB.FES), FEB.offset, FEB.last_index, entries)
+end
+
 get_ncomponents(FB::FEVectorBlock) = get_ncomponents(get_FEType(FB.FES))
 
 function Base.show(io::IO, FEB::FEVectorBlock; tol = 1e-14)
@@ -47,6 +51,11 @@ struct FEVector{T, Tv, Ti} #<: AbstractVector{T}
 	tags::Vector{Any}
 end
 
+function Base.copy(FEV::FEVector{T,Tv,Ti}) where {T,Tv,Ti}
+    entries = deepcopy(FEV.entries)
+    return FEVector{T,Tv,Ti}([copy(B, entries) for B in FEV.FEVectorBlocks], entries, [t for t in FEV.tags])
+end
+
 # overload stuff for AbstractArray{T,1} behaviour
 Base.getindex(FEF::FEVector{T, Tv, Ti}, tag) where {T, Tv, Ti} = FEF.FEVectorBlocks[findfirst(==(tag), FEF.tags)]
 Base.getindex(FEF::FEVector, i::Int) = FEF.FEVectorBlocks[i]

src/finiteelements.jl

@@ -44,10 +44,15 @@ struct FESpace{Tv, Ti, FEType <: AbstractFiniteElement, AT <: AssemblyType}
 	broken::Bool                          # if true, broken dofmaps are generated
 	ndofs::Int64                          # total number of dofs
 	coffset::Int                          # offset for component dofs
-	xgrid::ExtendableGrid{Tv, Ti}          # link to xgrid 
+	xgrid::ExtendableGrid{Tv, Ti}         # link to (master/parent) grid 
+	dofgrid::ExtendableGrid{Tv,Ti}	      # link to (sub) grid used for dof numbering (expected to be equal to or child grid of xgrid)
 	dofmaps::Dict{Type{<:AbstractGridComponent}, Any} # backpack with dofmaps
 end
 
+function Base.copy(FES::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT}
+	return FESpace{Tv, Ti, FEType, AT}(deepcopy(FES.name), FES.broken, FES.ndofs, FES.coffset, FES.xgrid, FES.dofgrid, FES.dofmaps)
+end
+
 get_AT(::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT} = AT
 get_FEType(::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT} = FEType
 
@@ -75,6 +80,7 @@ If no AT is provided, the space is generated ON_CELLS.
 function FESpace{FEType, AT}(
 	xgrid::ExtendableGrid{Tv, Ti};
 	name = "",
+	regions = nothing,
 	broken::Bool = false) where {Tv, Ti, FEType <: AbstractFiniteElement, AT <: AssemblyType}
 
 	# piecewise constants are always broken
@@ -102,12 +108,22 @@ function FESpace{FEType, AT}(
 		xgrid[BFaceVolumes] = zeros(Tv, 0)
 	end
 
+	if isnothing(regions)
+		regions = unique(xgrid[CellRegions])
+	end
+
+	if regions != unique(xgrid[CellRegions])
+		dofgrid = subgrid(xgrid, regions)
+	else
+		dofgrid = xgrid
+	end
+
 	# first generate some empty FESpace
 	if name == ""
 		name = broken ? "$FEType (broken)" : "$FEType"
 	end
-	ndofs, coffset = count_ndofs(xgrid, FEType, broken)
-	FES = FESpace{Tv, Ti, FEType, AT}(name, broken, ndofs, coffset, xgrid, Dict{Type{<:AbstractGridComponent}, Any}())
+	ndofs, coffset = count_ndofs(dofgrid, FEType, broken)
+	FES = FESpace{Tv, Ti, FEType, AT}(name, broken, ndofs, coffset, xgrid, dofgrid, Dict{Type{<:AbstractGridComponent}, Any}())
 
 	@debug "Generated FESpace $name ($AT, ndofs=$ndofs)"
 
@@ -116,9 +132,8 @@ end
 
 function FESpace{FEType}(
 	xgrid::ExtendableGrid{Tv, Ti};
-	name = "",
-	broken::Bool = false) where {Tv, Ti, FEType <: AbstractFiniteElement}
-	return FESpace{FEType, ON_CELLS}(xgrid; name = name, broken = broken)
+	kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement}
+	return FESpace{FEType, ON_CELLS}(xgrid; kwargs...)
 end
 
 """