new subgrid functionalities from respective updates in ExtendableGrid…

…s/ExtendableFEMBase; new Example206 demonstrates the new features; replaced certain deepcopies by new copy functions

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
 authors = ["Christian Merdon <[email protected]>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
@@ -30,8 +30,8 @@ ExtendableFEMDiffEQExt = "DifferentialEquations"
 CommonSolve = "0.2"
 DiffResults = "1"
 DocStringExtensions = "0.8,0.9"
-ExtendableFEMBase = "^0.2, 0.3"
-ExtendableGrids = "1.0"
+ExtendableFEMBase = "^0.3.1, 0.4"
+ExtendableGrids = "1.3"
 ExtendableSparse = "1.2"
 ForwardDiff = "^0.10.35"
 GridVisualize = "1.5"

docs/make.jl

@@ -22,6 +22,7 @@ function make_all(; with_examples::Bool = true, run_examples::Bool = true, run_n
 			"Example203_PoissonProblemDG.jl",
 			"Example204_LaplaceEVProblem.jl",
 			"Example205_HeatEquation.jl",
+			"Example206_CoupledSubGridProblems.jl",
 			"Example210_LshapeAdaptivePoissonProblem.jl",
 			"Example211_LshapeAdaptiveEQPoissonProblem.jl",
 			"Example220_ReactionConvectionDiffusion.jl",

examples/Example206_CoupledSubGridProblems.jl

@@ -0,0 +1,113 @@
+#=
+
+# 206 : CoupledSubGridProblems
+([source code](@__SOURCE_URL__))
+
+This example demonstrates how to solve a coupled problem where two
+variables only live on a sub-domain and are coupled through an interface condition.
+Consider the unit square domain cut in half through on of its diagonals.
+On each subdomain a solutiong ``u_j`` of the two-dimensional Poisson problem
+```math
+\begin{aligned}
+-\Delta u & = f \quad \text{in } \Omega
+\end{aligned}
+```
+with inhomogeneous boundary conditions on the former boundaries of the full square is searched.
+Along the common boundary between the two subdomains a new interface region is assigned (appended to BFaceNodes)
+and an interface condition is assembled that couples the two solutions ``u_1`` and ``u_2``
+to each other.
+In this toy example, this interface conditions penalizes the jump between the two solutions on each side
+of the diagonal. Oberserve, that if the penalization factor ``\tau`` is large, the two solutions are
+almost equal along the interface.
+
+The computed solution(s) looks like this:
+
+![](example206.svg)
+
+Each column of the plot shows the solution, the subgrid it lives on. The last row shows the full grid.
+
+=#
+
+module Example206_CoupledSubGridProblems
+
+using ExtendableFEM
+using ExtendableFEMBase
+using ExtendableGrids
+using ExtendableSparse
+using GridVisualize
+using UnicodePlots
+using Test #
+
+
+function boundary_conditions!(result, qpinfo)
+    result[1] = 1 - qpinfo.x[1] - qpinfo.x[2] # used for both subsolutions
+end
+
+function interface_condition!(result, u, qpinfo)
+    result[1] = -(u[2] - u[1])
+    result[2] = -result[1]
+end
+
+
+function main(; μ = [1.0,1.0], f = [10,-10], τ = 1, nref = 4, order = 2, Plotter = nothing, kwargs...)
+
+	## Finite element type
+	FEType = H1Pk{1, 2, order}
+
+	## generate mesh
+	xgrid = grid_unitsquare(Triangle2D)
+
+    ## define regions
+    xgrid[CellRegions] = Int32[1,2,2,1]
+
+    ## add an interface between region 1 and 2
+    ## (one can use the BFace storages for that)
+    xgrid[BFaceNodes] = Int32[xgrid[BFaceNodes] [2 5; 5 4]]
+    append!(xgrid[BFaceRegions], [5,5])
+    xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 6)
+
+    ## refine
+    xgrid = uniform_refine(xgrid, nref)
+
+    ## define an FESpace just on region 1 and one just on region 2
+    FES1 = FESpace{FEType}(xgrid; regions = [1])
+    FES2 = FESpace{FEType}(xgrid; regions = [2])
+
+    ## define variables
+    u1 = Unknown("u1"; name = "potential in region 1")
+    u2 = Unknown("u2"; name = "potential in region 2")
+
+    ## problem description
+	PD = ProblemDescription()
+	assign_unknown!(PD, u1)
+	assign_unknown!(PD, u2)
+	assign_operator!(PD, BilinearOperator([grad(u1)]; regions = [1], factor = μ[1], kwargs...))
+	assign_operator!(PD, BilinearOperator([grad(u2)]; regions = [2], factor = μ[2], kwargs...))
+    assign_operator!(PD, LinearOperator([id(u1)]; regions = [1], factor = f[1]))
+    assign_operator!(PD, LinearOperator([id(u2)]; regions = [2], factor = f[2]))
+	assign_operator!(PD, BilinearOperator(interface_condition!, [id(u1), id(u2)]; regions = [5], factor = τ, entities = ON_FACES, kwargs...))
+	assign_operator!(PD, InterpolateBoundaryData(u1, boundary_conditions!; regions = 1:4))
+	assign_operator!(PD, InterpolateBoundaryData(u2, boundary_conditions!; regions = 1:4))
+
+    sol = solve(PD, [FES1, FES2])
+
+    plt = plot([id(u1), id(u2), dofgrid(u1), dofgrid(u2), grid(u1)], sol; Plotter = Plotter)
+
+	return sol, plt
+end
+
+generateplots = default_generateplots(Example206_CoupledSubGridProblems, "example206.svg") #hide
+
+
+function jump_l2norm!(result, u, qpinfo)
+    result[1] = (u[1] - u[2])^2
+end
+function runtests() #hide
+    ## test if jump at interface vanishes for large penalty
+	sol, plt = main(; τ = 1e9, nrefs = 2, order = 2) #hide
+    jump_integrator = ItemIntegrator(jump_l2norm!, [id(1), id(2)]; entities = ON_BFACES, regions = [5], resultdim = 1, quadorder = 4)
+    jump_error = sqrt(sum(evaluate(jump_integrator, sol)))
+    @info "||[u_1 - u_2]|| = $(jump_error)"
+	@test jump_error < 1e-8 #hide
+end #hide
+end #module

src/ExtendableFEM.jl

@@ -51,7 +51,7 @@ export print_table
 
 include("unknowns.jl")
 export Unknown
-export grid
+export grid, dofgrid
 export id, grad, hessian, div, normalflux, Δ, apply, curl1, curl2, curl3, laplace
 export id_jump, grad_jump, normalflux_jump
 

src/common_operators/bilinear_operator.jl

@@ -294,9 +294,11 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz, FE_arg
 			@info ".... building assembler for $(O.parameters[:name])"
 		end
 
+		## determine grid
+		xgrid = determine_assembly_grid(FES_test, FES_ansatz, FES_args)
+
 		## prepare assembly
 		AT = O.parameters[:entities]
-		xgrid = FES_test[1].xgrid
 		Ti = typeof(xgrid).parameters[2]
 		gridAT = ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[1]), AT)
 		itemassemblygroups = xgrid[GridComponentAssemblyGroups4AssemblyType(AT)]
