From 6064415c98c17a5d0599e0291dca31cd170d6023 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 25 Jan 2021 20:36:14 -0500 Subject: [PATCH 01/46] Allow evaluating FE functions at arbitrary points This is a first cut only, still very inefficient and restrictive. --- Project.toml | 4 ++++ src/FESpaces/FESpaces.jl | 3 +++ 2 files changed, 7 insertions(+) diff --git a/Project.toml b/Project.toml index 1963d7e5e..9b6948549 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.15.1" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" +Bernstein = "e9b0fb4c-9cb7-4f61-9c14-701a41c684d7" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -18,6 +19,7 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -27,6 +29,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] AbstractTrees = "0.3.3" BSON = "0.2.5" +Bernstein = "1.0.1" BlockArrays = "0.12.12, 0.13" Combinatorics = "1.0.0" DocStringExtensions = "0.8.1" @@ -38,6 +41,7 @@ JLD2 = "0.1.11, 0.3" JSON = "0.21.0" LineSearches = "7.0.1" NLsolve = "4.3.0" +NearestNeighbors = "0.4.8" QuadGK = "2.3.1, 2.4" StaticArrays = "0.12.1, 1.0" WriteVTK = "1.7, 1.8" diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 98ed37d29..5a03c639a 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -10,7 +10,10 @@ using Test using FillArrays using BlockArrays using SparseArrays +using StaticArrays using LinearAlgebra +using NearestNeighbors +using Bernstein using Gridap.Helpers using Gridap.Algebra From f123693b5e3a15b406b7440da1cadea7646bdcd5 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 26 Jan 2021 08:35:27 -0500 Subject: [PATCH 02/46] Add function to evaluate FE functions --- src/FESpaces/SingleFieldFESpaces.jl | 76 +++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/src/FESpaces/SingleFieldFESpaces.jl b/src/FESpaces/SingleFieldFESpaces.jl index 0a29e73f2..e1f14560f 100644 --- a/src/FESpaces/SingleFieldFESpaces.jl +++ b/src/FESpaces/SingleFieldFESpaces.jl @@ -148,6 +148,82 @@ get_free_values(f::SingleFieldFEFunction) = f.free_values get_cell_dof_values(f::SingleFieldFEFunction) = f.cell_dof_values get_fe_space(f::SingleFieldFEFunction) = f.fe_space +function (f::SingleFieldFEFunction)(x) + # Create k-d tree for mesh + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + @assert !isempty(node_coordinates) + sample_coord = node_coordinates[1] + D = length(sample_coord) + T = eltype(sample_coord) + @assert D > 0 + kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) + kdtree::KDTree{SVector{D,T},Euclidean,T} + + # Find neighbouring cells + id,dist = nn(kdtree, SVector(x.data)) + cells = Int[] + # TODO: Use better algorithm + # Francesc Verdugo @fverdugo 14:07 + # You have to start with a DiscreteModel, then + # topo=get_grid_topology(model) + # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + for (c,ids) in enumerate(trian.cell_node_ids) + id ∈ ids && push!(cells, c) + end + @assert !isempty(cells) + + # Find polytope + reffes = trian.reffes + @assert length(reffes) == 1 + reffe = reffes[1].reffe + polytope = reffe.polytope + extrusion = polytope.extrusion + @assert length(extrusion) == D + + # Loop over neighbouring cells and calculate distances + dists = T[] + + if all(e == TET_AXIS for e in extrusion) + # Calculate barycentric coordinates for a point x in the given cell + function bary(cell, x) + local ids = trian.cell_node_ids[cell] + @assert length(ids) == D+1 + local ys = zero(SVector{D+1,SVector{D,T}}) + for (n,id) in enumerate(ids) + local y = SVector(node_coordinates[id].data) + ys = setindex(ys, y, n) + end + local λ = cartesian2barycentric(ys, SVector(x.data)) + return λ + end + + for c in cells + λ = bary(c, x) + # Negative distances are outside, positive distance are + # inside the simplex + dist = min(minimum(λ), 1 - maximum(λ)) + push!(dists, dist) + end + + dist, m = findmax(dists) + cell = cells[m] + λ = bary(cell, x) + + # Ensure the point is inside one of the cells, up to round-off errors + @assert dist ≥ -1000 * eps(T) + + # TODO: Actually evaluate basis functions + values = f.cell_dof_values[cell] + fx = sum(λ[n] * values[n] for n in 1:D+1) + + else + @assert false "Unsupported polytope" + end + + return fx +end + """ FEFunction( fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector) From 7e785be4087d93b9b787e6b2b8d0721fcb2aa092 Mon Sep 17 00:00:00 2001 From: eneiva Date: Mon, 4 Jan 2021 13:45:48 +0100 Subject: [PATCH 03/46] New function inverse_map --- src/Fields/AffineMaps.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Fields/AffineMaps.jl b/src/Fields/AffineMaps.jl index b1233d74d..6f156be19 100644 --- a/src/Fields/AffineMaps.jl +++ b/src/Fields/AffineMaps.jl @@ -69,3 +69,11 @@ function lazy_map( cell_∇∇bt = lazy_map(transpose,cell_∇∇b) cell_∇∇bt end + +function inverse_map(f::AffineMap) + Jt = f.gradient + y0 = f.origin + invJt = inv(Jt) + x0 = -y0⋅invJt + AffineMap(invJt,x0) +end From 027fa14e45854e14c14f5958e15f4e6df639c8e9 Mon Sep 17 00:00:00 2001 From: eneiva Date: Tue, 19 Jan 2021 12:00:53 +0100 Subject: [PATCH 04/46] bugfixed in CellField change_domain --- src/CellData/CellFields.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 1175d8301..57f6fac2c 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -113,7 +113,7 @@ function change_domain(a::CellField,::ReferenceDomain,::PhysicalDomain) trian = get_triangulation(a) cell_map = get_cell_map(trian) cell_invmap = lazy_map(inverse_map,cell_map) - cell_field_ref = get_data(cell_field) + cell_field_ref = get_data(a) cell_field_phys = lazy_map(Broadcasting(∘),cell_field_ref,cell_invmap) GenericCellField(cell_field_phys,trian,PhysicalDomain()) end From 81bf213c4b7d7dd688862e9b6e22cdd385e5f66e Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 26 Jan 2021 10:03:00 -0500 Subject: [PATCH 05/46] Move code to the right location --- src/CellData/CellData.jl | 3 + src/CellData/CellFields.jl | 86 +++++++++++++++++++++++++---- src/FESpaces/FESpaces.jl | 3 - src/FESpaces/SingleFieldFESpaces.jl | 76 ------------------------- 4 files changed, 79 insertions(+), 89 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 4335b9b58..8fb527e84 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -8,6 +8,9 @@ module CellData using Test using DocStringExtensions using FillArrays +using Bernstein +using NearestNeighbors +using StaticArrays using Gridap.Helpers using Gridap.Arrays diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 57f6fac2c..c3a86745b 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -205,13 +205,79 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() (a::CellField)(x) = evaluate(a,x) function evaluate!(cache,f::CellField,x::Point) - @notimplemented """\n - Evaluation of a CellField at a given Point is not implemented yet. + # Create k-d tree for mesh + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + @assert !isempty(node_coordinates) + sample_coord = node_coordinates[1] + D = length(sample_coord) + T = eltype(sample_coord) + @assert D > 0 + kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) + kdtree::KDTree{SVector{D,T},Euclidean,T} + + # Find neighbouring cells + id,dist = nn(kdtree, SVector(x.data)) + cells = Int[] + # TODO: Use better algorithm + # Francesc Verdugo @fverdugo 14:07 + # You have to start with a DiscreteModel, then + # topo=get_grid_topology(model) + # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + for (c,ids) in enumerate(trian.cell_node_ids) + id ∈ ids && push!(cells, c) + end + @assert !isempty(cells) + + # Find polytope + reffes = trian.reffes + @assert length(reffes) == 1 + reffe = reffes[1].reffe + polytope = reffe.polytope + extrusion = polytope.extrusion + @assert length(extrusion) == D + + # Loop over neighbouring cells and calculate distances + dists = T[] + + if all(e == TET_AXIS for e in extrusion) + # Calculate barycentric coordinates for a point x in the given cell + function bary(cell, x) + local ids = trian.cell_node_ids[cell] + @assert length(ids) == D+1 "simplex has correct number of vertices" + local ys = zero(SVector{D+1,SVector{D,T}}) + for (n,id) in enumerate(ids) + local y = SVector(node_coordinates[id].data) + ys = setindex(ys, y, n) + end + local λ = cartesian2barycentric(ys, SVector(x.data)) + return λ + end + + for c in cells + λ = bary(c, x) + # Negative distances are outside, positive distance are + # inside the simplex + dist = min(minimum(λ), 1 - maximum(λ)) + push!(dists, dist) + end + + dist, m = findmax(dists) + cell = cells[m] + λ = bary(cell, x) + + # Ensure the point is inside one of the cells, up to round-off errors + @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" + + # TODO: Actually evaluate basis functions + values = f.cell_dof_values[cell] + fx = sum(λ[n] * values[n] for n in 1:D+1) - This is a feature that we want to have at some point in Gridap. - If you are ready to help with this implementation, please contact the - Gridap administrators. - """ + else + @assert false "Unsupported polytope" + end + + return fx end function evaluate!(cache,f::CellField,x::CellPoint) @@ -324,11 +390,11 @@ struct OperationCellField{DS} <: CellField r = Fields.BroadcastingFieldOpMap(op.op)(axi...) catch @unreachable """\n - It is not possible to perform operation $(op.op) on the given cell fields. + It is not possible to perform the operation "$(op.op)" on the given cell fields. - See the catched error for more information. (If you are using Visual - Studio Code REPL you might not see the catched error, please use the - command-line REPL). + See the caught error for more information. (If you are using the Visual + Studio Code REPL you might not see the caught error, please use the + command-line REPL instead). """ end end diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 5a03c639a..98ed37d29 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -10,10 +10,7 @@ using Test using FillArrays using BlockArrays using SparseArrays -using StaticArrays using LinearAlgebra -using NearestNeighbors -using Bernstein using Gridap.Helpers using Gridap.Algebra diff --git a/src/FESpaces/SingleFieldFESpaces.jl b/src/FESpaces/SingleFieldFESpaces.jl index e1f14560f..0a29e73f2 100644 --- a/src/FESpaces/SingleFieldFESpaces.jl +++ b/src/FESpaces/SingleFieldFESpaces.jl @@ -148,82 +148,6 @@ get_free_values(f::SingleFieldFEFunction) = f.free_values get_cell_dof_values(f::SingleFieldFEFunction) = f.cell_dof_values get_fe_space(f::SingleFieldFEFunction) = f.fe_space -function (f::SingleFieldFEFunction)(x) - # Create k-d tree for mesh - trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - @assert !isempty(node_coordinates) - sample_coord = node_coordinates[1] - D = length(sample_coord) - T = eltype(sample_coord) - @assert D > 0 - kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) - kdtree::KDTree{SVector{D,T},Euclidean,T} - - # Find neighbouring cells - id,dist = nn(kdtree, SVector(x.data)) - cells = Int[] - # TODO: Use better algorithm - # Francesc Verdugo @fverdugo 14:07 - # You have to start with a DiscreteModel, then - # topo=get_grid_topology(model) - # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) - for (c,ids) in enumerate(trian.cell_node_ids) - id ∈ ids && push!(cells, c) - end - @assert !isempty(cells) - - # Find polytope - reffes = trian.reffes - @assert length(reffes) == 1 - reffe = reffes[1].reffe - polytope = reffe.polytope - extrusion = polytope.extrusion - @assert length(extrusion) == D - - # Loop over neighbouring cells and calculate distances - dists = T[] - - if all(e == TET_AXIS for e in extrusion) - # Calculate barycentric coordinates for a point x in the given cell - function bary(cell, x) - local ids = trian.cell_node_ids[cell] - @assert length(ids) == D+1 - local ys = zero(SVector{D+1,SVector{D,T}}) - for (n,id) in enumerate(ids) - local y = SVector(node_coordinates[id].data) - ys = setindex(ys, y, n) - end - local λ = cartesian2barycentric(ys, SVector(x.data)) - return λ - end - - for c in cells - λ = bary(c, x) - # Negative distances are outside, positive distance are - # inside the simplex - dist = min(minimum(λ), 1 - maximum(λ)) - push!(dists, dist) - end - - dist, m = findmax(dists) - cell = cells[m] - λ = bary(cell, x) - - # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≥ -1000 * eps(T) - - # TODO: Actually evaluate basis functions - values = f.cell_dof_values[cell] - fx = sum(λ[n] * values[n] for n in 1:D+1) - - else - @assert false "Unsupported polytope" - end - - return fx -end - """ FEFunction( fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector) From 0fd60ffd79e926372870b73a603bdbd9836393c2 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 26 Jan 2021 11:15:17 -0500 Subject: [PATCH 06/46] Reformat code --- src/CellData/CellFields.jl | 155 ++++++++++++++++++++----------------- 1 file changed, 86 insertions(+), 69 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index c3a86745b..602a4a85b 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -204,80 +204,97 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() (a::CellField)(x) = evaluate(a,x) +struct CellFieldCache + kdtree::KDTree{SVector{D,T},Euclidean,T} where {D,T} +end + +function return_cache(f::CellField,x::Point) + # Create k-d tree for mesh + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + @assert !isempty(node_coordinates) + sample_coord = node_coordinates[1] + D = length(sample_coord) + T = eltype(sample_coord) + @assert D > 0 + kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) + kdtree::KDTree{SVector{D,T},Euclidean,T} + return CellFieldCache(kdtree) +end + function evaluate!(cache,f::CellField,x::Point) - # Create k-d tree for mesh - trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - @assert !isempty(node_coordinates) - sample_coord = node_coordinates[1] - D = length(sample_coord) - T = eltype(sample_coord) - @assert D > 0 - kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) - kdtree::KDTree{SVector{D,T},Euclidean,T} - - # Find neighbouring cells - id,dist = nn(kdtree, SVector(x.data)) - cells = Int[] - # TODO: Use better algorithm - # Francesc Verdugo @fverdugo 14:07 - # You have to start with a DiscreteModel, then - # topo=get_grid_topology(model) - # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) - for (c,ids) in enumerate(trian.cell_node_ids) - id ∈ ids && push!(cells, c) + cache::CellFieldCache + + # Get k-d tree for mesh from cache + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + @assert !isempty(node_coordinates) + sample_coord = node_coordinates[1] + D = length(sample_coord) + T = eltype(sample_coord) + @assert D > 0 + kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} + + # Find neighbouring cells + id,dist = nn(kdtree, SVector(x.data)) + cells = Int[] + # TODO: Use better algorithm + # Francesc Verdugo @fverdugo 14:07 + # You have to start with a DiscreteModel, then + # topo=get_grid_topology(model) + # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + for (c,ids) in enumerate(trian.cell_node_ids) + id ∈ ids && push!(cells, c) + end + @assert !isempty(cells) + + # Find polytope + reffes = trian.reffes + @assert length(reffes) == 1 + reffe = reffes[1].reffe + polytope = reffe.polytope + extrusion = polytope.extrusion + @assert length(extrusion) == D + + # Loop over neighbouring cells and calculate distances + dists = T[] + + if all(e == TET_AXIS for e in extrusion) + # Calculate barycentric coordinates for a point x in the given cell + function bary(cell, x) + local ids = trian.cell_node_ids[cell] + @assert length(ids) == D+1 "simplex has correct number of vertices" + local ys = zero(SVector{D+1,SVector{D,T}}) + for (n,id) in enumerate(ids) + local y = SVector(node_coordinates[id].data) + ys = setindex(ys, y, n) + end + local λ = cartesian2barycentric(ys, SVector(x.data)) + return λ end - @assert !isempty(cells) - - # Find polytope - reffes = trian.reffes - @assert length(reffes) == 1 - reffe = reffes[1].reffe - polytope = reffe.polytope - extrusion = polytope.extrusion - @assert length(extrusion) == D - - # Loop over neighbouring cells and calculate distances - dists = T[] - - if all(e == TET_AXIS for e in extrusion) - # Calculate barycentric coordinates for a point x in the given cell - function bary(cell, x) - local ids = trian.cell_node_ids[cell] - @assert length(ids) == D+1 "simplex has correct number of vertices" - local ys = zero(SVector{D+1,SVector{D,T}}) - for (n,id) in enumerate(ids) - local y = SVector(node_coordinates[id].data) - ys = setindex(ys, y, n) - end - local λ = cartesian2barycentric(ys, SVector(x.data)) - return λ - end - - for c in cells - λ = bary(c, x) - # Negative distances are outside, positive distance are - # inside the simplex - dist = min(minimum(λ), 1 - maximum(λ)) - push!(dists, dist) - end - - dist, m = findmax(dists) - cell = cells[m] - λ = bary(cell, x) - - # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" - - # TODO: Actually evaluate basis functions - values = f.cell_dof_values[cell] - fx = sum(λ[n] * values[n] for n in 1:D+1) - else - @assert false "Unsupported polytope" + for c in cells + λ = bary(c, x) + # Negative distances are outside, positive distance are inside the simplex + dist = min(minimum(λ), 1 - maximum(λ)) + push!(dists, dist) end - return fx + dist, m = findmax(dists) + cell = cells[m] + λ = bary(cell, x) + + # Ensure the point is inside one of the cells, up to round-off errors + @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" + + # TODO: Actually evaluate basis functions + values = f.cell_dof_values[cell] + fx = sum(λ[n] * values[n] for n in 1:D+1) + else + @assert false "Unsupported polytope" + end + + return fx end function evaluate!(cache,f::CellField,x::CellPoint) From b6a18d4ed8141c26f27734e21c6c177e0a274570 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 26 Jan 2021 11:44:51 -0500 Subject: [PATCH 07/46] Import return_cache --- src/CellData/CellData.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 8fb527e84..291513177 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -24,6 +24,7 @@ using Gridap.Integration import Gridap.Arrays: lazy_append import Gridap.Arrays: get_array import Gridap.Arrays: evaluate! +import Gridap.Arrays: return_cache import Gridap.Fields: gradient import Gridap.Fields: ∇∇ import Gridap.Fields: integrate From 06b3524630b2761db68b69810d82a92c4a32c7ba Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 27 Jan 2021 09:08:21 -0500 Subject: [PATCH 08/46] Remove dependency on Bernstein.jl --- Project.toml | 2 -- src/CellData/CellData.jl | 1 - src/CellData/CellFields.jl | 35 +++++++++++++++++++++++++++++++++++ 3 files changed, 35 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 9b6948549..e01f3ef41 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ version = "0.15.1" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" -Bernstein = "e9b0fb4c-9cb7-4f61-9c14-701a41c684d7" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -29,7 +28,6 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] AbstractTrees = "0.3.3" BSON = "0.2.5" -Bernstein = "1.0.1" BlockArrays = "0.12.12, 0.13" Combinatorics = "1.0.0" DocStringExtensions = "0.8.1" diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 291513177..4fc9d11a3 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -8,7 +8,6 @@ module CellData using Test using DocStringExtensions using FillArrays -using Bernstein using NearestNeighbors using StaticArrays diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 602a4a85b..50a29d02f 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -297,6 +297,41 @@ function evaluate!(cache,f::CellField,x::Point) return fx end +""" + λ = cartesian2barycentric(s, p) + +Convert Cartesian to barycentric coordinates. + +# Arguments +- `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` + vertices in `D` dimensions. +- `p`: Point in Cartesian coordinates + +# Result +- `λ`: Point in barycentric coordinates +""" +function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} + @assert N ≤ D + 1 + # Algorithm as described on + # , + # section "Conversion between barycentric and Cartesian coordinates" + A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) + b = SVector{D+1}(p..., T(1)) + return A \ b +end + +function cartesian2barycentric(s::SVector{N,SVector{D,T}}, + p::SVector{D,T}) where {D,N,T} + return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), + p) +end + +function barycentric2cartesian(s::SVector{N,SVector{D,T}}, + λ::SVector{N,T}) where {D,N,T} + @assert N ≤ D + 1 + return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) +end + function evaluate!(cache,f::CellField,x::CellPoint) _f, _x = _to_common_domain(f,x) cell_field = get_data(_f) From c363d2dd20381f997660728ef261b237da293aae Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 27 Jan 2021 09:09:25 -0500 Subject: [PATCH 09/46] Change sign convention for distance from simplex --- src/CellData/CellFields.jl | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 50a29d02f..68d6e1c0a 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -235,7 +235,7 @@ function evaluate!(cache,f::CellField,x::Point) @assert D > 0 kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} - # Find neighbouring cells + # Find all neighbouring cells id,dist = nn(kdtree, SVector(x.data)) cells = Int[] # TODO: Use better algorithm @@ -273,19 +273,27 @@ function evaluate!(cache,f::CellField,x::Point) return λ end + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # negative for exactly one cell and positive for all other ones. + # Due to round-off, the distance can be slightly negative or + # slightly positive for points near cell boundaries, in particular + # near vertices. In this case, choose the cell with the smallest + # distance, and check that the distance (if positive) is at most + # at round-off level. for c in cells λ = bary(c, x) - # Negative distances are outside, positive distance are inside the simplex - dist = min(minimum(λ), 1 - maximum(λ)) + # Positive distances are outside, negative distance are inside the simplex + dist = min(-minimum(λ), maximum(λ) - 1) push!(dists, dist) end - dist, m = findmax(dists) + dist, m = findmin(dists) cell = cells[m] λ = bary(cell, x) # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" + @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" # TODO: Actually evaluate basis functions values = f.cell_dof_values[cell] From 489a377f7f3ddbd3753ff4c0268c2755ab575851 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 27 Jan 2021 09:19:51 -0500 Subject: [PATCH 10/46] Simplify finding minimum distance from simplices --- src/CellData/CellFields.jl | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index c80eb604c..17de65d2e 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -257,8 +257,6 @@ function evaluate!(cache,f::CellField,x::Point) @assert length(extrusion) == D # Loop over neighbouring cells and calculate distances - dists = T[] - if all(e == TET_AXIS for e in extrusion) # Calculate barycentric coordinates for a point x in the given cell function bary(cell, x) @@ -281,16 +279,22 @@ function evaluate!(cache,f::CellField,x::Point) # near vertices. In this case, choose the cell with the smallest # distance, and check that the distance (if positive) is at most # at round-off level. + mindist = T(Inf) + minc = 0 + minλ = zero(SVector{D+1,T}) for c in cells λ = bary(c, x) - # Positive distances are outside, negative distance are inside the simplex - dist = min(-minimum(λ), maximum(λ) - 1) - push!(dists, dist) + dist = max(-minimum(λ), maximum(λ) - 1) + if dist < mindist + minc = c + minλ = λ + mindist = dist + end end - - dist, m = findmin(dists) - cell = cells[m] - λ = bary(cell, x) + @assert 1 ≤ minc + cell = minc + λ = minλ + dist = mindist # Ensure the point is inside one of the cells, up to round-off errors @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" From 01ef2f18d616ac314a09ef7c6654c6fc645c57b2 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 27 Jan 2021 11:10:53 -0500 Subject: [PATCH 11/46] Calculate k-d tree only when needed --- src/CellData/CellFields.jl | 63 +++++++++++++++++++------------------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 17de65d2e..297e56aa3 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -204,28 +204,16 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() (a::CellField)(x) = evaluate(a,x) -struct CellFieldCache - kdtree::KDTree{SVector{D,T},Euclidean,T} where {D,T} +mutable struct CellFieldCache + kdtree::Union{Nothing, KDTree{SVector{D,T},Euclidean,T} where {D,T}} end -function return_cache(f::CellField,x::Point) - # Create k-d tree for mesh - trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - @assert !isempty(node_coordinates) - sample_coord = node_coordinates[1] - D = length(sample_coord) - T = eltype(sample_coord) - @assert D > 0 - kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) - kdtree::KDTree{SVector{D,T},Euclidean,T} - return CellFieldCache(kdtree) -end +return_cache(f::CellField,x::Point) = CellFieldCache(nothing) function evaluate!(cache,f::CellField,x::Point) cache::CellFieldCache - # Get k-d tree for mesh from cache + # Examine triangulation trian = get_triangulation(f) node_coordinates = get_node_coordinates(trian) @assert !isempty(node_coordinates) @@ -233,20 +221,6 @@ function evaluate!(cache,f::CellField,x::Point) D = length(sample_coord) T = eltype(sample_coord) @assert D > 0 - kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} - - # Find all neighbouring cells - id,dist = nn(kdtree, SVector(x.data)) - cells = Int[] - # TODO: Use better algorithm - # Francesc Verdugo @fverdugo 14:07 - # You have to start with a DiscreteModel, then - # topo=get_grid_topology(model) - # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) - for (c,ids) in enumerate(trian.cell_node_ids) - id ∈ ids && push!(cells, c) - end - @assert !isempty(cells) # Find polytope reffes = trian.reffes @@ -256,8 +230,27 @@ function evaluate!(cache,f::CellField,x::Point) extrusion = polytope.extrusion @assert length(extrusion) == D - # Loop over neighbouring cells and calculate distances if all(e == TET_AXIS for e in extrusion) + + # Get k-d tree for mesh from cache + if cache.kdtree === nothing + cache.kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) + end + kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} + + # Find all neighbouring cells + id,dist = nn(kdtree, SVector(x.data)) + cells = Int[] + # TODO: Use better algorithm + # Francesc Verdugo @fverdugo 14:07 + # You have to start with a DiscreteModel, then + # topo=get_grid_topology(model) + # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + for (c,ids) in enumerate(trian.cell_node_ids) + id ∈ ids && push!(cells, c) + end + @assert !isempty(cells) + # Calculate barycentric coordinates for a point x in the given cell function bary(cell, x) local ids = trian.cell_node_ids[cell] @@ -302,8 +295,16 @@ function evaluate!(cache,f::CellField,x::Point) # TODO: Actually evaluate basis functions values = f.cell_dof_values[cell] fx = sum(λ[n] * values[n] for n in 1:D+1) + + elseif all(e == HEX_AXIS for e in extrusion) + + # Find cell + @assert false "n-cube not implemented" + else + @assert false "Unsupported polytope" + end return fx From 5103e07e1f5fc9830bd707fb0198642a950f5185 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 28 Jan 2021 11:36:41 -0500 Subject: [PATCH 12/46] simplexify: return tet simplex instead of hex simplex in 1d --- src/ReferenceFEs/ExtrusionPolytopes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ReferenceFEs/ExtrusionPolytopes.jl b/src/ReferenceFEs/ExtrusionPolytopes.jl index 43e089583..0d8f0a965 100644 --- a/src/ReferenceFEs/ExtrusionPolytopes.jl +++ b/src/ReferenceFEs/ExtrusionPolytopes.jl @@ -230,7 +230,7 @@ function simplexify(p::ExtrusionPolytope{0}) end function simplexify(p::ExtrusionPolytope{1}) - [[1,2],], SEGMENT + [[1,2],], Polytope(TET_AXIS) end function simplexify(p::ExtrusionPolytope{2}) From bc6d04956b5a660eeb338a0e0c8e04bcf111b48d Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 28 Jan 2021 11:38:03 -0500 Subject: [PATCH 13/46] Move evaluation function to MultiFields module --- Project.toml | 1 + src/CellData/CellData.jl | 6 +- src/CellData/CellFields.jl | 141 ---------------- src/MultiField/MultiField.jl | 8 +- src/MultiField/MultiFieldCellFields.jl | 156 ++++++++++++++++++ test/CellDataTests/CellFieldsTests.jl | 6 +- .../MultiFieldCellFieldsTests.jl | 35 ++++ 7 files changed, 204 insertions(+), 149 deletions(-) diff --git a/Project.toml b/Project.toml index e01f3ef41..eba2be63e 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 4fc9d11a3..c2365d933 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -8,17 +8,15 @@ module CellData using Test using DocStringExtensions using FillArrays -using NearestNeighbors -using StaticArrays using Gridap.Helpers -using Gridap.Arrays using Gridap.Algebra +using Gridap.Arrays using Gridap.TensorValues using Gridap.Fields +using Gridap.Integration using Gridap.ReferenceFEs using Gridap.Geometry -using Gridap.Integration import Gridap.Arrays: lazy_append import Gridap.Arrays: get_array diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 297e56aa3..c866b2697 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -204,147 +204,6 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() (a::CellField)(x) = evaluate(a,x) -mutable struct CellFieldCache - kdtree::Union{Nothing, KDTree{SVector{D,T},Euclidean,T} where {D,T}} -end - -return_cache(f::CellField,x::Point) = CellFieldCache(nothing) - -function evaluate!(cache,f::CellField,x::Point) - cache::CellFieldCache - - # Examine triangulation - trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - @assert !isempty(node_coordinates) - sample_coord = node_coordinates[1] - D = length(sample_coord) - T = eltype(sample_coord) - @assert D > 0 - - # Find polytope - reffes = trian.reffes - @assert length(reffes) == 1 - reffe = reffes[1].reffe - polytope = reffe.polytope - extrusion = polytope.extrusion - @assert length(extrusion) == D - - if all(e == TET_AXIS for e in extrusion) - - # Get k-d tree for mesh from cache - if cache.kdtree === nothing - cache.kdtree = KDTree(map(nc -> SVector(nc.data), node_coordinates)) - end - kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} - - # Find all neighbouring cells - id,dist = nn(kdtree, SVector(x.data)) - cells = Int[] - # TODO: Use better algorithm - # Francesc Verdugo @fverdugo 14:07 - # You have to start with a DiscreteModel, then - # topo=get_grid_topology(model) - # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) - for (c,ids) in enumerate(trian.cell_node_ids) - id ∈ ids && push!(cells, c) - end - @assert !isempty(cells) - - # Calculate barycentric coordinates for a point x in the given cell - function bary(cell, x) - local ids = trian.cell_node_ids[cell] - @assert length(ids) == D+1 "simplex has correct number of vertices" - local ys = zero(SVector{D+1,SVector{D,T}}) - for (n,id) in enumerate(ids) - local y = SVector(node_coordinates[id].data) - ys = setindex(ys, y, n) - end - local λ = cartesian2barycentric(ys, SVector(x.data)) - return λ - end - - # Calculate the distance from the point to all the cells. Without - # round-off, and with non-overlapping cells, the distance would be - # negative for exactly one cell and positive for all other ones. - # Due to round-off, the distance can be slightly negative or - # slightly positive for points near cell boundaries, in particular - # near vertices. In this case, choose the cell with the smallest - # distance, and check that the distance (if positive) is at most - # at round-off level. - mindist = T(Inf) - minc = 0 - minλ = zero(SVector{D+1,T}) - for c in cells - λ = bary(c, x) - dist = max(-minimum(λ), maximum(λ) - 1) - if dist < mindist - minc = c - minλ = λ - mindist = dist - end - end - @assert 1 ≤ minc - cell = minc - λ = minλ - dist = mindist - - # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" - - # TODO: Actually evaluate basis functions - values = f.cell_dof_values[cell] - fx = sum(λ[n] * values[n] for n in 1:D+1) - - elseif all(e == HEX_AXIS for e in extrusion) - - # Find cell - @assert false "n-cube not implemented" - - else - - @assert false "Unsupported polytope" - - end - - return fx -end - -""" - λ = cartesian2barycentric(s, p) - -Convert Cartesian to barycentric coordinates. - -# Arguments -- `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` - vertices in `D` dimensions. -- `p`: Point in Cartesian coordinates - -# Result -- `λ`: Point in barycentric coordinates -""" -function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} - @assert N ≤ D + 1 - # Algorithm as described on - # , - # section "Conversion between barycentric and Cartesian coordinates" - A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) - b = SVector{D+1}(p..., T(1)) - return A \ b -end - -function cartesian2barycentric(s::SVector{N,SVector{D,T}}, - p::SVector{D,T}) where {D,N,T} - return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), - p) -end - -function barycentric2cartesian(s::SVector{N,SVector{D,T}}, - λ::SVector{N,T}) where {D,N,T} - @assert N ≤ D + 1 - return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) -end - function evaluate!(cache,f::CellField,x::CellPoint) _f, _x = _to_common_domain(f,x) cell_field = get_data(_f) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index a399c175b..67f5ad0c4 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -16,13 +16,19 @@ using Gridap.TensorValues using Gridap.CellData using Gridap.Fields -using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis +using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.Arrays: BlockArrayCooMap +using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, get_node_coordinates using FillArrays using SparseArrays using LinearAlgebra using BlockArrays +using NearestNeighbors +using StaticArrays + +import Gridap.Arrays: evaluate! +import Gridap.Arrays: return_cache export num_fields export compute_field_offsets diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 92750b220..13c9c7bed 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -39,3 +39,159 @@ num_fields(a::MultiFieldCellField) = length(a.single_fields) Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i] Base.iterate(a::MultiFieldCellField) = iterate(a.single_fields) Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state) + +# Evaluation of CellFields + +# This code lives here (instead of in the CellData module) because we +# need access to MultiField functionality. + +mutable struct CellFieldCache + kdtree::Union{Nothing, KDTree{SVector{D,T},Euclidean,T} where {D,T}} +end + +return_cache(f::CellField,x::Point) = CellFieldCache(nothing) + +function evaluate!(cache,f::CellField,x::Point) + cache::CellFieldCache + + # Examine triangulation + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + @assert !isempty(node_coordinates) + sample_coord = node_coordinates[1] + D = length(sample_coord) + T = eltype(sample_coord) + @assert D > 0 + + # Find polytope + reffes = trian.reffes + @assert length(reffes) == 1 + reffe = reffes[1].reffe + polytope = reffe.polytope + extrusion = polytope.extrusion + @assert length(extrusion) == D + + if all(e == TET_AXIS for e in extrusion) + + # Get k-d tree for mesh from cache + if cache.kdtree === nothing + cache.kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) + end + kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} + + # Find all neighbouring cells + id,dist = nn(kdtree, SVector(Tuple(x))) + cells = Int[] + # TODO: Use better algorithm + # Francesc Verdugo @fverdugo 14:07 + # You have to start with a DiscreteModel, then + # topo=get_grid_topology(model) + # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + for (c,ids) in enumerate(trian.cell_node_ids) + id ∈ ids && push!(cells, c) + end + @assert !isempty(cells) + + # Calculate barycentric coordinates for a point x in the given cell + function bary(cell, x) + local ids = trian.cell_node_ids[cell] + @assert length(ids) == D+1 "simplex has correct number of vertices" + local ys = zero(SVector{D+1,SVector{D,T}}) + for (n,id) in enumerate(ids) + local y = SVector(Tuple(node_coordinates[id])) + ys = setindex(ys, y, n) + end + local λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) + return λ + end + + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # negative for exactly one cell and positive for all other ones. + # Due to round-off, the distance can be slightly negative or + # slightly positive for points near cell boundaries, in particular + # near vertices. In this case, choose the cell with the smallest + # distance, and check that the distance (if positive) is at most + # at round-off level. + mindist = T(Inf) + minc = 0 + minλ = zero(SVector{D+1,T}) + for c in cells + λ = bary(c, x) + dist = max(-minimum(λ), maximum(λ) - 1) + if dist < mindist + minc = c + minλ = λ + mindist = dist + end + end + @assert 1 ≤ minc + cell = minc + λ = minλ + dist = mindist + + # Ensure the point is inside one of the cells, up to round-off errors + @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" + + # TODO: Actually evaluate basis functions + # if f isa FESpaces.SingleFieldFEFunction + # values = f.cell_dof_values[cell] + # elseif f isa CellData.GenericCellField + # f.cell_field::AbstractArray{<:Field} + # values = f.cell_field[cell] + # else + # @assert false "unkonwn CellField type" + # end + values = get_cell_dof_values(f)[cell] + fx = sum(λ[n] * values[n] for n in 1:D+1) + + elseif all(e == HEX_AXIS for e in extrusion) + + # Find cell + # x - + + @assert false "n-cube not implemented" + + else + + @assert false "Unsupported polytope" + + end + + return fx +end + +""" + λ = cartesian2barycentric(s, p) + +Convert Cartesian to barycentric coordinates. + +# Arguments +- `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` + vertices in `D` dimensions. +- `p`: Point in Cartesian coordinates + +# Result +- `λ`: Point in barycentric coordinates +""" +function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} + @assert N ≤ D + 1 + # Algorithm as described on + # , + # section "Conversion between barycentric and Cartesian coordinates" + A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) + b = SVector{D+1}(p..., T(1)) + return A \ b +end + +function cartesian2barycentric(s::SVector{N,SVector{D,T}}, + p::SVector{D,T}) where {D,N,T} + return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), + p) +end + +function barycentric2cartesian(s::SVector{N,SVector{D,T}}, + λ::SVector{N,T}) where {D,N,T} + @assert N ≤ D + 1 + return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) +end diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index f5a35ec1d..cf79a273e 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -4,10 +4,12 @@ using Test using BlockArrays using FillArrays using Gridap.Arrays -using Gridap.Fields using Gridap.TensorValues +using Gridap.Fields +using Gridap.ReferenceFEs using Gridap.Geometry using Gridap.CellData +using Gridap.FESpaces domain = (0,1,0,1) cells = (3,3) @@ -145,8 +147,6 @@ test_array(hx_N,collect(hx_N)) - - #np = 3 #ndofs = 4 # diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index d7f21b07e..472e82d97 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -11,6 +11,8 @@ using Gridap.Integration using Gridap.CellData using Gridap.MultiField using Gridap.TensorValues +using Random +using StaticArrays using Test using Gridap.Arrays: BlockArrayCooMap @@ -128,6 +130,39 @@ cellmat1_Γ = integrate(((n⋅dv.⁺)-dq.⁻)*((n⋅du.⁺)+dp.⁻),quad_Γ) cellmat2_Γ = integrate((n⋅dv.⁺)*(n⋅du.⁺)+(n⋅dv.⁺)*dp.⁻-dq.⁻*(n⋅du.⁺)-dq.⁻*dp.⁻,quad_Γ) test_array(cellmat1_Γ,cellmat2_Γ,≈) +# Test function evaluation + +# Set reproducible random number seed +Random.seed!(0) +@testset "evaluating functions" for D in 1:5 + xmin = 0 + xmax = 1 + domain = repeat([xmin, xmax], D) + ncells = 3 + partition = repeat([ncells], D) + model = CartesianDiscreteModel(domain, partition) + # TODO: test both with and without this + model = simplexify(model) + + order = 2 + reffe = ReferenceFE(lagrangian, Float64, order) + V = FESpace(model, reffe) + + coeff0 = rand(Float64) + coeffs = rand(SVector{D,Float64}) + f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 + fh = interpolate_everywhere(f, V) + fhcache = return_cache(fh, VectorValue(zeros(D)...)) + + for i in 1:10 + x = VectorValue(rand(D)...) + fx = f(x) + fhx = evaluate!(fhcache, fh, x) + @test fhx ≈ fx + end +end + + #a = cellmat_Γ #using BenchmarkTools From d1548f1b2c8dfe7c5677465350311736ec28bd89 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 28 Jan 2021 12:05:54 -0500 Subject: [PATCH 14/46] Comment code --- src/MultiField/MultiFieldCellFields.jl | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 13c9c7bed..b07863a99 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -131,25 +131,16 @@ function evaluate!(cache,f::CellField,x::Point) dist = mindist # Ensure the point is inside one of the cells, up to round-off errors + # TODO: Add proper run-time check @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" # TODO: Actually evaluate basis functions - # if f isa FESpaces.SingleFieldFEFunction - # values = f.cell_dof_values[cell] - # elseif f isa CellData.GenericCellField - # f.cell_field::AbstractArray{<:Field} - # values = f.cell_field[cell] - # else - # @assert false "unkonwn CellField type" - # end values = get_cell_dof_values(f)[cell] fx = sum(λ[n] * values[n] for n in 1:D+1) elseif all(e == HEX_AXIS for e in extrusion) - # Find cell - # x - - + # TOOD @assert false "n-cube not implemented" else From 9bd3a22eabd13014fce0328afbfac4f1eb00c9a7 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 28 Feb 2021 18:20:40 -0500 Subject: [PATCH 15/46] Add various comments with ideas or questions --- src/MultiField/MultiFieldCellFields.jl | 27 +++++++++++++++++++ .../MultiFieldCellFieldsTests.jl | 7 +++++ 2 files changed, 34 insertions(+) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index b07863a99..4ae41b525 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -49,6 +49,8 @@ mutable struct CellFieldCache kdtree::Union{Nothing, KDTree{SVector{D,T},Euclidean,T} where {D,T}} end +# TODO: Should this return the cache already in `f` instead of +# creating a new one? return_cache(f::CellField,x::Point) = CellFieldCache(nothing) function evaluate!(cache,f::CellField,x::Point) @@ -59,6 +61,7 @@ function evaluate!(cache,f::CellField,x::Point) node_coordinates = get_node_coordinates(trian) @assert !isempty(node_coordinates) sample_coord = node_coordinates[1] + # TODO: call Gridap.ReferenceFEs.num_cell_dims instead? D = length(sample_coord) T = eltype(sample_coord) @assert D > 0 @@ -83,6 +86,12 @@ function evaluate!(cache,f::CellField,x::Point) id,dist = nn(kdtree, SVector(Tuple(x))) cells = Int[] # TODO: Use better algorithm + # + # Francesc Verdugo @fverdugo Dec 13 2020 12:46 + # vertex_to_cells = get_faces(get_grid_topology(model),0,num_cell_dims(model)) + # This will provide cells around each vertex (vertex == node for + # linear meshes) + # # Francesc Verdugo @fverdugo 14:07 # You have to start with a DiscreteModel, then # topo=get_grid_topology(model) @@ -135,9 +144,27 @@ function evaluate!(cache,f::CellField,x::Point) @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" # TODO: Actually evaluate basis functions + # Francesc Verdugo @fverdugo Dec 09 2020 13:06 + # for 2., take a look at the functions in this file https://github.com/gridap/Gridap.jl/blob/master/src/Geometry/FaceLabelings.jl + # In particular this one: + # https://github.com/gridap/Gridap.jl/blob/4342137dfbe635867f39e64ecb01e8dbbc5fa6e7/src/Geometry/FaceLabelings.jl#L289 values = get_cell_dof_values(f)[cell] fx = sum(λ[n] * values[n] for n in 1:D+1) + # you can just evaluate the (finite element) solution at this + # point at this cell. + + # TODO: If you define the finite element space in the physical + # space, i.e., DomainStyle being PhysicalDomain(), that is all. + # However, if DomainStyle equal to ReferenceDomain() one needs the + # inverse geometrical map that takes the physical cell and maps it + # to the reference one (in which we evaluate shape functions) + + # TODO: are these lines useful? + # v = get_cell_shapefuns(V) + # u = get_cell_shapefuns_trial(U) + # cellmat = get_array( ∫( ∇(u)⋅∇(v) )dΩ ) + elseif all(e == HEX_AXIS for e in extrusion) # TOOD diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index 472e82d97..e54f3157e 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -151,6 +151,13 @@ Random.seed!(0) coeff0 = rand(Float64) coeffs = rand(SVector{D,Float64}) f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 + # TODO: use this mechanism instead to project + # Francesc Verdugo @fverdugo 13:11 + # a(u,v) = ∫( u*v )dΩ + # l(v) = a(f,v) + # Solve a fe problem with this weak form + # See also tutorial 10, "Isotropic damage model", section "L2 + # projection", function "project" fh = interpolate_everywhere(f, V) fhcache = return_cache(fh, VectorValue(zeros(D)...)) From c6ebff7b0a43cee4374ee393f8241fb13cdddd62 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 28 Feb 2021 18:20:54 -0500 Subject: [PATCH 16/46] Add progress output while running test case --- test/MultiFieldTests/MultiFieldCellFieldsTests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index e54f3157e..3448cf777 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -135,6 +135,7 @@ test_array(cellmat1_Γ,cellmat2_Γ,≈) # Set reproducible random number seed Random.seed!(0) @testset "evaluating functions" for D in 1:5 + @show "evaluating functions" D xmin = 0 xmax = 1 domain = repeat([xmin, xmax], D) From af56037db6678002a1940e9c9f4df0162ff7632d Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 28 Feb 2021 18:21:09 -0500 Subject: [PATCH 17/46] Shorten code --- src/MultiField/MultiFieldCellFields.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 4ae41b525..5ec67ffc3 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -123,20 +123,17 @@ function evaluate!(cache,f::CellField,x::Point) # distance, and check that the distance (if positive) is at most # at round-off level. mindist = T(Inf) - minc = 0 - minλ = zero(SVector{D+1,T}) + minc,minλ = 0, zero(SVector{D+1,T}) for c in cells λ = bary(c, x) dist = max(-minimum(λ), maximum(λ) - 1) if dist < mindist - minc = c - minλ = λ + minc,minλ = c,λ mindist = dist end end @assert 1 ≤ minc - cell = minc - λ = minλ + cell,λ = minc,minλ dist = mindist # Ensure the point is inside one of the cells, up to round-off errors From 6b081cf0bc0b9ec2818340940c9b804e37dc1e3f Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 3 Mar 2021 20:05:09 -0500 Subject: [PATCH 18/46] Rewrite code, making use of existing Gridap features --- src/MultiField/MultiFieldCellFields.jl | 178 ++++++++----------------- 1 file changed, 59 insertions(+), 119 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 5ec67ffc3..dfa23ed53 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -45,134 +45,74 @@ Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state) # This code lives here (instead of in the CellData module) because we # need access to MultiField functionality. -mutable struct CellFieldCache - kdtree::Union{Nothing, KDTree{SVector{D,T},Euclidean,T} where {D,T}} +function return_cache(f::CellField,x::Point) + trian = get_triangulation(f) + node_coordinates = get_node_coordinates(trian) + node_coordinates = reshape(node_coordinates, length(node_coordinates)) + kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) + cell_f = get_array(f) + c1 = array_cache(cell_f) + f = testitem(cell_f) + @show typeof(f) typeof(x) typeof(cell_f) + c2 = return_cache(f,x) + return kdtree, c1, c2, cell_f, f end -# TODO: Should this return the cache already in `f` instead of -# creating a new one? -return_cache(f::CellField,x::Point) = CellFieldCache(nothing) +# function evaluate!(cache,f::CellField,x::Point) +function evaluate!(cache,f::SingleFieldFEFunction,x::Point) + kdtree, c1, c2, cell_g, g = cache + if f === g + cell_f = cell_g + else + cell_f = get_array(f) + end -function evaluate!(cache,f::CellField,x::Point) - cache::CellFieldCache + # Find nearest vertex + id,dist = nn(kdtree, SVector(Tuple(x))) - # Examine triangulation + # Find all neighbouring cells trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - @assert !isempty(node_coordinates) - sample_coord = node_coordinates[1] - # TODO: call Gridap.ReferenceFEs.num_cell_dims instead? - D = length(sample_coord) - T = eltype(sample_coord) - @assert D > 0 - - # Find polytope - reffes = trian.reffes - @assert length(reffes) == 1 - reffe = reffes[1].reffe - polytope = reffe.polytope - extrusion = polytope.extrusion - @assert length(extrusion) == D - - if all(e == TET_AXIS for e in extrusion) - - # Get k-d tree for mesh from cache - if cache.kdtree === nothing - cache.kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) + # topo = get_grid_topology(trian) + topo = GridTopology(trian) + vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + cells = vertex_to_cells[id] + + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # negative for exactly one cell and positive for all other ones. Due + # to round-off, the distance can be slightly negative or slightly + # positive for points near cell boundaries, in particular near + # vertices. In this case, choose the cell with the smallest + # distance, and check that the distance (if positive) is at most at + # round-off level. + mindist = T(Inf) + minc,minλ = 0,zero(SVector{D+1,T}) + for c in cells + # Calculate barycentric coordinates λ for x in cell c + ids = trian.cell_node_ids[c] + @assert length(ids) == D+1 "simplex has correct number of vertices" + ys = zero(SVector{D+1,SVector{D,T}}) + for (n,id) in enumerate(ids) + y = SVector(Tuple(node_coordinates[id])) + ys = setindex(ys, y, n) end - kdtree = cache.kdtree::KDTree{SVector{D,T},Euclidean,T} - - # Find all neighbouring cells - id,dist = nn(kdtree, SVector(Tuple(x))) - cells = Int[] - # TODO: Use better algorithm - # - # Francesc Verdugo @fverdugo Dec 13 2020 12:46 - # vertex_to_cells = get_faces(get_grid_topology(model),0,num_cell_dims(model)) - # This will provide cells around each vertex (vertex == node for - # linear meshes) - # - # Francesc Verdugo @fverdugo 14:07 - # You have to start with a DiscreteModel, then - # topo=get_grid_topology(model) - # vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) - for (c,ids) in enumerate(trian.cell_node_ids) - id ∈ ids && push!(cells, c) + λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) + dist = max(-minimum(λ), maximum(λ) - 1) + if dist < mindist + minc,minλ = c,λ + mindist = dist end - @assert !isempty(cells) - - # Calculate barycentric coordinates for a point x in the given cell - function bary(cell, x) - local ids = trian.cell_node_ids[cell] - @assert length(ids) == D+1 "simplex has correct number of vertices" - local ys = zero(SVector{D+1,SVector{D,T}}) - for (n,id) in enumerate(ids) - local y = SVector(Tuple(node_coordinates[id])) - ys = setindex(ys, y, n) - end - local λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) - return λ - end - - # Calculate the distance from the point to all the cells. Without - # round-off, and with non-overlapping cells, the distance would be - # negative for exactly one cell and positive for all other ones. - # Due to round-off, the distance can be slightly negative or - # slightly positive for points near cell boundaries, in particular - # near vertices. In this case, choose the cell with the smallest - # distance, and check that the distance (if positive) is at most - # at round-off level. - mindist = T(Inf) - minc,minλ = 0, zero(SVector{D+1,T}) - for c in cells - λ = bary(c, x) - dist = max(-minimum(λ), maximum(λ) - 1) - if dist < mindist - minc,minλ = c,λ - mindist = dist - end - end - @assert 1 ≤ minc - cell,λ = minc,minλ - dist = mindist - - # Ensure the point is inside one of the cells, up to round-off errors - # TODO: Add proper run-time check - @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" - - # TODO: Actually evaluate basis functions - # Francesc Verdugo @fverdugo Dec 09 2020 13:06 - # for 2., take a look at the functions in this file https://github.com/gridap/Gridap.jl/blob/master/src/Geometry/FaceLabelings.jl - # In particular this one: - # https://github.com/gridap/Gridap.jl/blob/4342137dfbe635867f39e64ecb01e8dbbc5fa6e7/src/Geometry/FaceLabelings.jl#L289 - values = get_cell_dof_values(f)[cell] - fx = sum(λ[n] * values[n] for n in 1:D+1) - - # you can just evaluate the (finite element) solution at this - # point at this cell. - - # TODO: If you define the finite element space in the physical - # space, i.e., DomainStyle being PhysicalDomain(), that is all. - # However, if DomainStyle equal to ReferenceDomain() one needs the - # inverse geometrical map that takes the physical cell and maps it - # to the reference one (in which we evaluate shape functions) - - # TODO: are these lines useful? - # v = get_cell_shapefuns(V) - # u = get_cell_shapefuns_trial(U) - # cellmat = get_array( ∫( ∇(u)⋅∇(v) )dΩ ) - - elseif all(e == HEX_AXIS for e in extrusion) - - # TOOD - @assert false "n-cube not implemented" - - else - - @assert false "Unsupported polytope" - end + @assert 1 ≤ minc + cell,λ = minc,minλ + dist = mindist + + # Ensure the point is inside one of the cells, up to round-off errors + # TODO: Add proper run-time check + @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" + f = getindex!(c1,cell_f,cell) + fx = evaluate!(c2,f,x) return fx end From d6187d1edf6583048c5eb2a60dcbe6c717fcfd0a Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 3 Mar 2021 23:04:18 -0500 Subject: [PATCH 19/46] Commit current state --- src/MultiField/MultiField.jl | 4 +- src/MultiField/MultiFieldCellFields.jl | 141 ++++++++++++++++++------- 2 files changed, 107 insertions(+), 38 deletions(-) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 67f5ad0c4..73915629c 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -18,7 +18,9 @@ using Gridap.Fields using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.Arrays: BlockArrayCooMap -using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, get_node_coordinates +using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, + get_faces, get_grid_topology, get_node_coordinates, + num_cell_dims using FillArrays using SparseArrays diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index dfa23ed53..c81d05297 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -50,16 +50,15 @@ function return_cache(f::CellField,x::Point) node_coordinates = get_node_coordinates(trian) node_coordinates = reshape(node_coordinates, length(node_coordinates)) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) + # TODO: This does not work for SingleFieldFEFunction cell_f = get_array(f) c1 = array_cache(cell_f) f = testitem(cell_f) - @show typeof(f) typeof(x) typeof(cell_f) c2 = return_cache(f,x) return kdtree, c1, c2, cell_f, f end -# function evaluate!(cache,f::CellField,x::Point) -function evaluate!(cache,f::SingleFieldFEFunction,x::Point) +function evaluate!(cache,f::CellField,x::Point) kdtree, c1, c2, cell_g, g = cache if f === g cell_f = cell_g @@ -69,47 +68,115 @@ function evaluate!(cache,f::SingleFieldFEFunction,x::Point) # Find nearest vertex id,dist = nn(kdtree, SVector(Tuple(x))) + @show id + T = typeof(dist) # Find all neighbouring cells trian = get_triangulation(f) - # topo = get_grid_topology(trian) - topo = GridTopology(trian) - vertex_to_cells = get_faces(topo,0,num_cell_dims(model)) + @show typeof(trian) + node_coordinates = get_node_coordinates(trian) + node_coordinates = reshape(node_coordinates, length(node_coordinates)) + model = trian # this is wrong + # topo = get_grid_topology(model) + topo = GridTopology(model) + D = num_cell_dims(trian) + vertex_to_cells = get_faces(topo,0,D) cells = vertex_to_cells[id] - - # Calculate the distance from the point to all the cells. Without - # round-off, and with non-overlapping cells, the distance would be - # negative for exactly one cell and positive for all other ones. Due - # to round-off, the distance can be slightly negative or slightly - # positive for points near cell boundaries, in particular near - # vertices. In this case, choose the cell with the smallest - # distance, and check that the distance (if positive) is at most at - # round-off level. - mindist = T(Inf) - minc,minλ = 0,zero(SVector{D+1,T}) - for c in cells - # Calculate barycentric coordinates λ for x in cell c - ids = trian.cell_node_ids[c] - @assert length(ids) == D+1 "simplex has correct number of vertices" - ys = zero(SVector{D+1,SVector{D,T}}) - for (n,id) in enumerate(ids) - y = SVector(Tuple(node_coordinates[id])) - ys = setindex(ys, y, n) + @show cells + + # Find polytope + # TODO: this only works for simplices, not for hypercubes + # type CartesianGrid has no field reffes + @show typeof(trian) + reffes = trian.reffes + @assert length(reffes) == 1 + reffe = reffes[1].reffe + polytope = reffe.polytope + extrusion = polytope.extrusion + @assert length(extrusion) == D + + if all(e == TET_AXIS for e in extrusion) + + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # positive for exactly one cell and positive for all other ones. + # Due to round-off, the distance can be slightly negative or + # slightly positive for points near cell boundaries, in particular + # near vertices. In this case, choose the cell with the largest + # distance, and check that the distance (if negative) is at most + # at round-off level. + maxdist = -T(Inf) + maxcell = 0 + for cell in cells + @show "for" cell + # Calculate barycentric coordinates λ for x in cell c + ids = trian.cell_node_ids[cell] + @show "for" ids + N = length(ids) + # @show N D+1 + @assert N == D+1 "simplex has correct number of vertices" + ys = zero(SVector{N,SVector{D,T}}) + for (n,id) in enumerate(ids) + y = SVector(Tuple(node_coordinates[id])) + ys = setindex(ys, y, n) + end + @show "for" ys + λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) + @show "for" λ + # dist is the distance to the closest face. If positive, the point + # is inside this cell. + dist = minimum(λ) + @show "for" dist + if dist > maxdist + maxcell = cell + maxdist = dist + end end - λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) - dist = max(-minimum(λ), maximum(λ) - 1) - if dist < mindist - minc,minλ = c,λ - mindist = dist + @assert maxcell > 0 + cell = maxcell + dist = maxdist + + # Ensure the point is inside one of the cells, up to round-off errors + # TODO: Add proper run-time check + @show cell dist + @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" + + elseif all(e == HEX_AXIS for e in extrusion) + + maxdist = -T(Inf) + maxcell = 0 + for cell in cells + ids = trian.cell_node_ids[cell] + N = length(ids) + @assert N == 2^D "hypercube has correct number of vertices" + dist = T(Inf) + for d in 1:D + for f in 1:2 + face = [i+1 for i in 0:2^D-1 if (i & (1<<(D-d)) == 0) == (f==1)] + yc = ids[face[1]][d] + @assert all(ids[j][d] == y for j in face) + dist = min(dist, f == 1 ? x[d] - y : y - x[d]) + end + end + if dist > maxdist + maxcell = cell + maxdist = dist + end end - end - @assert 1 ≤ minc - cell,λ = minc,minλ - dist = mindist + @assert maxcell > 0 + cell = maxcell + dist = maxdist + + # Ensure the point is inside one of the cells, up to round-off errors + # TODO: Add proper run-time check + @show cell dist + @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" - # Ensure the point is inside one of the cells, up to round-off errors - # TODO: Add proper run-time check - @assert dist ≤ 1000 * eps(T) "Point is not inside a cell" + else + + @notimplemented "Need either tet or hex mesh" + + end f = getindex!(c1,cell_f,cell) fx = evaluate!(c2,f,x) From 7544e8cedee042a80b5dfa0740993d946fe71993 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 14:06:02 -0500 Subject: [PATCH 20/46] Things seem to work now. --- src/MultiField/MultiField.jl | 4 +- src/MultiField/MultiFieldCellFields.jl | 271 ++++++++++++------------- 2 files changed, 133 insertions(+), 142 deletions(-) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 73915629c..47e5165d8 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -18,8 +18,8 @@ using Gridap.Fields using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.Arrays: BlockArrayCooMap -using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, - get_faces, get_grid_topology, get_node_coordinates, +using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, ExtrusionPolytope, + get_faces, get_node_coordinates, get_polytope, num_cell_dims using FillArrays diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index c81d05297..b6e73d374 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -45,12 +45,42 @@ Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state) # This code lives here (instead of in the CellData module) because we # need access to MultiField functionality. +""" + dist = distance(polytope::ExtrusionPolytope, + inv_cmap::Field, + x::Point) + +Calculate distance from point `x` to the polytope. The polytope is +given by its type and by the inverse cell map, i.e. by the map from +the physical to the reference space. + +Positive distances are outside the polytope, negative distances are +inside the polytope. + +The distance is measured in an unspecified norm, currently the L∞ +norm. +""" +function distance(polytope::ExtrusionPolytope, inv_cmap::Field, x::Point) + extrusion = polytope.extrusion + isempty(extrusion) && return zero(eltype(x)) + p = inv_cmap(x) + if all(e == HEX_AXIS for e in extrusion) + # Boundaries are at `a=0` and `a=1` in each direction + return maximum(max(0 - a, a - 1) for a in p) + elseif all(e == TET_AXIS for e in extrusion) + # Calculate barycentric coordinates + λ = Point(p..., 1 - sum(p)) + return maximum(-λ) + else + @notimplemented "Only hypercubes and simplices are implemented so far" + end +end + function return_cache(f::CellField,x::Point) trian = get_triangulation(f) node_coordinates = get_node_coordinates(trian) node_coordinates = reshape(node_coordinates, length(node_coordinates)) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) - # TODO: This does not work for SingleFieldFEFunction cell_f = get_array(f) c1 = array_cache(cell_f) f = testitem(cell_f) @@ -66,154 +96,115 @@ function evaluate!(cache,f::CellField,x::Point) cell_f = get_array(f) end + # Left-over comments from discussing with @fverdugo: + + # TODO: check contained-ness in a helper function; keep this + # function here independent of the topology type + # e.g. is_inside(reffe::LagrangianRefFE,xs,x) + # xs: physical nodes; x: point to check + # + # maybe better this instead: + # pass cell map instead, convert to reference space, check there + # e.g. is_inside(polytope,cell_map,x) + # or pass inverse cell map instead? maybe it is already computed? + # e.g. is_inside(polytope,invcell_map,x) + # + # cell_map = get_cell_map(trian) + # cell_map[cell] + # inverse_map(cell_map[cell]) + # cell_to_ctype = get_cell_type(trian) + # ctype = cell_to_ctype[cell] + # poly = ctype_to_poly[ctype] + # + # Find polytope + # TODO: this only works for simplices, not for hypercubes + # type CartesianGrid has no field reffes + # + # TODO: ctype_to_reffe = get_reffes(trian) + # TODO: ctype_to_poly = map(get_polytope, ctype_to_reffe) + # TODO: @assert length(ctype_to_poly) == 1 + # TODO: poly = first(ctype_to_poly) + # + # End of left-over comments. + # Find nearest vertex id,dist = nn(kdtree, SVector(Tuple(x))) - @show id - T = typeof(dist) # Find all neighbouring cells trian = get_triangulation(f) - @show typeof(trian) - node_coordinates = get_node_coordinates(trian) - node_coordinates = reshape(node_coordinates, length(node_coordinates)) - model = trian # this is wrong - # topo = get_grid_topology(model) - topo = GridTopology(model) + topo = GridTopology(trian) # TODO: this is expensive D = num_cell_dims(trian) vertex_to_cells = get_faces(topo,0,D) cells = vertex_to_cells[id] - @show cells - - # Find polytope - # TODO: this only works for simplices, not for hypercubes - # type CartesianGrid has no field reffes - @show typeof(trian) - reffes = trian.reffes - @assert length(reffes) == 1 - reffe = reffes[1].reffe - polytope = reffe.polytope - extrusion = polytope.extrusion - @assert length(extrusion) == D - - if all(e == TET_AXIS for e in extrusion) - - # Calculate the distance from the point to all the cells. Without - # round-off, and with non-overlapping cells, the distance would be - # positive for exactly one cell and positive for all other ones. - # Due to round-off, the distance can be slightly negative or - # slightly positive for points near cell boundaries, in particular - # near vertices. In this case, choose the cell with the largest - # distance, and check that the distance (if negative) is at most - # at round-off level. - maxdist = -T(Inf) - maxcell = 0 - for cell in cells - @show "for" cell - # Calculate barycentric coordinates λ for x in cell c - ids = trian.cell_node_ids[cell] - @show "for" ids - N = length(ids) - # @show N D+1 - @assert N == D+1 "simplex has correct number of vertices" - ys = zero(SVector{N,SVector{D,T}}) - for (n,id) in enumerate(ids) - y = SVector(Tuple(node_coordinates[id])) - ys = setindex(ys, y, n) - end - @show "for" ys - λ = cartesian2barycentric(ys, SVector{D,T}(Tuple(x))) - @show "for" λ - # dist is the distance to the closest face. If positive, the point - # is inside this cell. - dist = minimum(λ) - @show "for" dist - if dist > maxdist - maxcell = cell - maxdist = dist - end - end - @assert maxcell > 0 - cell = maxcell - dist = maxdist - - # Ensure the point is inside one of the cells, up to round-off errors - # TODO: Add proper run-time check - @show cell dist - @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" - - elseif all(e == HEX_AXIS for e in extrusion) - - maxdist = -T(Inf) - maxcell = 0 - for cell in cells - ids = trian.cell_node_ids[cell] - N = length(ids) - @assert N == 2^D "hypercube has correct number of vertices" - dist = T(Inf) - for d in 1:D - for f in 1:2 - face = [i+1 for i in 0:2^D-1 if (i & (1<<(D-d)) == 0) == (f==1)] - yc = ids[face[1]][d] - @assert all(ids[j][d] == y for j in face) - dist = min(dist, f == 1 ? x[d] - y : y - x[d]) - end - end - if dist > maxdist - maxcell = cell - maxdist = dist - end - end - @assert maxcell > 0 - cell = maxcell - dist = maxdist - - # Ensure the point is inside one of the cells, up to round-off errors - # TODO: Add proper run-time check - @show cell dist - @assert dist ≥ -1000 * eps(T) "Point is not inside a cell" - - else - - @notimplemented "Need either tet or hex mesh" - + @assert !isempty(cells) + + cell_to_ctype = get_cell_type(trian) + ctype_to_reffe = get_reffes(trian) + ctype_to_polytope = map(get_polytope, ctype_to_reffe) + + cell_map = get_cell_map(trian) + + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # negative for exactly one cell and positive for all other ones. Due + # to round-off, the distance can be slightly negative or slightly + # positive for points near cell boundaries, in particular near + # vertices. In this case, choose the cell with the smallest + # distance, and check that the distance (if positive) is at most at + # round-off level. + T = eltype(dist) + cell_distances = similar(cells, T) + for (i, cell) in enumerate(cells) + ctype = cell_to_ctype[cell] + polytope = ctype_to_polytope[ctype] + cmap = cell_map[cell] + inv_cmap = inverse_map(cmap) + cell_distances[i] = distance(polytope, inv_cmap, x) end - - f = getindex!(c1,cell_f,cell) - fx = evaluate!(c2,f,x) + # Ensure at most one cell contains the point + @assert count(<(1000 * eps(T)), cell_distances) ≤ 1 + i = argmin(cell_distances) + dist = cell_distances[i] + # Ensure the point is inside one of the cells, up to round-off errors + @assert dist ≤ -1000 * eps(T) "Point is not inside any cell" + cell = cells[i] + + f = getindex!(c1, cell_f, cell) + fx = evaluate!(c2, f, x) return fx end -""" - λ = cartesian2barycentric(s, p) - -Convert Cartesian to barycentric coordinates. - -# Arguments -- `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` - vertices in `D` dimensions. -- `p`: Point in Cartesian coordinates - -# Result -- `λ`: Point in barycentric coordinates -""" -function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} - @assert N ≤ D + 1 - # Algorithm as described on - # , - # section "Conversion between barycentric and Cartesian coordinates" - A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) - b = SVector{D+1}(p..., T(1)) - return A \ b -end - -function cartesian2barycentric(s::SVector{N,SVector{D,T}}, - p::SVector{D,T}) where {D,N,T} - return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), - p) -end - -function barycentric2cartesian(s::SVector{N,SVector{D,T}}, - λ::SVector{N,T}) where {D,N,T} - @assert N ≤ D + 1 - return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) -end +# """ +# λ = cartesian2barycentric(s, p) +# +# Convert Cartesian to barycentric coordinates. +# +# # Arguments +# - `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` +# vertices in `D` dimensions. +# - `p`: Point in Cartesian coordinates +# +# # Result +# - `λ`: Point in barycentric coordinates +# """ +# function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} +# @assert N ≤ D + 1 +# # Algorithm as described on +# # , +# # section "Conversion between barycentric and Cartesian coordinates" +# A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) +# b = SVector{D+1}(p..., T(1)) +# return A \ b +# end +# +# function cartesian2barycentric(s::SVector{N,SVector{D,T}}, +# p::SVector{D,T}) where {D,N,T} +# return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), +# p) +# end +# +# function barycentric2cartesian(s::SVector{N,SVector{D,T}}, +# λ::SVector{N,T}) where {D,N,T} +# @assert N ≤ D + 1 +# return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) +# end From 0d765f134c162dfd1c3a52ebffdf744f38f39507 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 14:21:38 -0500 Subject: [PATCH 21/46] Remove left-over comments --- src/MultiField/MultiFieldCellFields.jl | 66 ------------------- .../MultiFieldCellFieldsTests.jl | 1 - 2 files changed, 67 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index b6e73d374..ddbf685e7 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -96,37 +96,6 @@ function evaluate!(cache,f::CellField,x::Point) cell_f = get_array(f) end - # Left-over comments from discussing with @fverdugo: - - # TODO: check contained-ness in a helper function; keep this - # function here independent of the topology type - # e.g. is_inside(reffe::LagrangianRefFE,xs,x) - # xs: physical nodes; x: point to check - # - # maybe better this instead: - # pass cell map instead, convert to reference space, check there - # e.g. is_inside(polytope,cell_map,x) - # or pass inverse cell map instead? maybe it is already computed? - # e.g. is_inside(polytope,invcell_map,x) - # - # cell_map = get_cell_map(trian) - # cell_map[cell] - # inverse_map(cell_map[cell]) - # cell_to_ctype = get_cell_type(trian) - # ctype = cell_to_ctype[cell] - # poly = ctype_to_poly[ctype] - # - # Find polytope - # TODO: this only works for simplices, not for hypercubes - # type CartesianGrid has no field reffes - # - # TODO: ctype_to_reffe = get_reffes(trian) - # TODO: ctype_to_poly = map(get_polytope, ctype_to_reffe) - # TODO: @assert length(ctype_to_poly) == 1 - # TODO: poly = first(ctype_to_poly) - # - # End of left-over comments. - # Find nearest vertex id,dist = nn(kdtree, SVector(Tuple(x))) @@ -173,38 +142,3 @@ function evaluate!(cache,f::CellField,x::Point) fx = evaluate!(c2, f, x) return fx end - -# """ -# λ = cartesian2barycentric(s, p) -# -# Convert Cartesian to barycentric coordinates. -# -# # Arguments -# - `s`: Simplex vertices in Cartesian coordinates. `s` has `N ≤ D + 1` -# vertices in `D` dimensions. -# - `p`: Point in Cartesian coordinates -# -# # Result -# - `λ`: Point in barycentric coordinates -# """ -# function cartesian2barycentric(s::SMatrix{N,D,T}, p::SVector{D,T}) where {N,D,T} -# @assert N ≤ D + 1 -# # Algorithm as described on -# # , -# # section "Conversion between barycentric and Cartesian coordinates" -# A = SMatrix{D+1,N}(i == D+1 ? T(1) : s[j, i] for i in 1:D+1, j in 1:N) -# b = SVector{D+1}(p..., T(1)) -# return A \ b -# end -# -# function cartesian2barycentric(s::SVector{N,SVector{D,T}}, -# p::SVector{D,T}) where {D,N,T} -# return cartesian2barycentric(SMatrix{N,D,T}(s[i][j] for i in 1:N, j in 1:D), -# p) -# end -# -# function barycentric2cartesian(s::SVector{N,SVector{D,T}}, -# λ::SVector{N,T}) where {D,N,T} -# @assert N ≤ D + 1 -# return SVector{D,T}(sum(s[i][j] * λ[i] for i in 1:N) for j in 1:D) -# end diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index 3448cf777..e54f3157e 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -135,7 +135,6 @@ test_array(cellmat1_Γ,cellmat2_Γ,≈) # Set reproducible random number seed Random.seed!(0) @testset "evaluating functions" for D in 1:5 - @show "evaluating functions" D xmin = 0 xmax = 1 domain = repeat([xmin, xmax], D) From 06ec40331515b51a6d368bb18c4758a392fc01d4 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 14:58:21 -0500 Subject: [PATCH 22/46] Correct error bounds in self check --- src/MultiField/MultiFieldCellFields.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index ddbf685e7..71c1e1850 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -131,11 +131,11 @@ function evaluate!(cache,f::CellField,x::Point) cell_distances[i] = distance(polytope, inv_cmap, x) end # Ensure at most one cell contains the point - @assert count(<(1000 * eps(T)), cell_distances) ≤ 1 + @assert count(<(-1000 * eps(T)), cell_distances) ≤ 1 i = argmin(cell_distances) dist = cell_distances[i] # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≤ -1000 * eps(T) "Point is not inside any cell" + @assert dist ≤ 1000 * eps(T) "Point is not inside any cell" cell = cells[i] f = getindex!(c1, cell_f, cell) From 73672a3d108c395c08f3cf276f74ee4167d91ba4 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 16:35:22 -0500 Subject: [PATCH 23/46] Avoid temporary array --- src/MultiField/MultiFieldCellFields.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 71c1e1850..e8d00c66d 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -122,21 +122,25 @@ function evaluate!(cache,f::CellField,x::Point) # distance, and check that the distance (if positive) is at most at # round-off level. T = eltype(dist) - cell_distances = similar(cells, T) - for (i, cell) in enumerate(cells) + function cell_distance(cell::Int32) ctype = cell_to_ctype[cell] polytope = ctype_to_polytope[ctype] cmap = cell_map[cell] inv_cmap = inverse_map(cmap) - cell_distances[i] = distance(polytope, inv_cmap, x) + return distance(polytope, inv_cmap, x) + end + # findmin, without allocating an array + cell = Int32(0) + dist = T(Inf) + for jcell in cells + jdist = cell_distance(jcell) + if jdist < dist + cell = jcell + dist = jdist + end end - # Ensure at most one cell contains the point - @assert count(<(-1000 * eps(T)), cell_distances) ≤ 1 - i = argmin(cell_distances) - dist = cell_distances[i] # Ensure the point is inside one of the cells, up to round-off errors @assert dist ≤ 1000 * eps(T) "Point is not inside any cell" - cell = cells[i] f = getindex!(c1, cell_f, cell) fx = evaluate!(c2, f, x) From 43afc113d8aad1571a5852cae852ccc50efbd864 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 16:57:40 -0500 Subject: [PATCH 24/46] Allow more general cell index types --- src/MultiField/MultiFieldCellFields.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index e8d00c66d..af96d8f20 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -122,7 +122,7 @@ function evaluate!(cache,f::CellField,x::Point) # distance, and check that the distance (if positive) is at most at # round-off level. T = eltype(dist) - function cell_distance(cell::Int32) + function cell_distance(cell::Integer) ctype = cell_to_ctype[cell] polytope = ctype_to_polytope[ctype] cmap = cell_map[cell] @@ -130,7 +130,7 @@ function evaluate!(cache,f::CellField,x::Point) return distance(polytope, inv_cmap, x) end # findmin, without allocating an array - cell = Int32(0) + cell = zero(eltype(cells)) dist = T(Inf) for jcell in cells jdist = cell_distance(jcell) From f8998e3100da7e256cb6dd3a25c40dca514d8c79 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 16:57:51 -0500 Subject: [PATCH 25/46] Correct indentation --- src/MultiField/MultiFieldCellFields.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index af96d8f20..68b865ff4 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -133,11 +133,11 @@ function evaluate!(cache,f::CellField,x::Point) cell = zero(eltype(cells)) dist = T(Inf) for jcell in cells - jdist = cell_distance(jcell) - if jdist < dist - cell = jcell - dist = jdist - end + jdist = cell_distance(jcell) + if jdist < dist + cell = jcell + dist = jdist + end end # Ensure the point is inside one of the cells, up to round-off errors @assert dist ≤ 1000 * eps(T) "Point is not inside any cell" From 2f2e82f298e80072c1e8f3cddb50c2badde60e16 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 5 Mar 2021 23:43:25 -0500 Subject: [PATCH 26/46] Use vertex coordinates, not node coordinates to locate cells --- src/MultiField/MultiFieldCellFields.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 68b865ff4..eaacb9468 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -78,9 +78,9 @@ end function return_cache(f::CellField,x::Point) trian = get_triangulation(f) - node_coordinates = get_node_coordinates(trian) - node_coordinates = reshape(node_coordinates, length(node_coordinates)) - kdtree = KDTree(map(nc -> SVector(Tuple(nc)), node_coordinates)) + topo = GridTopology(trian) # TODO: this is expensive + vertex_coordinates = Geometry.get_vertex_coordinates(topo) + kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) cell_f = get_array(f) c1 = array_cache(cell_f) f = testitem(cell_f) From 8d11d5cd9fd18a3159a47e34c9b205ffde2c0526 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 26 Apr 2021 11:07:56 -0400 Subject: [PATCH 27/46] Commit all changes --- src/CellData/CellData.jl | 1 - src/CellData/CellFields.jl | 2 +- src/FESpaces/ConformingFESpaces.jl | 2 +- src/Fields/FieldArrays.jl | 34 ++++++++++++++++++++++++++ src/Fields/Fields.jl | 1 + src/Geometry/SkeletonTriangulations.jl | 2 +- src/Geometry/Triangulations.jl | 2 +- 7 files changed, 39 insertions(+), 5 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index c2365d933..271e5c45d 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -14,7 +14,6 @@ using Gridap.Algebra using Gridap.Arrays using Gridap.TensorValues using Gridap.Fields -using Gridap.Integration using Gridap.ReferenceFEs using Gridap.Geometry diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index cd220703b..d7392b410 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -519,7 +519,7 @@ function Base.getproperty(x::CellField, sym::Symbol) end end -function Base.propertynames(x::CellField, private=false) +function Base.propertynames(x::CellField, private::Bool=false) (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus) end diff --git a/src/FESpaces/ConformingFESpaces.jl b/src/FESpaces/ConformingFESpaces.jl index c4e1714b7..1d6d53814 100644 --- a/src/FESpaces/ConformingFESpaces.jl +++ b/src/FESpaces/ConformingFESpaces.jl @@ -21,7 +21,7 @@ function Base.getproperty(a::CellConformity, sym::Symbol) end end -function Base.propertynames(x::CellConformity, private=false) +function Base.propertynames(x::CellConformity, private::Bool=false) (fieldnames(typeof(x))...,:d_ctype_offset,:d_ctype_ldface_own_ldofs) end diff --git a/src/Fields/FieldArrays.jl b/src/Fields/FieldArrays.jl index 2f4771b64..ed7584578 100644 --- a/src/Fields/FieldArrays.jl +++ b/src/Fields/FieldArrays.jl @@ -664,3 +664,37 @@ for op in (:*,:⋅,:⊙,:⊗) end end end + +# Inverse fields + +struct InverseField{F} <: Field + original::F +end + +inverse_map(a::Field) = InverseField(a) +inverse_map(a::InverseField) = a.original + +function return_cache(a::InverseField,x::Point) + return return_cache(a.original,x), return_cache(∇(a.original),x) +end + +function evaluate!(caches,a::InverseField,x::Point) + P = typeof(x) + D = length(x) + cache,∇cache = caches + # Initial guess. (Can this be improved?) + y₀ = zero([x...]) + # Function and its derivative + f!(F,y) = F .= SVector{D}(Tuple(evaluate!(cache,a.original,P(y...)) - x)) + j!(J,y) = J .= SMatrix{D,D}(Tuple(evaluate!(∇cache,∇(a.original),P(y...)))) + # Solve + res = nlsolve(f!,j!,y₀) + @assert converged(res) + # Extract solution + y = res.zero + return P(y...) +end + +function evaluate!(cache,a::InverseField,xs::AbstractVector{<:Point}) + map(x->evaluate!(cache,a,x), xs) +end diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 95aa7aead..81d0afc0b 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -20,6 +20,7 @@ using LinearAlgebra: mul!, Transpose using ForwardDiff using FillArrays +using NLsolve using Test using BlockArrays using StaticArrays diff --git a/src/Geometry/SkeletonTriangulations.jl b/src/Geometry/SkeletonTriangulations.jl index df0c71296..36250650e 100644 --- a/src/Geometry/SkeletonTriangulations.jl +++ b/src/Geometry/SkeletonTriangulations.jl @@ -26,7 +26,7 @@ function Base.getproperty(x::SkeletonTriangulation, sym::Symbol) end end -function Base.propertynames(x::SkeletonTriangulation, private=false) +function Base.propertynames(x::SkeletonTriangulation, private::Bool=false) (fieldnames(typeof(x))...,:⁺,:⁻) end diff --git a/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl index 821ff5075..2c449ac55 100644 --- a/src/Geometry/Triangulations.jl +++ b/src/Geometry/Triangulations.jl @@ -204,7 +204,7 @@ function Base.getproperty(x::SkeletonPair, sym::Symbol) end end -function Base.propertynames(x::SkeletonPair, private=false) +function Base.propertynames(x::SkeletonPair, private::Bool=false) (fieldnames(typeof(x))...,:⁺,:⁻) end From 68736738b9f014d9ea540eef3d153d0c27726208 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 28 Apr 2021 15:14:39 -0400 Subject: [PATCH 28/46] Add tests for inverse maps --- src/Fields/FieldArrays.jl | 34 ----------------- src/Fields/Fields.jl | 2 + src/Fields/InverseFields.jl | 39 +++++++++++++++++++ test/FieldsTests/InverseFieldsTests.jl | 53 ++++++++++++++++++++++++++ test/FieldsTests/runtests.jl | 2 + 5 files changed, 96 insertions(+), 34 deletions(-) create mode 100644 src/Fields/InverseFields.jl create mode 100644 test/FieldsTests/InverseFieldsTests.jl diff --git a/src/Fields/FieldArrays.jl b/src/Fields/FieldArrays.jl index ed7584578..2f4771b64 100644 --- a/src/Fields/FieldArrays.jl +++ b/src/Fields/FieldArrays.jl @@ -664,37 +664,3 @@ for op in (:*,:⋅,:⊙,:⊗) end end end - -# Inverse fields - -struct InverseField{F} <: Field - original::F -end - -inverse_map(a::Field) = InverseField(a) -inverse_map(a::InverseField) = a.original - -function return_cache(a::InverseField,x::Point) - return return_cache(a.original,x), return_cache(∇(a.original),x) -end - -function evaluate!(caches,a::InverseField,x::Point) - P = typeof(x) - D = length(x) - cache,∇cache = caches - # Initial guess. (Can this be improved?) - y₀ = zero([x...]) - # Function and its derivative - f!(F,y) = F .= SVector{D}(Tuple(evaluate!(cache,a.original,P(y...)) - x)) - j!(J,y) = J .= SMatrix{D,D}(Tuple(evaluate!(∇cache,∇(a.original),P(y...)))) - # Solve - res = nlsolve(f!,j!,y₀) - @assert converged(res) - # Extract solution - y = res.zero - return P(y...) -end - -function evaluate!(cache,a::InverseField,xs::AbstractVector{<:Point}) - map(x->evaluate!(cache,a,x), xs) -end diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 81d0afc0b..789d5ad57 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -100,4 +100,6 @@ include("AutoDiff.jl") include("BlockFieldArrays.jl") +include("InverseFields.jl") + end diff --git a/src/Fields/InverseFields.jl b/src/Fields/InverseFields.jl new file mode 100644 index 000000000..8cd3de12d --- /dev/null +++ b/src/Fields/InverseFields.jl @@ -0,0 +1,39 @@ + +# Inverse fields + +struct InverseField{F} <: Field + original::F +end + +inverse_map(a::Field) = InverseField(a) +inverse_map(a::InverseField) = a.original + +function return_cache(a::InverseField,x::Point) + return return_cache(a.original,x), return_cache(∇(a.original),x) +end + +function evaluate!(caches,a::InverseField,x::Point) + P = typeof(x) + D = length(x) + cache,∇cache = caches + # Initial guess. (Can this be improved?) + # We cannot use an SVector here; nlsolve expects a mutable object. + # MVector also does not work (nlsolve then encounters nans). + y₀ = [zero(P)...] + # Function and its derivative + f!(F,y) = F .= SVector{D}(Tuple(evaluate!(cache,a.original,P(y)) - x)) + j!(J,y) = J .= SMatrix{D,D}(Tuple(evaluate!(∇cache,∇(a.original),P(y)))) + # Solve + res = nlsolve(f!,j!,y₀) + @assert converged(res) + # Extract solution + y = res.zero + return P(y) +end + +function return_cache(a::InverseField,xs::AbstractVector{<:Point}) + return return_cache(a,testitem(xs)) +end +function evaluate!(cache,a::InverseField,xs::AbstractVector{<:Point}) + map(x->evaluate!(cache,a,x), xs) +end diff --git a/test/FieldsTests/InverseFieldsTests.jl b/test/FieldsTests/InverseFieldsTests.jl new file mode 100644 index 000000000..8cf7bb5d2 --- /dev/null +++ b/test/FieldsTests/InverseFieldsTests.jl @@ -0,0 +1,53 @@ +module InverseFieldsTests + +using Gridap.Fields +using Gridap.Arrays +using Gridap.TensorValues + +using Test + +b0 = Point(0,0) +m0 = TensorValue(1,0,0,1) +id = AffineMap(m0,b0) + +b1 = Point(1,1) +m1 = TensorValue(2,0,0,2) +h1 = AffineMap(m1,b1) + +b2 = Point(3,3) +m2 = TensorValue(4,0,0,4) +h2 = AffineMap(m2,b2) + +h = h2 ∘ h1 + +# (4x+3) ∘ (2x+1) = 8x+7 +b3 = Point(7,7) +m3 = TensorValue(8,0,0,8) +h3 = AffineMap(m3,b3) + +h1inv = inverse_map(h1) +h2inv = inverse_map(h2) +hinv = inverse_map(h) +hinvinv = inverse_map(hinv) + +function fields_are_equal(f::Field, g::Field) + local xs = Point{2,Float64}[Point(0,0), Point(1,0), Point(0,1)] + # all(f(x) ≈ g(x) for x in xs) + all(f(xs) ≈ g(xs)) +end + +@test fields_are_equal(id ∘ id, id) + +@test fields_are_equal(h, h3) + +@test fields_are_equal(id ∘ h, h) +@test fields_are_equal(h ∘ id, h) + +@test fields_are_equal(hinv ∘ h, id) +@test fields_are_equal(h ∘ hinv, id) + +@test fields_are_equal(hinvinv, h) + +@test fields_are_equal(hinv, h1inv ∘ h2inv) + +end # module diff --git a/test/FieldsTests/runtests.jl b/test/FieldsTests/runtests.jl index 4cae4f66d..a9a0c6d9f 100644 --- a/test/FieldsTests/runtests.jl +++ b/test/FieldsTests/runtests.jl @@ -18,4 +18,6 @@ using Test @testset "BlockFieldArrays" begin include("BlockFieldArraysTests.jl") end +@testset "InverseFields" begin include("InverseFieldsTests.jl") end + end From 7197cecc23da49b76f79dec0d80291934dfa44cb Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 29 Apr 2021 15:50:27 -0400 Subject: [PATCH 29/46] Correct module imports --- src/MultiField/MultiField.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 00eab8da4..966eeb5fc 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -17,7 +17,6 @@ using Gridap.Fields using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis -using Gridap.Arrays: BlockArrayCooMap using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, ExtrusionPolytope, get_faces, get_node_coordinates, get_polytope, num_cell_dims From 785876af9c8cf57c1b01a90a78a2e7f7d33a134f Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 29 Apr 2021 20:10:50 -0400 Subject: [PATCH 30/46] Make evaluate! for CellField more efficient --- src/MultiField/MultiFieldCellFields.jl | 43 ++++++++++++-------------- 1 file changed, 19 insertions(+), 24 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index d7e1bdbb0..64365cf9b 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -78,41 +78,36 @@ end function return_cache(f::CellField,x::Point) trian = get_triangulation(f) - topo = GridTopology(trian) # TODO: this is expensive + topo = GridTopology(trian) # Note: this is expensive vertex_coordinates = Geometry.get_vertex_coordinates(topo) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) + + D = num_cell_dims(trian) + vertex_to_cells = get_faces(topo,0,D) + cell_to_ctype = get_cell_type(trian) + ctype_to_reffe = get_reffes(trian) + ctype_to_polytope = map(get_polytope, ctype_to_reffe) + cell_map = get_cell_map(trian) + cell_f = get_array(f) - c1 = array_cache(cell_f) - f = testitem(cell_f) - c2 = return_cache(f,x) - return kdtree, c1, c2, cell_f, f + cell_f_cache = array_cache(cell_f) + cf = testitem(cell_f) + f_cache = return_cache(cf,x) + + return kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, cell_f_cache, f_cache, cell_f, f end function evaluate!(cache,f::CellField,x::Point) - kdtree, c1, c2, cell_g, g = cache - if f === g - cell_f = cell_g - else - cell_f = get_array(f) - end + kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, cell_f_cache, f_cache, cell_f, f₀ = cache + @check f === f₀ "Wrong cache" # Find nearest vertex id,dist = nn(kdtree, SVector(Tuple(x))) # Find all neighbouring cells - trian = get_triangulation(f) - topo = GridTopology(trian) # TODO: this is expensive - D = num_cell_dims(trian) - vertex_to_cells = get_faces(topo,0,D) cells = vertex_to_cells[id] @assert !isempty(cells) - cell_to_ctype = get_cell_type(trian) - ctype_to_reffe = get_reffes(trian) - ctype_to_polytope = map(get_polytope, ctype_to_reffe) - - cell_map = get_cell_map(trian) - # Calculate the distance from the point to all the cells. Without # round-off, and with non-overlapping cells, the distance would be # negative for exactly one cell and positive for all other ones. Due @@ -140,9 +135,9 @@ function evaluate!(cache,f::CellField,x::Point) end end # Ensure the point is inside one of the cells, up to round-off errors - @assert dist ≤ 1000 * eps(T) "Point is not inside any cell" + @check dist ≤ 1000eps(T) "Point is not inside any cell" - f = getindex!(c1, cell_f, cell) - fx = evaluate!(c2, f, x) + cf = getindex!(cell_f_cache, cell_f, cell) + fx = evaluate!(f_cache, cf, x) return fx end From b65f9f697ce682c7165ce127030eca8da2a85a85 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 29 Apr 2021 20:11:28 -0400 Subject: [PATCH 31/46] Test evaluating cell fields on arrays of points --- src/MultiField/MultiFieldCellFields.jl | 6 ++++++ test/MultiFieldTests/MultiFieldCellFieldsTests.jl | 5 ++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 64365cf9b..6b46c8e72 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -141,3 +141,9 @@ function evaluate!(cache,f::CellField,x::Point) fx = evaluate!(f_cache, cf, x) return fx end + +# Simple version: +function evaluate!(::Nothing,f::CellField,xs::AbstractVector{<:Point}) + cache = return_cache(f,testitem(xs)) + return map(x->evaluate!(cache,f,x), xs) +end diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index 5a78e1eb9..506fe774a 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -149,12 +149,15 @@ Random.seed!(0) fh = interpolate_everywhere(f, V) fhcache = return_cache(fh, VectorValue(zeros(D)...)) - for i in 1:10 + xs = [VectorValue(rand(D)...) for i in 1:10] + for x in xs x = VectorValue(rand(D)...) fx = f(x) fhx = evaluate!(fhcache, fh, x) @test fhx ≈ fx end + fhxs = fh(xs) + @test fhxs ≈ f.(xs) end From 57ac76a8b7c70f299826c0abb92e9742baa7bb68 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 30 Apr 2021 14:51:45 -0400 Subject: [PATCH 32/46] Remove duplicate include statement --- src/Fields/Fields.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 316bb4cb4..82c2f111a 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -105,6 +105,4 @@ include("ArrayBlocks.jl") include("InverseFields.jl") -include("InverseFields.jl") - end From 4fef3f18023d8a5800a2d8f018ce8cc3e09eae1a Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 30 Apr 2021 14:52:41 -0400 Subject: [PATCH 33/46] Implement efficient vectorized evaluation --- src/MultiField/MultiFieldCellFields.jl | 74 ++++++++++++++++++++++---- 1 file changed, 65 insertions(+), 9 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 6b46c8e72..633d21af9 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -81,25 +81,25 @@ function return_cache(f::CellField,x::Point) topo = GridTopology(trian) # Note: this is expensive vertex_coordinates = Geometry.get_vertex_coordinates(topo) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) - D = num_cell_dims(trian) vertex_to_cells = get_faces(topo,0,D) cell_to_ctype = get_cell_type(trian) ctype_to_reffe = get_reffes(trian) ctype_to_polytope = map(get_polytope, ctype_to_reffe) cell_map = get_cell_map(trian) + cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map cell_f = get_array(f) cell_f_cache = array_cache(cell_f) cf = testitem(cell_f) f_cache = return_cache(cf,x) + cache2 = cell_f_cache, f_cache, cell_f, f - return kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, cell_f_cache, f_cache, cell_f, f + return cache1,cache2 end -function evaluate!(cache,f::CellField,x::Point) - kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, cell_f_cache, f_cache, cell_f, f₀ = cache - @check f === f₀ "Wrong cache" +function point_to_cell!(cache,f::CellField,x::Point) + kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache # Find nearest vertex id,dist = nn(kdtree, SVector(Tuple(x))) @@ -137,13 +137,69 @@ function evaluate!(cache,f::CellField,x::Point) # Ensure the point is inside one of the cells, up to round-off errors @check dist ≤ 1000eps(T) "Point is not inside any cell" + return cell +end + +function evaluate!(cache,f::CellField,x::Point) + cache1,cache2 = cache + cell_f_cache, f_cache, cell_f, f₀ = cache2 + @check f === f₀ "Wrong cache" + + cell = point_to_cell!(cache1,f,x) cf = getindex!(cell_f_cache, cell_f, cell) fx = evaluate!(f_cache, cf, x) return fx end -# Simple version: -function evaluate!(::Nothing,f::CellField,xs::AbstractVector{<:Point}) - cache = return_cache(f,testitem(xs)) - return map(x->evaluate!(cache,f,x), xs) +return_cache(f::CellField,xs::AbstractVector{<:Point}) = return_cache(f,testitem(xs)) + +# # Simple version: +# function evaluate!(cache,f::CellField,xs::AbstractVector{<:Point}) +# return map(x->evaluate!(cache,f,x), xs) +# end + +function make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int) + ni = length(i2j) + @assert nj≥0 + + p = sortperm(i2j) + # i2j[p] is a sorted array of js + + data = p + ptrs = Array{Int}(undef, nj+1) + i2lis = Array{Int}(undef, ni) + jprev = zero(eltype(i2j)) + liprev = 0 + for (n,i) in enumerate(p) + j = i2j[i] + @assert jprev≤j≤nj + li = (j==jprev ? liprev : 0) + 1 + ptrs[jprev+1:j] .= n + i2lis[i] = li + jprev = j + liprev = li + end + ptrs[jprev+1:nj+1] .= ni+1 + j2is = Table(data,ptrs) + + return j2is,i2lis +end + +# Efficient version: +function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) + cache1,cache2 = return_cache(f,testitem(point_to_x)) + kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1 + cell_f_cache, f_cache, cell_f, f₀ = cache2 + @check f === f₀ "Wrong cache" + + ncells = length(cell_map) + x_to_cell(x) = point_to_cell!(cache1,f,x) + point_to_cell = map(x_to_cell,point_to_x) + cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells) + cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points) + cell_to_f = get_array(f) + cell_to_fxs = lazy_map(evaluate,cell_to_f,cell_to_xs) + point_to_fxs = lazy_map(Reindex(cell_to_fxs),point_to_cell) + point_to_fx = lazy_map(getindex,point_to_fxs,point_to_lpoint) + collect(point_to_fx) # Collect into a plain array end From 7548314b4c88e9a1112f494abfe50a4ac8ad854e Mon Sep 17 00:00:00 2001 From: Balaje K Date: Mon, 31 May 2021 17:55:47 +1000 Subject: [PATCH 34/46] Move code to src/CellData from src/MultiField --- src/CellData/CellData.jl | 3 + src/CellData/CellFields.jl | 164 +++++++++++++++++- src/MultiField/MultiFieldCellFields.jl | 163 ----------------- test/CellDataTests/CellFieldsTests.jl | 44 ++++- .../MultiFieldCellFieldsTests.jl | 43 ----- 5 files changed, 209 insertions(+), 208 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 271e5c45d..005d9c76b 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -17,6 +17,9 @@ using Gridap.Fields using Gridap.ReferenceFEs using Gridap.Geometry +using NearestNeighbors +using StaticArrays + import Gridap.Arrays: lazy_append import Gridap.Arrays: get_array import Gridap.Arrays: evaluate! diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 4966a8d02..f75885d91 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -1,4 +1,3 @@ - """ A single point or an array of points on the cells of a Triangulation CellField objects can be evaluated efficiently at CellPoint instances. @@ -210,6 +209,169 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() # Evaluation of CellFields +# This code lives here (instead of in the CellData module) because we +# need access to MultiField functionality. + +""" + dist = distance(polytope::ExtrusionPolytope, + inv_cmap::Field, + x::Point) + +Calculate distance from point `x` to the polytope. The polytope is +given by its type and by the inverse cell map, i.e. by the map from +the physical to the reference space. + +Positive distances are outside the polytope, negative distances are +inside the polytope. + +The distance is measured in an unspecified norm, currently the L∞ +norm. +""" +function distance(polytope::ExtrusionPolytope, inv_cmap::Field, x::Point) + extrusion = polytope.extrusion + isempty(extrusion) && return zero(eltype(x)) + p = inv_cmap(x) + if all(e == HEX_AXIS for e in extrusion) + # Boundaries are at `a=0` and `a=1` in each direction + return maximum(max(0 - a, a - 1) for a in p) + elseif all(e == TET_AXIS for e in extrusion) + # Calculate barycentric coordinates + λ = Point(p..., 1 - sum(p)) + return maximum(-λ) + else + @notimplemented "Only hypercubes and simplices are implemented so far" + end +end + +function return_cache(f::CellField,x::Point) + trian = get_triangulation(f) + topo = GridTopology(trian) # Note: this is expensive + vertex_coordinates = Geometry.get_vertex_coordinates(topo) + kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) + D = num_cell_dims(trian) + vertex_to_cells = get_faces(topo,0,D) + cell_to_ctype = get_cell_type(trian) + ctype_to_reffe = get_reffes(trian) + ctype_to_polytope = map(get_polytope, ctype_to_reffe) + cell_map = get_cell_map(trian) + cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map + + cell_f = get_array(f) + cell_f_cache = array_cache(cell_f) + cf = testitem(cell_f) + f_cache = return_cache(cf,x) + cache2 = cell_f_cache, f_cache, cell_f, f + + return cache1,cache2 +end + +function point_to_cell!(cache,f::CellField,x::Point) + kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache + + # Find nearest vertex + id,dist = nn(kdtree, SVector(Tuple(x))) + + # Find all neighbouring cells + cells = vertex_to_cells[id] + @assert !isempty(cells) + + # Calculate the distance from the point to all the cells. Without + # round-off, and with non-overlapping cells, the distance would be + # negative for exactly one cell and positive for all other ones. Due + # to round-off, the distance can be slightly negative or slightly + # positive for points near cell boundaries, in particular near + # vertices. In this case, choose the cell with the smallest + # distance, and check that the distance (if positive) is at most at + # round-off level. + T = eltype(dist) + function cell_distance(cell::Integer) + ctype = cell_to_ctype[cell] + polytope = ctype_to_polytope[ctype] + cmap = cell_map[cell] + inv_cmap = inverse_map(cmap) + return distance(polytope, inv_cmap, x) + end + # findmin, without allocating an array + cell = zero(eltype(cells)) + dist = T(Inf) + for jcell in cells + jdist = cell_distance(jcell) + if jdist < dist + cell = jcell + dist = jdist + end + end + # Ensure the point is inside one of the cells, up to round-off errors + @check dist ≤ 1000eps(T) "Point is not inside any cell" + + return cell +end + +function evaluate!(cache,f::CellField,x::Point) + cache1,cache2 = cache + cell_f_cache, f_cache, cell_f, f₀ = cache2 + @check f === f₀ "Wrong cache" + + cell = point_to_cell!(cache1,f,x) + cf = getindex!(cell_f_cache, cell_f, cell) + fx = evaluate!(f_cache, cf, x) + return fx +end + +return_cache(f::CellField,xs::AbstractVector{<:Point}) = return_cache(f,testitem(xs)) + +# # Simple version: +# function evaluate!(cache,f::CellField,xs::AbstractVector{<:Point}) +# return map(x->evaluate!(cache,f,x), xs) +# end + +function make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int) + ni = length(i2j) + @assert nj≥0 + + p = sortperm(i2j) + # i2j[p] is a sorted array of js + + data = p + ptrs = Array{Int}(undef, nj+1) + i2lis = Array{Int}(undef, ni) + jprev = zero(eltype(i2j)) + liprev = 0 + for (n,i) in enumerate(p) + j = i2j[i] + @assert jprev≤j≤nj + li = (j==jprev ? liprev : 0) + 1 + ptrs[jprev+1:j] .= n + i2lis[i] = li + jprev = j + liprev = li + end + ptrs[jprev+1:nj+1] .= ni+1 + j2is = Table(data,ptrs) + + return j2is,i2lis +end + +# Efficient version: +function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) + cache1,cache2 = return_cache(f,testitem(point_to_x)) + kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1 + cell_f_cache, f_cache, cell_f, f₀ = cache2 + @check f === f₀ "Wrong cache" + + ncells = length(cell_map) + x_to_cell(x) = point_to_cell!(cache1,f,x) + point_to_cell = map(x_to_cell,point_to_x) + cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells) + cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points) + cell_to_f = get_array(f) + cell_to_fxs = lazy_map(evaluate,cell_to_f,cell_to_xs) + point_to_fxs = lazy_map(Reindex(cell_to_fxs),point_to_cell) + point_to_fx = lazy_map(getindex,point_to_fxs,point_to_lpoint) + collect(point_to_fx) # Collect into a plain array +end + + (a::CellField)(x) = evaluate(a,x) function evaluate!(cache,f::CellField,x::CellPoint) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 633d21af9..bc8dbca04 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -1,4 +1,3 @@ - struct MultiFieldCellField{DS<:DomainStyle} <: CellField single_fields::Vector{<:CellField} trian::Triangulation @@ -41,165 +40,3 @@ Base.iterate(a::MultiFieldCellField) = iterate(a.single_fields) Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state) # Evaluation of CellFields - -# This code lives here (instead of in the CellData module) because we -# need access to MultiField functionality. - -""" - dist = distance(polytope::ExtrusionPolytope, - inv_cmap::Field, - x::Point) - -Calculate distance from point `x` to the polytope. The polytope is -given by its type and by the inverse cell map, i.e. by the map from -the physical to the reference space. - -Positive distances are outside the polytope, negative distances are -inside the polytope. - -The distance is measured in an unspecified norm, currently the L∞ -norm. -""" -function distance(polytope::ExtrusionPolytope, inv_cmap::Field, x::Point) - extrusion = polytope.extrusion - isempty(extrusion) && return zero(eltype(x)) - p = inv_cmap(x) - if all(e == HEX_AXIS for e in extrusion) - # Boundaries are at `a=0` and `a=1` in each direction - return maximum(max(0 - a, a - 1) for a in p) - elseif all(e == TET_AXIS for e in extrusion) - # Calculate barycentric coordinates - λ = Point(p..., 1 - sum(p)) - return maximum(-λ) - else - @notimplemented "Only hypercubes and simplices are implemented so far" - end -end - -function return_cache(f::CellField,x::Point) - trian = get_triangulation(f) - topo = GridTopology(trian) # Note: this is expensive - vertex_coordinates = Geometry.get_vertex_coordinates(topo) - kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) - D = num_cell_dims(trian) - vertex_to_cells = get_faces(topo,0,D) - cell_to_ctype = get_cell_type(trian) - ctype_to_reffe = get_reffes(trian) - ctype_to_polytope = map(get_polytope, ctype_to_reffe) - cell_map = get_cell_map(trian) - cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map - - cell_f = get_array(f) - cell_f_cache = array_cache(cell_f) - cf = testitem(cell_f) - f_cache = return_cache(cf,x) - cache2 = cell_f_cache, f_cache, cell_f, f - - return cache1,cache2 -end - -function point_to_cell!(cache,f::CellField,x::Point) - kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache - - # Find nearest vertex - id,dist = nn(kdtree, SVector(Tuple(x))) - - # Find all neighbouring cells - cells = vertex_to_cells[id] - @assert !isempty(cells) - - # Calculate the distance from the point to all the cells. Without - # round-off, and with non-overlapping cells, the distance would be - # negative for exactly one cell and positive for all other ones. Due - # to round-off, the distance can be slightly negative or slightly - # positive for points near cell boundaries, in particular near - # vertices. In this case, choose the cell with the smallest - # distance, and check that the distance (if positive) is at most at - # round-off level. - T = eltype(dist) - function cell_distance(cell::Integer) - ctype = cell_to_ctype[cell] - polytope = ctype_to_polytope[ctype] - cmap = cell_map[cell] - inv_cmap = inverse_map(cmap) - return distance(polytope, inv_cmap, x) - end - # findmin, without allocating an array - cell = zero(eltype(cells)) - dist = T(Inf) - for jcell in cells - jdist = cell_distance(jcell) - if jdist < dist - cell = jcell - dist = jdist - end - end - # Ensure the point is inside one of the cells, up to round-off errors - @check dist ≤ 1000eps(T) "Point is not inside any cell" - - return cell -end - -function evaluate!(cache,f::CellField,x::Point) - cache1,cache2 = cache - cell_f_cache, f_cache, cell_f, f₀ = cache2 - @check f === f₀ "Wrong cache" - - cell = point_to_cell!(cache1,f,x) - cf = getindex!(cell_f_cache, cell_f, cell) - fx = evaluate!(f_cache, cf, x) - return fx -end - -return_cache(f::CellField,xs::AbstractVector{<:Point}) = return_cache(f,testitem(xs)) - -# # Simple version: -# function evaluate!(cache,f::CellField,xs::AbstractVector{<:Point}) -# return map(x->evaluate!(cache,f,x), xs) -# end - -function make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int) - ni = length(i2j) - @assert nj≥0 - - p = sortperm(i2j) - # i2j[p] is a sorted array of js - - data = p - ptrs = Array{Int}(undef, nj+1) - i2lis = Array{Int}(undef, ni) - jprev = zero(eltype(i2j)) - liprev = 0 - for (n,i) in enumerate(p) - j = i2j[i] - @assert jprev≤j≤nj - li = (j==jprev ? liprev : 0) + 1 - ptrs[jprev+1:j] .= n - i2lis[i] = li - jprev = j - liprev = li - end - ptrs[jprev+1:nj+1] .= ni+1 - j2is = Table(data,ptrs) - - return j2is,i2lis -end - -# Efficient version: -function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) - cache1,cache2 = return_cache(f,testitem(point_to_x)) - kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1 - cell_f_cache, f_cache, cell_f, f₀ = cache2 - @check f === f₀ "Wrong cache" - - ncells = length(cell_map) - x_to_cell(x) = point_to_cell!(cache1,f,x) - point_to_cell = map(x_to_cell,point_to_x) - cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells) - cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points) - cell_to_f = get_array(f) - cell_to_fxs = lazy_map(evaluate,cell_to_f,cell_to_xs) - point_to_fxs = lazy_map(Reindex(cell_to_fxs),point_to_cell) - point_to_fx = lazy_map(getindex,point_to_fxs,point_to_lpoint) - collect(point_to_fx) # Collect into a plain array -end diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index f56c08f56..5cce58845 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -10,6 +10,8 @@ using Gridap.ReferenceFEs using Gridap.Geometry using Gridap.CellData using Gridap.FESpaces +using Random +using StaticArrays domain = (0,1,0,1) cells = (3,3) @@ -192,7 +194,47 @@ h_N = (2*f_N+g)⋅g hx_N = h_N(x_N) test_array(hx_N,collect(hx_N)) - +# Test function evaluation + +# Set reproducible random number seed +Random.seed!(0) +@testset "evaluating functions" for D in 1:3 + xmin = 0 + xmax = 1 + domain = repeat([xmin, xmax], D) + ncells = 3 + partition = repeat([ncells], D) + model = CartesianDiscreteModel(domain, partition) + # TODO: test both with and without this + model = simplexify(model) + + order = 2 + reffe = ReferenceFE(lagrangian, Float64, order) + V = FESpace(model, reffe) + + coeff0 = rand(Float64) + coeffs = rand(SVector{D,Float64}) + f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 + # TODO: use this mechanism instead to project + # Francesc Verdugo @fverdugo 13:11 + # a(u,v) = ∫( u*v )dΩ + # l(v) = a(f,v) + # Solve a fe problem with this weak form + # See also tutorial 10, "Isotropic damage model", section "L2 + # projection", function "project" + fh = interpolate_everywhere(f, V) + fhcache = return_cache(fh, VectorValue(zeros(D)...)) + + xs = [VectorValue(rand(D)...) for i in 1:10] + for x in xs + x = VectorValue(rand(D)...) + fx = f(x) + fhx = evaluate!(fhcache, fh, x) + @test fhx ≈ fx + end + fhxs = fh(xs) + @test fhxs ≈ f.(xs) +end #np = 3 #ndofs = 4 diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index 8ab29e672..e059e4057 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -118,49 +118,6 @@ cellmat1_Γ = integrate(((n⋅dv.⁺)-dq.⁻)*((n⋅du.⁺)+dp.⁻),quad_Γ) cellmat2_Γ = integrate((n⋅dv.⁺)*(n⋅du.⁺)+(n⋅dv.⁺)*dp.⁻-dq.⁻*(n⋅du.⁺)-dq.⁻*dp.⁻,quad_Γ) test_array(cellmat1_Γ,cellmat2_Γ,≈) -# Test function evaluation - -# Set reproducible random number seed -Random.seed!(0) -@testset "evaluating functions" for D in 1:5 - xmin = 0 - xmax = 1 - domain = repeat([xmin, xmax], D) - ncells = 3 - partition = repeat([ncells], D) - model = CartesianDiscreteModel(domain, partition) - # TODO: test both with and without this - model = simplexify(model) - - order = 2 - reffe = ReferenceFE(lagrangian, Float64, order) - V = FESpace(model, reffe) - - coeff0 = rand(Float64) - coeffs = rand(SVector{D,Float64}) - f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 - # TODO: use this mechanism instead to project - # Francesc Verdugo @fverdugo 13:11 - # a(u,v) = ∫( u*v )dΩ - # l(v) = a(f,v) - # Solve a fe problem with this weak form - # See also tutorial 10, "Isotropic damage model", section "L2 - # projection", function "project" - fh = interpolate_everywhere(f, V) - fhcache = return_cache(fh, VectorValue(zeros(D)...)) - - xs = [VectorValue(rand(D)...) for i in 1:10] - for x in xs - x = VectorValue(rand(D)...) - fx = f(x) - fhx = evaluate!(fhcache, fh, x) - @test fhx ≈ fx - end - fhxs = fh(xs) - @test fhxs ≈ f.(xs) -end - - #a = cellmat_Γ #using BenchmarkTools From 5823b9a615db0d837c3557be4b5f63b9df666e78 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Mon, 31 May 2021 19:34:59 +1000 Subject: [PATCH 35/46] Fixing comment in CellFields and remove repeating InverseFields Tests --- src/CellData/CellFields.jl | 3 --- test/FieldsTests/runtests.jl | 2 -- 2 files changed, 5 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index f75885d91..940c005e0 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -209,9 +209,6 @@ DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() # Evaluation of CellFields -# This code lives here (instead of in the CellData module) because we -# need access to MultiField functionality. - """ dist = distance(polytope::ExtrusionPolytope, inv_cmap::Field, diff --git a/test/FieldsTests/runtests.jl b/test/FieldsTests/runtests.jl index 2c87f69f2..51378f23d 100644 --- a/test/FieldsTests/runtests.jl +++ b/test/FieldsTests/runtests.jl @@ -20,6 +20,4 @@ using Test @testset "InverseFields" begin include("InverseFieldsTests.jl") end -@testset "InverseFields" begin include("InverseFieldsTests.jl") end - end From 0f6663f2d5b34a211fe97136e2e31abe590695e1 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 12:58:32 +1000 Subject: [PATCH 36/46] Add more tests for evaluate --- test/CellDataTests/CellFieldsTests.jl | 79 ++++++++++++++++++--------- 1 file changed, 52 insertions(+), 27 deletions(-) diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 5cce58845..29f014bae 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -204,36 +204,61 @@ Random.seed!(0) domain = repeat([xmin, xmax], D) ncells = 3 partition = repeat([ncells], D) - model = CartesianDiscreteModel(domain, partition) - # TODO: test both with and without this - model = simplexify(model) - + base_model = CartesianDiscreteModel(domain, partition) order = 2 reffe = ReferenceFE(lagrangian, Float64, order) - V = FESpace(model, reffe) - - coeff0 = rand(Float64) - coeffs = rand(SVector{D,Float64}) - f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 - # TODO: use this mechanism instead to project - # Francesc Verdugo @fverdugo 13:11 - # a(u,v) = ∫( u*v )dΩ - # l(v) = a(f,v) - # Solve a fe problem with this weak form - # See also tutorial 10, "Isotropic damage model", section "L2 - # projection", function "project" - fh = interpolate_everywhere(f, V) - fhcache = return_cache(fh, VectorValue(zeros(D)...)) - - xs = [VectorValue(rand(D)...) for i in 1:10] - for x in xs - x = VectorValue(rand(D)...) - fx = f(x) - fhx = evaluate!(fhcache, fh, x) - @test fhx ≈ fx + + for M ∈ ["Hypercubes", "Simplices"] + model = (M == "Hypercubes") ? base_model : simplexify(base_model) + + V = FESpace(model, reffe) + + coeff0 = rand(Float64) + coeffs = rand(SVector{D,Float64}) + f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0 + # TODO: use this mechanism instead to project + # Francesc Verdugo @fverdugo 13:11 + # a(u,v) = ∫( u*v )dΩ + # l(v) = a(f,v) + # Solve a fe problem with this weak form + # See also tutorial 10, "Isotropic damage model", section "L2 + # projection", function "project" + fh = interpolate_everywhere(f, V) + fhcache = return_cache(fh, VectorValue(zeros(D)...)) + + # Test Random points + xs = [VectorValue(rand(D)...) for i in 1:10] + for x in xs + x = VectorValue(rand(D)...) + fx = f(x) + fhx = evaluate!(fhcache, fh, x) + @test fhx ≈ fx + end + fhxs = fh(xs) + @test fhxs ≈ f.(xs) + + nv = num_vertices(model) # Number of vertices + nf = num_faces(model,D-1) # Number of faces + + topo = GridTopology(model) + + pts = get_vertex_coordinates(topo) # Vertex coordinates + face_nodes = get_face_nodes(model, D-1) # face-to-node numbering + face_coords = lazy_map(Broadcasting(Reindex(pts)), face_nodes) # Get LazyArray of coordinates of face + + # Test a random vertex from the triangulation + pt = pts[rand(1:nv)] + fhpt = evaluate!(fhcache, fh, pt) + @test fhpt .≈ f(pt) + + # Test a random point lying on the face of the polytope + face_coord = face_coords[rand(1:nf)] + λ = rand(length(face_coord)); + λ = (D > 1) ? λ./sum(λ) : λ + pt = face_coord ⋅ λ # Point on the face + fhpt = evaluate!(fhcache, fh, pt) + @test fhpt .≈ f(pt) end - fhxs = fh(xs) - @test fhxs ≈ f.(xs) end #np = 3 From 008a46a682704be14168d038fdf230d4dd271687 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 19:13:20 +1000 Subject: [PATCH 37/46] Remove unused deps in MultiField.jl, more test for evaluate --- src/CellData/CellData.jl | 4 ++++ src/CellData/CellFields.jl | 3 +-- src/MultiField/MultiField.jl | 11 +---------- test/CellDataTests/CellFieldsTests.jl | 17 ++++++++++++++++- 4 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 005d9c76b..93329acdf 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -61,6 +61,10 @@ export CellDof export get_normal_vector export get_cell_measure +export make_inverse_table +export point_to_cell! +export distance + export DomainContribution export num_domains export get_domains diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 940c005e0..ca70d3c29 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -207,7 +207,6 @@ get_data(f::GenericCellField) = f.cell_field get_triangulation(f::GenericCellField) = f.trian DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() -# Evaluation of CellFields """ dist = distance(polytope::ExtrusionPolytope, @@ -351,7 +350,7 @@ end # Efficient version: function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) - cache1,cache2 = return_cache(f,testitem(point_to_x)) + cache1,cache2 = cache kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1 cell_f_cache, f_cache, cell_f, f₀ = cache2 @check f === f₀ "Wrong cache" diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 966eeb5fc..cc7ae08fd 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -1,4 +1,4 @@ -""" +.""" The exported names are $(EXPORTS) @@ -17,19 +17,10 @@ using Gridap.Fields using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis -using Gridap.ReferenceFEs: HEX_AXIS, TET_AXIS, ExtrusionPolytope, - get_faces, get_node_coordinates, get_polytope, - num_cell_dims using FillArrays using SparseArrays using LinearAlgebra -using BlockArrays -using NearestNeighbors -using StaticArrays - -import Gridap.Arrays: evaluate! -import Gridap.Arrays: return_cache export num_fields export compute_field_offsets diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 29f014bae..0c1e08656 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -239,7 +239,7 @@ Random.seed!(0) nv = num_vertices(model) # Number of vertices nf = num_faces(model,D-1) # Number of faces - + trian = Triangulation(model) topo = GridTopology(model) pts = get_vertex_coordinates(topo) # Vertex coordinates @@ -258,6 +258,21 @@ Random.seed!(0) pt = face_coord ⋅ λ # Point on the face fhpt = evaluate!(fhcache, fh, pt) @test fhpt .≈ f(pt) + + # Test with CellPoint + cache1,cache2 = fhcache + ncells = num_cells(model) + x_to_cell(x) = point_to_cell!(cache1, fh, x) + point_to_cell = map(x_to_cell, xs) + cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells) + cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points) + cell_to_f = get_array(fh) + cell_to_fxs = lazy_map(evaluate, cell_to_f, cell_to_xs) + # Now build CellPoint with cell_to_xs map + cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain()) + cell_point_fxs = evaluate(fh, cell_point_xs) + @test cell_point_fxs ≈ cell_to_fxs + end end From e034de13998cf1aaa349d98fe9f1f3dfa7b430dc Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 19:22:53 +1000 Subject: [PATCH 38/46] Fix an error --- src/MultiField/MultiField.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index cc7ae08fd..f9b64173e 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -1,4 +1,4 @@ -.""" +""" The exported names are $(EXPORTS) From 54962d9f91207c42911c09e8de5bd8c11bbb9b27 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 19:44:50 +1000 Subject: [PATCH 39/46] Fix error in MultiField.jl --- src/MultiField/MultiField.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index f9b64173e..fd80bbf68 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -21,6 +21,7 @@ using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis using FillArrays using SparseArrays using LinearAlgebra +using BlockArrays export num_fields export compute_field_offsets From 5c445d1001f6a47b47278dfb3894101881db5fdc Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 20:37:59 +1000 Subject: [PATCH 40/46] Remove comment --- src/MultiField/MultiFieldCellFields.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index bc8dbca04..db41c07b3 100644 --- a/src/MultiField/MultiFieldCellFields.jl +++ b/src/MultiField/MultiFieldCellFields.jl @@ -38,5 +38,3 @@ num_fields(a::MultiFieldCellField) = length(a.single_fields) Base.getindex(a::MultiFieldCellField,i::Integer) = a.single_fields[i] Base.iterate(a::MultiFieldCellField) = iterate(a.single_fields) Base.iterate(a::MultiFieldCellField,state) = iterate(a.single_fields,state) - -# Evaluation of CellFields From 8ade674b15048ce6ab06ea60da44d9cbabd4e444 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 21:22:07 +1000 Subject: [PATCH 41/46] Change to not export distance --- src/CellData/CellData.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 93329acdf..529478fd5 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -63,7 +63,6 @@ export get_cell_measure export make_inverse_table export point_to_cell! -export distance export DomainContribution export num_domains From 0b4d2d957f27c92ea0f7e5258b26243b91595493 Mon Sep 17 00:00:00 2001 From: Balaje K Date: Tue, 1 Jun 2021 21:32:16 +1000 Subject: [PATCH 42/46] Remove unused argument from point_to_cell --- src/CellData/CellFields.jl | 6 +++--- test/CellDataTests/CellFieldsTests.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index ca70d3c29..1495bc24e 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -261,7 +261,7 @@ function return_cache(f::CellField,x::Point) return cache1,cache2 end -function point_to_cell!(cache,f::CellField,x::Point) +function point_to_cell!(cache, x::Point) kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache # Find nearest vertex @@ -308,7 +308,7 @@ function evaluate!(cache,f::CellField,x::Point) cell_f_cache, f_cache, cell_f, f₀ = cache2 @check f === f₀ "Wrong cache" - cell = point_to_cell!(cache1,f,x) + cell = point_to_cell!(cache1, x) cf = getindex!(cell_f_cache, cell_f, cell) fx = evaluate!(f_cache, cf, x) return fx @@ -356,7 +356,7 @@ function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) @check f === f₀ "Wrong cache" ncells = length(cell_map) - x_to_cell(x) = point_to_cell!(cache1,f,x) + x_to_cell(x) = point_to_cell!(cache1,x) point_to_cell = map(x_to_cell,point_to_x) cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells) cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points) diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 0c1e08656..e8320c0de 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -262,7 +262,7 @@ Random.seed!(0) # Test with CellPoint cache1,cache2 = fhcache ncells = num_cells(model) - x_to_cell(x) = point_to_cell!(cache1, fh, x) + x_to_cell(x) = point_to_cell!(cache1, x) point_to_cell = map(x_to_cell, xs) cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells) cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points) From 7a79a6eccf0e5a794b8d9e09c633da7e6c26a50f Mon Sep 17 00:00:00 2001 From: Balaje K Date: Wed, 2 Jun 2021 13:00:34 +1000 Subject: [PATCH 43/46] Add interface for AbstractVector{<:Point} --- src/CellData/CellFields.jl | 36 +++++++++++++++++++-------- test/CellDataTests/CellFieldsTests.jl | 6 +++-- 2 files changed, 30 insertions(+), 12 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 1495bc24e..a59937172 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -241,16 +241,7 @@ end function return_cache(f::CellField,x::Point) trian = get_triangulation(f) - topo = GridTopology(trian) # Note: this is expensive - vertex_coordinates = Geometry.get_vertex_coordinates(topo) - kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) - D = num_cell_dims(trian) - vertex_to_cells = get_faces(topo,0,D) - cell_to_ctype = get_cell_type(trian) - ctype_to_reffe = get_reffes(trian) - ctype_to_polytope = map(get_polytope, ctype_to_reffe) - cell_map = get_cell_map(trian) - cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map + cache1 = return_cache(trian) cell_f = get_array(f) cell_f_cache = array_cache(cell_f) @@ -261,6 +252,19 @@ function return_cache(f::CellField,x::Point) return cache1,cache2 end +function return_cache(trian::Triangulation) + topo = GridTopology(trian) + vertex_coordinates = Geometry.get_vertex_coordinates(topo) + kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) + D = num_cell_dims(trian) + vertex_to_cells = get_faces(topo, 0, D) + cell_to_ctype = get_cell_type(trian) + ctype_to_reffe = get_reffes(trian) + ctype_to_polytope = map(get_polytope, ctype_to_reffe) + cell_map = get_cell_map(trian) + cache1 = kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map +end + function point_to_cell!(cache, x::Point) kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache @@ -367,6 +371,18 @@ function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) collect(point_to_fx) # Collect into a plain array end +# New CellPoint implementation +function CellPoint(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain) + #Find the location of the point in the triangulation + cache1 = return_cache(trian) + x_to_cell(x) = point_to_cell!(cache1, x) + point_to_cell = map(x_to_cell, xs) + ncells = num_cells(trian) + cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells) + cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points) + #Create CellPoint + cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain()) +end (a::CellField)(x) = evaluate(a,x) diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index e8320c0de..5b172e987 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -260,6 +260,7 @@ Random.seed!(0) @test fhpt .≈ f(pt) # Test with CellPoint + # Build cell_to_fxs manually cache1,cache2 = fhcache ncells = num_cells(model) x_to_cell(x) = point_to_cell!(cache1, x) @@ -268,8 +269,9 @@ Random.seed!(0) cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points) cell_to_f = get_array(fh) cell_to_fxs = lazy_map(evaluate, cell_to_f, cell_to_xs) - # Now build CellPoint with cell_to_xs map - cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain()) + + # Now build CellPoint with xs instead of building cell_to_xs + cell_point_xs = CellPoint(xs, trian, PhysicalDomain()) cell_point_fxs = evaluate(fh, cell_point_xs) @test cell_point_fxs ≈ cell_to_fxs From 620ee7f5ec28facb31543b26bb08491712c2c20f Mon Sep 17 00:00:00 2001 From: Balaje K Date: Wed, 2 Jun 2021 16:19:52 +1000 Subject: [PATCH 44/46] Fix function names --- src/CellData/CellFields.jl | 6 +++--- test/CellDataTests/CellFieldsTests.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index a59937172..22f947bed 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -241,7 +241,7 @@ end function return_cache(f::CellField,x::Point) trian = get_triangulation(f) - cache1 = return_cache(trian) + cache1 = my_trian_data(trian) cell_f = get_array(f) cell_f_cache = array_cache(cell_f) @@ -252,7 +252,7 @@ function return_cache(f::CellField,x::Point) return cache1,cache2 end -function return_cache(trian::Triangulation) +function my_trian_data(trian::Triangulation) topo = GridTopology(trian) vertex_coordinates = Geometry.get_vertex_coordinates(topo) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) @@ -374,7 +374,7 @@ end # New CellPoint implementation function CellPoint(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain) #Find the location of the point in the triangulation - cache1 = return_cache(trian) + cache1 = my_trian_data(trian) x_to_cell(x) = point_to_cell!(cache1, x) point_to_cell = map(x_to_cell, xs) ncells = num_cells(trian) diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 5b172e987..3d2ee1c4c 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -208,8 +208,8 @@ Random.seed!(0) order = 2 reffe = ReferenceFE(lagrangian, Float64, order) - for M ∈ ["Hypercubes", "Simplices"] - model = (M == "Hypercubes") ? base_model : simplexify(base_model) + for M ∈ [:hypercubes, :simplices] + model = (M == :hypercubes) ? base_model : simplexify(base_model) V = FESpace(model, reffe) From 0af55ff6133f7f6bf35745c5d3db3bdb67fb216c Mon Sep 17 00:00:00 2001 From: Balaje K Date: Wed, 2 Jun 2021 18:00:22 +1000 Subject: [PATCH 45/46] Change from CellPoint to compute_cell_points_from_vector_of_points --- src/CellData/CellData.jl | 1 + src/CellData/CellFields.jl | 11 ++++------- test/CellDataTests/CellFieldsTests.jl | 2 +- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 529478fd5..c51e8d9ae 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -63,6 +63,7 @@ export get_cell_measure export make_inverse_table export point_to_cell! +export compute_cell_points_from_vector_of_points export DomainContribution export num_domains diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index 22f947bed..93a845855 100644 --- a/src/CellData/CellFields.jl +++ b/src/CellData/CellFields.jl @@ -241,7 +241,7 @@ end function return_cache(f::CellField,x::Point) trian = get_triangulation(f) - cache1 = my_trian_data(trian) + cache1 = _point_to_cell_cache(trian) cell_f = get_array(f) cell_f_cache = array_cache(cell_f) @@ -252,7 +252,7 @@ function return_cache(f::CellField,x::Point) return cache1,cache2 end -function my_trian_data(trian::Triangulation) +function _point_to_cell_cache(trian::Triangulation) topo = GridTopology(trian) vertex_coordinates = Geometry.get_vertex_coordinates(topo) kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates)) @@ -371,16 +371,13 @@ function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) collect(point_to_fx) # Collect into a plain array end -# New CellPoint implementation -function CellPoint(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain) - #Find the location of the point in the triangulation - cache1 = my_trian_data(trian) +function compute_cell_points_from_vector_of_points(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain) + cache1 = _point_to_cell_cache(trian) x_to_cell(x) = point_to_cell!(cache1, x) point_to_cell = map(x_to_cell, xs) ncells = num_cells(trian) cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells) cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points) - #Create CellPoint cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain()) end diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 3d2ee1c4c..14561521a 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -271,7 +271,7 @@ Random.seed!(0) cell_to_fxs = lazy_map(evaluate, cell_to_f, cell_to_xs) # Now build CellPoint with xs instead of building cell_to_xs - cell_point_xs = CellPoint(xs, trian, PhysicalDomain()) + cell_point_xs = compute_cell_points_from_vector_of_points(xs, trian, PhysicalDomain()) cell_point_fxs = evaluate(fh, cell_point_xs) @test cell_point_fxs ≈ cell_to_fxs From b7f83ba59032e53812309f52e35bb562fd83aa2a Mon Sep 17 00:00:00 2001 From: Balaje K Date: Wed, 2 Jun 2021 18:10:29 +1000 Subject: [PATCH 46/46] Update NEWS.md --- NEWS.md | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index fc80f6fd6..bc91be8b2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Experimental support for mixed dimensional PDEs. Since PR [#567](https://github.com/gridap/Gridap.jl/pull/567). - Added `get_cell_dof_basis(model,cell_reffes,::Conformity)` and `get_cell_shapefuns(model,cell_reffes,::Conformity)`. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579). - Implemented `get_cell_dof_basis` and `get_cell_shapefuns` for global RT FE spaces in a new file `DivConformingFESpaces.jl`. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579). +- Added support to allow evaluation of FE functions at arbitrary points. Since PR [#523](https://github.com/gridap/Gridap.jl/pull/523). +- Implemented `compute_cell_points_from_vector_of_points` to build `CellPoint` from a vector of points. Since PR [#523](https://github.com/gridap/Gridap.jl/pull/523). ### Changed - Major refactoring in the handling of blocks (e.g. in multi-field and skeleton terms). The new code follows a much more simple approach based in the new type `ArrayBlock`. Since PR [#583](https://github.com/gridap/Gridap.jl/pull/583). @@ -28,7 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - All code associated with with `BlockArrayCoo`. Since PR [#583](https://github.com/gridap/Gridap.jl/pull/583). - Module `Gridap.Integration` has been deleted and its contents have been merged into `Gridap.ReferenceFEs` module. - Types `SparseMatrixCSR` and `SymSparseMatrixCSR` have been moved to the registered package [`SparseMatricesCSR`](https://github.com/gridap/SparseMatricesCSR.jl). To use them simply add `SparseMatricesCSR` into your environment and type `using SparseMatricesCSR`. Since Since PR [#568](https://github.com/gridap/Gridap.jl/pull/568). -- Removed `PushForwardMap` and all code depending upon it. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579). +- Removed `PushForwardMap` and all code depending upon it. Since PR [#579](https://github.com/gridap/Gridap.jl/pull/579). ## [0.15.5] - 2021-05-25 @@ -43,7 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.15.3] - 2021-03-16 -### Added +### Added - `get_cell_map` now returns array of `AffineMap` for linear grids of simplices. Needed to compute Laplacian operator, inverse maps etc. Since PR [#553](https://github.com/gridap/Gridap.jl/pull/553). ### Fixed @@ -83,7 +85,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic - Bug-fix for 32-bit Julia: Replace all occurences of Int64 by Int. Since PR [#445](https://github.com/gridap/Gridap.jl/pull/445). - Bug-fix for 32-bit Julia. Using inttype=Int keyword argument for JSON parsing. Since PR [#456](https://github.com/gridap/Gridap.jl/pull/456). -## [0.14.1] - 2020-09-17 +## [0.14.1] - 2020-09-17 ### Added - Added VectorWithEntryInserted and VectorWithEntryRemoved. Since PR [#401](https://github.com/gridap/Gridap.jl/pull/401/). @@ -102,7 +104,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic - First and second argument switch in `update_state_variables!` in order to have function-first style. Since PR [#376](https://github.com/gridap/Gridap.jl/pull/376/). - Table struct has been generalized such that data and ptrs arrays can be of an arbitrary type extending AbstractArray. Since PR [#310](https://github.com/gridap/Gridap.jl/pull/310/) - `interpolate, interpolate!, interpolate_dirichlet...` switched argument order to function first style. For instance `interpolate(u, V)` instead of `interpolate(V, u)` - + ### Added - Allowing the construction of an `HomogeneousTrialFESpace` from a `TrialFESpace`. Since PR [#384](https://github.com/gridap/Gridap.jl/pull/384). - Support for automatic differentiation of residuals and Jacobians in multi-field computations since PR [#383](https://github.com/gridap/Gridap.jl/pull/383/). @@ -113,7 +115,7 @@ This version is a major (backwards-incompatible) refactoring of the project whic ## [0.13.4] - 2020-08-23 -### Added +### Added - New `FilteredCellArray` since PR [#372](https://github.com/gridap/Gridap.jl/pull/372/). ## [0.13.3] - 2020-08-12