fixed hardcoded Int32 type in PointEvaluator, evalute! now evaluates …

…in global coordinates and evaluate_bary! in local coordinates
WIAS-PDELib · Nov 17, 2023 · 41c5140 · 41c5140
1 parent fadb326
commit 41c5140
Show file tree

Hide file tree

Showing 5 changed files with 81 additions and 16 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 authors = ["Christian Merdon <[email protected]>"]
-version = "0.0.13"
+version = "0.0.14"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl
@@ -125,7 +125,7 @@ include("segment_integrator.jl")
 export SegmentIntegrator, initialize!, integrate_segment!
 
 include("point_evaluator.jl")
-export PointEvaluator, evaluate!, eval_func
+export PointEvaluator, evaluate!, evaluate_bary!, eval_func
 
 include("lazy_interpolate.jl")
 export lazy_interpolate!

diff --git a/src/dofmaps.jl b/src/dofmaps.jl
@@ -387,7 +387,7 @@ function init_broken_dofmap!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Union{Type{B
 	dofmap4cellEG::Array{ParsedDofMap, 1} = Array{ParsedDofMap, 1}(undef, length(cellEG))
 	dofmap4EG::Array{ParsedDofMap, 1} = Array{ParsedDofMap, 1}(undef, length(EG))
 	for j ∈ 1:length(cellEG)
-		pattern = get_dofmap_pattern(FEType, DM, cellEG[j])
+		pattern = get_dofmap_pattern(FEType, CellDofs, cellEG[j])
 		dofmap4cellEG[j] = ParsedDofMap(pattern, ncomponents, cellEG[j])
 	end
 	for j ∈ 1:length(EG)

diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl
@@ -4,8 +4,11 @@ mutable struct PointEvaluator{Tv <: Real, UT, KFT <: Function}
 	kernel::KFT
 	BE_args::Any
 	L2G::Any
+	CF::Any
 	eval_selector::Any
+	evaluator_bary::Any
 	evaluator::Any
+	xref::Vector{Tv}
 	parameters::Dict{Symbol, Any}
 end
 
@@ -47,7 +50,7 @@ function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, k
 	parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_peval_kwargs())
 	_update_params!(parameters, kwargs)
 	@assert length(u_args) == length(ops_args)
-	PE = PointEvaluator{Tv, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, nothing, parameters)
+	PE = PointEvaluator{Tv, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, nothing, nothing, nothing, zeros(Tv, 2), parameters)
 	if sol !== nothing
 		initialize!(PE, sol)
 	end
@@ -87,13 +90,16 @@ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {
 	FES_args = [sol[j].FES for j in ind_args]
 	nargs = length(FES_args)
 	xgrid = FES_args[1].xgrid
+	Ti = eltype(xgrid[CellNodes])
 	EGs = xgrid[UniqueCellGeometries]
 	AT = ON_CELLS
 	gridAT = ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_args[1]), AT)
 	xgrid = FES_args[1].xgrid
 	itemregions = xgrid[CellRegions]
 	itemgeometries = xgrid[CellGeometries]
 
+	O.CF = CellFinder(xgrid)
+	O.xref = zeros(T, size(xgrid[Coordinates],1))
 
 	O.BE_args = Array{Array{<:FEEvaluator, 1}, 1}([])
 	O.L2G = []