@@ -404,16 +406,16 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz, FE_arg
 		if O.parameters[:parallel_groups]
 			Aj = Array{typeof(A), 1}(undef, length(EGs))
 			for j ∈ 1:length(EGs)
-				Aj[j] = deepcopy(A)
+				Aj[j] = copy(A)
 			end
 		end
 
 		FEATs_test = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[j]), AT) for j ∈ 1:ntest]
 		FEATs_ansatz = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_ansatz[j]), AT) for j ∈ 1:nansatz]
 		FEATs_args = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_args[j]), AT) for j ∈ 1:nargs]
-		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_test[j][Dofmap4AssemblyType(FEATs_test[j])] for j ∈ 1:ntest]
-		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_ansatz[j][Dofmap4AssemblyType(FEATs_ansatz[j])] for j ∈ 1:nansatz]
-		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j ∈ 1:nargs]
+		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_test[j], xgrid, FEATs_test[j]) for j = 1 : ntest]
+		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_ansatz[j], xgrid, FEATs_ansatz[j]) for j = 1 : nansatz]
+		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_args[j], xgrid, FEATs_args[j]) for j = 1 : nargs]
 		factor = O.parameters[:factor]
 		transposed_copy = O.parameters[:transposed_copy]
 		entry_tol = O.parameters[:entry_tolerance]
