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 diff --git a/Project.toml b/Project.toml index 57bd4ff35..a7fca05d8 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,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" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -40,6 +41,7 @@ JLD2 = "0.1.11, 0.3, 0.4" JSON = "0.21.0" LineSearches = "7.0.1" NLsolve = "4.3.0" +NearestNeighbors = "0.4.8" QuadGK = "2.3.1, 2.4" SparseMatricesCSR = "0.6" StaticArrays = "0.12.1, 1.0" diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index 28b85c190..c51e8d9ae 100644 --- a/src/CellData/CellData.jl +++ b/src/CellData/CellData.jl @@ -10,16 +10,20 @@ using DocStringExtensions using FillArrays using Gridap.Helpers -using Gridap.Arrays using Gridap.Algebra +using Gridap.Arrays using Gridap.TensorValues 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! +import Gridap.Arrays: return_cache import Gridap.Fields: gradient import Gridap.Fields: ∇∇ import Gridap.Fields: integrate @@ -57,6 +61,10 @@ export CellDof export get_normal_vector 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 export get_domains diff --git a/src/CellData/CellFields.jl b/src/CellData/CellFields.jl index c01331fbe..93a845855 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. @@ -208,20 +207,182 @@ get_data(f::GenericCellField) = f.cell_field get_triangulation(f::GenericCellField) = f.trian DomainStyle(::Type{GenericCellField{DS}}) where DS = DS() -# Evaluation of CellFields -(a::CellField)(x) = evaluate(a,x) +""" + 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) + cache1 = _point_to_cell_cache(trian) + + 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(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 + + # 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) - @notimplemented """\n - Evaluation of a CellField at a given Point is not implemented yet. + cache1,cache2 = cache + cell_f_cache, f_cache, cell_f, f₀ = cache2 + @check f === f₀ "Wrong cache" + + cell = point_to_cell!(cache1, 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) - 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. - """ + return j2is,i2lis end +# Efficient version: +function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point}) + 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" + + ncells = length(cell_map) + 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) + 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 + +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) + cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain()) +end + +(a::CellField)(x) = evaluate(a,x) + function evaluate!(cache,f::CellField,x::CellPoint) _f, _x = _to_common_domain(f,x) cell_field = get_data(_f) @@ -334,11 +495,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 caught error for more information. (If you are using Visual + 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). + command-line REPL instead). """ end end diff --git a/src/MultiField/MultiField.jl b/src/MultiField/MultiField.jl index 3ab239136..fd80bbf68 100644 --- a/src/MultiField/MultiField.jl +++ b/src/MultiField/MultiField.jl @@ -15,12 +15,13 @@ using Gridap.TensorValues using Gridap.CellData using Gridap.Fields +using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis using FillArrays using SparseArrays using LinearAlgebra -import BlockArrays +using BlockArrays export num_fields export compute_field_offsets diff --git a/src/MultiField/MultiFieldCellFields.jl b/src/MultiField/MultiFieldCellFields.jl index 5c704d6ff..db41c07b3 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 diff --git a/test/CellDataTests/CellFieldsTests.jl b/test/CellDataTests/CellFieldsTests.jl index 6bea038da..14561521a 100644 --- a/test/CellDataTests/CellFieldsTests.jl +++ b/test/CellDataTests/CellFieldsTests.jl @@ -4,10 +4,14 @@ 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 +using Random +using StaticArrays domain = (0,1,0,1) cells = (3,3) @@ -190,10 +194,89 @@ 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) + base_model = CartesianDiscreteModel(domain, partition) + order = 2 + reffe = ReferenceFE(lagrangian, Float64, order) + + 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 + trian = Triangulation(model) + 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) + + # Test with CellPoint + # Build cell_to_fxs manually + cache1,cache2 = fhcache + ncells = num_cells(model) + 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) + cell_to_f = get_array(fh) + 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 = 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 + + end +end #np = 3 #ndofs = 4 diff --git a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl index d74f6d301..e059e4057 100644 --- a/test/MultiFieldTests/MultiFieldCellFieldsTests.jl +++ b/test/MultiFieldTests/MultiFieldCellFieldsTests.jl @@ -9,6 +9,8 @@ using Gridap.Fields using Gridap.CellData using Gridap.MultiField using Gridap.TensorValues +using Random +using StaticArrays using Test domain = (0,1,0,1) @@ -116,6 +118,7 @@ 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_Γ,≈) + #a = cellmat_Γ #using BenchmarkTools #cache = array_cache(a)