diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl index 43c9860cc3..b6512c3b61 100644 --- a/src/FEValues/FacetValues.jl +++ b/src/FEValues/FacetValues.jl @@ -51,9 +51,9 @@ function FacetValues( ) where {T, sdim, FunDiffOrder, GeoDiffOrder} # max(GeoDiffOrder, 1) ensures that we get the jacobian needed to calculate the normal. - geo_mapping = map(qr -> GeometryMapping{max(GeoDiffOrder, 1)}(T, ip_geo.ip, qr), fqr.face_rules) - fun_values = map(qr -> FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo), fqr.face_rules) - max_nquadpoints = maximum(qr -> length(getweights(qr)), fqr.face_rules) + geo_mapping = map(qr -> GeometryMapping{max(GeoDiffOrder, 1)}(T, ip_geo.ip, qr), fqr.facet_rules) + fun_values = map(qr -> FunctionValues{FunDiffOrder}(T, ip_fun, qr, ip_geo), fqr.facet_rules) + max_nquadpoints = maximum(qr -> length(getweights(qr)), fqr.facet_rules) # detJdV always calculated, since we needed to calculate the jacobian anyways for the normal. detJdV = fill(T(NaN), max_nquadpoints) normals = fill(zero(Vec{sdim, T}) * T(NaN), max_nquadpoints) @@ -112,21 +112,21 @@ getnormal(fv::FacetValues, qp::Int) = fv.normals[qp] nfacets(fv::FacetValues) = length(fv.geo_mapping) -function set_current_facet!(fv::FacetValues, face_nr::Int) - # Checking face_nr before setting current_facet allows us to use @inbounds +function set_current_facet!(fv::FacetValues, facet_nr::Int) + # Checking facet_nr before setting current_facet allows us to use @inbounds # when indexing by getcurrentfacet(fv) in other places! - checkbounds(Bool, 1:nfacets(fv), face_nr) || throw(ArgumentError("Face index out of range.")) - fv.current_facet = face_nr + checkbounds(Bool, 1:nfacets(fv), facet_nr) || throw(ArgumentError("Facet nr is out of range.")) + fv.current_facet = facet_nr return end -@inline function reinit!(fv::FacetValues, x::AbstractVector, face_nr::Int) - return reinit!(fv, nothing, x, face_nr) +@inline function reinit!(fv::FacetValues, x::AbstractVector, facet_nr::Int) + return reinit!(fv, nothing, x, facet_nr) end -function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::AbstractVector{Vec{dim, T}}, face_nr::Int) where {dim, T} +function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::AbstractVector{Vec{dim, T}}, facet_nr::Int) where {dim, T} check_reinit_sdim_consistency(:FacetValues, shape_gradient_type(fv), eltype(x)) - set_current_facet!(fv, face_nr) + set_current_facet!(fv, facet_nr) n_geom_basefuncs = getngeobasefunctions(fv) if !checkbounds(Bool, x, 1:n_geom_basefuncs) || length(x) != n_geom_basefuncs throw_incompatible_coord_length(length(x), n_geom_basefuncs) @@ -139,11 +139,11 @@ function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::Abstrac throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings")) end - @inbounds for (q_point, w) in pairs(getweights(fv.fqr, face_nr)) + @inbounds for (q_point, w) in pairs(getweights(fv.fqr, facet_nr)) mapping = calculate_mapping(geo_mapping, q_point, x) J = getjacobian(mapping) # See the `Ferrite.embedding_det` docstring for more background - weight_norm = weighted_normal(J, getrefshape(geo_mapping.ip), face_nr) + weight_norm = weighted_normal(J, getrefshape(geo_mapping.ip), facet_nr) detJ = norm(weight_norm) detJ > 0.0 || throw_detJ_not_pos(detJ) @inbounds fv.detJdV[q_point] = detJ * w @@ -161,11 +161,11 @@ function Base.show(io::IO, d::MIME"text/plain", fv::FacetValues) sdim = GradT === nothing ? nothing : sdim_from_gradtype(GradT) vstr = vdim == 0 ? "scalar" : "vdim=$vdim" print(io, "FacetValues(", vstr, ", rdim=$rdim, sdim=$sdim): ") - nqp = getnquadpoints.(fv.fqr.face_rules) + nqp = getnquadpoints.(fv.fqr.facet_rules) if all(n == first(nqp) for n in nqp) - println(io, first(nqp), " quadrature points per face") + println(io, first(nqp), " quadrature points per facet") else - println(io, tuple(nqp...), " quadrature points on each face") + println(io, tuple(nqp...), " quadrature points on each facet") end print(io, " Function interpolation: "); show(io, d, function_interpolation(fv)) print(io, "\nGeometric interpolation: ") @@ -225,9 +225,9 @@ function spatial_coordinate(bcv::BCValues, q_point::Int, xh::AbstractVector{Vec{ 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_entity[] + facet = 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] + x += bcv.M[i, q_point, facet] * xh[i] # geometric_value(fe_v, q_point, i) * xh[i] end return x end diff --git a/src/FEValues/InterfaceValues.jl b/src/FEValues/InterfaceValues.jl index d452a28242..64ea4dc3fd 100644 --- a/src/FEValues/InterfaceValues.jl +++ b/src/FEValues/InterfaceValues.jl @@ -17,10 +17,10 @@ The first element of the interface is denoted "here" and the second element "the geometric interpolation. * `InterfaceValues(qr_here::FacetQuadratureRule, ip_here::Interpolation, ip_geo_here::Interpolation, qr_there::FacetQuadratureRule, ip_there::Interpolation, ip_geo_there::Interpolation)`: same as above but with given geometric interpolation. -* `InterfaceValues(fv::FacetValues)`: quadrature rule and interpolations from face values +* `InterfaceValues(fv::FacetValues)`: quadrature rule and interpolations from facet values (same on both sides). * `InterfaceValues(fv_here::FacetValues, fv_there::FacetValues)`: quadrature rule and - interpolations from the face values. + interpolations from the facet values. **Associated methods:** * [`shape_value_average`](@ref) @@ -396,12 +396,12 @@ end Relative orientation information for 1D and 2D interfaces in 2D and 3D elements respectively. This information is used to construct the transformation matrix to transform the quadrature points from facet_a to facet_b achieving synced -spatial coordinates. Face B's orientation relative to Face A's can +spatial coordinates. Facet B's orientation relative to Facet A's can possibly be flipped (i.e. the vertices indices order is reversed) and the vertices can be rotated against each other. -The reference orientation of face B is such that the first node +The reference orientation of facet B is such that the first node has the lowest vertex index. Thus, this structure also stores the -shift of the lowest vertex index which is used to reorient the face in +shift of the lowest vertex index which is used to reorient the facet in case of flipping [`transform_interface_points!`](@ref). """ struct InterfaceOrientationInfo{RefShapeA, RefShapeB} @@ -415,7 +415,7 @@ end """ InterfaceOrientationInfo(cell_a::AbstractCell, cell_b::AbstractCell, facet_a::Int, facet_b::Int) -Return the relative orientation info for face B with regards to face A. +Return the relative orientation info for facet B with regards to facet A. Relative orientation is computed using a [`OrientationInfo`](@ref) for each side of the interface. """ function InterfaceOrientationInfo(cell_a::AbstractCell{RefShapeA}, cell_b::AbstractCell{RefShapeB}, facet_a::Int, facet_b::Int) where {RefShapeA <: AbstractRefShape, RefShapeB <: AbstractRefShape} @@ -435,8 +435,8 @@ end Returns the transformation matrix corresponding to the interface orientation information stored in `InterfaceOrientationInfo`. The transformation matrix is constructed using a combination of affine transformations defined for each interface reference shape. -The transformation for a flipped face is a function of both relative orientation and the orientation of the second face. -If the face is not flipped then the transformation is a function of relative orientation only. +The transformation for a flipped facet is a function of both relative orientation and the orientation of the second facet. +If the facet is not flipped then the transformation is a function of relative orientation only. """ get_transformation_matrix @@ -485,9 +485,9 @@ end @doc raw""" transform_interface_points!(dst::AbstractVector{Vec{3, Float64}}, points::AbstractVector{Vec{3, Float64}}, interface_transformation::InterfaceOrientationInfo) -Transform the points from face A to face B using the orientation information of the interface and store it in the vector dst. -For 3D, the faces are transformed into regular polygons such that the rotation angle is the shift in reference node index × 2π ÷ number of edges in face. -If the face is flipped then the flipping is about the axis that preserves the position of the first node (which is the reference node after being rotated to be in the first position, +Transform the points from facet A to facet B using the orientation information of the interface and store it in the vector dst. +For 3D, the facets are transformed into regular polygons such that the rotation angle is the shift in reference node index × 2π ÷ number of edges in facet. +If the facet is flipped then the flipping is about the axis that preserves the position of the first node (which is the reference node after being rotated to be in the first position, it's rotated back in the opposite direction after flipping). Take for example the interface ``` diff --git a/src/FEValues/face_integrals.jl b/src/FEValues/facet_integrals.jl similarity index 57% rename from src/FEValues/face_integrals.jl rename to src/FEValues/facet_integrals.jl index 2fbd0506d9..6c172cdf10 100644 --- a/src/FEValues/face_integrals.jl +++ b/src/FEValues/facet_integrals.jl @@ -15,8 +15,8 @@ This is the inverse of `facet_to_element_transformation`. element_to_facet_transformation """ - weighted_normal(J::AbstractTensor, fv::FacetValues, face::Int) - weighted_normal(J::AbstractTensor, ::Type{<:AbstractRefShape}, face::Int) + weighted_normal(J::AbstractTensor, fv::FacetValues, facet::Int) + weighted_normal(J::AbstractTensor, ::Type{<:AbstractRefShape}, facet::Int) Compute the vector normal to the facet weighted by the area ratio between the facet and the reference facet. This is computed by taking the cross product of the Jacobian components that @@ -68,23 +68,23 @@ end ################## # Mapping from to 0D node to 1D line vertex. -function facet_to_element_transformation(::Union{Vec{0, T}, Vec{1, T}}, ::Type{RefLine}, face::Int) where {T} - face == 1 && return Vec{1, T}((-one(T),)) - face == 2 && return Vec{1, T}((one(T),)) +function facet_to_element_transformation(::Union{Vec{0, T}, Vec{1, T}}, ::Type{RefLine}, facet::Int) where {T} + facet == 1 && return Vec{1, T}((-one(T),)) + facet == 2 && return Vec{1, T}((one(T),)) throw(ArgumentError("unknown facet number")) end # Mapping from 1D line to point. -function element_to_facet_transformation(point::Vec{1, T}, ::Type{RefLine}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{1, T}, ::Type{RefLine}, facet::Int) where {T} x = point[] - face == 1 && return Vec(-x) - face == 2 && return Vec(x) + facet == 1 && return Vec(-x) + facet == 2 && return Vec(x) throw(ArgumentError("unknown facet number")) end -function weighted_normal(::Tensor{2, 1, T}, ::Type{RefLine}, face::Int) where {T} - face == 1 && return Vec{1, T}((-one(T),)) - face == 2 && return Vec{1, T}((one(T),)) +function weighted_normal(::Tensor{2, 1, T}, ::Type{RefLine}, facet::Int) where {T} + facet == 1 && return Vec{1, T}((-one(T),)) + facet == 2 && return Vec{1, T}((one(T),)) throw(ArgumentError("unknown facet number")) end @@ -93,31 +93,31 @@ end ########################### # Mapping from 1D line to 2D face of a quadrilateral. -function facet_to_element_transformation(point::Vec{1, T}, ::Type{RefQuadrilateral}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{1, T}, ::Type{RefQuadrilateral}, facet::Int) where {T} x = point[1] - face == 1 && return Vec{2, T}((x, -one(T))) - face == 2 && return Vec{2, T}((one(T), x)) - face == 3 && return Vec{2, T}((-x, one(T))) - face == 4 && return Vec{2, T}((-one(T), -x)) + facet == 1 && return Vec{2, T}((x, -one(T))) + facet == 2 && return Vec{2, T}((one(T), x)) + facet == 3 && return Vec{2, T}((-x, one(T))) + facet == 4 && return Vec{2, T}((-one(T), -x)) throw(ArgumentError("unknown facet number")) end # Mapping from 2D face of a quadrilateral to 1D line. -function element_to_facet_transformation(point::Vec{2, T}, ::Type{RefQuadrilateral}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{2, T}, ::Type{RefQuadrilateral}, facet::Int) where {T} x, y = point - face == 1 && return Vec(x) - face == 2 && return Vec(y) - face == 3 && return Vec(-x) - face == 4 && return Vec(-y) + facet == 1 && return Vec(x) + facet == 2 && return Vec(y) + facet == 3 && return Vec(-x) + facet == 4 && return Vec(-y) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 2}, ::Type{RefQuadrilateral}, face::Int) +function weighted_normal(J::Tensor{2, 2}, ::Type{RefQuadrilateral}, facet::Int) @inbounds begin - face == 1 && return Vec{2}((J[2, 1], -J[1, 1])) - face == 2 && return Vec{2}((J[2, 2], -J[1, 2])) - face == 3 && return Vec{2}((-J[2, 1], J[1, 1])) - face == 4 && return Vec{2}((-J[2, 2], J[1, 2])) + facet == 1 && return Vec{2}((J[2, 1], -J[1, 1])) + facet == 2 && return Vec{2}((J[2, 2], -J[1, 2])) + facet == 3 && return Vec{2}((-J[2, 1], J[1, 1])) + facet == 4 && return Vec{2}((-J[2, 2], J[1, 2])) end throw(ArgumentError("unknown facet number")) end @@ -127,28 +127,28 @@ end ###################### # Mapping from 1D line to 2D face of a triangle. -function facet_to_element_transformation(point::Vec{1, T}, ::Type{RefTriangle}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{1, T}, ::Type{RefTriangle}, facet::Int) where {T} x = (point[1] + one(T)) / 2 - face == 1 && return Vec{2, T}((one(T) - x, x)) - face == 2 && return Vec{2, T}((zero(T), one(T) - x)) - face == 3 && return Vec{2, T}((x, zero(T))) + facet == 1 && return Vec{2, T}((one(T) - x, x)) + facet == 2 && return Vec{2, T}((zero(T), one(T) - x)) + facet == 3 && return Vec{2, T}((x, zero(T))) throw(ArgumentError("unknown facet number")) end # Mapping from 2D face of a triangle to 1D line. -function element_to_facet_transformation(point::Vec{2, T}, ::Type{RefTriangle}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{2, T}, ::Type{RefTriangle}, facet::Int) where {T} x, y = point - face == 1 && return Vec(one(T) - x * 2) - face == 2 && return Vec(one(T) - y * 2) - face == 3 && return Vec(x * 2 - one(T)) + facet == 1 && return Vec(one(T) - x * 2) + facet == 2 && return Vec(one(T) - y * 2) + facet == 3 && return Vec(x * 2 - one(T)) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 2}, ::Type{RefTriangle}, face::Int) +function weighted_normal(J::Tensor{2, 2}, ::Type{RefTriangle}, facet::Int) @inbounds begin - face == 1 && return Vec{2}((-(J[2, 1] - J[2, 2]), J[1, 1] - J[1, 2])) - face == 2 && return Vec{2}((-J[2, 2], J[1, 2])) - face == 3 && return Vec{2}((J[2, 1], -J[1, 1])) + facet == 1 && return Vec{2}((-(J[2, 1] - J[2, 2]), J[1, 1] - J[1, 2])) + facet == 2 && return Vec{2}((-J[2, 2], J[1, 2])) + facet == 3 && return Vec{2}((J[2, 1], -J[1, 1])) end throw(ArgumentError("unknown facet number")) end @@ -158,37 +158,37 @@ end ######################## # Mapping from 2D quadrilateral to 3D face of a hexahedron. -function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefHexahedron}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefHexahedron}, facet::Int) where {T} x, y = point - face == 1 && return Vec{3, T}((y, x, -one(T))) - face == 2 && return Vec{3, T}((x, -one(T), y)) - face == 3 && return Vec{3, T}((one(T), x, y)) - face == 4 && return Vec{3, T}((-x, one(T), y)) - face == 5 && return Vec{3, T}((-one(T), y, x)) - face == 6 && return Vec{3, T}((x, y, one(T))) + facet == 1 && return Vec{3, T}((y, x, -one(T))) + facet == 2 && return Vec{3, T}((x, -one(T), y)) + facet == 3 && return Vec{3, T}((one(T), x, y)) + facet == 4 && return Vec{3, T}((-x, one(T), y)) + facet == 5 && return Vec{3, T}((-one(T), y, x)) + facet == 6 && return Vec{3, T}((x, y, one(T))) throw(ArgumentError("unknown facet number")) end # Mapping from 3D face of a hexahedron to 2D quadrilateral. -function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefHexahedron}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefHexahedron}, facet::Int) where {T} x, y, z = point - face == 1 && return Vec(y, x) - face == 2 && return Vec(x, z) - face == 3 && return Vec(y, z) - face == 4 && return Vec(-x, z) - face == 5 && return Vec(z, y) - face == 6 && return Vec(x, y) + facet == 1 && return Vec(y, x) + facet == 2 && return Vec(x, z) + facet == 3 && return Vec(y, z) + facet == 4 && return Vec(-x, z) + facet == 5 && return Vec(z, y) + facet == 6 && return Vec(x, y) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 3}, ::Type{RefHexahedron}, face::Int) +function weighted_normal(J::Tensor{2, 3}, ::Type{RefHexahedron}, facet::Int) @inbounds begin - face == 1 && return J[:, 2] × J[:, 1] - face == 2 && return J[:, 1] × J[:, 3] - face == 3 && return J[:, 2] × J[:, 3] - face == 4 && return J[:, 3] × J[:, 1] - face == 5 && return J[:, 3] × J[:, 2] - face == 6 && return J[:, 1] × J[:, 2] + facet == 1 && return J[:, 2] × J[:, 1] + facet == 2 && return J[:, 1] × J[:, 3] + facet == 3 && return J[:, 2] × J[:, 3] + facet == 4 && return J[:, 3] × J[:, 1] + facet == 5 && return J[:, 3] × J[:, 2] + facet == 6 && return J[:, 1] × J[:, 2] end throw(ArgumentError("unknown facet number")) end @@ -198,31 +198,31 @@ end ######################### # Mapping from 2D triangle to 3D face of a tetrahedon. -function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefTetrahedron}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefTetrahedron}, facet::Int) where {T} x, y = point - face == 1 && return Vec{3, T}((one(T) - x - y, y, zero(T))) - face == 2 && return Vec{3, T}((y, zero(T), one(T) - x - y)) - face == 3 && return Vec{3, T}((x, y, one(T) - x - y)) - face == 4 && return Vec{3, T}((zero(T), one(T) - x - y, y)) + facet == 1 && return Vec{3, T}((one(T) - x - y, y, zero(T))) + facet == 2 && return Vec{3, T}((y, zero(T), one(T) - x - y)) + facet == 3 && return Vec{3, T}((x, y, one(T) - x - y)) + facet == 4 && return Vec{3, T}((zero(T), one(T) - x - y, y)) throw(ArgumentError("unknown facet number")) end # Mapping from 3D face of a tetrahedon to 2D triangle. -function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefTetrahedron}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefTetrahedron}, facet::Int) where {T} x, y, z = point - face == 1 && return Vec(one(T) - x - y, y) - face == 2 && return Vec(one(T) - z - x, x) - face == 3 && return Vec(x, y) - face == 4 && return Vec(one(T) - y - z, z) + facet == 1 && return Vec(one(T) - x - y, y) + facet == 2 && return Vec(one(T) - z - x, x) + facet == 3 && return Vec(x, y) + facet == 4 && return Vec(one(T) - y - z, z) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 3}, ::Type{RefTetrahedron}, face::Int) +function weighted_normal(J::Tensor{2, 3}, ::Type{RefTetrahedron}, facet::Int) @inbounds begin - face == 1 && return J[:, 2] × J[:, 1] - face == 2 && return J[:, 1] × J[:, 3] - face == 3 && return (J[:, 1] - J[:, 3]) × (J[:, 2] - J[:, 3]) - face == 4 && return J[:, 3] × J[:, 2] + facet == 1 && return J[:, 2] × J[:, 1] + facet == 2 && return J[:, 1] × J[:, 3] + facet == 3 && return (J[:, 1] - J[:, 3]) × (J[:, 2] - J[:, 3]) + facet == 4 && return J[:, 3] × J[:, 2] end throw(ArgumentError("unknown facet number")) end @@ -232,35 +232,35 @@ end ################### # Mapping from 2D quadrilateral/triangle to 3D face of a wedge. -function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefPrism}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefPrism}, facet::Int) where {T} # Note that for quadrilaterals the domain is [-1, 1]² but for triangles it is [0, 1]² x, y = point - face == 1 && return Vec{3, T}((one(T) - x - y, y, zero(T))) - face == 2 && return Vec{3, T}(((one(T) + x) / 2, zero(T), (one(T) + y) / 2)) - face == 3 && return Vec{3, T}((zero(T), one(T) - (one(T) + x) / 2, (one(T) + y) / 2)) - face == 4 && return Vec{3, T}((one(T) - (one(T) + x) / 2, (one(T) + x) / 2, (one(T) + y) / 2)) - face == 5 && return Vec{3, T}((y, one(T) - x - y, one(T))) + facet == 1 && return Vec{3, T}((one(T) - x - y, y, zero(T))) + facet == 2 && return Vec{3, T}(((one(T) + x) / 2, zero(T), (one(T) + y) / 2)) + facet == 3 && return Vec{3, T}((zero(T), one(T) - (one(T) + x) / 2, (one(T) + y) / 2)) + facet == 4 && return Vec{3, T}((one(T) - (one(T) + x) / 2, (one(T) + x) / 2, (one(T) + y) / 2)) + facet == 5 && return Vec{3, T}((y, one(T) - x - y, one(T))) throw(ArgumentError("unknown facet number")) end # Mapping from 3D face of a wedge to 2D triangle or 2D quadrilateral. -function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefPrism}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefPrism}, facet::Int) where {T} x, y, z = point - face == 1 && return Vec(one(T) - x - y, y) - face == 2 && return Vec(2 * x - one(T), 2 * z - one(T)) - face == 3 && return Vec(2 * (one(T) - y) - one(T), 2 * z - one(T)) - face == 4 && return Vec(2 * y - one(T), 2 * z - one(T)) - face == 5 && return Vec(one(T) - x - y, x) + facet == 1 && return Vec(one(T) - x - y, y) + facet == 2 && return Vec(2 * x - one(T), 2 * z - one(T)) + facet == 3 && return Vec(2 * (one(T) - y) - one(T), 2 * z - one(T)) + facet == 4 && return Vec(2 * y - one(T), 2 * z - one(T)) + facet == 5 && return Vec(one(T) - x - y, x) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 3}, ::Type{RefPrism}, face::Int) +function weighted_normal(J::Tensor{2, 3}, ::Type{RefPrism}, facet::Int) @inbounds begin - face == 1 && return J[:, 2] × J[:, 1] - face == 2 && return J[:, 1] × J[:, 3] - face == 3 && return J[:, 3] × J[:, 2] - face == 4 && return (J[:, 2] - J[:, 1]) × J[:, 3] - face == 5 && return J[:, 1] × J[:, 2] + facet == 1 && return J[:, 2] × J[:, 1] + facet == 2 && return J[:, 1] × J[:, 3] + facet == 3 && return J[:, 3] × J[:, 2] + facet == 4 && return (J[:, 2] - J[:, 1]) × J[:, 3] + facet == 5 && return J[:, 1] × J[:, 2] end throw(ArgumentError("unknown facet number")) end @@ -270,34 +270,34 @@ end ##################### # Mapping from 2D face to 3D face of a pyramid. -function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefPyramid}, face::Int) where {T} +function facet_to_element_transformation(point::Vec{2, T}, ::Type{RefPyramid}, facet::Int) where {T} x, y = point - face == 1 && return Vec{3, T}(((y + one(T)) / 2, (x + one(T)) / 2, zero(T))) - face == 2 && return Vec{3, T}((y, zero(T), one(T) - x - y)) - face == 3 && return Vec{3, T}((zero(T), one(T) - x - y, y)) - face == 4 && return Vec{3, T}((x + y, y, one(T) - x - y)) - face == 5 && return Vec{3, T}((one(T) - x - y, one(T) - y, y)) + facet == 1 && return Vec{3, T}(((y + one(T)) / 2, (x + one(T)) / 2, zero(T))) + facet == 2 && return Vec{3, T}((y, zero(T), one(T) - x - y)) + facet == 3 && return Vec{3, T}((zero(T), one(T) - x - y, y)) + facet == 4 && return Vec{3, T}((x + y, y, one(T) - x - y)) + facet == 5 && return Vec{3, T}((one(T) - x - y, one(T) - y, y)) throw(ArgumentError("unknown facet number")) end # Mapping from 3D face of a pyramid to 2D triangle or 2D quadrilateral. -function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefPyramid}, face::Int) where {T} +function element_to_facet_transformation(point::Vec{3, T}, ::Type{RefPyramid}, facet::Int) where {T} x, y, z = point - face == 1 && return Vec(2 * y - one(T), 2 * x - one(T)) - face == 2 && return Vec(one(T) - z - x, x) - face == 3 && return Vec(one(T) - y - z, z) - face == 4 && return Vec(x - y, y) - face == 5 && return Vec(one(T) - x - z, z) + facet == 1 && return Vec(2 * y - one(T), 2 * x - one(T)) + facet == 2 && return Vec(one(T) - z - x, x) + facet == 3 && return Vec(one(T) - y - z, z) + facet == 4 && return Vec(x - y, y) + facet == 5 && return Vec(one(T) - x - z, z) throw(ArgumentError("unknown facet number")) end -function weighted_normal(J::Tensor{2, 3}, ::Type{RefPyramid}, face::Int) +function weighted_normal(J::Tensor{2, 3}, ::Type{RefPyramid}, facet::Int) @inbounds begin - face == 1 && return J[:, 2] × J[:, 1] - face == 2 && return J[:, 1] × J[:, 3] - face == 3 && return J[:, 3] × J[:, 2] - face == 4 && return J[:, 2] × (J[:, 3] - J[:, 1]) - face == 5 && return (J[:, 3] - J[:, 2]) × J[:, 1] + facet == 1 && return J[:, 2] × J[:, 1] + facet == 2 && return J[:, 1] × J[:, 3] + facet == 3 && return J[:, 3] × J[:, 2] + facet == 4 && return J[:, 2] × (J[:, 3] - J[:, 1]) + facet == 5 && return (J[:, 3] - J[:, 2]) × J[:, 1] end throw(ArgumentError("unknown facet number")) end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index cd2b2131b5..9dd051c642 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -132,7 +132,7 @@ include("FEValues/EdgeValues.jl") include("FEValues/InterfaceValues.jl") include("FEValues/PointValues.jl") include("FEValues/common_values.jl") -#include("FEValues/face_integrals.jl") +#include("FEValues/facet_integrals.jl") include("FEValues/boundary_integrals.jl") # Grid diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 3d27431bd1..c7aa5250bd 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -528,7 +528,7 @@ Returns all nodesets of the `grid`. """ getfacetset(grid::AbstractGrid, setname::String) -Returns all faces as `FacetIndex` in the set with name `setname`. +Returns all facets as `FacetIndex` in the set with name `setname`. """ @inline getfacetset(grid::AbstractGrid, setname::String) = grid.facetsets[setname] """ @@ -648,10 +648,10 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex) Base.getindex(I::($INDEX), i::Int) = I.idx[i] - #To be able to do a,b = faceidx + #To be able to do, e.g., `a,b = facetidx` Base.iterate(I::($INDEX), state::Int = 1) = (state == 3) ? nothing : (I[state], state + 1) - # Necessary to check if, e.g. `(cellid, faceidx) in faceset` + # Necessary to check if, e.g. `(cellid, facetnr) in facetset` Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx Base.isequal(x::Tuple{Int, Int}, y::$INDEX) = x[1] == y.idx[1] && x[2] == y.idx[2] Base.isequal(y::$INDEX, x::Tuple{Int, Int}) = x[1] == y.idx[1] && x[2] == y.idx[2] diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl index 1b47c7ef50..a54839517f 100644 --- a/src/Grid/grid_generators.jl +++ b/src/Grid/grid_generators.jl @@ -27,13 +27,13 @@ function generate_grid(::Type{Line}, nel::NTuple{1, Int}, left::Vec{1, T} = Vec{ end - # Cell faces + # Cell facets boundary = [ FacetIndex(1, 1), FacetIndex(nel_x, 2), ] - # Cell face sets + # Cell facet sets facetsets = Dict( "left" => OrderedSet{FacetIndex}([boundary[1]]), "right" => OrderedSet{FacetIndex}([boundary[2]]) @@ -60,13 +60,13 @@ function generate_grid(::Type{QuadraticLine}, nel::NTuple{1, Int}, left::Vec{1, push!(cells, QuadraticLine((2 * i - 1, 2 * i + 1, 2 * i))) end - # Cell faces + # Cell facets boundary = FacetIndex[ FacetIndex(1, 1), FacetIndex(nel_x, 2), ] - # Cell face sets + # Cell facet sets facetsets = Dict( "left" => OrderedSet{FacetIndex}([boundary[1]]), "right" => OrderedSet{FacetIndex}([boundary[2]]) @@ -125,7 +125,7 @@ function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2, Int}, LL::Vec{2, T push!(cells, Quadrilateral((node_array[i, j], node_array[i + 1, j], node_array[i + 1, j + 1], node_array[i, j + 1]))) end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (nel_x, nel_y)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[:, 1]]; @@ -134,7 +134,7 @@ function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2, Int}, LL::Vec{2, T [FacetIndex(cl, 4) for cl in cell_array[1, :]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[:, 1])) .+ offset]); offset += length(cell_array[:, 1]) @@ -170,7 +170,7 @@ function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2, Int}, LL:: push!(cells, cell) end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (nel_x, nel_y)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[:, 1]]; @@ -179,7 +179,7 @@ function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2, Int}, LL:: [FacetIndex(cl, 4) for cl in cell_array[1, :]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[:, 1])) .+ offset]); offset += length(cell_array[:, 1]) @@ -219,7 +219,7 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3, Int}, left::Vec{3, T} push!(cells, cell) end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (nel_x, nel_y, nel_z)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[:, :, 1][:]]; @@ -230,7 +230,7 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3, Int}, left::Vec{3, T} [FacetIndex(cl, 6) for cl in cell_array[:, :, end][:]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[:, :, 1][:])) .+ offset]); offset += length(cell_array[:, :, 1][:]) @@ -403,7 +403,7 @@ function generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTuple{3, In push!(cells, cell) end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (nel_x, nel_y, nel_z)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[:, :, 1][:]]; @@ -414,7 +414,7 @@ function generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTuple{3, In [FacetIndex(cl, 6) for cl in cell_array[:, :, end][:]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[:, :, 1][:])) .+ offset]); offset += length(cell_array[:, :, 1][:]) @@ -446,7 +446,7 @@ function generate_grid(::Type{Triangle}, nel::NTuple{2, Int}, LL::Vec{2, T}, LR: push!(cells, Triangle((node_array[i + 1, j], node_array[i + 1, j + 1], node_array[i, j + 1]))) # ◹ end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (2, nel_x, nel_y)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[1, :, 1]]; @@ -455,7 +455,7 @@ function generate_grid(::Type{Triangle}, nel::NTuple{2, Int}, LL::Vec{2, T}, LR: [FacetIndex(cl, 3) for cl in cell_array[1, 1, :]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[1, :, 1])) .+ offset]); offset += length(cell_array[1, :, 1]) @@ -499,7 +499,7 @@ function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2, Int}, LL::Vec{2 push!(cells, triangle2) end - # Cell faces + # Cell facets cell_array = reshape(collect(1:nel_tot), (2, nel_x, nel_y)) boundary = FacetIndex[ [FacetIndex(cl, 1) for cl in cell_array[1, :, 1]]; @@ -508,7 +508,7 @@ function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2, Int}, LL::Vec{2 [FacetIndex(cl, 3) for cl in cell_array[1, 1, :]] ] - # Cell face sets + # Cell facet sets offset = 0 facetsets = Dict{String, OrderedSet{FacetIndex}}() facetsets["bottom"] = OrderedSet{FacetIndex}(boundary[(1:length(cell_array[1, :, 1])) .+ offset]); offset += length(cell_array[1, :, 1]) diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl index ce5937010e..984c8a00e2 100644 --- a/src/Grid/utils.jl +++ b/src/Grid/utils.jl @@ -93,7 +93,7 @@ addboundaryvertexset!(grid::AbstractGrid, topology::ExclusiveTopology, name::Str Adds a boundary vertexset to the grid with key `name`. A vertexset maps a `String` key to an `OrderedSet` of tuples corresponding to `(global_cell_id, local_vertex_id)`. `all=true` implies that `f(x)` must return `true` for all nodal -coordinates `x` on the face if the face should be added to the set, otherwise it suffices +coordinates `x` on the facet if the facet should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. """ function addboundaryvertexset!(grid::AbstractGrid, top::ExclusiveTopology, name::String, f::Function; kwargs...) diff --git a/src/Quadrature/quadrature.jl b/src/Quadrature/quadrature.jl index b4368a15f3..21be7cf16e 100644 --- a/src/Quadrature/quadrature.jl +++ b/src/Quadrature/quadrature.jl @@ -173,29 +173,29 @@ end """ FacetQuadratureRule{shape}([::Type{T},] [quad_rule_type::Symbol,] order::Int) - FacetQuadratureRule{shape}(face_rules::NTuple{<:Any, <:QuadratureRule{shape}}) - FacetQuadratureRule{shape}(face_rules::AbstractVector{<:QuadratureRule{shape}}) + FacetQuadratureRule{shape}(facet_rules::NTuple{<:Any, <:QuadratureRule{shape}}) + FacetQuadratureRule{shape}(facet_rules::AbstractVector{<:QuadratureRule{shape}}) -Create a `FacetQuadratureRule` used for integration of the faces of the refshape `shape` (of +Create a `FacetQuadratureRule` used for integration of the facets of the refshape `shape` (of type [`AbstractRefShape`](@ref)). `order` is the order of the quadrature rule. If no symbol is provided, the default `quad_rule_type` for each facet's reference shape is used (see [QuadratureRule](@ref)). For non-default `quad_rule_type`s on cells with mixed facet types (e.g. `RefPrism` and `RefPyramid`), the -`face_rules` must be provided explicitly. +`facet_rules` must be provided explicitly. `FacetQuadratureRule` is used as one of the components to create [`FacetValues`](@ref). """ struct FacetQuadratureRule{shape, FacetRulesType} - face_rules::FacetRulesType # E.g. Tuple{QuadratureRule{RefLine,...}, QuadratureRule{RefLine,...}} - function FacetQuadratureRule{shape}(face_rules::Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}) where {shape, QRType <: QuadratureRule{shape}} - if length(face_rules) != nfacets(shape) - throw(ArgumentError("number of quadrature rules does not not match number of facets (#rules=$(length(face_rules)) != #facets=$(nfacets(shape)))")) + facet_rules::FacetRulesType # E.g. Tuple{QuadratureRule{RefLine,...}, QuadratureRule{RefLine,...}} + function FacetQuadratureRule{shape}(facet_rules::Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}) where {shape, QRType <: QuadratureRule{shape}} + if length(facet_rules) != nfacets(shape) + throw(ArgumentError("number of quadrature rules does not not match number of facets (#rules=$(length(facet_rules)) != #facets=$(nfacets(shape)))")) end - return new{shape, typeof(face_rules)}(face_rules) + return new{shape, typeof(facet_rules)}(facet_rules) end end -function FacetQuadratureRule(face_rules::Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}) where {shape, QRType <: QuadratureRule{shape}} - return FacetQuadratureRule{shape}(face_rules) +function FacetQuadratureRule(facet_rules::Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}) where {shape, QRType <: QuadratureRule{shape}} + return FacetQuadratureRule{shape}(facet_rules) end # Fill in defaults T=Float64 @@ -206,8 +206,8 @@ function FacetQuadratureRule{shape}(quad_type::Symbol, order::Int) where {shape return FacetQuadratureRule{shape}(Float64, quad_type, order) end -# For RefShapes with equal face-shapes: generate quad rule for the face shape -# and expand to each face +# For RefShapes with equal facet shapes: generate quad rule for the facet shape +# and expand to each facet function FacetQuadratureRule{RefLine}(::Type{T}, ::Int) where {T} w, p = T[1], Vec{0, T}[Vec{0, T}(())] return create_facet_quad_rule(RefLine, w, p) @@ -306,16 +306,16 @@ Return the number of quadrature points in `qr`. getnquadpoints(qr::QuadratureRule) = length(getweights(qr)) """ - getnquadpoints(qr::FacetQuadratureRule, face::Int) + getnquadpoints(qr::FacetQuadratureRule, facet::Int) -Return the number of quadrature points in `qr` for local face index `face`. +Return the number of quadrature points in `qr` for local facet index `facet`. """ -getnquadpoints(qr::FacetQuadratureRule, face::Int) = getnquadpoints(qr.face_rules[face]) +getnquadpoints(qr::FacetQuadratureRule, facet::Int) = getnquadpoints(qr.facet_rules[facet]) getnquadpoints(qr::EdgeQuadratureRule, edge::Int) = getnquadpoints(qr.edge_rules[edge]) """ getweights(qr::QuadratureRule) - getweights(qr::FacetQuadratureRule, face::Int) + getweights(qr::FacetQuadratureRule, facet::Int) Return the weights of the quadrature rule. @@ -331,12 +331,12 @@ julia> getweights(qr) ``` """ getweights(qr::QuadratureRule) = qr.weights -getweights(qr::FacetQuadratureRule, face::Int) = getweights(qr.face_rules[face]) +getweights(qr::FacetQuadratureRule, facet::Int) = getweights(qr.facet_rules[facet]) getweights(qr::EdgeQuadratureRule, edge::Int) = getweights(qr.edge_rules[edge]) """ getpoints(qr::QuadratureRule) - getpoints(qr::FacetQuadratureRule, face::Int) + getpoints(qr::FacetQuadratureRule, facet::Int) Return the points of the quadrature rule. @@ -352,10 +352,10 @@ julia> getpoints(qr) ``` """ getpoints(qr::QuadratureRule) = qr.points -getpoints(qr::FacetQuadratureRule, face::Int) = getpoints(qr.face_rules[face]) +getpoints(qr::FacetQuadratureRule, facet::Int) = getpoints(qr.facet_rules[facet]) getpoints(qr::EdgeQuadratureRule, edge::Int) = getpoints(qr.edge_rules[edge]) getrefshape(::QuadratureRule{RefShape}) where {RefShape} = RefShape -# TODO: This is used in copy(::(Cell|Face)Values), but it it useful to get an actual copy? +# TODO: This is used in copy(::(Cell|Facet)Values), but it it useful to get an actual copy? Base.copy(qr::Union{QuadratureRule, FacetQuadratureRule, EdgeQuadratureRule}) = qr diff --git a/src/iterators.jl b/src/iterators.jl index 2b9613accf..64b38f2e0e 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -101,11 +101,11 @@ nfacets(cc::CellCache) = nfacets(getcells(cc.grid, cc.cellid)) FacetCache(dh::AbstractDofHandler) Create a cache object with pre-allocated memory for the nodes, coordinates, and dofs of a -cell suitable for looping over *faces* in a grid. The cache is updated for a new face by +cell suitable for looping over *facets* in a grid. The cache is updated for a new facet by calling `reinit!(cache, fi::FacetIndex)`. **Methods with `fc::FacetCache`** - - `reinit!(fc, fi)`: reinitialize the cache for face `fi::FacetIndex` + - `reinit!(fc, fi)`: reinitialize the cache for facet `fi::FacetIndex` - `cellid(fc)`: get the current cellid - `getnodes(fc)`: get the global node ids of the *cell* - `getcoordinates(fc)`: get the coordinates of the *cell* @@ -150,11 +150,11 @@ end Create a cache object with pre-allocated memory for the nodes, coordinates, and dofs of an interface. The cache is updated for a new cell by calling `reinit!(cache, facet_a, facet_b)` where -`facet_a::FacetIndex` and `facet_b::FacetIndex` are the two interface faces. +`facet_a::FacetIndex` and `facet_b::FacetIndex` are the two interface facets. **Struct fields of `InterfaceCache`** - - `ic.a :: FacetCache`: face cache for the first face of the interface - - `ic.b :: FacetCache`: face cache for the second face of the interface + - `ic.a :: FacetCache`: facet cache for the first facet of the interface + - `ic.b :: FacetCache`: facet cache for the second facet of the interface - `ic.dofs :: Vector{Int}`: global dof ids for the interface (union of `ic.a.dofs` and `ic.b.dofs`) **Methods with `InterfaceCache`** @@ -327,12 +327,13 @@ end ``` is thus simply convenience for the following equivalent snippet for grids of dimensions > 1: ```julia -ic = InterfaceCache(grid, topology) -for face in topology.face_skeleton - neighborhood = topology.face_face_neighbor[face[1], face[2]] - isempty(neighborhood) && continue - neighbor_face = neighborhood[1] - reinit!(ic, face, neighbor_face) +ic = InterfaceCache(grid) +neighborhood = Ferrite.get_facet_facet_neighborhood(topology, grid) +for facet in facetskeleton(topology, grid) + neighbors = neighborhood[facet[1], facet[2]] + isempty(neighbors) && continue + neighbor_facet = neighbors[1] + reinit!(ic, facet, neighbor_facet) # ... end ``` diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 672713b315..cd77c95024 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -185,7 +185,7 @@ @test contains(showstring, "Function interpolation: Lagrange{RefQuadrilateral, 2}()") @test contains(showstring, "Geometric interpolation: Lagrange{RefQuadrilateral, 1}()^2") fv2 = copy(fv) - push!(Ferrite.getweights(fv2.fqr.face_rules[1]), 1) + push!(Ferrite.getweights(fv2.fqr.facet_rules[1]), 1) showstring = sprint(show, MIME"text/plain"(), fv2) @test startswith(showstring, "FacetValues(scalar, rdim=2, sdim=2): (3, 2, 2, 2) quadrature points on each face") end