@@ -613,9 +615,12 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz; time =
 				@info "...... ANSATZ : $([(get_FEType(FES_ansatz[j]), O.ops_ansatz[j]) for j = 1 : nansatz])"
 			end
 		end
+
+		## determine grid
+		xgrid = determine_assembly_grid(FES_test, FES_ansatz)
+
 		## prepare assembly
 		AT = O.parameters[:entities]
-		xgrid = FES_test[1].xgrid
 		Ti = typeof(xgrid).parameters[2]
 		gridAT = ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[1]), AT)
 		itemassemblygroups = xgrid[GridComponentAssemblyGroups4AssemblyType(gridAT)]
@@ -712,14 +717,14 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz; time =
 		if O.parameters[:parallel_groups]
 			Aj = Array{typeof(A), 1}(undef, length(EGs))
 			for j ∈ 1:length(EGs)
-				Aj[j] = deepcopy(A)
+				Aj[j] = copy(A)
 			end
 		end
 
 		FEATs_test = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[j]), AT) for j ∈ 1:ntest]
 		FEATs_ansatz = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_ansatz[j]), AT) for j ∈ 1:nansatz]
-		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_test[j][Dofmap4AssemblyType(FEATs_test[j])] for j ∈ 1:ntest]
-		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_ansatz[j][Dofmap4AssemblyType(FEATs_ansatz[j])] for j ∈ 1:nansatz]
+		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_test[j], xgrid, FEATs_test[j]) for j = 1 : ntest]
+		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_ansatz[j], xgrid, FEATs_ansatz[j]) for j = 1 : nansatz]
 		factor = O.parameters[:factor]
 		transposed_copy = O.parameters[:transposed_copy]
 		entry_tol = O.parameters[:entry_tolerance]
@@ -757,6 +762,10 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz; time =
 					if !(visit_region[itemregions[item]]) || AT == ON_IFACES
 						continue
 					end
+				else
+					if length(regions) > 0
+						continue
+					end
 				end
 				QPinfos.region = itemregions[item]
 				QPinfos.item = item
@@ -820,7 +829,6 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz; time =
 				for id ∈ 1:nansatz, idt in couples_with[id]
 					Aloc[idt, id] .*= itemvolumes[item]
 					for j ∈ 1:ndofs_test[idt]
-
 						dof_j = itemdofs_test[idt][j, item] + offsets_test[idt]
 						for k ∈ 1:ndofs_ansatz[id]
 							dof_k = itemdofs_ansatz[id][k, item] + offsets_ansatz[id]

src/common_operators/bilinear_operator_dg.jl

@@ -237,10 +237,12 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz, FE_a
 			@info ".... building assembler for $(O.parameters[:name])"
 		end
 
+		## determine grid
+		xgrid = determine_assembly_grid(FES_test, FES_ansatz, FES_args)
+
 		## prepare assembly
 		AT = O.parameters[:entities]
 		@assert AT <: ON_FACES "only works for entities <: ON_FACES"