@@ -115,14 +121,14 @@ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {
 	input_args = zeros(T, op_offsets_args[end])
 
 	FEATs_args = [ExtendableFEMBase.EffAT4AssemblyType(get_AT(FES_args[j]), AT) for j ∈ 1:nargs]
-	itemdofs_args::Array{Union{Adjacency{Int32}, SerialVariableTargetAdjacency{Int32}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j ∈ 1:nargs]
+	itemdofs_args::Array{Union{Adjacency{Ti}, SerialVariableTargetAdjacency{Ti}}, 1} = [FES_args[j][Dofmap4AssemblyType(FEATs_args[j])] for j ∈ 1:nargs]
 	kernel = O.kernel
 
 	function eval_selector(item)
 		return findfirst(==(itemgeometries[item]), EGs)
 	end
 
-	function _evaluate!(
+	function _evaluate_bary!(
 		result,
 		BE_args::Array{<:FEEvaluator, 1},
 		L2G::L2GTransformer,
@@ -160,7 +166,25 @@ function initialize!(O::PointEvaluator{T, UT}, sol; time = 0, kwargs...) where {
 
 		return nothing
 	end
+
+	## initialize cell finder
+	CF = CellFinder(xgrid)
+	xref = zeros(T, size(xgrid[Coordinates],1))
+	function _evaluate!(
+		result,
+		BE_args::Array{<:FEEvaluator, 1},
+		L2G::L2GTransformer,
+		x
+	)
+
+
+		# evaluate in barycentric coordinates
+		_evaluate_bary!(result, BE_args, L2G, xref, item)
+
+		return nothing
+	end
 	O.evaluator = _evaluate!
+	O.evaluator_bary = _evaluate_bary!
 	O.eval_selector = eval_selector
 
 	return nothing
@@ -169,7 +193,7 @@ end
 
 """
 ````
-function evaluate!(
+function evaluate_bary!(
 	result,
 	PE::PointEvaluator,
 	xref, 
@@ -179,7 +203,7 @@ function evaluate!(
 
 Evaluates the PointEvaluator at the specified reference coordinates in the cell with the specified item number.
 """
-function evaluate!(
+function evaluate_bary!(
 	result,
 	PE::PointEvaluator,
 	xref,
@@ -190,19 +214,58 @@ function evaluate!(
 	j = PE.eval_selector(item)
 
 	## evaluate
-	PE.evaluator(result, PE.BE_args[j], PE.L2G[j], xref, item)
+	PE.evaluator_bary(result, PE.BE_args[j], PE.L2G[j], xref, item)
+end
+
+"""
+````
+function evaluate!(
+	result,
+	PE::PointEvaluator,
+	x
+	)
+````
+
+Evaluates the PointEvaluator at the specified coordinates x.
+(To do so it internally calls CellFinder to find the cell and the barycentric
+coordinates of x and calls evaluate_bary!.)
+"""
+function evaluate!(
+	result,
+	PE::PointEvaluator,
+	x
+)
+
+	# find cell
+	item = gFindLocal!(PE.xref, PE.CF, x)
+
+	## find cell geometry id
+	j = PE.eval_selector(item)
+
+	## evaluate
+	PE.evaluator_bary(result, PE.BE_args[j], PE.L2G[j], PE.xref, item)
 end
 
 
 
+"""
+````
+function eval_func_bary(PE::PointEvaluator)
+````
+
+Yields the function (result, xref, item) -> evaluate_bary!(result,PE,xref,item).
+"""
+function eval_func_bary(PE::PointEvaluator)
+	return (result, xref, item) -> evaluate_bary!(result, PE, xref, item)
+end
+
 """
 ````
 function eval_func(PE::PointEvaluator)
 ````
 
-Yields the function (result, xref, item) -> evaluate!(result,PE,xref,item)
-that can be, e.g., used in a vectorplot.
+Yields the function (result, x) -> evaluate!(result,PE,x).
 """
 function eval_func(PE::PointEvaluator)
-	return (result, xref, item) -> evaluate!(result, PE, xref, item)
+	return (result, x) -> evaluate!(result, PE, x)
 end
diff --git a/test/test_pointevaluator.jl b/test/test_pointevaluator.jl
@@ -6,8 +6,7 @@ function run_pointevaluator_tests()
         println("======================")
         println("Testing PointEvaluator")
         println("======================")
-        error = test_pointevaluation()
-        @test error < 1e-15
+        test_pointevaluation()
     end
 end
 
@@ -25,6 +24,9 @@ function test_pointevaluation()
     x = [0.15234,0.2234] # point inside the cell
     xref = [0.0,0.0]
     cell = gFindLocal!(xref, CF, x; icellstart = 1)
-    evaluate!(eval, PE, xref, cell)
-    return abs(eval[1] - sum(x))
+    evaluate_bary!(eval, PE, xref, cell)
+    @test abs(eval[1] - sum(x)) < 1e-15
+
+    evaluate!(eval, PE, x)
+    @test abs(eval[1] - sum(x)) < 1e-15
 end