From 1c0176147b653135ff6a0355e85a0a864dd86fd2 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Fri, 17 Mar 2023 17:11:01 +0100 Subject: [PATCH] Implement support for Wedge/Triangular prism elements with corresponding linear and quadratic interpolation (#581) --- CHANGELOG.md | 12 + docs/src/devdocs/interpolations.md | 22 +- src/Dofs/ConstraintHandler.jl | 27 +- src/Dofs/DofHandler.jl | 125 +++--- src/Dofs/MixedDofHandler.jl | 98 ++--- src/Export/VTK.jl | 1 + src/FEValues/face_values.jl | 32 +- src/Ferrite.jl | 29 ++ src/Grid/grid.jl | 38 +- src/Grid/grid_generators.jl | 49 +++ src/Quadrature/gaussquad_prism_table.jl | 350 ++++++++++++++++ src/Quadrature/quadrature.jl | 19 + src/deprecations.jl | 5 + src/exports.jl | 2 + src/interpolations.jl | 518 +++++++++++++++++------- src/iterators.jl | 4 +- test/test_facevalues.jl | 4 +- test/test_interpolations.jl | 89 +++- test/test_utils.jl | 14 +- 19 files changed, 1102 insertions(+), 336 deletions(-) create mode 100644 src/Quadrature/gaussquad_prism_table.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index fa565f8e60..480cbc9b5e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Added + - Support for classical trilinear and triquadratic wedge elements. ([#581][github-581]) + - Symmetric quadrature rules up to order 10 for prismatic elements. ([#581][github-581]) + - Finer granulation of dof distribution, allowing to distribute different amounts of dofs + per entity. ([#581][github-581]) +### Fixed + - Dof distribution for embedded elements. ([#581][github-581]) +### Other improvements + - To clarify the dof management `vertices(ip)`, `edges(ip)` and `faces(ip)` has been + deprecated in favor of `vertexdof_indices(ip)`, `edgedof_indices(ip)` and + `facedof_indices(ip)`. ([#578][github-578]) ## [0.3.12] - 2023-02-28 ### Added @@ -324,6 +335,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [github-574]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/574 [github-575]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/575 [github-578]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/578 +[github-581]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/581 [github-583]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/583 [github-588]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/588 [github-591]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/591 diff --git a/docs/src/devdocs/interpolations.md b/docs/src/devdocs/interpolations.md index 5dfe0bfe5f..08fbd570d0 100644 --- a/docs/src/devdocs/interpolations.md +++ b/docs/src/devdocs/interpolations.md @@ -13,19 +13,27 @@ Ferrite.getrefshape(::Interpolation) Ferrite.getorder(::Interpolation) Ferrite.value(::Interpolation{dim}, ::Vec{dim,T}) where {dim,T} Ferrite.derivative(::Interpolation{dim}, ::Vec{dim}) where {dim} +Ferrite.boundarydof_indices ``` ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element +Depending on the dimension of the reference element the following functions have to be implemented + ```@docs Ferrite.value(::Interpolation, ::Int, ::Vec) -Ferrite.vertices(::Interpolation) -Ferrite.nvertexdofs(::Interpolation) -Ferrite.faces(::Interpolation) -Ferrite.nfacedofs(::Interpolation) -Ferrite.edges(::Interpolation) -Ferrite.nedgedofs(::Interpolation) -Ferrite.ncelldofs(::Interpolation) +Ferrite.vertexdof_indices(::Interpolation) +Ferrite.facedof_indices(::Interpolation) +Ferrite.facedof_interior_indices(::Interpolation) +Ferrite.edgedof_indices(::Interpolation) +Ferrite.edgedof_interior_indices(::Interpolation) +Ferrite.celldof_interior_indices(::Interpolation) Ferrite.getnbasefunctions(::Interpolation) Ferrite.reference_coordinates(::Interpolation) ``` + +for all entities which exist on that reference element. The dof functions default to having no +dofs defined on a specific entity. Hence, not overloading of the dof functions will result in an +element with zero dofs. Also, it should always be double checked that everything is consistent as +specified in the docstring of the corresponding function, as inconsistent implementations can +lead to bugs which are really difficult to track down. diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index c95452b115..635e24e7ca 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -308,7 +308,7 @@ end function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) where {Index<:BoundaryIndex} local_face_dofs, local_face_dofs_offset = - _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundaryfunction(eltype(bcfaces))) + _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundarydof_indices(eltype(bcfaces))) copy!(dbc.local_face_dofs, local_face_dofs) copy!(dbc.local_face_dofs_offset, local_face_dofs_offset) @@ -337,7 +337,7 @@ end # Calculate which local dof index live on each face: # face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1] -function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=faces) where F +function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=facedof_indices) where F @assert issorted(components) local_face_dofs = Int[] local_face_dofs_offset = Int[1] @@ -442,24 +442,23 @@ function update!(ch::ConstraintHandler, time::Real=0.0) return nothing end -# for faces -function _update!(inhomogeneities::Vector{Float64}, f::Function, faces::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int}, - components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues, +# for vertices, faces and edges +function _update!(inhomogeneities::Vector{Float64}, f::Function, boundary_entities::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int}, + components::Vector{Int}, dh::AbstractDofHandler, boundaryvalues::BCValues, dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where {T} cc = CellCache(dh, UpdateFlags(; nodes=false, coords=true, dofs=true)) - for (cellidx, faceidx) in faces + for (cellidx, entityidx) in boundary_entities reinit!(cc, cellidx) - # no need to reinit!, enough to update current_face since we only need geometric shape functions M - facevalues.current_face[] = faceidx + # no need to reinit!, enough to update current_entity since we only need geometric shape functions M + boundaryvalues.current_entity[] = entityidx # local dof-range for this face - r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1) + r = local_face_dofs_offset[entityidx]:(local_face_dofs_offset[entityidx+1]-1) counter = 1 - - for location in 1:getnquadpoints(facevalues) - x = spatial_coordinate(facevalues, location, cc.coords) + for location in 1:getnquadpoints(boundaryvalues) + x = spatial_coordinate(boundaryvalues, location, cc.coords) bc_value = f(x, time) @assert length(bc_value) == length(components) @@ -1302,7 +1301,7 @@ end function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{2,<:Union{RefCube,RefTetrahedron}}, n::Int) # For 2D we always permute since Ferrite defines dofs counter-clockwise ret = collect(1:length(local_face_dofs)) - for (i, f) in enumerate(faces(ip)) + for (i, f) in enumerate(facedof_indices(ip)) this_offset = local_face_dofs_offset[i] other_offset = this_offset + n for d in 1:n @@ -1323,7 +1322,7 @@ function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange ret = collect(1:length(local_face_dofs)) # Mirror by changing from counter-clockwise to clockwise - for (i, f) in enumerate(faces(ip)) + for (i, f) in enumerate(facedof_indices(ip)) r = local_face_dofs_offset[i]:(local_face_dofs_offset[i+1] - 1) # 1. Rotate the corners vertex_range = r[1:(N*n)] diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 2cc46d4de7..aafd6f255a 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -9,13 +9,12 @@ Operates slightly faster than [`MixedDofHandler`](@ref). Supports: - `Grid`s with a single concrete cell type. - One or several fields on the whole domaine. """ -struct DofHandler{dim,T,G<:AbstractGrid{dim}} <: AbstractDofHandler +struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler field_names::Vector{Symbol} field_dims::Vector{Int} # TODO: field_interpolations can probably be better typed: We should at least require # all the interpolations to have the same dimension and reference shape field_interpolations::Vector{Interpolation} - bc_values::Vector{BCValues{T}} # TODO: BcValues is created/handeld by the constrainthandler, so this can be removed cell_dofs::Vector{Int} cell_dofs_offset::Vector{Int} closed::ScalarWrapper{Bool} @@ -25,7 +24,7 @@ end function DofHandler(grid::AbstractGrid) isconcretetype(getcelltype(grid)) || error("Grid includes different celltypes. Use MixedDofHandler instead of DofHandler") - DofHandler(Symbol[], Int[], Interpolation[], BCValues{Float64}[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1)) + DofHandler(Symbol[], Int[], Interpolation[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1)) end function Base.show(io::IO, ::MIME"text/plain", dh::DofHandler) @@ -70,7 +69,6 @@ end getfieldinterpolation(dh::DofHandler, field_idx::Int) = dh.field_interpolations[field_idx] getfielddim(dh::DofHandler, field_idx::Int) = dh.field_dims[field_idx] -getbcvalue(dh::DofHandler, field_idx::Int) = dh.bc_values[field_idx] function getfielddim(dh::DofHandler, name::Symbol) field_pos = findfirst(i->i == name, getfieldnames(dh)) @@ -108,9 +106,9 @@ end Add a `dim`-dimensional `Field` called `name` which is approximated by `ip` to `dh`. -The field is added to all cells of the underlying grid. In case no interpolation `ip` is given, -the default interpolation of the grid's celltype is used. -If the grid uses several celltypes, [`add!(dh::MixedDofHandler, fh::FieldHandler)`](@ref) must be used instead. +The field is added to all cells of the underlying grid. In case no interpolation `ip` is +given, the default interpolation of the grid's celltype is used. If the grid uses several +celltypes, [`add!(dh::MixedDofHandler, fh::FieldHandler)`](@ref) must be used instead. """ function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(dh.grid))) @assert !isclosed(dh) @@ -118,10 +116,14 @@ function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_ push!(dh.field_names, name) push!(dh.field_dims, dim) push!(dh.field_interpolations, ip) - push!(dh.bc_values, BCValues(ip, default_interpolation(getcelltype(dh.grid)))) return dh end +# Method for supporting dim=1 default +function add!(dh::DofHandler, name::Symbol, ip::Interpolation=default_interpolation(getcelltype(dh.grid))) + return add!(dh, name, 1, ip) +end + """ sortedge(edge::Tuple{Int,Int}) @@ -203,7 +205,7 @@ function __close!(dh::DofHandler{dim}) where {dim} end # not implemented yet: more than one facedof per face in 3D - dim == 3 && @assert(!any(x->x.nfacedofs > 1, interpolation_infos)) + dim == 3 && @assert(!any(x->any(y->y > 1, x.nfacedofs), interpolation_infos)) nextdof = 1 # next free dof to distribute push!(dh.cell_dofs_offset, 1) # dofs for the first cell start at 1 @@ -211,89 +213,112 @@ function __close!(dh::DofHandler{dim}) where {dim} # loop over all the cells, and distribute dofs for all the fields for (ci, cell) in enumerate(getcells(dh.grid)) @debug println("cell #$ci") - for fi in 1:nfields(dh) - interpolation_info = interpolation_infos[fi] - @debug println(" field: $(dh.field_names[fi])") - if interpolation_info.nvertexdofs > 0 - for vertex in vertices(cell) + for field_idx in 1:nfields(dh) + interpolation_info = interpolation_infos[field_idx] + @debug println(" field: $(dh.field_names[field_idx])") + for (vi, vertex) in enumerate(vertices(cell)) + if interpolation_info.nvertexdofs[vi] > 0 @debug println(" vertex#$vertex") - token = Base.ht_keyindex2!(vertexdicts[fi], vertex) + token = Base.ht_keyindex2!(vertexdicts[field_idx], vertex) if token > 0 # haskey(vertexdicts[fi], vertex) # reuse dofs - reuse_dof = vertexdicts[fi].vals[token] # vertexdicts[fi][vertex] - for d in 1:dh.field_dims[fi] + reuse_dof = vertexdicts[field_idx].vals[token] # vertexdicts[fi][vertex] + for d in 1:dh.field_dims[field_idx] @debug println(" reusing dof #$(reuse_dof + (d-1))") push!(dh.cell_dofs, reuse_dof + (d-1)) end else # token <= 0, distribute new dofs - for vertexdof in 1:interpolation_info.nvertexdofs - Base._setindex!(vertexdicts[fi], nextdof, vertex, -token) # vertexdicts[fi][vertex] = nextdof - for d in 1:dh.field_dims[fi] + for vertexdof in 1:interpolation_info.nvertexdofs[vi] + Base._setindex!(vertexdicts[field_idx], nextdof, vertex, -token) # vertexdicts[field_idx][vertex] = nextdof + for d in 1:dh.field_dims[field_idx] @debug println(" adding dof#$nextdof") push!(dh.cell_dofs, nextdof) nextdof += 1 end end end - end # vertex loop - end + end + end # vertex loop if dim == 3 # edges only in 3D - if interpolation_info.nedgedofs > 0 - for edge in edges(cell) + for (ei, edge) in enumerate(edges(cell)) + if interpolation_info.dim == 3 && interpolation_info.nedgedofs[ei] > 0 sedge, dir = sortedge(edge) @debug println(" edge#$sedge dir: $(dir)") - token = Base.ht_keyindex2!(edgedicts[fi], sedge) - if token > 0 # haskey(edgedicts[fi], sedge), reuse dofs - startdof, olddir = edgedicts[fi].vals[token] # edgedicts[fi][sedge] # first dof for this edge (if dir == true) - for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs) : (interpolation_info.nedgedofs:-1:1)) - for d in 1:dh.field_dims[fi] - reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[fi] + token = Base.ht_keyindex2!(edgedicts[field_idx], sedge) + if token > 0 # haskey(edgedicts[field_idx], sedge), reuse dofs + startdof, olddir = edgedicts[field_idx].vals[token] # edgedicts[field_idx][sedge] # first dof for this edge (if dir == true) + for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs[ei]) : (interpolation_info.nedgedofs[ei]:-1:1)) + for d in 1:dh.field_dims[field_idx] + reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[field_idx] @debug println(" reusing dof#$(reuse_dof)") push!(dh.cell_dofs, reuse_dof) end end else # token <= 0, distribute new dofs - Base._setindex!(edgedicts[fi], (nextdof, dir), sedge, -token) # edgedicts[fi][sedge] = (nextdof, dir), store only the first dof for the edge - for edgedof in 1:interpolation_info.nedgedofs - for d in 1:dh.field_dims[fi] + Base._setindex!(edgedicts[field_idx], (nextdof, dir), sedge, -token) # edgedicts[field_idx][sedge] = (nextdof, dir), store only the first dof for the edge + for edgedof in 1:interpolation_info.nedgedofs[ei] + for d in 1:dh.field_dims[field_idx] @debug println(" adding dof#$nextdof") push!(dh.cell_dofs, nextdof) nextdof += 1 end end end - end # edge loop - end + elseif interpolation_info.dim == 2 && interpolation_info.nfacedofs[ei] > 0 # hotfix for current embedded elements + sedge, dir = sortedge(edge) + @debug println(" edge#$sedge dir: $(dir)") + token = Base.ht_keyindex2!(edgedicts[field_idx], sedge) + if token > 0 # haskey(edgedicts[field_idx], sedge), reuse dofs + startdof, olddir = edgedicts[field_idx].vals[token] # edgedicts[field_idx][sedge] # first dof for this edge (if dir == true) + for edgedof in (dir == olddir ? (1:interpolation_info.nfacedofs[ei]) : (interpolation_info.nfacedofs[ei]:-1:1)) + for d in 1:dh.field_dims[field_idx] + reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[field_idx] + @debug println(" reusing dof#$(reuse_dof)") + push!(dh.cell_dofs, reuse_dof) + end + end + else # token <= 0, distribute new dofs + Base._setindex!(edgedicts[field_idx], (nextdof, dir), sedge, -token) # edgedicts[field_idx][sedge] = (nextdof, dir), store only the first dof for the edge + for edgedof in 1:interpolation_info.nfacedofs[ei] + for d in 1:dh.field_dims[field_idx] + @debug println(" adding dof#$nextdof") + push!(dh.cell_dofs, nextdof) + nextdof += 1 + end + end + end + end + end # edge loop end - if interpolation_info.nfacedofs > 0 && (interpolation_info.dim == dim) - for face in faces(cell) - sface = sortface(face) # TODO: faces(cell) may as well just return the sorted list + for (fi, face) in enumerate(faces(cell)) + if interpolation_info.nfacedofs[fi] > 0 && (interpolation_info.dim == dim) + sface = sortface(face) @debug println(" face#$sface") - token = Base.ht_keyindex2!(facedicts[fi], sface) - if token > 0 # haskey(facedicts[fi], sface), reuse dofs - startdof = facedicts[fi].vals[token] # facedicts[fi][sface] - for facedof in interpolation_info.nfacedofs:-1:1 # always reverse (YOLO) - for d in 1:dh.field_dims[fi] - reuse_dof = startdof + (d-1) + (facedof-1)*dh.field_dims[fi] + token = Base.ht_keyindex2!(facedicts[field_idx], sface) + if token > 0 # haskey(facedicts[field_idx], sface), reuse dofs + startdof = facedicts[field_idx].vals[token] # facedicts[field_idx][sface] + for facedof in interpolation_info.nfacedofs[fi]:-1:1 # always reverse (YOLO) + for d in 1:dh.field_dims[field_idx] + reuse_dof = startdof + (d-1) + (facedof-1)*dh.field_dims[field_idx] @debug println(" reusing dof#$(reuse_dof)") push!(dh.cell_dofs, reuse_dof) end end else # distribute new dofs - Base._setindex!(facedicts[fi], nextdof, sface, -token)# facedicts[fi][sface] = nextdof, store the first dof for this face - for facedof in 1:interpolation_info.nfacedofs - for d in 1:dh.field_dims[fi] + Base._setindex!(facedicts[field_idx], nextdof, sface, -token)# facedicts[field_idx][sface] = nextdof, store the first dof for this face + for facedof in 1:interpolation_info.nfacedofs[fi] + for d in 1:dh.field_dims[field_idx] @debug println(" adding dof#$nextdof") push!(dh.cell_dofs, nextdof) nextdof += 1 end end end - end # face loop - end + end + end # face loop if interpolation_info.ncelldofs > 0 # always distribute new dofs for cell @debug println(" cell#$ci") for celldof in 1:interpolation_info.ncelldofs - for d in 1:dh.field_dims[fi] + for d in 1:dh.field_dims[field_idx] @debug println(" adding dof#$nextdof") push!(dh.cell_dofs, nextdof) nextdof += 1 diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index 1bc5973286..d6468a0d39 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -58,7 +58,7 @@ end function MixedDofHandler(grid::Grid{dim,C,T}) where {dim,C,T} ncells = getncells(grid) - MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], CellVector(Int[],zeros(Int,ncells),zeros(Int,ncells)), CellVector(Int[],Int[],Int[]), CellVector(Vec{dim,T}[],Int[],Int[]), Ferrite.ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1)) + MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], CellVector(Int[],zeros(Int,ncells),zeros(Int,ncells)), CellVector(Int[],Int[],Int[]), CellVector(Vec{dim,T}[],Int[],Int[]), ScalarWrapper(false), grid, ScalarWrapper(-1)) end getfieldnames(fh::FieldHandler) = [field.name for field in fh.fields] @@ -205,9 +205,8 @@ function close!(dh::MixedDofHandler) end function __close!(dh::MixedDofHandler{dim}) where {dim} - - @assert !Ferrite.isclosed(dh) - field_names = Ferrite.getfieldnames(dh) # all the fields in the problem + @assert !isclosed(dh) + field_names = getfieldnames(dh) # all the fields in the problem numfields = length(field_names) # Create dicts that store created dofs @@ -228,9 +227,9 @@ function __close!(dh::MixedDofHandler{dim}) where {dim} dh, cellnumbers, field_names, - Ferrite.getfieldnames(fh), - Ferrite.getfielddims(fh), - Ferrite.getfieldinterpolations(fh), + getfieldnames(fh), + getfielddims(fh), + getfieldinterpolations(fh), nextdof, vertexdicts, edgedicts, @@ -259,15 +258,12 @@ function __close!(dh::MixedDofHandler{dim}) where {dim} end function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, field_names, field_dims, field_interpolations, nextdof, vertexdicts, edgedicts, facedicts, celldicts) where {dim} - - ip_infos = Ferrite.InterpolationInfo[] + ip_infos = InterpolationInfo[] for interpolation in field_interpolations - ip_info = Ferrite.InterpolationInfo(interpolation) + ip_info = InterpolationInfo(interpolation) # these are not implemented yet (or have not been tested) - @assert(ip_info.nvertexdofs <= 1) - @assert(ip_info.nedgedofs <= 1) - @assert(ip_info.nfacedofs <= 1) - @assert(ip_info.ncelldofs <= 1) # not tested but probably works + @assert(all(ip_info.nedgedofs .<= 1)) + @assert(all(ip_info.nfacedofs .<= 1)) push!(ip_infos, ip_info) end @@ -285,22 +281,25 @@ function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fiel @debug "\tfield: $(field_name)" ip_info = ip_infos[local_num] - if ip_info.nvertexdofs > 0 - nextdof = add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[local_num], ip_info.nvertexdofs, nextdof) - end + # We first distribute the vertex dofs + nextdof = add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[local_num], ip_info.nvertexdofs, nextdof) - if ip_info.nedgedofs > 0 && dim == 3 #Edges only in 3d - nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nedgedofs, nextdof) + # Then the edge dofs + if dim == 3 + if ip_info.dim == 3 # Regular 3D element + nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nedgedofs, nextdof) + elseif ip_info.dim == 2 # 2D embedded element in 3D + nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nfacedofs, nextdof) + end end - if ip_info.nfacedofs > 0 && (ip_info.dim == dim) + # Then the face dofs + if ip_info.dim == dim # Regular 3D element nextdof = add_face_dofs(cell_dofs, cell, facedicts[fi], field_dims[local_num], ip_info.nfacedofs, nextdof) end - if ip_info.ncelldofs > 0 - nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof) - end - + # And finally the celldofs + nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof) end # after done creating dofs for the cell, push them to the global list append!(dh.cell_dofs.values, cell_dofs) @@ -316,13 +315,11 @@ Returns the next global dof number and an array of dofs. If dofs have already been created for the object (vertex, face) then simply return those, otherwise create new dofs. """ function get_or_create_dofs!(nextdof, field_dim; dict, key) - token = Base.ht_keyindex2!(dict, key) if token > 0 # vertex, face etc. visited before # reuse stored dofs (TODO unless field is discontinuous) @debug "\t\tkey: $key dofs: $(dict[key]) (reused dofs)" return nextdof, dict[key] - else # create new dofs dofs = nextdof : (nextdof + field_dim-1) @debug "\t\tkey: $key dofs: $dofs" @@ -333,33 +330,39 @@ function get_or_create_dofs!(nextdof, field_dim; dict, key) end function add_vertex_dofs(cell_dofs, cell, vertexdict, field_dim, nvertexdofs, nextdof) - for vertex in Ferrite.vertices(cell) - @debug "\tvertex #$vertex" - nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex) - append!(cell_dofs, dofs) + for (vi, vertex) in enumerate(vertices(cell)) + if nvertexdofs[vi] > 0 + @debug "\tvertex #$vertex" + nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex) + append!(cell_dofs, dofs) + end end return nextdof end function add_face_dofs(cell_dofs, cell, facedict, field_dim, nfacedofs, nextdof) - @assert nfacedofs == 1 "Currently only supports interpolations with nfacedofs = 1" - - for face in Ferrite.faces(cell) - sface = Ferrite.sortface(face) - @debug "\tface #$sface" - nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface) - append!(cell_dofs, dofs) + @debug @assert all(nfacedofs .<= 1) "Currently only supports interpolations with less that 2 dofs per face" + + for (fi,face) in enumerate(faces(cell)) + if nfacedofs[fi] > 0 + sface = sortface(face) + @debug "\tface #$sface" + nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface) + # TODO permutate dofs according to face orientation + append!(cell_dofs, dofs) + end end return nextdof end function add_edge_dofs(cell_dofs, cell, edgedict, field_dim, nedgedofs, nextdof) - @assert nedgedofs == 1 "Currently only supports interpolations with nedgedofs = 1" - for edge in Ferrite.edges(cell) - sedge, dir = Ferrite.sortedge(edge) - @debug "\tedge #$sedge" - nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge) - append!(cell_dofs, dofs) + for (ei,edge) in enumerate(edges(cell)) + if nedgedofs[ei] > 0 + sedge, dir = sortedge(edge) + @debug "\tedge #$sedge" + nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge) + append!(cell_dofs, dofs) + end end return nextdof end @@ -389,9 +392,9 @@ function field_offset(fh::FieldHandler, field_name::Symbol) return offset end -function Ferrite.dof_range(fh::FieldHandler, field_name::Symbol) - f = Ferrite.find_field(fh, field_name) - offset = Ferrite.field_offset(fh, field_name) +function dof_range(fh::FieldHandler, field_name::Symbol) + f = find_field(fh, field_name) + offset = field_offset(fh, field_name) field_interpolation = fh.fields[f].interpolation field_dim = fh.fields[f].dim n_field_dofs = getnbasefunctions(field_interpolation)::Int * field_dim @@ -406,9 +409,8 @@ getfielddim(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].dim getfielddim(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].dim function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) where T - # make sure the field exists - fieldname ∈ Ferrite.getfieldnames(dh) || error("Field $fieldname not found.") + fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.") field_dim = getfielddim(dh, fieldname) space_dim = field_dim == 2 ? 3 : field_dim diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 1977d8a692..2d54ad99c4 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -14,6 +14,7 @@ cell_to_vtkcell(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON cell_to_vtkcell(::Type{Cell{3,20,6}}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON cell_to_vtkcell(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA cell_to_vtkcell(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA +cell_to_vtkcell(::Type{Wedge}) = VTKCellTypes.VTK_WEDGE nodes_to_vtkorder(cell::AbstractCell) = collect(cell.nodes) diff --git a/src/FEValues/face_values.jl b/src/FEValues/face_values.jl index 25e941b18d..6b65180627 100644 --- a/src/FEValues/face_values.jl +++ b/src/FEValues/face_values.jl @@ -232,48 +232,50 @@ for each dof-position determined by the `func_interpol`. Used mainly by the `Con """ struct BCValues{T} M::Array{T,3} - current_face::ScalarWrapper{Int} + nqp::Array{Int} + current_entity::ScalarWrapper{Int} end BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) = BCValues(Float64, func_interpol, geom_interpol, boundary_type) function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_interpol::Interpolation{dim,refshape}, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) where {T,dim,refshape} - # set up quadrature rules for each face with dof-positions + # set up quadrature rules for each boundary entity with dof-positions # (determined by func_interpol) as the quadrature points interpolation_coords = reference_coordinates(func_interpol) qrs = QuadratureRule{dim,refshape,T}[] - faces = boundaryfunction(boundary_type) # faces, edges or vertices - for face in faces(func_interpol) + for boundarydofs in boundarydof_indices(boundary_type)(func_interpol) dofcoords = Vec{dim,T}[] - for facedof in face - push!(dofcoords, interpolation_coords[facedof]) + for boundarydof in boundarydofs + push!(dofcoords, interpolation_coords[boundarydof]) end qrf = QuadratureRule{dim,refshape,T}(fill(T(NaN), length(dofcoords)), dofcoords) # weights will not be used push!(qrs, qrf) end - n_faces = length(qrs) - n_qpoints = n_faces == 0 ? 0 : length(getweights(qrs[1])) # assume same in all + n_boundary_entities = length(qrs) + n_qpoints = n_boundary_entities == 0 ? 0 : maximum(qr->length(getweights(qr)), qrs) # Bound number of qps correctly. n_geom_basefuncs = getnbasefunctions(geom_interpol) - M = fill(zero(T) * T(NaN), n_geom_basefuncs, n_qpoints, n_faces) + M = fill(zero(T) * T(NaN), n_geom_basefuncs, n_qpoints, n_boundary_entities) + nqp = zeros(Int,n_boundary_entities) - for face in 1:n_faces, (qp, ξ) in enumerate(qrs[face].points) - for i in 1:n_geom_basefuncs - M[i, qp, face] = value(geom_interpol, i, ξ) + for n_boundary_entity in 1:n_boundary_entities + for (qp, ξ) in enumerate(qrs[n_boundary_entity].points), i in 1:n_geom_basefuncs + M[i, qp, n_boundary_entity] = value(geom_interpol, i, ξ) end + nqp[n_boundary_entity] = length(qrs[n_boundary_entity].points) end - BCValues{T}(M, ScalarWrapper(0)) + BCValues{T}(M, nqp, ScalarWrapper(0)) end -getnquadpoints(bcv::BCValues) = size(bcv.M, 2) +getnquadpoints(bcv::BCValues) = bcv.nqp[bcv.current_entity.x] function spatial_coordinate(bcv::BCValues, q_point::Int, xh::AbstractVector{Vec{dim,T}}) where {dim,T} n_base_funcs = size(bcv.M, 1) length(xh) == n_base_funcs || throw_incompatible_coord_length(length(xh), n_base_funcs) x = zero(Vec{dim,T}) - face = bcv.current_face[] + face = bcv.current_entity[] @inbounds for i in 1:n_base_funcs x += bcv.M[i,q_point,face] * xh[i] # geometric_value(fe_v, q_point, i) * xh[i] end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 63f26d8f46..3ee75b6e95 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -20,6 +20,7 @@ abstract type AbstractRefShape end struct RefTetrahedron <: AbstractRefShape end struct RefCube <: AbstractRefShape end +struct RefPrism <: AbstractRefShape end """ Abstract type which has `CellValues` and `FaceValues` as subtypes @@ -33,6 +34,34 @@ Abstract type which is used as identifier for faces, edges and verices """ abstract type BoundaryIndex end +""" +A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh +""" +struct CellIndex + idx::Int +end + +""" +A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face). +""" +struct FaceIndex <: BoundaryIndex + idx::Tuple{Int,Int} # cell and side +end + +""" +A `EdgeIndex` wraps an (Int, Int) and defines a local edge by pointing to a (cell, edge). +""" +struct EdgeIndex <: BoundaryIndex + idx::Tuple{Int,Int} # cell and side +end + +""" +A `VertexIndex` wraps an (Int, Int) and defines a local vertex by pointing to a (cell, vert). +""" +struct VertexIndex <: BoundaryIndex + idx::Tuple{Int,Int} # cell and side +end + include("utils.jl") # Matrix/Vector utilities diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index b1d24a395f..ee862abab1 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -100,36 +100,10 @@ const implemented_celltypes = ( (const QuadraticTetrahedron = Cell{3,10,4}), (const Hexahedron = Cell{3,8,6}), - (Cell{2,20,6}) -) - -""" -A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh -""" -struct CellIndex - idx::Int -end - -""" -A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face). -""" -struct FaceIndex <: BoundaryIndex - idx::Tuple{Int,Int} # cell and side -end + (Cell{2,20,6}), -""" -A `EdgeIndex` wraps an (Int, Int) and defines a local edge by pointing to a (cell, edge). -""" -struct EdgeIndex <: BoundaryIndex - idx::Tuple{Int,Int} # cell and side -end - -""" -A `VertexIndex` wraps an (Int, Int) and defines a local vertex by pointing to a (cell, vert). -""" -struct VertexIndex <: BoundaryIndex - idx::Tuple{Int,Int} # cell and side -end + (const Wedge = Cell{3,6,5}) +) struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}} neighbor_info::Vector{T} @@ -863,6 +837,11 @@ faces(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c. edges(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1])) faces(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2],c.nodes[3],c.nodes[4]),) +vertices(c::Wedge) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6]) +edges(c::Wedge) = ((c.nodes[2],c.nodes[1]), (c.nodes[1],c.nodes[3]), (c.nodes[1],c.nodes[4]), (c.nodes[3],c.nodes[2]), (c.nodes[2],c.nodes[5]), (c.nodes[3],c.nodes[6]), (c.nodes[4],c.nodes[5]), (c.nodes[4],c.nodes[6]), (c.nodes[6],c.nodes[5])) +faces(c::Wedge) = ((c.nodes[1],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[5],c.nodes[4]), (c.nodes[3],c.nodes[1],c.nodes[4],c.nodes[6]), (c.nodes[2],c.nodes[3],c.nodes[6],c.nodes[5]), (c.nodes[4],c.nodes[5],c.nodes[6])) + +# random stuff default_interpolation(::Union{Type{Line},Type{Line2D},Type{Line3D}}) = Lagrange{1,RefCube,1}() default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}() default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}() @@ -873,6 +852,7 @@ default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}() default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3,RefTetrahedron,2}() default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}() default_interpolation(::Type{Cell{3,20,6}}) = Serendipity{3,RefCube,2}() +default_interpolation(::Type{Wedge}) = Lagrange{3,RefPrism,1}() """ boundaryfunction(::Type{<:BoundaryIndex}) diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl index 7b2b0e424a..d5ca1c8f36 100644 --- a/src/Grid/grid_generators.jl +++ b/src/Grid/grid_generators.jl @@ -239,6 +239,55 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix) end +# Wedge +function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T} + nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z + n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1 + n_nodes = n_nodes_x * n_nodes_y * n_nodes_z + + # Generate nodes + coords_x = range(left[1], stop=right[1], length=n_nodes_x) + coords_y = range(left[2], stop=right[2], length=n_nodes_y) + coords_z = range(left[3], stop=right[3], length=n_nodes_z) + nodes = Node{3,T}[] + for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x + push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k]))) + end + + # Generate cells + node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z)) + cells = Wedge[] + for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x + push!(cells, Wedge((node_array[i,j,k], node_array[i+1,j,k], node_array[i,j+1,k], + node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i,j+1,k+1]))) # ◺ + push!(cells, Wedge((node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k], + node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1]))) # ◹ + end + + # Order the cells as c_nxyz[2, x, y, z] such that we can look up boundary cells + c_nxyz = reshape(1:length(cells), (2, nel...)) + + @views le = map(x -> FaceIndex(x,3), c_nxyz[1, 1, :, :][:]) + @views ri = map(x -> FaceIndex(x,2), c_nxyz[2, end, :, :][:]) + @views fr = map(x -> FaceIndex(x,2), c_nxyz[1, :, 1, :][:]) + @views ba = map(x -> FaceIndex(x,4), c_nxyz[2, :, end, :][:]) + @views bo = [map(x -> FaceIndex(x,1), c_nxyz[1, :, :, 1][:]) ; map(x -> FaceIndex(x,1), c_nxyz[2, :, :, 1][:])] + @views to = [map(x -> FaceIndex(x,5), c_nxyz[1, :, :, end][:]) ; map(x -> FaceIndex(x,5), c_nxyz[2, :, :, end][:])] + + boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba]) + + facesets = Dict( + "left" => Set(le), + "right" => Set(ri), + "front" => Set(fr), + "back" => Set(ba), + "bottom" => Set(bo), + "top" => Set(to), + ) + + return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix) +end + function Ferrite.generate_grid(::Type{Cell{3,20,6}}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T} nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z nnode_x = 2nel_x + 1; nnode_y = 2nel_y + 1; nnode_z = 2nel_z + 1 #Note: not the actually number of nodes in x/y/z, just a temporary variables diff --git a/src/Quadrature/gaussquad_prism_table.jl b/src/Quadrature/gaussquad_prism_table.jl new file mode 100644 index 0000000000..8fe0c88724 --- /dev/null +++ b/src/Quadrature/gaussquad_prism_table.jl @@ -0,0 +1,350 @@ +# Symmetric quadrature rules takes from +# Witherden, Freddie D., and Peter E. Vincent. "On the identification of +# symmetric quadrature rules for finite element methods." Computers & +# Mathematics with Applications 69.10 (2015): 1232-1241. +# Note that the original rule is defined on [-1,1]^3 while our reference prism is defined on [0,1]^3, hence we transform in the end. +function _get_gauss_prismdata_polyquad(n::Int) + if n == 1 + return reshape([1/3 1/3 1/2 1/2], (1,4)) + elseif n == 2 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.99999999999999985119303811076522004183 0.66666666666666686507594918564641756458 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.99999999999999985119303811076522004183 0.66666666666666686507594918564641756458 + -0.74158162379719638007464124598498266439 0.48316324759439276014928249196996532878 0 0.88888888888888875661603387623572162361 + 0.48316324759439276014928249196996532878 -0.74158162379719638007464124598498266439 0 0.88888888888888875661603387623572162361 + -0.74158162379719638007464124598498266439 -0.74158162379719638007464124598498266439 0 0.88888888888888875661603387623572162361 + ] + elseif n == 3 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.86066920446477793249332885107223920359 0.89998695257960196814387903509173005101 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.86066920446477793249332885107223920359 0.89998695257960196814387903509173005101 + -0.52223312808684327863699399822724711768 0.52222829574978345631488370728601281383 0 0.366671015806799343952040321636089983 + 0.52222829574978345631488370728601281383 -0.52223312808684327863699399822724711768 0 0.366671015806799343952040321636089983 + -0.99999516766294017767788970905876569615 0.52222829574978345631488370728601281383 0 0.366671015806799343952040321636089983 + 0.52222829574978345631488370728601281383 -0.99999516766294017767788970905876569615 0 0.366671015806799343952040321636089983 + -0.99999516766294017767788970905876569615 -0.52223312808684327863699399822724711768 0 0.366671015806799343952040321636089983 + -0.52223312808684327863699399822724711768 -0.99999516766294017767788970905876569615 0 0.366671015806799343952040321636089983 + ] + elseif n == 4 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.86686197400903479841583223937151297546 0.43164789926214201915560571539043719684 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.86686197400903479841583223937151297546 0.43164789926214201915560571539043719684 + -0.062688380276009575455822705856210492672 -0.87462323944798084908835458828757901466 0 0.54565845042191057226153838655080393041 + -0.87462323944798084908835458828757901466 -0.062688380276009575455822705856210492672 0 0.54565845042191057226153838655080393041 + -0.062688380276009575455822705856210492672 -0.062688380276009575455822705856210492672 0 0.54565845042191057226153838655080393041 + -0.79851918840217872745314255427215462609 0.59703837680435745490628510854430925218 -0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + 0.59703837680435745490628510854430925218 -0.79851918840217872745314255427215462609 -0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + -0.79851918840217872745314255427215462609 -0.79851918840217872745314255427215462609 -0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + -0.79851918840217872745314255427215462609 0.59703837680435745490628510854430925218 0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + 0.59703837680435745490628510854430925218 -0.79851918840217872745314255427215462609 0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + -0.79851918840217872745314255427215462609 -0.79851918840217872745314255427215462609 0.67563982368225971739661759355428029031 0.24995480836833070748402890159445230251 + ] + elseif n == 5 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0 0.71145555593148684996197343022099556338 + -0.025400070899508667767082370124706051859 -0.94919985820098266446583525975058789628 0 0.22471006722826685202639587426081419504 + -0.94919985820098266446583525975058789628 -0.025400070899508667767082370124706051859 0 0.22471006722826685202639587426081419504 + -0.025400070899508667767082370124706051859 -0.025400070899508667767082370124706051859 0 0.22471006722826685202639587426081419504 + -0.10880379065925599247444934872713650457 -0.78239241868148801505110130254572699087 -0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.78239241868148801505110130254572699087 -0.10880379065925599247444934872713650457 -0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.10880379065925599247444934872713650457 -0.10880379065925599247444934872713650457 -0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.10880379065925599247444934872713650457 -0.78239241868148801505110130254572699087 0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.78239241868148801505110130254572699087 -0.10880379065925599247444934872713650457 0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.10880379065925599247444934872713650457 -0.10880379065925599247444934872713650457 0.8710029348654440527292590830738210619 0.18566142131615813987257487985503458339 + -0.79828210803458293816870337538129563058 0.59656421606916587633740675076259126116 -0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + 0.59656421606916587633740675076259126116 -0.79828210803458293816870337538129563058 -0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + -0.79828210803458293816870337538129563058 -0.79828210803458293816870337538129563058 -0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + -0.79828210803458293816870337538129563058 0.59656421606916587633740675076259126116 0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + 0.59656421606916587633740675076259126116 -0.79828210803458293816870337538129563058 0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + -0.79828210803458293816870337538129563058 -0.79828210803458293816870337538129563058 0.57042698070515927206196786842334325646 0.25007428574779395912056494464439239186 + ] + elseif n == 6 + wx = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.52662270480497475869798007748379244322 0.22269313409222502250038629788629976598 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.52662270480497475869798007748379244322 0.22269313409222502250038629788629976598 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.99083630081924474869286718300205167763 0.13839610937412451172575590081432159364 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.99083630081924474869286718300205167763 0.13839610937412451172575590081432159364 + -0.66959346994381506531292392577889361708 0.33918693988763013062584785155778723416 0 0.31182801481333532900066338577135735737 + 0.33918693988763013062584785155778723416 -0.66959346994381506531292392577889361708 0 0.31182801481333532900066338577135735737 + -0.66959346994381506531292392577889361708 -0.66959346994381506531292392577889361708 0 0.31182801481333532900066338577135735737 + -0.95905346077548616195214584157965381704 0.91810692155097232390429168315930763408 0 0.050146535238696554002545078611329875622 + 0.91810692155097232390429168315930763408 -0.95905346077548616195214584157965381704 0 0.050146535238696554002545078611329875622 + -0.95905346077548616195214584157965381704 -0.95905346077548616195214584157965381704 0 0.050146535238696554002545078611329875622 + -0.067776668819680981202280000599855057023 -0.86444666236063803759543999880028988595 -0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.86444666236063803759543999880028988595 -0.067776668819680981202280000599855057023 -0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.067776668819680981202280000599855057023 -0.067776668819680981202280000599855057023 -0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.067776668819680981202280000599855057023 -0.86444666236063803759543999880028988595 0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.86444666236063803759543999880028988595 -0.067776668819680981202280000599855057023 0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.067776668819680981202280000599855057023 -0.067776668819680981202280000599855057023 0.48167413488438751267079566597126764678 0.19605056020608216687021916179149498019 + -0.59485220654953844609181984939099585377 0.52549348166744454603903998678937679713 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + 0.52549348166744454603903998678937679713 -0.59485220654953844609181984939099585377 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.93064127511790609994722013739838094335 0.52549348166744454603903998678937679713 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + 0.52549348166744454603903998678937679713 -0.93064127511790609994722013739838094335 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.93064127511790609994722013739838094335 -0.59485220654953844609181984939099585377 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.59485220654953844609181984939099585377 -0.93064127511790609994722013739838094335 -0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.59485220654953844609181984939099585377 0.52549348166744454603903998678937679713 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + 0.52549348166744454603903998678937679713 -0.59485220654953844609181984939099585377 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.93064127511790609994722013739838094335 0.52549348166744454603903998678937679713 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + 0.52549348166744454603903998678937679713 -0.93064127511790609994722013739838094335 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.93064127511790609994722013739838094335 -0.59485220654953844609181984939099585377 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + -0.59485220654953844609181984939099585377 -0.93064127511790609994722013739838094335 0.80928579325583275231645432806703865787 0.084632875139559356776397936558477141722 + ] + elseif n == 7 + wx = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.98022806959089160171504882914128008603 0.11378272445075715563392222208600913026 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.98022806959089160171504882914128008603 0.11378272445075715563392222208600913026 + -0.98332480906795705560590831565479581293 0.96664961813591411121181663130959162585 0 0.024472968692380158977782543188506640362 + 0.96664961813591411121181663130959162585 -0.98332480906795705560590831565479581293 0 0.024472968692380158977782543188506640362 + -0.98332480906795705560590831565479581293 -0.98332480906795705560590831565479581293 0 0.024472968692380158977782543188506640362 + -0.036956049341726798726440754485546402356 -0.92608790131654640254711849102890719529 -0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.92608790131654640254711849102890719529 -0.036956049341726798726440754485546402356 -0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.036956049341726798726440754485546402356 -0.036956049341726798726440754485546402356 -0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.036956049341726798726440754485546402356 -0.92608790131654640254711849102890719529 0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.92608790131654640254711849102890719529 -0.036956049341726798726440754485546402356 0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.036956049341726798726440754485546402356 -0.036956049341726798726440754485546402356 0.84138735420652599295096060100359489758 0.086203456344905850875965416526967283444 + -0.80903350325702126768802840491559333013 0.61806700651404253537605680983118666027 -0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + 0.61806700651404253537605680983118666027 -0.80903350325702126768802840491559333013 -0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + -0.80903350325702126768802840491559333013 -0.80903350325702126768802840491559333013 -0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + -0.80903350325702126768802840491559333013 0.61806700651404253537605680983118666027 0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + 0.61806700651404253537605680983118666027 -0.80903350325702126768802840491559333013 0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + -0.80903350325702126768802840491559333013 -0.80903350325702126768802840491559333013 0.79584909058698306626909009833843112408 0.11671409960839383433248991439345238969 + -0.97570173680324342035479366133901540544 0.48599336414579127066992364784746385989 0 0.10205942534059700976786329599159184945 + 0.48599336414579127066992364784746385989 -0.97570173680324342035479366133901540544 0 0.10205942534059700976786329599159184945 + -0.51029162734254785031512998650844845445 0.48599336414579127066992364784746385989 0 0.10205942534059700976786329599159184945 + 0.48599336414579127066992364784746385989 -0.51029162734254785031512998650844845445 0 0.10205942534059700976786329599159184945 + -0.51029162734254785031512998650844845445 -0.97570173680324342035479366133901540544 0 0.10205942534059700976786329599159184945 + -0.97570173680324342035479366133901540544 -0.51029162734254785031512998650844845445 0 0.10205942534059700976786329599159184945 + -0.38968756713554786569859973469754143293 -0.69403080315040476305982696739689230156 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.69403080315040476305982696739689230156 -0.38968756713554786569859973469754143293 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + 0.083718370285952628758426702094433734489 -0.69403080315040476305982696739689230156 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.69403080315040476305982696739689230156 0.083718370285952628758426702094433734489 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + 0.083718370285952628758426702094433734489 -0.38968756713554786569859973469754143293 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.38968756713554786569859973469754143293 0.083718370285952628758426702094433734489 -0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.38968756713554786569859973469754143293 -0.69403080315040476305982696739689230156 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.69403080315040476305982696739689230156 -0.38968756713554786569859973469754143293 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + 0.083718370285952628758426702094433734489 -0.69403080315040476305982696739689230156 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.69403080315040476305982696739689230156 0.083718370285952628758426702094433734489 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + 0.083718370285952628758426702094433734489 -0.38968756713554786569859973469754143293 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + -0.38968756713554786569859973469754143293 0.083718370285952628758426702094433734489 0.4039345605321099116998238341362959742 0.15576281310483042016174134706586605691 + ] + elseif n == 8 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.22443386610599057286592461466529935585 0.082939746571421834033180571412028667484 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.22443386610599057286592461466529935585 0.082939746571421834033180571412028667484 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.68182544157086581170215168002374039522 0.20435770844872619128158195584935122282 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.68182544157086581170215168002374039522 0.20435770844872619128158195584935122282 + -0.89317529812618572284389086719897089711 0.78635059625237144568778173439794179421 0 0.072608995136598867700988758862181257017 + 0.78635059625237144568778173439794179421 -0.89317529812618572284389086719897089711 0 0.072608995136598867700988758862181257017 + -0.89317529812618572284389086719897089711 -0.89317529812618572284389086719897089711 0 0.072608995136598867700988758862181257017 + -0.079822074372578761440178266060788079463 -0.84035585125484247711964346787842384107 0 0.2107989796260286741645220068699562506 + -0.84035585125484247711964346787842384107 -0.079822074372578761440178266060788079463 0 0.2107989796260286741645220068699562506 + -0.079822074372578761440178266060788079463 -0.079822074372578761440178266060788079463 0 0.2107989796260286741645220068699562506 + -0.082861862418097326120000870246452008956 -0.83427627516380534775999825950709598209 -0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.83427627516380534775999825950709598209 -0.082861862418097326120000870246452008956 -0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.082861862418097326120000870246452008956 -0.082861862418097326120000870246452008956 -0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.082861862418097326120000870246452008956 -0.83427627516380534775999825950709598209 0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.83427627516380534775999825950709598209 -0.082861862418097326120000870246452008956 0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.082861862418097326120000870246452008956 -0.082861862418097326120000870246452008956 0.87619591000025424474025327234232345955 0.084346037174345969123924723988899397111 + -0.65187678415125922476703025933216051413 0.30375356830251844953406051866432102826 -0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + 0.30375356830251844953406051866432102826 -0.65187678415125922476703025933216051413 -0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + -0.65187678415125922476703025933216051413 -0.65187678415125922476703025933216051413 -0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + -0.65187678415125922476703025933216051413 0.30375356830251844953406051866432102826 0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + 0.30375356830251844953406051866432102826 -0.65187678415125922476703025933216051413 0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + -0.65187678415125922476703025933216051413 -0.65187678415125922476703025933216051413 0.40608422935459024220472507696440370377 0.16251705063595730345650668876189645828 + -0.90552242832046124932977560924577476395 0.8110448566409224986595512184915495279 -0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + 0.8110448566409224986595512184915495279 -0.90552242832046124932977560924577476395 -0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + -0.90552242832046124932977560924577476395 -0.90552242832046124932977560924577476395 -0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + -0.90552242832046124932977560924577476395 0.8110448566409224986595512184915495279 0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + 0.8110448566409224986595512184915495279 -0.90552242832046124932977560924577476395 0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + -0.90552242832046124932977560924577476395 -0.90552242832046124932977560924577476395 0.82881544305868020493486275746014093078 0.028581431014154463281085418946775972209 + -0.68050147211482206696479393170526142918 0.36100294422964413392958786341052285835 -0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + 0.36100294422964413392958786341052285835 -0.68050147211482206696479393170526142918 -0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + -0.68050147211482206696479393170526142918 -0.68050147211482206696479393170526142918 -0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + -0.68050147211482206696479393170526142918 0.36100294422964413392958786341052285835 0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + 0.36100294422964413392958786341052285835 -0.68050147211482206696479393170526142918 0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + -0.68050147211482206696479393170526142918 -0.68050147211482206696479393170526142918 0.96586516654065435448682371732369762245 0.044572109393861233332081245114974812183 + -0.4743723986175179830847346899909858658 0.45719614360199995718682043629988661246 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + 0.45719614360199995718682043629988661246 -0.4743723986175179830847346899909858658 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.98282374498448197410208574630890074666 0.45719614360199995718682043629988661246 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + 0.45719614360199995718682043629988661246 -0.98282374498448197410208574630890074666 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.98282374498448197410208574630890074666 -0.4743723986175179830847346899909858658 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.4743723986175179830847346899909858658 -0.98282374498448197410208574630890074666 -0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.4743723986175179830847346899909858658 0.45719614360199995718682043629988661246 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + 0.45719614360199995718682043629988661246 -0.4743723986175179830847346899909858658 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.98282374498448197410208574630890074666 0.45719614360199995718682043629988661246 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + 0.45719614360199995718682043629988661246 -0.98282374498448197410208574630890074666 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.98282374498448197410208574630890074666 -0.4743723986175179830847346899909858658 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + -0.4743723986175179830847346899909858658 -0.98282374498448197410208574630890074666 0.57735026918962576450914878050195745565 0.054590116363492292384362848950462321489 + ] + elseif n == 9 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0 0.19114798844620108508537046584425003736 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.87583915879987761660772171906042408394 0.1060969195893273614653564550346218398 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.87583915879987761660772171906042408394 0.1060969195893273614653564550346218398 + -0.0010699296586281549894456411507705393926 -0.99786014068274369002110871769845892121 0 0.049462405441895223670279761487537762472 + -0.99786014068274369002110871769845892121 -0.0010699296586281549894456411507705393926 0 0.049462405441895223670279761487537762472 + -0.0010699296586281549894456411507705393926 -0.0010699296586281549894456411507705393926 0 0.049462405441895223670279761487537762472 + -0.10986499991499267966311447997302837166 -0.78027000017001464067377104005394325668 -0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.78027000017001464067377104005394325668 -0.10986499991499267966311447997302837166 -0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.10986499991499267966311447997302837166 -0.10986499991499267966311447997302837166 -0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.10986499991499267966311447997302837166 -0.78027000017001464067377104005394325668 0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.78027000017001464067377104005394325668 -0.10986499991499267966311447997302837166 0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.10986499991499267966311447997302837166 -0.10986499991499267966311447997302837166 0.52963355520607675696300739352733925666 0.15080793261275860374457825128943921743 + -0.57769113257897895683607913964943464903 0.15538226515795791367215827929886929806 -0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + 0.15538226515795791367215827929886929806 -0.57769113257897895683607913964943464903 -0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + -0.57769113257897895683607913964943464903 -0.57769113257897895683607913964943464903 -0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + -0.57769113257897895683607913964943464903 0.15538226515795791367215827929886929806 0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + 0.15538226515795791367215827929886929806 -0.57769113257897895683607913964943464903 0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + -0.57769113257897895683607913964943464903 -0.57769113257897895683607913964943464903 0.32397688216877005165250824998597314548 0.09338246974553840864616210022748459774 + -0.91231503711701435232450671034134809014 0.82463007423402870464901342068269618029 -0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + 0.82463007423402870464901342068269618029 -0.91231503711701435232450671034134809014 -0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + -0.91231503711701435232450671034134809014 -0.91231503711701435232450671034134809014 -0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + -0.91231503711701435232450671034134809014 0.82463007423402870464901342068269618029 0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + 0.82463007423402870464901342068269618029 -0.91231503711701435232450671034134809014 0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + -0.91231503711701435232450671034134809014 -0.91231503711701435232450671034134809014 0.40165166857078872104499428459378049449 0.03838828638520642757038156568015552951 + -0.03450013665826621486161695118676603904 -0.93099972668346757027676609762646792192 -0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.93099972668346757027676609762646792192 -0.03450013665826621486161695118676603904 -0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.03450013665826621486161695118676603904 -0.03450013665826621486161695118676603904 -0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.03450013665826621486161695118676603904 -0.93099972668346757027676609762646792192 0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.93099972668346757027676609762646792192 -0.03450013665826621486161695118676603904 0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.03450013665826621486161695118676603904 -0.03450013665826621486161695118676603904 0.9677921254642105035410819146400981344 0.026528431202320433845464036710085263608 + -0.6456435117766454460710959796587760834 0.29128702355329089214219195931755216681 -0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + 0.29128702355329089214219195931755216681 -0.6456435117766454460710959796587760834 -0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + -0.6456435117766454460710959796587760834 -0.6456435117766454460710959796587760834 -0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + -0.6456435117766454460710959796587760834 0.29128702355329089214219195931755216681 0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + 0.29128702355329089214219195931755216681 -0.6456435117766454460710959796587760834 0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + -0.6456435117766454460710959796587760834 -0.6456435117766454460710959796587760834 0.89078174045697066578383528302064143989 0.066719942573448407007300296577585702464 + -0.90753763917811234690103494985844043039 0.81507527835622469380206989971688086078 -0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + 0.81507527835622469380206989971688086078 -0.90753763917811234690103494985844043039 -0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + -0.90753763917811234690103494985844043039 -0.90753763917811234690103494985844043039 -0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + -0.90753763917811234690103494985844043039 0.81507527835622469380206989971688086078 0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + 0.81507527835622469380206989971688086078 -0.90753763917811234690103494985844043039 0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + -0.90753763917811234690103494985844043039 -0.90753763917811234690103494985844043039 0.9577062038929269606887790227035291939 0.012131173836407148944270419030858660558 + -0.88126598051326876444819587362986475675 0.44939053910599895722040631866985064464 0 0.087691161312839084087580790508946261895 + 0.44939053910599895722040631866985064464 -0.88126598051326876444819587362986475675 0 0.087691161312839084087580790508946261895 + -0.56812455859273019277221044503998588789 0.44939053910599895722040631866985064464 0 0.087691161312839084087580790508946261895 + 0.44939053910599895722040631866985064464 -0.56812455859273019277221044503998588789 0 0.087691161312839084087580790508946261895 + -0.56812455859273019277221044503998588789 -0.88126598051326876444819587362986475675 0 0.087691161312839084087580790508946261895 + -0.88126598051326876444819587362986475675 -0.56812455859273019277221044503998588789 0 0.087691161312839084087580790508946261895 + -0.5517449084932878730050536869300298641 0.49681091913207739922639584227322548803 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + 0.49681091913207739922639584227322548803 -0.5517449084932878730050536869300298641 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.94506601063878952622134215534319562393 0.49681091913207739922639584227322548803 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + 0.49681091913207739922639584227322548803 -0.94506601063878952622134215534319562393 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.94506601063878952622134215534319562393 -0.5517449084932878730050536869300298641 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.5517449084932878730050536869300298641 -0.94506601063878952622134215534319562393 -0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.5517449084932878730050536869300298641 0.49681091913207739922639584227322548803 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + 0.49681091913207739922639584227322548803 -0.5517449084932878730050536869300298641 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.94506601063878952622134215534319562393 0.49681091913207739922639584227322548803 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + 0.49681091913207739922639584227322548803 -0.94506601063878952622134215534319562393 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.94506601063878952622134215534319562393 -0.5517449084932878730050536869300298641 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + -0.5517449084932878730050536869300298641 -0.94506601063878952622134215534319562393 0.69522055387163478698981196744187399375 0.0495312141698622864915543816230467997 + ] + elseif n == 10 + xw = [ + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -1 0.038043108084240763009572271119163370293 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 1 0.038043108084240763009572271119163370293 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 -0.2838364013985715684657060577249329797 0.12804920643255787761445727849451714328 + -0.33333333333333333333333333333333333333 -0.33333333333333333333333333333333333333 0.2838364013985715684657060577249329797 0.12804920643255787761445727849451714328 + -0.58509925002156576804484663776694671178 0.17019850004313153608969327553389342357 0 0.045423158480136973856219994221831113886 + 0.17019850004313153608969327553389342357 -0.58509925002156576804484663776694671178 0 0.045423158480136973856219994221831113886 + -0.58509925002156576804484663776694671178 -0.58509925002156576804484663776694671178 0 0.045423158480136973856219994221831113886 + -0.087745274263836468012344092631834458624 -0.82450945147232706397531181473633108275 0 0.12078661810034629566292168130821741978 + -0.82450945147232706397531181473633108275 -0.087745274263836468012344092631834458624 0 0.12078661810034629566292168130821741978 + -0.087745274263836468012344092631834458624 -0.087745274263836468012344092631834458624 0 0.12078661810034629566292168130821741978 + -0.83336149569989450879622615827585222698 0.66672299139978901759245231655170445397 0 0.070523638001611290833549221721909796384 + 0.66672299139978901759245231655170445397 -0.83336149569989450879622615827585222698 0 0.070523638001611290833549221721909796384 + -0.83336149569989450879622615827585222698 -0.83336149569989450879622615827585222698 0 0.070523638001611290833549221721909796384 + -0.58566729572033783230931886738991607412 0.17133459144067566461863773477983214824 -0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + 0.17133459144067566461863773477983214824 -0.58566729572033783230931886738991607412 -0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + -0.58566729572033783230931886738991607412 -0.58566729572033783230931886738991607412 -0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + -0.58566729572033783230931886738991607412 0.17133459144067566461863773477983214824 0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + 0.17133459144067566461863773477983214824 -0.58566729572033783230931886738991607412 0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + -0.58566729572033783230931886738991607412 -0.58566729572033783230931886738991607412 0.47829026090036813970247985116366265217 0.10435318130434033235265655019864541098 + -0.71401760691431835996629386931253079667 0.42803521382863671993258773862506159333 -0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + 0.42803521382863671993258773862506159333 -0.71401760691431835996629386931253079667 -0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + -0.71401760691431835996629386931253079667 -0.71401760691431835996629386931253079667 -0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + -0.71401760691431835996629386931253079667 0.42803521382863671993258773862506159333 0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + 0.42803521382863671993258773862506159333 -0.71401760691431835996629386931253079667 0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + -0.71401760691431835996629386931253079667 -0.71401760691431835996629386931253079667 0.8520666639991396363392553699935906628 0.045353049770314160210925207718221960406 + -0.0059182255205082649355625415633530872236 -0.98816354895898347012887491687329382555 -0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.98816354895898347012887491687329382555 -0.0059182255205082649355625415633530872236 -0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.0059182255205082649355625415633530872236 -0.0059182255205082649355625415633530872236 -0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.0059182255205082649355625415633530872236 -0.98816354895898347012887491687329382555 0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.98816354895898347012887491687329382555 -0.0059182255205082649355625415633530872236 0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.0059182255205082649355625415633530872236 -0.0059182255205082649355625415633530872236 0.62749011164795345095028595959683589121 0.028126573654953943549801873601641305154 + -0.14443508029301579714242731927948333445 -0.71112983941396840571514536144103333111 -0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.71112983941396840571514536144103333111 -0.14443508029301579714242731927948333445 -0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.14443508029301579714242731927948333445 -0.14443508029301579714242731927948333445 -0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.14443508029301579714242731927948333445 -0.71112983941396840571514536144103333111 0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.71112983941396840571514536144103333111 -0.14443508029301579714242731927948333445 0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.14443508029301579714242731927948333445 -0.14443508029301579714242731927948333445 0.76897188255294104672722599221030791768 0.09730006751474251906683045094709204535 + -0.92655062106089238633040808286062465799 0.85310124212178477266081616572124931598 -0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + 0.85310124212178477266081616572124931598 -0.92655062106089238633040808286062465799 -0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + -0.92655062106089238633040808286062465799 -0.92655062106089238633040808286062465799 -0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + -0.92655062106089238633040808286062465799 0.85310124212178477266081616572124931598 0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + 0.85310124212178477266081616572124931598 -0.92655062106089238633040808286062465799 0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + -0.92655062106089238633040808286062465799 -0.92655062106089238633040808286062465799 0.94536185592175700771401046143327478453 0.0075697428403356175236518069263715176348 + -0.99621476080190120325860648054065438352 0.99582008190465265332555390082623777504 0 0.0031208636401827010315915455956052660089 + 0.99582008190465265332555390082623777504 -0.99621476080190120325860648054065438352 0 0.0031208636401827010315915455956052660089 + -0.99960532110275145006694742028558339152 0.99582008190465265332555390082623777504 0 0.0031208636401827010315915455956052660089 + 0.99582008190465265332555390082623777504 -0.99960532110275145006694742028558339152 0 0.0031208636401827010315915455956052660089 + -0.99960532110275145006694742028558339152 -0.99621476080190120325860648054065438352 0 0.0031208636401827010315915455956052660089 + -0.99621476080190120325860648054065438352 -0.99960532110275145006694742028558339152 0 0.0031208636401827010315915455956052660089 + -0.91342063244101618001407082339238449491 0.3255064864791910018042153093720481789 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + 0.3255064864791910018042153093720481789 -0.91342063244101618001407082339238449491 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.41208585403817482179014448597966368399 0.3255064864791910018042153093720481789 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + 0.3255064864791910018042153093720481789 -0.41208585403817482179014448597966368399 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.41208585403817482179014448597966368399 -0.91342063244101618001407082339238449491 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.91342063244101618001407082339238449491 -0.41208585403817482179014448597966368399 -0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.91342063244101618001407082339238449491 0.3255064864791910018042153093720481789 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + 0.3255064864791910018042153093720481789 -0.91342063244101618001407082339238449491 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.41208585403817482179014448597966368399 0.3255064864791910018042153093720481789 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + 0.3255064864791910018042153093720481789 -0.41208585403817482179014448597966368399 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.41208585403817482179014448597966368399 -0.91342063244101618001407082339238449491 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.91342063244101618001407082339238449491 -0.41208585403817482179014448597966368399 0.94236346745083749482627300953387466221 0.024183540649806648348002441154357646858 + -0.9234169616555928413056660388640430162 0.37513146835671456389900781390965761438 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + 0.37513146835671456389900781390965761438 -0.9234169616555928413056660388640430162 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.45171450670112172259334177504561459817 0.37513146835671456389900781390965761438 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + 0.37513146835671456389900781390965761438 -0.45171450670112172259334177504561459817 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.45171450670112172259334177504561459817 -0.9234169616555928413056660388640430162 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.9234169616555928413056660388640430162 -0.45171450670112172259334177504561459817 -0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.9234169616555928413056660388640430162 0.37513146835671456389900781390965761438 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + 0.37513146835671456389900781390965761438 -0.9234169616555928413056660388640430162 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.45171450670112172259334177504561459817 0.37513146835671456389900781390965761438 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + 0.37513146835671456389900781390965761438 -0.45171450670112172259334177504561459817 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.45171450670112172259334177504561459817 -0.9234169616555928413056660388640430162 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.9234169616555928413056660388640430162 -0.45171450670112172259334177504561459817 0.33533610866739435446116498131265141599 0.05460310873669147926531122821390556603 + -0.96792121274553310722108842071500064645 0.74052721727873797554079834546410281777 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + 0.74052721727873797554079834546410281777 -0.96792121274553310722108842071500064645 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.77260600453320486831970992474910217133 0.74052721727873797554079834546410281777 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + 0.74052721727873797554079834546410281777 -0.77260600453320486831970992474910217133 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.77260600453320486831970992474910217133 -0.96792121274553310722108842071500064645 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.96792121274553310722108842071500064645 -0.77260600453320486831970992474910217133 -0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.96792121274553310722108842071500064645 0.74052721727873797554079834546410281777 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + 0.74052721727873797554079834546410281777 -0.96792121274553310722108842071500064645 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.77260600453320486831970992474910217133 0.74052721727873797554079834546410281777 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + 0.74052721727873797554079834546410281777 -0.77260600453320486831970992474910217133 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.77260600453320486831970992474910217133 -0.96792121274553310722108842071500064645 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + -0.96792121274553310722108842071500064645 -0.77260600453320486831970992474910217133 0.61432350187707741015732547519719343903 0.024769538519410488660113297222678366236 + ] + else + throw(ArgumentError("unsupported order for prism polyquad integration")) + end + + # Transform from [-1,1] × [-1,1] × [-1,1] to [0,1] × [0,1] × [0,1] + #x + xw[:,1] .+= 1.0 + xw[:,1] ./= 2.0 + #y + xw[:,2] .+= 1.0 + xw[:,2] ./= 2.0 + #z + xw[:,3] .+= 1.0 + xw[:,3] ./= 2.0 + #w + xw[:,4] ./= 2^3 + + @show xw[:,3] + + return xw +end diff --git a/src/Quadrature/quadrature.jl b/src/Quadrature/quadrature.jl index a2c0de7375..69d9dbbbf3 100644 --- a/src/Quadrature/quadrature.jl +++ b/src/Quadrature/quadrature.jl @@ -1,5 +1,6 @@ include("gaussquad_tri_table.jl") include("gaussquad_tet_table.jl") +include("gaussquad_prism_table.jl") include("generate_quadrature.jl") import Base.Cartesian: @nloops, @nref, @ntuple, @nexprs @@ -85,6 +86,7 @@ julia> getpoints(qr) getpoints(qr::QuadratureRule) = qr.points QuadratureRule{dim,shape}(order::Int) where {dim,shape} = QuadratureRule{dim,shape}(:legendre, order) +QuadratureRule{3,RefPrism}(order::Int) = QuadratureRule{3,RefPrism}(:polyquad, order) # No legendre rule available yet for prism... # Special case for face integration of 1D problems function (::Type{QuadratureRule{0, RefCube}})(quad_type::Symbol, order::Int) @@ -162,3 +164,20 @@ function (::Type{QuadratureRule{1,RefTetrahedron}})(quad_type::Symbol, order::In end return QuadratureRule{1,RefTetrahedron,Float64}(weights, points) end + +# Grab prism quadrature rule from table +function (::Type{QuadratureRule{3, RefPrism}})(quad_type::Symbol, order::Int) + if quad_type == :polyquad + data = _get_gauss_prismdata_polyquad(order) + else + throw(ArgumentError("unsupported quadrature rule")) + end + n_points = size(data,1) + points = Vector{Vec{3,Float64}}(undef, n_points) + + for p in 1:size(data, 1) + points[p] = Vec{3,Float64}(@ntuple 3 i -> data[p, i]) + end + weights = data[:, 4] + QuadratureRule{3,RefPrism,Float64}(weights, points) +end diff --git a/src/deprecations.jl b/src/deprecations.jl index da4358d54b..11b2a1436d 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -1,4 +1,9 @@ Base.@deprecate_binding DirichletBoundaryConditions ConstraintHandler Base.@deprecate_binding DirichletBoundaryCondition Dirichlet + import Base: push! @deprecate push!(dh::AbstractDofHandler, args...) add!(dh, args...) + +@deprecate vertices(ip::Interpolation) vertexdof_indices(ip) false +@deprecate faces(ip::Interpolation) facedof_indices(ip) false +@deprecate edges(ip::Interpolation) edgedof_indices(ip) false diff --git a/src/exports.jl b/src/exports.jl index 83e3d9d6ad..b735672237 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -3,6 +3,7 @@ export Interpolation, RefCube, RefTetrahedron, + RefPrism, BubbleEnrichedLagrange, CrouzeixRaviart, Lagrange, @@ -55,6 +56,7 @@ export QuadraticTetrahedron, Hexahedron, #QuadraticHexahedron, + Wedge, CellIndex, FaceIndex, EdgeIndex, diff --git a/src/interpolations.jl b/src/interpolations.jl index 8a10b7c301..6082a37fbf 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -24,6 +24,8 @@ The following interpolations are implemented: * `Lagrange{3,RefCube,2}` * `Lagrange{3,RefTetrahedron,1}` * `Lagrange{3,RefTetrahedron,2}` +* `Lagrange{3,RefPrism,1}` +* `Lagrange{3,RefPrism,2}` * `Serendipity{2,RefCube,2}` * `Serendipity{3,RefCube,2}` @@ -38,6 +40,72 @@ julia> getnbasefunctions(ip) """ abstract type Interpolation{dim,shape,order} end +# TODO: Remove: this is a hotfix to apply constraints to embedded elements. +edges(ip::Interpolation{2}) = faces(ip) +edgedof_indices(ip::Interpolation{2}) = facedof_indices(ip) +edgedof_interior_indices(ip::Interpolation{2}) = facedof_interior_indices(ip) +facedof_indices(ip::Interpolation{1}) = vertexdof_indices(ip) + +""" + InterpolationInfo + +Gathers all the information needed to distribute dofs for a given interpolation. Note that +this cache is of the same type no matter the interpolation: the purpose is to make +dof-distribution type-stable. +""" +struct InterpolationInfo + nvertexdofs::Vector{Int} + nedgedofs::Vector{Int} + nfacedofs::Vector{Int} + ncelldofs::Int + dim::Int + function InterpolationInfo(interpolation::Interpolation{3}) + new( + [length(i) for i ∈ vertexdof_indices(interpolation)], + [length(i) for i ∈ edgedof_interior_indices(interpolation)], + [length(i) for i ∈ facedof_interior_indices(interpolation)], + length(celldof_interior_indices(interpolation)), + 3, + ) + end + function InterpolationInfo(interpolation::Interpolation{2}) + new( + [length(i) for i ∈ vertexdof_indices(interpolation)], + Int[], + [length(i) for i ∈ facedof_interior_indices(interpolation)], + length(celldof_interior_indices(interpolation)), + 2, + ) + end + function InterpolationInfo(interpolation::Interpolation{1}) + new( + [length(i) for i ∈ vertexdof_indices(interpolation)], + Int[], + [0, 0], # Well. Apparently vertices are also faces. :) + length(celldof_interior_indices(interpolation)), + 1, + ) + end +end + +# Some redundant information about the geometry of the reference cells. +nfaces(::Interpolation{dim, RefCube}) where {dim}= 2*dim +nfaces(::Interpolation{2, RefTetrahedron}) = 3 +nfaces(::Interpolation{3, RefTetrahedron}) = 4 +nfaces(::Interpolation{3, RefPrism}) = 5 + +nedges(::Interpolation{1, RefCube}) = 0 +nedges(::Interpolation{2, RefCube}) = 0 +nedges(::Interpolation{3, RefCube}) = 12 +nedges(::Interpolation{2, RefTetrahedron}) = 0 +nedges(::Interpolation{3, RefTetrahedron}) = 6 +nedges(::Interpolation{3, RefPrism}) = 9 + +nvertices(::Interpolation{dim, RefCube}) where {dim} = 2^dim +nvertices(::Interpolation{2, RefTetrahedron}) = 3 +nvertices(::Interpolation{3, RefTetrahedron}) = 4 +nvertices(::Interpolation{3, RefPrism}) = 6 + Base.copy(ip::Interpolation) = ip """ @@ -94,52 +162,12 @@ Return the number of base functions for the interpolation `ip`. """ getnbasefunctions(::Interpolation) -# struct that gathers all the information needed to distribute -# dofs for a given interpolation. -struct InterpolationInfo - # TODO: Can be smaller than `Int` if that matters... - nvertexdofs::Int - nedgedofs::Int - nfacedofs::Int - ncelldofs::Int - dim::Int - InterpolationInfo(interpolation::Interpolation{dim}) where {dim} = - new(nvertexdofs(interpolation), nedgedofs(interpolation), - nfacedofs(interpolation), ncelldofs(interpolation), dim) -end - # The following functions are used to distribute the dofs. Definitions: # vertexdof: dof on a "corner" of the reference shape # facedof: dof in the dim-1 dimension (line in 2D, surface in 3D) # edgedof: dof on a line between 2 vertices (i.e. "corners") (3D only) # celldof: dof that is local to the element -# Necessary for correct distribution of dofs. -""" - Ferrite.nvertexdofs(::Interpolation) - -Number of dofs per vertex. -""" -nvertexdofs(::Interpolation) = 0 -""" - Ferrite.nedgedofs(::Interpolation) - -Number of dofs on the interior of a single edge. -""" -nedgedofs(::Interpolation) = 0 -""" - Ferrite.nfacedofs(::Interpolation) - -Number of dofs on the interior of a single face. -""" -nfacedofs(::Interpolation) = 0 -""" - Ferrite.ncelldofs(::Interpolation) - -Total number of dofs in the interior of the element. -""" -ncelldofs(::Interpolation) = 0 - """ value(ip::Interpolation, i::Int, ξ::Vec) @@ -154,9 +182,9 @@ indices of [`reference_coordinates(::Interpolation)`](@ref). value(ip::Interpolation, i::Int, ξ::Vec) """ - reference_coordinates(::Interpolation) + reference_coordinates(ip::Interpolation) -Returns a vector of coordinates with length [`getnbasefunctions(::Interpolation)`](@ref) +Returns a vector of coordinates with length [`getnbasefunctions(::Interpolation)`](@ref) and indices corresponding to the indices of a dof in [`vertices`](@ref), [`faces`](@ref) and [`edges`](@ref). @@ -164,41 +192,78 @@ and indices corresponding to the indices of a dof in [`vertices`](@ref), [`faces TODO: Separate nodal and non-nodal interpolations. """ -reference_coordinates(ip::Interpolation) +reference_coordinates(::Interpolation) """ - vertices(::Interpolation) + vertexdof_indices(ip::Interpolation) -A tuple containing tuples of local dof indices for the respective -vertex in local enumeration on a cell defined by [`vertices(::Cell)`](@ref). The vertex enumeration must +A tuple containing tuples of local dof indices for the respective vertex in local +enumeration on a cell defined by [`vertices(::Cell)`](@ref). The vertex enumeration must match the vertex enumeration of the corresponding geometrical cell. -The length of each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref). """ -vertices(::Interpolation) +vertexdof_indices(ip::Interpolation) = ntuple(_ -> (), nvertices(ip)) + +""" + edgedof_indices(ip::Interpolation) + +A tuple containing tuples of local dof indices for the respective edge in local enumeration +on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must match the edge +enumeration of the corresponding geometrical cell. """ - edges(::Interpolation) +edgedof_indices(::Interpolation) -A tuple containing tuples of local dof indices for the respective -edge in local enumeration on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must -match the edge enumeration of the corresponding geometrical cell. -Note that the vertex dofs are included here. The length of -each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref)+[`nedgedofs(::Interpolation)`](@ref). """ -edges(::Interpolation) + edgedof_interior_indices(ip::Interpolation) + +A tuple containing tuples of the local dof indices on the interior of the respective edge in +local enumeration on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must +match the edge enumeration of the corresponding geometrical cell. Note that the vertex dofs +are included here. """ - faces(::Interpolation) +edgedof_interior_indices(::Interpolation) -A tuple containing tuples of all local dof indices for the respective -face in local enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must -match the face enumeration of the corresponding geometrical cell. -Note that the vertex and edge dofs are included here. The -length of each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref)+[`nedgedofs(::Interpolation)`](@ref)+[`nfacedofs(::Interpolation)`](@ref). """ -faces(::Interpolation) + facedof_indices(ip::Interpolation) -# Needed for distributing dofs on shells correctly (face in 2d is edge in 3d) -edges(ip::Interpolation{2}) = faces(ip) -nedgedofs(ip::Interpolation{2}) = nfacedofs(ip) +A tuple containing tuples of all local dof indices for the respective face in local +enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must match +the face enumeration of the corresponding geometrical cell. +""" +facedof_indices(::Interpolation) + +""" + facedof_interior_indices(ip::Interpolation) + +A tuple containing tuples of the local dof indices on the interior of the respective face in +local enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must +match the face enumeration of the corresponding geometrical cell. Note that the vertex and +edge dofs are included here. +""" +facedof_interior_indices(::Interpolation) + +""" + celldof_interior_indices(ip::Interpolation) + +Tuple containing the dof indices associated with the interior of the cell. +""" +celldof_interior_indices(::Interpolation) = () + +# Some helpers to skip boilerplate +edgedof_indices(ip::Interpolation{3}) = ntuple(_ -> (), nedges(ip)) +edgedof_interior_indices(ip::Interpolation{3}) = ntuple(_ -> (), nedges(ip)) +facedof_indices(ip::Union{Interpolation{2}, Interpolation{3}}) = ntuple(_ -> (), nfaces(ip)) +facedof_interior_indices(ip::Union{Interpolation{2}, Interpolation{3}}) = ntuple(_ -> (), nfaces(ip)) + +""" + boundarydof_indices(::Type{<:BoundaryIndex}) + +Helper function to generically dispatch on the correct dof sets of a boundary entity. +""" +boundarydof_indices(::Type{<:BoundaryIndex}) + +boundarydof_indices(::Type{FaceIndex}) = Ferrite.facedof_indices +boundarydof_indices(::Type{EdgeIndex}) = Ferrite.edgedof_indices +boundarydof_indices(::Type{VertexIndex}) = Ferrite.vertexdof_indices ######################### # DiscontinuousLagrange # @@ -212,17 +277,11 @@ struct DiscontinuousLagrange{dim,shape,order} <: Interpolation{dim,shape,order} getlowerdim(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = DiscontinuousLagrange{dim-1,shape,order}() getlowerorder(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = DiscontinuousLagrange{dim,shape,order-1}() -# TODO generalize properly -# ncelldofs(::DiscontinuousLagrange{dim,RefCube,order}) where {dim,order} = (order+1)^dim -# ncelldofs(::DiscontinuousLagrange{2,RefTetrahedron,order}) where {order} = (order+1)*(order+2)/2 -# ncelldofs(::DiscontinuousLagrange{3,RefTetrahedron,order}) where {order} = (order+1)*(order+2)/2 getnbasefunctions(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = getnbasefunctions(Lagrange{dim,shape,order}()) -ncelldofs(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = getnbasefunctions(DiscontinuousLagrange{dim,shape,order}()) - getnbasefunctions(::DiscontinuousLagrange{dim,shape,0}) where {dim,shape} = 1 -ncelldofs(::DiscontinuousLagrange{dim,shape,0}) where {dim,shape} = 1 -faces(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = () +# This just moves all dofs into the interior of the element. +celldof_interior_indices(ip::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = ntuple(i->i, getnbasefunctions(ip)) # Mirror the Lagrange element for now. function reference_coordinates(ip::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} @@ -246,6 +305,7 @@ function reference_coordinates(ip::DiscontinuousLagrange{3,RefTetrahedron,0}) end function value(ip::DiscontinuousLagrange{dim,shape,0}, i::Int, ξ::Vec{dim}) where {dim,shape} + i > 1 && throw(ArgumentError("no shape function $i for interpolation $ip")) return 1.0 end @@ -255,11 +315,12 @@ end struct Lagrange{dim,shape,order} <: Interpolation{dim,shape,order} end # Vertices for all Lagrange interpolations are the same -vertices(::Lagrange{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,)) -vertices(::Lagrange{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,)) -vertices(::Lagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,)) -vertices(::Lagrange{3,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,),(4,)) -nvertexdofs(::Lagrange{dim,shape,order}) where {dim,shape,order} = 1 +vertexdof_indices(::Lagrange{1,RefCube,order}) where {order} = ((1,),(2,)) +vertexdof_indices(::Lagrange{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,)) +vertexdof_indices(::Lagrange{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,)) +vertexdof_indices(::Lagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,)) +vertexdof_indices(::Lagrange{3,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,),(4,)) +vertexdof_indices(::Lagrange{3,RefPrism,order}) where {order} = ((1,), (2,), (3,), (4,), (5,), (6,)) getlowerdim(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim-1,shape,order}() getlowerorder(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim,shape,order-1}() @@ -270,8 +331,6 @@ getlowerorder(::Lagrange{dim,shape,1}) where {dim,shape} = DiscontinuousLagrange ################################## getnbasefunctions(::Lagrange{1,RefCube,1}) = 2 -faces(::Lagrange{1,RefCube,1}) = ((1,), (2,)) - function reference_coordinates(::Lagrange{1,RefCube,1}) return [Vec{1, Float64}((-1.0,)), Vec{1, Float64}(( 1.0,))] @@ -288,9 +347,9 @@ end # Lagrange dim 1 RefCube order 2 # ################################## getnbasefunctions(::Lagrange{1,RefCube,2}) = 3 -ncelldofs(::Lagrange{1,RefCube,2}) = 1 -faces(::Lagrange{1,RefCube,2}) = ((1,), (2,)) +facedof_indices(::Lagrange{1,RefCube,2}) = ((1,), (2,)) +celldof_interior_indices(::Lagrange{1,RefCube,2}) = (3,) function reference_coordinates(::Lagrange{1,RefCube,2}) return [Vec{1, Float64}((-1.0,)), @@ -311,7 +370,7 @@ end ################################## getnbasefunctions(::Lagrange{2,RefCube,1}) = 4 -faces(::Lagrange{2,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1)) +facedof_indices(::Lagrange{2,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1)) function reference_coordinates(::Lagrange{2,RefCube,1}) return [Vec{2, Float64}((-1.0, -1.0)), @@ -334,10 +393,10 @@ end # Lagrange dim 2 RefCube order 2 # ################################## getnbasefunctions(::Lagrange{2,RefCube,2}) = 9 -nfacedofs(::Lagrange{2,RefCube,2}) = 1 -ncelldofs(::Lagrange{2,RefCube,2}) = 1 -faces(::Lagrange{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8)) +facedof_indices(::Lagrange{2,RefCube,2}) = ((1,2, 5), (2,3, 6), (3,4, 7), (4,1, 8)) +facedof_interior_indices(::Lagrange{2,RefCube,2}) = ((5,), (6,), (7,), (8,)) +celldof_interior_indices(::Lagrange{2,RefCube,2}) = (9,) function reference_coordinates(::Lagrange{2,RefCube,2}) return [Vec{2, Float64}((-1.0, -1.0)), @@ -372,8 +431,7 @@ end getnbasefunctions(::Lagrange{2,RefTetrahedron,1}) = 3 getlowerdim(::Lagrange{2, RefTetrahedron, order}) where {order} = Lagrange{1, RefCube, order}() -vertices(::Lagrange{2,RefTetrahedron,1}) = (1,2,3) -faces(::Lagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1)) +facedof_indices(::Lagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1)) function reference_coordinates(::Lagrange{2,RefTetrahedron,1}) return [Vec{2, Float64}((1.0, 0.0)), @@ -394,10 +452,9 @@ end # Lagrange dim 2 RefTetrahedron order 2 # ######################################### getnbasefunctions(::Lagrange{2,RefTetrahedron,2}) = 6 -nfacedofs(::Lagrange{2,RefTetrahedron,2}) = 1 -vertices(::Lagrange{2,RefTetrahedron,2}) = (1,2,3) -faces(::Lagrange{2,RefTetrahedron,2}) = ((1,2,4), (2,3,5), (3,1,6)) +facedof_indices(::Lagrange{2,RefTetrahedron,2}) = ((1,2,4), (2,3,5), (3,1,6)) +facedof_interior_indices(::Lagrange{2,RefTetrahedron,2}) = ((4,), (5,), (6,)) function reference_coordinates(::Lagrange{2,RefTetrahedron,2}) return [Vec{2, Float64}((1.0, 0.0)), @@ -436,16 +493,8 @@ function getnbasefunctions(ip::Lagrange2Tri345) return (order + 1) * (order + 2) ÷ 2 end -nvertexdofs(::Lagrange2Tri345) = 1 -nedgedofs(ip::Lagrange2Tri345) = getorder(ip) - 1 -nfacedofs(ip::Lagrange2Tri345) = getorder(ip) - 1 -function ncelldofs(ip::Lagrange2Tri345) - order = getorder(ip) - return (order + 1) * (order + 2) ÷ 2 - 3 * order -end - # Permutation to switch numbering to Ferrite ordering -const permdof2D = Dict{Int,Vector{Int}}( +const permdof2DLagrange2Tri345 = Dict{Int,Vector{Int}}( 1 => [1, 2, 3], 2 => [3, 6, 1, 5, 4, 2], 3 => [4, 10, 1, 7, 9, 8, 5, 2, 3, 6], @@ -453,23 +502,35 @@ const permdof2D = Dict{Int,Vector{Int}}( 5 => [6, 21, 1, 11, 15, 18, 20, 19, 16, 12, 7, 2, 3, 4, 5, 8, 9, 10, 13, 14, 17], ) -vertices(::Lagrange2Tri345) = (1, 2, 3) +vertexdof_indices(::Lagrange2Tri345) = (1, 2, 3) -function faces(ip::Lagrange2Tri345) +function facedof_indices(ip::Lagrange2Tri345) order = getorder(ip) - if order == 1 - return ((1,2), (2,3), (3,1)) - elseif order == 2 - return ((1,2,4), (2,3,5), (3,1,6)) - elseif order == 3 - return ((1,2,4,5), (2,3,6,7), (3,1,8,9)) - elseif order == 4 - return ((1,2,4,5,6), (2,3,7,8,9), (3,1,10,11,12)) - elseif order == 5 - return ((1,2,4,5,6,7), (2,3,8,9,10,11), (3,1,12,13,14,15)) - else - error("unreachable") - end + order == 1 && return ((1,2), (2,3), (3,1)) + order == 2 && return ((1,2,4), (2,3,5), (3,1,6)) + order == 3 && return ((1,2,4,5), (2,3,6,7), (3,1,8,9)) + order == 4 && return ((1,2,4,5,6), (2,3,7,8,9), (3,1,10,11,12)) + order == 5 && return ((1,2,4,5,6,7), (2,3,8,9,10,11), (3,1,12,13,14,15)) + + throw(ArgumentError("Unsupported order $order for Lagrange on triangles.")) +end + +function facedof_interior_indices(ip::Lagrange2Tri345) + order = getorder(ip) + order == 1 && return ((), (), ()) + order == 2 && return ((4,), (5,), (6,)) + order == 3 && return ((4,5), (6,7), (8,9)) + order == 4 && return ((4,5,6), (7,8,9), (10,11,12)) + order == 5 && return ((4,5,6,7), (8,9,10,11), (12,13,14,15)) + + throw(ArgumentError("Unsupported order $order for Lagrange on triangles.")) +end + +function celldof_interior_indices(ip::Lagrange2Tri345) + order = getorder(ip) + ncellintdofs = (order + 1) * (order + 2) ÷ 2 - 3 * order + totaldofs = getnbasefunctions(ip) + return ntuple(i->totaldofs-ncellintdofs+i, ncellintdofs) end function reference_coordinates(ip::Lagrange2Tri345) @@ -480,12 +541,15 @@ function reference_coordinates(ip::Lagrange2Tri345) push!(coordpts, Vec{2, Float64}((l / order, k / order))) end end - return permute!(coordpts, permdof2D[order]) + return permute!(coordpts, permdof2DLagrange2Tri345[order]) end function value(ip::Lagrange2Tri345, i::Int, ξ::Vec{2}) + if !(0 < i <= getnbasefunctions(ip)) + throw(ArgumentError("no shape function $i for interpolation $ip")) + end order = getorder(ip) - i = permdof2D[order][i] + i = permdof2DLagrange2Tri345[order][i] ξ_x = ξ[1] ξ_y = ξ[2] γ = 1. - ξ_x - ξ_y @@ -517,8 +581,8 @@ end ######################################### getnbasefunctions(::Lagrange{3,RefTetrahedron,1}) = 4 -faces(::Lagrange{3,RefTetrahedron,1}) = ((1,3,2), (1,2,4), (2,3,4), (1,4,3)) -edges(::Lagrange{3,RefTetrahedron,1}) = ((1,2), (2,3), (3,1), (1,4), (2,4), (3,4)) +facedof_indices(::Lagrange{3,RefTetrahedron,1}) = ((1,3,2), (1,2,4), (2,3,4), (1,4,3)) +edgedof_indices(::Lagrange{3,RefTetrahedron,1}) = ((1,2), (2,3), (3,1), (1,4), (2,4), (3,4)) function reference_coordinates(::Lagrange{3,RefTetrahedron,1}) return [Vec{3, Float64}((0.0, 0.0, 0.0)), @@ -542,10 +606,10 @@ end # Lagrange dim 3 RefTetrahedron order 2 # ######################################### getnbasefunctions(::Lagrange{3,RefTetrahedron,2}) = 10 -nedgedofs(::Lagrange{3,RefTetrahedron,2}) = 1 -faces(::Lagrange{3,RefTetrahedron,2}) = ((1,3,2,7,6,5), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7)) -edges(::Lagrange{3,RefTetrahedron,2}) = ((1,2,5), (2,3,6), (3,1,7), (1,4,8), (2,4,9), (3,4,10)) +facedof_indices(::Lagrange{3,RefTetrahedron,2}) = ((1,3,2,7,6,5), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7)) +edgedof_indices(::Lagrange{3,RefTetrahedron,2}) = ((1,2,5), (2,3,6), (3,1,7), (1,4,8), (2,4,9), (3,4,10)) +edgedof_interior_indices(::Lagrange{3,RefTetrahedron,2}) = ((5,), (6,), (7,), (8,), (9,), (10,)) function reference_coordinates(::Lagrange{3,RefTetrahedron,2}) return [Vec{3, Float64}((0.0, 0.0, 0.0)), @@ -584,8 +648,8 @@ end ################################## getnbasefunctions(::Lagrange{3,RefCube,1}) = 8 -faces(::Lagrange{3,RefCube,1}) = ((1,4,3,2), (1,2,6,5), (2,3,7,6), (3,4,8,7), (1,5,8,4), (5,6,7,8)) -edges(::Lagrange{3,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5)) +facedof_indices(::Lagrange{3,RefCube,1}) = ((1,4,3,2), (1,2,6,5), (2,3,7,6), (3,4,8,7), (1,5,8,4), (5,6,7,8)) +edgedof_indices(::Lagrange{3,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5)) function reference_coordinates(::Lagrange{3,RefCube,1}) return [Vec{3, Float64}((-1.0, -1.0, -1.0)), @@ -619,11 +683,8 @@ end ################################## # Based on vtkTriQuadraticHexahedron (see https://kitware.github.io/vtk-examples/site/Cxx/GeometricObjects/IsoparametricCellsDemo/) getnbasefunctions(::Lagrange{3,RefCube,2}) = 27 -nedgedofs(::Lagrange{3,RefCube,2}) = 1 -nfacedofs(::Lagrange{3,RefCube,2}) = 1 -ncelldofs(::Lagrange{3,RefCube,2}) = 1 -faces(::Lagrange{3,RefCube,2}) = ( +facedof_indices(::Lagrange{3,RefCube,2}) = ( (1,4,3,2, 12,11,10,9, 21), (1,2,6,5, 9,18,13,17, 22), (2,3,7,6, 10,19,14,18, 23), @@ -631,7 +692,29 @@ faces(::Lagrange{3,RefCube,2}) = ( (1,5,8,4, 17,16,20,12, 25), (5,6,7,8, 13,14,15,16, 26), ) -edges(::Lagrange{3,RefCube,2}) = ((1,2, 9), (2,3, 10), (3,4, 11), (4,1, 12), (1,5, 17), (2,6, 18), (3,7, 19), (4,8, 20), (5,6, 13), (6,7, 14), (7,8, 15), (8,5, 16)) +facedof_interior_indices(::Lagrange{3,RefCube,2}) = ( + (21,), (22,), (23,), (24,), (25,), (26,), +) + +edgedof_indices(::Lagrange{3,RefCube,2}) = ( + (1,2, 9), + (2,3, 10), + (3,4, 11), + (4,1, 12), + (1,5, 17), + (2,6, 18), + (3,7, 19), + (4,8, 20), + (5,6, 13), + (6,7, 14), + (7,8, 15), + (8,5, 16) +) +edgedof_interior_indices(::Lagrange{3,RefCube,2}) = ( + (9,), (10,), (11,), (12,), (17), (18,), (19,), (20,), (13,), (14,), (15,), (16,) +) + +celldof_interior_indices(::Lagrange{3,RefCube,2}) = (27,) function reference_coordinates(::Lagrange{3,RefCube,2}) # vertex @@ -707,6 +790,131 @@ function value(ip::Lagrange{3,RefCube,2}, i::Int, ξ::Vec{3, T}) where {T} throw(ArgumentError("no shape function $i for interpolation $ip")) end + +################################### +# Lagrange dim 3 RefPrism order 1 # +################################### +# Build on https://defelement.com/elements/examples/prism-Lagrange-1.html +getnbasefunctions(::Lagrange{3,RefPrism,1}) = 6 + +facedof_indices(::Lagrange{3,RefPrism,1}) = ((1,3,2), (1,2,5,4), (3,1,4,6), (2,3,6,5), (4,5,6)) +edgedof_indices(::Lagrange{3,RefPrism,1}) = ((2,1), (1,3), (1,4), (3,2), (2,5), (3,6), (4,5), (4,6), (6,5)) + +function reference_coordinates(::Lagrange{3,RefPrism,1}) + return [Vec{3, Float64}((0.0, 0.0, 0.0)), + Vec{3, Float64}((1.0, 0.0, 0.0)), + Vec{3, Float64}((0.0, 1.0, 0.0)), + Vec{3, Float64}((0.0, 0.0, 1.0)), + Vec{3, Float64}((1.0, 0.0, 1.0)), + Vec{3, Float64}((0.0, 1.0, 1.0))] +end + +function value(ip::Lagrange{3,RefPrism,1}, i::Int, ξ::Vec{3}) + (x,y,z) = ξ + i == 1 && return 1-x-y -z*(1-x-y) + i == 2 && return x*(1-z) + i == 3 && return y*(1-z) + i == 4 && return z*(1-x-y) + i == 5 && return x*z + i == 6 && return y*z + throw(ArgumentError("no shape function $i for interpolation $ip")) +end + +################################### +# Lagrange dim 3 RefPrism order 2 # +################################### +# Build on https://defelement.com/elements/examples/prism-Lagrange-2.html . +# This is simply the tensor-product of a quadratic triangle with a quadratic line. +getnbasefunctions(::Lagrange{3,RefPrism,2}) = 18 + +facedof_indices(::Lagrange{3,RefPrism,2}) = ( + #Vertices| Edges | Face + (1,3,2 , 8,10,7 ), + (1,2,5,4, 7,11,13,9, 16), + (3,1,4,6, 8,9,14,12, 17), + (2,3,6,5, 10,12,15,11, 18), + (4,5,6 , 13,15,14 ), +) +facedof_interior_indices(::Lagrange{3,RefPrism,2}) = ( + #Vertices| Edges | Face + (), + (16,), + (17,), + (18,), + (), +) +edgedof_indices(::Lagrange{3,RefPrism,2}) = ( + #Vert|Edge + (2,1, 7), + (1,3, 8), + (1,4, 9), + (3,2, 10), + (2,5, 11), + (3,6, 12), + (4,5, 13), + (4,6, 14), + (6,5, 15), +) +edgedof_interior_indices(::Lagrange{3,RefPrism,2}) = ( + #Vert|Edge + (7,), + (8,), + (9,), + (10,), + (11,), + (12,), + (13,), + (14,), + (15,), +) + +function reference_coordinates(::Lagrange{3,RefPrism,2}) + return [Vec{3, Float64}((0.0, 0.0, 0.0)), + Vec{3, Float64}((1.0, 0.0, 0.0)), + Vec{3, Float64}((0.0, 1.0, 0.0)), + Vec{3, Float64}((0.0, 0.0, 1.0)), + Vec{3, Float64}((1.0, 0.0, 1.0)), + Vec{3, Float64}((0.0, 1.0, 1.0)), + Vec{3, Float64}((1/2, 0.0, 0.0)), + Vec{3, Float64}((0.0, 1/2, 0.0)), + Vec{3, Float64}((0.0, 0.0, 1/2)), + Vec{3, Float64}((1/2, 1/2, 0.0)), + Vec{3, Float64}((1.0, 0.0, 1/2)), + Vec{3, Float64}((0.0, 1.0, 1/2)), + Vec{3, Float64}((1/2, 0.0, 1.0)), + Vec{3, Float64}((0.0, 1/2, 1.0)), + Vec{3, Float64}((1/2, 1/2, 1.0)), + Vec{3, Float64}((1/2, 0.0, 1/2)), + Vec{3, Float64}((0.0, 1/2, 1/2)), + Vec{3, Float64}((1/2, 1/2, 1/2)),] +end + +function value(ip::Lagrange{3,RefPrism,2}, i::Int, ξ::Vec{3}) + (x,y,z) = ξ + x² = x*x + y² = y*y + z² = z*z + i == 1 && return 4*x²*z² - 6x²*z +2x² +8x*y*z² -12x*y*z +4x*y -6x*z² +9x*z -3x +4y²*z² -6y²*z + 2y² -6y*z² +9y*z -3*y +2z² -3z +1 + i == 2 && return x*(4x*z² -6x*z +2x -2z² +3z -1) + i == 3 && return y*(4y*z² -6y*z +2y -2z² +3z -1) + i == 4 && return z*(4x²*z -2x² + 8x*y*z -4x*y -6x*z +3x +4y²*z -2y² -6y*z +3y +2z -1) + i == 5 && return x*z*(4x*z -2x -2z +1) + i == 6 && return y*z*(4y*z -2y -2z +1) + i == 7 && return 4x*(-2x*z² +3x*z -x -2*y*z² +3y*z -y +2z² -3z +1) + i == 8 && return 4y*(-2x*z² +3x*z -x -2*y*z² +3y*z -y +2z² -3z +1) + i == 9 && return 4z*(-2x²*z +2x² -4x*y*z +4x*y +3x*z -3x -2y²*z +2y² +3y*z -3y -z +1) + i == 10 && return 4x*y*(2z² -3z +1) + i == 11 && return 4x*z*(-2x*z +2x +z -1) + i == 12 && return 4y*z*(-2y*z +2y +z -1) + i == 13 && return 4x*z*(-2x*z +x -2y*z +y +2z -1) + i == 14 && return 4y*z*(-2x*z +x -2y*z +y +2z -1) + i == 15 && return 4x*y*z*(2z -1) + i == 16 && return 16x*z*(x*z -x +y*z -y -z +1) + i == 17 && return 16y*z*(x*z -x +y*z -y -z +1) + i == 18 && return 16x*y*z*(1 -z) + throw(ArgumentError("no shape function $i for interpolation $ip")) +end + ################### # Bubble elements # ################### @@ -716,18 +924,15 @@ Lagrange element with bubble stabilization. struct BubbleEnrichedLagrange{dim,ref_shape,order} <: Interpolation{dim,ref_shape,order} end getlowerdim(ip::BubbleEnrichedLagrange{dim,ref_shape,order}) where {dim,ref_shape,order} = Lagrange{dim-1,ref_shape,order}() -# Vertices for all Bubble interpolations are the same -vertices(::BubbleEnrichedLagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,)) -nvertexdofs(::BubbleEnrichedLagrange{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, order} = 1 - ################################################ # Lagrange-Bubble dim 2 RefTetrahedron order 1 # ################################################ # Taken from https://defelement.com/elements/bubble-enriched-lagrange.html getnbasefunctions(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = 4 -ncelldofs(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = 1 -faces(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1)) +vertexdof_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,), (2,), (3,)) +facedof_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1)) +celldof_interior_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = (4,) function reference_coordinates(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) return [Vec{2, Float64}((1.0, 0.0)), @@ -752,9 +957,8 @@ end struct Serendipity{dim,shape,order} <: Interpolation{dim,shape,order} end # Vertices for all Serendipity interpolations are the same -vertices(::Serendipity{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,)) -vertices(::Serendipity{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,)) -nvertexdofs(::Serendipity{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, order} = 1 +vertexdof_indices(::Serendipity{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,)) +vertexdof_indices(::Serendipity{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,)) ##################################### # Serendipity dim 2 RefCube order 2 # @@ -762,9 +966,9 @@ nvertexdofs(::Serendipity{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, or getnbasefunctions(::Serendipity{2,RefCube,2}) = 8 getlowerdim(::Serendipity{2,RefCube,2}) = Lagrange{1,RefCube,2}() getlowerorder(::Serendipity{2,RefCube,2}) = Lagrange{2,RefCube,1}() -nfacedofs(::Serendipity{2,RefCube,2}) = 1 -faces(::Serendipity{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8)) +facedof_indices(::Serendipity{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8)) +facedof_interior_indices(::Serendipity{2,RefCube,2}) = ((5,), (6,), (7,), (8,)) function reference_coordinates(::Serendipity{2,RefCube,2}) return [Vec{2, Float64}((-1.0, -1.0)), @@ -794,12 +998,37 @@ end ##################################### # Serendipity dim 3 RefCube order 2 # ##################################### +# Note that second order serendipity hex has no interior face indices. getnbasefunctions(::Serendipity{3,RefCube,2}) = 20 getlowerdim(::Serendipity{3,RefCube,2}) = Serendipity{2,RefCube,2}() getlowerorder(::Serendipity{3,RefCube,2}) = Lagrange{3,RefCube,1}() -nedgedofs(::Serendipity{3,RefCube,2}) = 1 -faces(::Serendipity{3,RefCube,2}) = ((1,4,3,2,12,11,10,9), (1,2,6,5,9,18,13,17), (2,3,7,6,10,19,14,18), (3,4,8,7,11,20,15,19), (1,5,8,4,17,16,20,12), (5,6,7,8,13,14,15,16)) +facedof_indices(::Serendipity{3,RefCube,2}) = ( + (1,4,3,2, 12,11,10,9), + (1,2,6,5, 9,18,13,17), + (2,3,7,6, 10,19,14,18), + (3,4,8,7, 11,20,15,19), + (1,5,8,4, 17,16,20,12), + (5,6,7,8, 13,14,15,16) +) +edgedof_indices(::Serendipity{3,RefCube,2}) = ( + (1,2, 9), + (2,3, 10), + (3,4, 11), + (4,1, 12), + (1,5, 17), + (2,6, 18), + (3,7, 19), + (4,8, 20), + (5,6, 13), + (6,7, 14), + (7,8, 15), + (8,5, 16) +) + +edgedof_interior_indices(::Serendipity{3,RefCube,2}) = ( + (9,), (10,), (11,), (12,), (17), (18,), (19,), (20,), (13,), (14,), (15,), (16,) +) function reference_coordinates(::Serendipity{3,RefCube,2}) return [Vec{3, Float64}((-1.0, -1.0, -1.0)), @@ -866,8 +1095,9 @@ and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 7.R3 (1 struct CrouzeixRaviart{dim,order} <: Interpolation{dim,RefTetrahedron,order} end getnbasefunctions(::CrouzeixRaviart{2,1}) = 3 -nfacedofs(::CrouzeixRaviart{2,1}) = 1 -faces(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,)) + +facedof_indices(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,)) +facedof_interior_indices(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,)) function reference_coordinates(::CrouzeixRaviart{2,1}) return [Vec{2, Float64}((0.5, 0.5)), diff --git a/src/iterators.jl b/src/iterators.jl index 0cd0c44b68..a37a9f28c5 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -57,10 +57,10 @@ function CellCache(grid::Grid{dim,C,T}, flags::UpdateFlags=UpdateFlags()) where return CellCache(flags, grid, ScalarWrapper(-1), nodes, coords, nothing, Int[]) end -function CellCache(dh::Union{DofHandler{dim,T},MixedDofHandler{dim,T}}, flags::UpdateFlags=UpdateFlags()) where {dim,T} +function CellCache(dh::Union{DofHandler{dim},MixedDofHandler{dim}}, flags::UpdateFlags=UpdateFlags()) where {dim} N = nnodes_per_cell(dh.grid) nodes = zeros(Int, N) - coords = zeros(Vec{dim,T}, N) + coords = zeros(Vec{dim, get_coordinate_eltype(dh.grid)}, N) n = ndofs_per_cell(dh) celldofs = zeros(Int, n) return CellCache(flags, dh.grid, ScalarWrapper(-1), nodes, coords, dh, celldofs) diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index eee76a9b0d..7bb0f65038 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -21,7 +21,7 @@ for (func_interpol, quad_rule) in ( fe_valtype == FaceVectorValues && @test getnbasefunctions(fv) == n_basefuncs * Ferrite.getdim(func_interpol) xs, n = valid_coordinates_and_normals(func_interpol) - for face in 1:getnfaces(func_interpol) + for face in 1:Ferrite.nfaces(func_interpol) reinit!(fv, xs, face) @test Ferrite.getcurrentface(fv) == face @@ -61,7 +61,7 @@ for (func_interpol, quad_rule) in ( for i in 1:getnquadpoints(fv) vol += getdetJdV(fv,i) end - x_face = xs[[Ferrite.faces(func_interpol)[face]...]] + x_face = xs[[Ferrite.facedof_indices(func_interpol)[face]...]] @test vol ≈ calculate_volume(Ferrite.getlowerdim(func_interpol), x_face) # Test quadrature rule after reinit! with ref. coords diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 278b87be9c..e33a0695b8 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -1,6 +1,6 @@ @testset "interpolations" begin -for interpolation in ( +@testset "$interpolation" for interpolation in ( Lagrange{1, RefCube, 1}(), Lagrange{1, RefCube, 2}(), Lagrange{2, RefCube, 1}(), @@ -9,13 +9,15 @@ for interpolation in ( Lagrange{2, RefTetrahedron, 2}(), Lagrange{2, RefTetrahedron, 3}(), Lagrange{2, RefTetrahedron, 4}(), - Lagrange{2, RefTetrahedron, 5}(), + #Lagrange{2, RefTetrahedron, 5}(), Lagrange{3, RefCube, 1}(), Lagrange{3, RefCube, 2}(), Serendipity{2, RefCube, 2}(), Serendipity{3, RefCube, 2}(), Lagrange{3, RefTetrahedron, 1}(), Lagrange{3, RefTetrahedron, 2}(), + Lagrange{3, RefPrism, 1}(), + Lagrange{3, RefPrism, 2}(), # DiscontinuousLagrange{1, RefCube, 0}(), DiscontinuousLagrange{2, RefCube, 0}(), @@ -34,35 +36,79 @@ for interpolation in ( ) # Test of utility functions - ndim = Ferrite.getdim(interpolation) - r_shape = Ferrite.getrefshape(interpolation) + ref_dim = Ferrite.getdim(interpolation) + ref_shape = Ferrite.getrefshape(interpolation) func_order = Ferrite.getorder(interpolation) - @test typeof(interpolation) <: Interpolation{ndim,r_shape,func_order} + @test typeof(interpolation) <: Interpolation{ref_dim,ref_shape,func_order} # Note that not every element formulation exists for every order and dimension. - applicable(Ferrite.getlowerdim, interpolation) && @test typeof(Ferrite.getlowerdim(interpolation)) <: Interpolation{ndim-1} - applicable(Ferrite.getlowerorder, interpolation) && @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ndim,r_shape,func_order-1} + applicable(Ferrite.getlowerdim, interpolation) && @test typeof(Ferrite.getlowerdim(interpolation)) <: Interpolation{ref_dim-1} + applicable(Ferrite.getlowerorder, interpolation) && @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ref_dim,ref_shape,func_order-1} + # Check partition of unity at random point. n_basefuncs = getnbasefunctions(interpolation) - x = rand(Tensor{1, ndim}) - f = (x) -> Ferrite.value(interpolation, Tensor{1, ndim}(x)) + x = rand(Tensor{1, ref_dim}) + f = (x) -> Ferrite.value(interpolation, Tensor{1, ref_dim}(x)) @test vec(ForwardDiff.jacobian(f, Array(x))') ≈ reinterpret(Float64, Ferrite.derivative(interpolation, x)) @test sum(Ferrite.value(interpolation, x)) ≈ 1.0 + # Check if the important functions are consistens coords = Ferrite.reference_coordinates(interpolation) - for node in 1:n_basefuncs - N_node = Ferrite.value(interpolation, coords[node]) - for k in 1:node - if k == node - @test N_node[k] ≈ 1.0 + @test length(coords) == n_basefuncs + @test Ferrite.value(interpolation, length(coords), x) == Ferrite.value(interpolation, length(coords), x) + @test_throws ArgumentError Ferrite.value(interpolation, length(coords)+1, x) + + # Test whether we have for each entity corresponding dof indices (possibly empty) + @test length(Ferrite.vertexdof_indices(interpolation)) == Ferrite.nvertices(interpolation) + if ref_dim > 1 + @test length(Ferrite.facedof_indices(interpolation)) == Ferrite.nfaces(interpolation) + @test length(Ferrite.facedof_interior_indices(interpolation)) == Ferrite.nfaces(interpolation) + elseif ref_dim > 2 + @test length(Ferrite.edgedof_indices(interpolation)) == Ferrite.nedges(interpolation) + @test length(Ferrite.edgedof_interior_indices(interpolation)) == Ferrite.nedges(interpolation) + end + # We have at least as many edge/face dofs as we have edge/face interior dofs + if ref_dim > 1 + @test all(length.(Ferrite.facedof_interior_indices(interpolation)) .<= length.(Ferrite.facedof_indices(interpolation))) + elseif ref_dim > 2 + @test all(length.(Ferrite.edgedof_interior_indices(interpolation)) .<= length.(Ferrite.edgedof_indices(interpolation))) + end + # The total number of dofs must match the number of base functions + totaldofs = sum(length.(Ferrite.vertexdof_indices(interpolation));init=0) + if ref_dim > 1 + totaldofs += sum(length.(Ferrite.facedof_interior_indices(interpolation));init=0) + end + if ref_dim > 2 + totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) + end + totaldofs += length(Ferrite.celldof_interior_indices(interpolation)) + @test totaldofs == n_basefuncs + + # The dof indices are valid. + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.vertexdof_indices(interpolation)]) + if ref_dim > 1 + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.facedof_indices(interpolation)]) + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.facedof_interior_indices(interpolation)]) + elseif ref_dim > 2 + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.edgedof_indices(interpolation)]) + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.edgedof_interior_indices(interpolation)]) + end + @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.celldof_interior_indices(interpolation)]) + + # Check for dirac delta property of interpolation + @testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs + N_dof = Ferrite.value(interpolation, coords[dof]) + for k in 1:n_basefuncs + if k == dof + @test N_dof[k] ≈ 1.0 else - @test N_node[k] ≈ 0.0 atol=4eps(Float64) + @test N_dof[k] ≈ 0.0 atol=4eps(Float64) #broken=typeof(interpolation)==Lagrange{2, RefTetrahedron, 5}&&dof==4&&k==18 end end end - # Test that faces(...) return in counter clockwise order (viewing from the outside) + # Test that facedof_indices(...) return in counter clockwise order (viewing from the outside) if interpolation isa Lagrange function __outward_normal(coords::Vector{<:Vec{1}}, nodes) n = coords[nodes[1]] @@ -81,21 +127,22 @@ for interpolation in ( n = (p3 - p2) × (p1 - p2) return n / norm(n) end - for (facenodes, normal) in zip(Ferrite.faces(interpolation), reference_normals(interpolation)) + for (facenodes, normal) in zip(Ferrite.facedof_indices(interpolation), reference_normals(interpolation)) @test __outward_normal(coords, facenodes) ≈ normal end end + # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/520 interpolation_type = typeof(interpolation).name.wrapper if func_order > 1 && interpolation_type != Ferrite.Serendipity - first_order = interpolation_type{ndim,r_shape,1}() - for (highorderface, firstorderface) in zip(Ferrite.faces(interpolation), Ferrite.faces(first_order)) + first_order = interpolation_type{ref_dim,ref_shape,1}() + for (highorderface, firstorderface) in zip(Ferrite.facedof_indices(interpolation), Ferrite.facedof_indices(first_order)) for (h_node, f_node) in zip(highorderface, firstorderface) @test h_node == f_node end end - if ndim > 2 - for (highorderedge, firstorderedge) in zip(Ferrite.edges(interpolation), Ferrite.edges(first_order)) + if ref_dim > 2 + for (highorderedge, firstorderedge) in zip(Ferrite.edgedof_indices(interpolation), Ferrite.edgedof_indices(first_order)) for (h_node, f_node) in zip(highorderedge, firstorderedge) @test h_node == f_node end diff --git a/test/test_utils.jl b/test/test_utils.jl index 00f1327f4f..b64f0f3247 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -5,6 +5,7 @@ ##################################### reference_volume(::Interpolation{dim, RefCube}) where {dim} = 2^dim reference_volume(::Interpolation{dim, RefTetrahedron}) where {dim} = 1 / factorial(dim) +reference_volume(::Interpolation{ 3, RefPrism}) = 1/2 # For faces reference_volume(fs::Interpolation, ::Int) = reference_volume(Ferrite.getlowerdim(fs)) reference_volume(fs::Interpolation{2, RefTetrahedron}, face::Int) = face == 1 ? sqrt(2) : 1.0 @@ -53,6 +54,15 @@ function reference_normals(::Lagrange{3, RefCube}) Vec{3, Float64}(( 0.0, 0.0, 1.0))] end +# Lagrange{3, Wedge} +function reference_normals(::Lagrange{3, RefPrism}) + return [Vec{3, Float64}(( 0.0, 0.0, -1.0)), + Vec{3, Float64}(( 0.0, -1.0, 0.0)), + Vec{3, Float64}((-1.0, 0.0, 0.0)), + Vec{3, Float64}((1/√2, 1/√2, 0.0)), + Vec{3, Float64}(( 0.0, 0.0, 1.0))] +end + # Serendipity{2, RefCube} reference_normals(::Serendipity{2, RefCube, 2}) = reference_normals(Lagrange{2, RefCube,1}()) @@ -161,10 +171,6 @@ function calculate_volume(::Lagrange{0, RefCube, order}, ::Vector{Vec{1, T}}) wh return one(T) end -getnfaces(::Interpolation{dim, RefCube}) where {dim} = 2*dim -getnfaces(::Interpolation{2, RefTetrahedron}) = 3 -getnfaces(::Interpolation{3, RefTetrahedron}) = 4 - coords_on_faces(x, ::Lagrange{1, RefCube, 1}) = ([x[1]], [x[2]]) coords_on_faces(x, ::Lagrange{1, RefCube, 2}) = ([x[1]], [x[2]]) coords_on_faces(x, ::Lagrange{2, RefCube, 1}) =