-		xgrid = FES_test[1].xgrid
 		Ti = typeof(xgrid).parameters[2]
 		gridAT = ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[1]), AT)
 		itemassemblygroups = xgrid[GridComponentAssemblyGroups4AssemblyType(AT)]
@@ -365,16 +367,16 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz, FE_a
 		if O.parameters[:parallel_groups]
 			Aj = Array{typeof(A), 1}(undef, length(EGs))
 			for j ∈ 1:length(EGs)
-				Aj[j] = deepcopy(A)
+				Aj[j] = copy(A)
 			end
 		end
 
 		FEATs_test = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[j]), ON_CELLS) for j ∈ 1:ntest]
 		FEATs_ansatz = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_ansatz[j]), ON_CELLS) for j ∈ 1:nansatz]
 		FEATs_args = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_args[j]), ON_CELLS) for j ∈ 1:nargs]
-		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_test[j][Dofmap4AssemblyType(FEATs_test[j])] for j ∈ 1:ntest]
-		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_ansatz[j][Dofmap4AssemblyType(FEATs_ansatz[j])] for j ∈ 1:nansatz]
-		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j ∈ 1:nargs]
+		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_test[j], xgrid, FEATs_test[j]) for j = 1 : ntest]
+		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_ansatz[j], xgrid, FEATs_ansatz[j]) for j = 1 : nansatz]
+		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_args[j], xgrid, FEATs_args[j]) for j = 1 : nargs]
 		factor = O.parameters[:factor]
 		transposed_copy = O.parameters[:transposed_copy]
 		entry_tol = O.parameters[:entry_tolerance]
@@ -600,10 +602,13 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz; time
 		if O.parameters[:verbosity] > 0
 			@info ".... building assembler for $(O.parameters[:name])"
 		end
+
+		## determine grid
+		xgrid = determine_assembly_grid(FES_test, FES_ansatz)
+
 		## prepare assembly
 		AT = O.parameters[:entities]
 		@assert AT <: ON_FACES "only works for entities <: ON_FACES"
-		xgrid = FES_test[1].xgrid
 		Ti = typeof(xgrid).parameters[2]
 		gridAT = ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[1]), AT)
 		itemassemblygroups = xgrid[GridComponentAssemblyGroups4AssemblyType(gridAT)]
@@ -714,14 +719,14 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz; time
 		if O.parameters[:parallel_groups]
 			Aj = Array{typeof(A), 1}(undef, length(EGs))
 			for j ∈ 1:length(EGs)
-				Aj[j] = deepcopy(A)
+				Aj[j] = copy(A)
 			end
 		end
 
 		FEATs_test = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_test[j]), ON_CELLS) for j ∈ 1:ntest]
 		FEATs_ansatz = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_ansatz[j]), ON_CELLS) for j ∈ 1:nansatz]
-		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_test[j][Dofmap4AssemblyType(FEATs_test[j])] for j ∈ 1:ntest]
-		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_ansatz[j][Dofmap4AssemblyType(FEATs_ansatz[j])] for j ∈ 1:nansatz]
+		itemdofs_test::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_test[j], xgrid, FEATs_test[j]) for j = 1 : ntest]
+		itemdofs_ansatz::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_ansatz[j], xgrid, FEATs_ansatz[j]) for j = 1 : nansatz]
 		factor = O.parameters[:factor]
 		transposed_copy = O.parameters[:transposed_copy]
 		entry_tol = O.parameters[:entry_tolerance]

src/common_operators/discface_interpolator.jl

