diff src/Grids/manifolds.jl @ 1584:d7483e8af705 feature/sbp_operators/laplace_curvilinear

Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 Apr 2024 08:45:54 +0200
parents f77c5309dd2b
children 84c3b9d71218
line wrap: on
line diff
--- a/src/Grids/manifolds.jl	Thu Apr 25 10:20:43 2024 +0200
+++ b/src/Grids/manifolds.jl	Fri Apr 26 08:45:54 2024 +0200
@@ -18,6 +18,7 @@
 [`Simplex`](@ref),
 """
 abstract type ParameterSpace{D} end
+Base.ndims(::ParameterSpace{D}) where D = D
 
 struct HyperBox{T,D} <: ParameterSpace{D}
     a::SVector{D,T}
@@ -43,53 +44,60 @@
 unithyperbox(D) = unithyperbox(Float64,D)
 
 
-struct Simplex{T,D} <: ParameterSpace{D}
-    verticies::NTuple{D,SVector{D,T}}
+struct Simplex{T,D,NV} <: ParameterSpace{D}
+    verticies::NTuple{NV,SVector{D,T}}
 end
 
+Simplex(verticies::Vararg{AbstractArray}) = Simplex(Tuple(SVector(v...) for v ∈ verticies))
+
+verticies(s::Simplex) = s.verticies
+
 Triangle{T} = Simplex{T,2}
 Tetrahedron{T} = Simplex{T,3}
 
-unittriangle(T) = unitsimplex(T,2)
-unittetrahedron(T) = unitsimplex(T,3)
+unittriangle(T=Float64) = unitsimplex(T,2)
+unittetrahedron(T=Float64) = unitsimplex(T,3)
 function unitsimplex(T,D)
     z = @SVector zeros(T,D)
     unitelement = one(eltype(z))
-    verticies = ntuple(i->setindex(z, unitelement, i), 4)
-    return Simplex(verticies)
+    verticies = ntuple(i->setindex(z, unitelement, i), D)
+    return Simplex((z,verticies...))
 end
-
+unitsimplex(D) = unitsimplex(Float64, D)
 
 """
+    Chart{D}
 
 A parametrized description of a manifold or part of a manifold.
-
-Should implement a methods for
-* `parameterspace`
-* `(::Chart)(ξs...)`
-
-There is a default implementation for `(::Chart{D})(::SVector{D})`
 """
-abstract type Chart{D} end
-# abstract type Chart{D,R} end
+struct Chart{D, PST<:ParameterSpace{D}, MT}
+    mapping::MT
+    parameterspace::PST
+end
 
 domain_dim(::Chart{D}) where D = D
-# range_dim(::Chart{D,R}) where {D,R} = R
+(c::Chart)(ξ) = c.mapping(ξ)
+parameterspace(c::Chart) = c.parameterspace
 
 """
-The parameterspace of a chart
-"""
-function parameterspace end
-
-(c::Chart{D})(x̄::SVector{D}) where D = c(x̄...)
-
+    jacobian(c::Chart, ξ)
 
-struct ConcereteChart{PST<:ParameterSpace, MT}
-    parameterspace::PST
-    mapping::MT
-end
-
-(c::Chart)(x̄) = c.mapping(x̄)
+The jacobian of the mapping evaluated at `ξ`. This defers to the
+implementation of `jacobian` for the mapping itself. If no implementation is
+available one can easily be specified for either the mapping function or the
+chart itself.
+```julia
+c = Chart(f, ps)
+jacobian(f::typeof(f), ξ) = f′(ξ)
+```
+or
+```julia
+c = Chart(f, ps)
+jacobian(c::typeof(c),ξ) = f′(ξ)
+```
+which will both allow calling `jacobian(c,ξ)`.
+"""
+jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ)
 
 
 """