Skip to content

Commit

Permalink
Merge master
Browse files Browse the repository at this point in the history
  • Loading branch information
KnutAM committed Jan 7, 2025
2 parents 368e277 + 951ec47 commit ac1d277
Show file tree
Hide file tree
Showing 10 changed files with 194 additions and 193 deletions.
36 changes: 18 additions & 18 deletions src/FEValues/FacetValues.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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: ")
Expand Down Expand Up @@ -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
22 changes: 11 additions & 11 deletions src/FEValues/InterfaceValues.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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}
Expand All @@ -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}
Expand All @@ -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

Expand Down Expand Up @@ -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
```
Expand Down
Loading

0 comments on commit ac1d277

Please sign in to comment.