Skip to content

Commit

Permalink
Many bugfixes
Browse files Browse the repository at this point in the history
  • Loading branch information
JordiManyer committed Jan 6, 2025
1 parent 186d26b commit f07fba1
Show file tree
Hide file tree
Showing 4 changed files with 15 additions and 12 deletions.
12 changes: 6 additions & 6 deletions src/FESpaces/Polytopal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -80,12 +80,12 @@ function centroid_map(poly::Polytope{D}) where D
v = get_vertex_coordinates(poly)
xc = mean(v)
h = maximum(vk -> norm(vk-xc), v)
origin = xc/h
gradient = TensorValues.diagonal_tensor(VectorValue([h for i in 1:D]))
origin = -xc/h
gradient = TensorValues.diagonal_tensor(VectorValue([1/h for i in 1:D]))
return affine_map(gradient,origin)
end

function Arrays.lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{typeof(centroid_map)}},x::AbstractVector{<:Point})
function Arrays.lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{typeof(centroid_map)}},x::AbstractVector)
polys = a.args[1]
lazy_map(CentroidCoordinateChangeMap(),polys,x)
end
Expand All @@ -97,12 +97,12 @@ function Arrays.evaluate!(cache,::CentroidCoordinateChangeMap,poly::Polytope{D},
return (x - xc) / h
end

function Arrays.return_cache(::CentroidCoordinateChangeMap,poly::Polytope{D},x::AbstractVector{Point{D}}) where D
function Arrays.return_cache(::CentroidCoordinateChangeMap,poly::Polytope{D},x::AbstractVector{<:Point{D}}) where D
return CachedArray(similar(x))
end

function Arrays.evaluate!(cache,::CentroidCoordinateChangeMap,poly::Polytope{D},x::AbstractVector{Point{D}}) where D
setsize!(cache,length(x))
function Arrays.evaluate!(cache,::CentroidCoordinateChangeMap,poly::Polytope{D},x::AbstractVector{<:Point{D}}) where D
setsize!(cache,size(x))
y = cache.array
v = get_vertex_coordinates(poly)
xc = mean(v)
Expand Down
2 changes: 1 addition & 1 deletion src/Fields/AffineMaps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,6 @@ end
# Constructor from a simplex given by D1+1 points
function affine_map(points::NTuple{D,Point{D2,T}}) where {D,D2,T}
origin, pk = first_and_tail(points)
gradient = TensorValues.tensor_from_columns(map(p -> p-origin, pk))
gradient = TensorValues.tensor_from_rows(map(p -> p-origin, pk))
AffineField(gradient,origin)
end
3 changes: 1 addition & 2 deletions src/ReferenceFEs/PolytopalQuadratures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,8 @@ end

function Arrays.evaluate!(cache,::PQuadCoordsMap, q::PolytopalQuadrature)
y_cache, cmap_cache = cache

conn = q.conn
coords = get_vertex_coordinates(q.poly)
conn = q.conn

x = get_coordinates(q.quad)
setsize!(y_cache,(length(x)*length(conn),))
Expand Down
10 changes: 7 additions & 3 deletions test/GeometryTests/PolytopalGridsTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,16 +14,20 @@ vmodel = Geometry.voronoi(simplexify(model))
polys = get_polytopes(vmodel)
# viz(vmodel;color=1:num_cells(vmodel),showpoints=true,showsegments=true)

order = 8
Ω = Triangulation(vmodel)
= Measure(Ω,2)
= Measure(Ω,2*order)
sum((1)dΩ)

V = FESpaces.PolytopalFESpace(Ω,1)
V = FESpaces.PolytopalFESpace(Ω,order)

u_exact(x) = x[1] + x[2]
u_exact(x) = x[1]^2 + x[2]^2
a(u,v) = (uv)dΩ
l(v) = (v*u_exact)dΩ

op = AffineFEOperator(a,l,V,V)
A = get_matrix(op)
Matrix(A)
uh = solve(op)

eh = uh - u_exact
Expand Down

0 comments on commit f07fba1

Please sign in to comment.