@@ -72,9 +72,11 @@ function build_assembler!(O::FaceInterpolator{Tv}, FE_args::Array{<:FEVectorBloc
 			@info ".... building assembler for $(O.parameters[:name])"
 		end
 
+		## determine grid
+		xgrid = determine_assembly_grid(FES_args)
+
 		## prepare assembly
 		AT = ON_CELLS
-		xgrid = FES_args[1].xgrid
 		Ti = typeof(xgrid).parameters[2]
 		itemassemblygroups = xgrid[CellAssemblyGroups]
 		itemgeometries = xgrid[CellGeometries]
@@ -151,7 +153,7 @@ function build_assembler!(O::FaceInterpolator{Tv}, FE_args::Array{<:FEVectorBloc
 		#A = FEMatrix(FES_target, FES_target)
 
 		FEATs_args = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_args[j]), AT) for j ∈ 1:nargs]
-		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j ∈ 1:nargs]
+		itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [get_dofmap(FES_args[j], xgrid, FEATs_args[j]) for j = 1 : nargs]
 		facedofs_target::Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}} = FES_target[CellDofs]
 
 		O.FES_args = FES_args

src/common_operators/homogeneousdata_operator.jl

@@ -68,18 +68,19 @@ end
 function assemble!(O::HomogeneousData{UT, AT}, FES = O.FES; offset = 0, kwargs...) where {UT,AT}
 	if O.FES !== FES
 		regions = O.parameters[:regions]
+		xgrid = FES.dofgrid
 		if AT <: ON_BFACES
 			itemdofs = FES[ExtendableFEMBase.BFaceDofs]
-			itemregions = FES.xgrid[BFaceRegions]
-			uniquegeometries = FES.xgrid[UniqueBFaceGeometries]
+			itemregions = xgrid[BFaceRegions]
+			uniquegeometries = xgrid[UniqueBFaceGeometries]
 		elseif AT <: ON_CELLS
 			itemdofs = FES[ExtendableFEMBase.CellDofs]
-			itemregions = FES.xgrid[CellRegions]
-			uniquegeometries = FES.xgrid[UniqueCellGeometries]
+			itemregions = xgrid[CellRegions]
+			uniquegeometries = xgrid[UniqueCellGeometries]
 		elseif AT <: ON_FACES
 			itemdofs = FES[ExtendableFEMBase.FaceDofs]
-			itemregions = FES.xgrid[FaceRegions]
-			uniquegeometries = FES.xgrid[UniqueFaceGeometries]
+			itemregions = xgrid[FaceRegions]
+			uniquegeometries = xgrid[UniqueFaceGeometries]
 		end
 		nitems = num_sources(itemdofs)
 		ndofs4item = max_num_targets_per_source(itemdofs)

src/common_operators/interpolateboundarydata_operator.jl

@@ -79,10 +79,11 @@ function assemble!(O::InterpolateBoundaryData, FES = O.FES; time = 0, offset = 0
 	bfaces::Array{Int,1} = O.bfaces
 	if O.FES !== FES
 		bddata = FEVector(FES)
-		Ti = eltype(FES.xgrid[CellNodes])
+		xgrid = FES.dofgrid
+		Ti = eltype(xgrid[CellNodes])
 		bfacedofs::Adjacency{Ti} = FES[BFaceDofs]
-		bfacefaces = FES.xgrid[BFaceFaces]
-		bfaceregions = FES.xgrid[BFaceRegions]
+		bfacefaces = xgrid[BFaceFaces]
+		bfaceregions = xgrid[BFaceRegions]
 		nbfaces = num_sources(bfacedofs)
 		ndofs4bface = max_num_targets_per_source(bfacedofs)
 		bdofs = []
@@ -107,7 +108,8 @@ function assemble!(O::InterpolateBoundaryData, FES = O.FES; time = 0, offset = 0
 		bfaces = O.bfaces
 		if FES.broken
 			FEType = eltype(FES)
-			FESc = FESpace{FEType}(FES.xgrid)
+			xgrid = FES.dofgrid
+			FESc = FESpace{FEType}(xgrid)
 			Targetc = FEVector(FESc)
 			interpolate!(Targetc[1], FESc, ON_FACES, data; items = bfaces, time = time, params = O.parameters[:params], bonus_quadorder = O.parameters[:bonus_quadorder])
 			bfacedofs = FES[BFaceDofs]