Skip to content

Commit

Permalink
Adopt proper codim/dim terminology for entities (#789, #914)
Browse files Browse the repository at this point in the history
This patch changes the naming and meaning of various entities in order
to be independent of the current dimension of the problem/grid/element
type.

After this patch, the following terminology is used to name entities
based on their dimension:

 - `vertex`: a zero-dimensional entity (point, node)
 - `edge`: a one-dimensional entity between two vertices
 - `face`: a two-dimensional entity between edges
 - `volume`: a three-dimensional entity between faces

and the following terminology is used to name entities based on their
codimension:
 - `cell`: codimension zero (same meaning as before this patch)
 - `facet`: codimension one

Note that the codimension is defined relative the *reference dimension*
of the cell. For example, in a problem in three spatial dimensions with
both three-dimensional hexahedrons and two-dimensional quadrilaterals,
the facets of the hexahedron are faces, and the facets of the
quadrilateral are edges.

These changes simplifies internal handling of degree-of-freedom
distribution, and simplifies handling of boundary conditions (which are
now consistently dealing with *facets* instead of *faces*, *edges*, and
*vertices for three-, two-, and one-dimensional problems respectively.

User-facing changes:

 - Face-sets in the `Grid` type are deprecated in favor of facet-sets:
   (`getfaceset` is replaced by `getfacetset`, `addfaceset!` is replaced
   by `addfacetset!`).
 - Facet-sets should be used instead of face-sets when defining boundary
   conditions (such as `Dirichlet`).
 - `FaceValues` have been renamed to `FacetValues`.
 - `FaceQuadratureRule` have been renamed to `FacetQuadratureRule`.
 - `nfaces` should be replaced with `nfacets` (when looping over the
   facets of an element to integrate Neumann boundary conditions, for
   example).

Fixes #899, fixes #901. Closes #789, closes #914.

Co-authored-by: Elias <elias.isak.borjesson@gmail.com>
Co-authored-by: Knut Andreas Meyer <knutam@gmail.com>
  • Loading branch information
2 people authored and fredrikekre committed May 19, 2024
1 parent 9fd9d34 commit 2ca4f94
Show file tree
Hide file tree
Showing 66 changed files with 2,068 additions and 2,292 deletions.
4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -33,13 +33,13 @@ WriteVTK = "1.13"
julia = "1.6"

[extras]
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464"
Expand All @@ -50,4 +50,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

[targets]
test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "JET", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"]
test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "Pkg", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"]
10 changes: 5 additions & 5 deletions benchmark/benchmarks-assembly.jl
Original file line number Diff line number Diff line change
Expand Up @@ -61,16 +61,16 @@ for spatial_dim ∈ 1:3
LAGRANGE_SUITE["petrov-galerkin"]["pressure-velocity"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $cvv, shape_divergence, $csv, shape_value, *)

if spatial_dim > 1
qr_face = FaceQuadratureRule{ref_type}(2*order-1)
fsv = FaceValues(qr_face, ip, ip_geo);
fsv2 = FaceValues(qr_face, ip, ip_geo);
qr_facet = FacetQuadratureRule{ref_type}(2*order-1)
fsv = FacetValues(qr_facet, ip, ip_geo);
fsv2 = FacetValues(qr_facet, ip, ip_geo);

LAGRANGE_SUITE["ritz-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, shape_value, *)
LAGRANGE_SUITE["petrov-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, $fsv2, shape_value, *)

ip = DiscontinuousLagrange{ref_type, order}()
isv = InterfaceValues(qr_face, ip, ip_geo);
isv2 = InterfaceValues(qr_face, ip, ip_geo);
isv = InterfaceValues(qr_facet, ip, ip_geo);
isv2 = InterfaceValues(qr_facet, ip, ip_geo);
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
Expand Down
2 changes: 1 addition & 1 deletion benchmark/benchmarks-boundary-conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ for spatial_dim ∈ [2]
close!(dh);

ch = ConstraintHandler(dh);
∂Ω = union(getfaceset.((grid, ), ["left"])...);
∂Ω = union(getfacetset.((grid, ), ["left"])...);
dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
add!(ch, dbc);
close!(ch);
Expand Down
40 changes: 20 additions & 20 deletions benchmark/helper.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,22 +40,22 @@ function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.Abstract
Ke
end

function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facevalues::FaceValues, f_shape, f_test, op)
n_basefuncs = getnbasefunctions(facevalues)
function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facetvalues::FacetValues, f_shape, f_test, op)
n_basefuncs = getnbasefunctions(facetvalues)

f = zeros(n_basefuncs)

X = getcoordinates(grid, 1)
for face in 1:nfaces(getcells(grid)[1])
reinit!(facevalues, X, face)
for facet in 1:nfacets(getcells(grid)[1])
reinit!(facetvalues, X, facet)

for q_point in 1:getnquadpoints(facevalues)
n = getnormal(facevalues, q_point)
= getdetJdV(facevalues, q_point)
for q_point in 1:getnquadpoints(facetvalues)
n = getnormal(facetvalues, q_point)
= getdetJdV(facetvalues, q_point)
for i in 1:n_basefuncs
test = f_test(facevalues, q_point, i)
test = f_test(facetvalues, q_point, i)
for j in 1:n_basefuncs
shape = f_shape(facevalues, q_point, j)
shape = f_shape(facetvalues, q_point, j)
f[i] += op(test, shape) n *
end
end
Expand Down Expand Up @@ -118,25 +118,25 @@ function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.Abstra
Ke
end

function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facevalues_shape::FaceValues, f_shape, facevalues_test::FaceValues, f_test, op)
n_basefuncs_shape = getnbasefunctions(facevalues_shape)
n_basefuncs_test = getnbasefunctions(facevalues_test)
function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facetvalues_shape::FacetValues, f_shape, facetvalues_test::FacetValues, f_test, op)
n_basefuncs_shape = getnbasefunctions(facetvalues_shape)
n_basefuncs_test = getnbasefunctions(facetvalues_test)

f = zeros(n_basefuncs_test)

X_shape = getcoordinates(grid, 1)
X_test = getcoordinates(grid, 1)
for face in 1:nfaces(getcells(grid)[1])
reinit!(facevalues_shape, X_shape, face)
reinit!(facevalues_test, X_test, face)
for facet in 1:nfacets(getcells(grid)[1])
reinit!(facetvalues_shape, X_shape, facet)
reinit!(facetvalues_test, X_test, facet)

for q_point in 1:getnquadpoints(facevalues_shape)
n = getnormal(facevalues_test, q_point)
= getdetJdV(facevalues_test, q_point)
for q_point in 1:getnquadpoints(facetvalues_shape)
n = getnormal(facetvalues_test, q_point)
= getdetJdV(facetvalues_test, q_point)
for i in 1:n_basefuncs_test
test = f_test(facevalues_test, q_point, i)
test = f_test(facetvalues_test, q_point, i)
for j in 1:n_basefuncs_shape
shape = f_shape(facevalues_shape, q_point, j)
shape = f_shape(facetvalues_shape, q_point, j)
f[i] += op(test, shape) n *
end
end
Expand Down
4 changes: 2 additions & 2 deletions docs/src/devdocs/FEValues.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
* `AbstractValues`
* `AbstractCellValues`
* [`CellValues`](@ref)
* `AbstractFaceValues`
* [`FaceValues`](@ref)
* `AbstractFacetValues`
* [`FacetValues`](@ref)
* [`BCValues`](@ref)
* [`PointValues`](@ref)

Expand Down
13 changes: 2 additions & 11 deletions docs/src/devdocs/elements.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,17 +26,8 @@ Ferrite.get_coordinate_eltype(::Node)
Ferrite.toglobal
Ferrite.sortface
Ferrite.sortedge
Ferrite.getfaceedges
Ferrite.getfacevertices
Ferrite.getedgevertices
Ferrite.getfaceinstances
Ferrite.getedgeinstances
Ferrite.getvertexinstances
Ferrite.filterfaces
Ferrite.filteredges
Ferrite.filtervertices
Ferrite.element_to_face_transformation
Ferrite.face_to_element_transformation
Ferrite.element_to_facet_transformation
Ferrite.facet_to_element_transformation
Ferrite.InterfaceOrientationInfo
Ferrite.transform_interface_points!
```
2 changes: 1 addition & 1 deletion docs/src/devdocs/interpolations.md
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ Ferrite.facedof_interior_indices(::Interpolation)
Ferrite.edgedof_indices(::Interpolation)
Ferrite.dirichlet_edgedof_indices(::Interpolation)
Ferrite.edgedof_interior_indices(::Interpolation)
Ferrite.celldof_interior_indices(::Interpolation)
Ferrite.volumedof_interior_indices(::Interpolation)
Ferrite.getnbasefunctions(::Interpolation)
Ferrite.reference_coordinates(::Interpolation)
Ferrite.is_discontinuous(::Interpolation)
Expand Down
30 changes: 15 additions & 15 deletions docs/src/literate-gallery/helmholtz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,9 @@ grid = generate_grid(Quadrilateral, (150, 150))

ip = Lagrange{RefQuadrilateral, 1}()
qr = QuadratureRule{RefQuadrilateral}(2)
qr_face = FaceQuadratureRule{RefQuadrilateral}(2)
qr_facet = FacetQuadratureRule{RefQuadrilateral}(2)
cellvalues = CellValues(qr, ip);
facevalues = FaceValues(qr_face, ip);
facetvalues = FacetValues(qr_facet, ip);

dh = DofHandler(grid)
add!(dh, :u, ip)
Expand All @@ -80,14 +80,14 @@ end;

dbcs = ConstraintHandler(dh)
# The (strong) Dirichlet boundary condition can be handled automatically by the Ferrite library.
dbc = Dirichlet(:u, union(getfaceset(grid, "top"), getfaceset(grid, "right")), (x,t) -> u_ana(x))
dbc = Dirichlet(:u, union(getfacetset(grid, "top"), getfacetset(grid, "right")), (x,t) -> u_ana(x))
add!(dbcs, dbc)
close!(dbcs)
update!(dbcs, 0.0)

K = create_sparsity_pattern(dh);

function doassemble(cellvalues::CellValues, facevalues::FaceValues,
function doassemble(cellvalues::CellValues, facetvalues::FacetValues,
K::SparseMatrixCSC, dh::DofHandler)
b = 1.0
f = zeros(ndofs(dh))
Expand Down Expand Up @@ -134,18 +134,18 @@ function doassemble(cellvalues::CellValues, facevalues::FaceValues,
# \int_{\Gamma_2} δu g_2 \, d\Gamma
# ```
#+
for face in 1:nfaces(cell)
if onboundary(cell, face) &&
((cellcount, face) getfaceset(grid, "left") ||
(cellcount, face) getfaceset(grid, "bottom"))
reinit!(facevalues, cell, face)
for q_point in 1:getnquadpoints(facevalues)
coords_qp = spatial_coordinate(facevalues, q_point, coords)
n = getnormal(facevalues, q_point)
for facet in 1:nfacets(cell)
if onboundary(cell, facet) &&
((cellcount, facet) getfacetset(grid, "left") ||
(cellcount, facet) getfacetset(grid, "bottom"))
reinit!(facetvalues, cell, facet)
for q_point in 1:getnquadpoints(facetvalues)
coords_qp = spatial_coordinate(facetvalues, q_point, coords)
n = getnormal(facetvalues, q_point)
g_2 = gradient(u_ana, coords_qp) n
= getdetJdV(facevalues, q_point)
= getdetJdV(facetvalues, q_point)
for i in 1:n_basefuncs
δu = shape_value(facevalues, q_point, i)
δu = shape_value(facetvalues, q_point, i)
fe[i] += (δu * g_2) *
end
end
Expand All @@ -158,7 +158,7 @@ function doassemble(cellvalues::CellValues, facevalues::FaceValues,
return K, f
end;

K, f = doassemble(cellvalues, facevalues, K, dh);
K, f = doassemble(cellvalues, facetvalues, K, dh);
apply!(K, f, dbcs)
u = Symmetric(K) \ f;

Expand Down
4 changes: 2 additions & 2 deletions docs/src/literate-gallery/landau.jl
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,8 @@ function LandauModel(α, G, gridsize, left::Vec{DIM, T}, right::Vec{DIM, T}, elp
startingconditions!(dofvector, dofhandler)
boundaryconds = ConstraintHandler(dofhandler)
#boundary conditions can be added but aren't necessary for optimization
#add!(boundaryconds, Dirichlet(:P, getfaceset(grid, "left"), (x, t) -> [0.0,0.0,0.53], [1,2,3]))
#add!(boundaryconds, Dirichlet(:P, getfaceset(grid, "right"), (x, t) -> [0.0,0.0,-0.53], [1,2,3]))
#add!(boundaryconds, Dirichlet(:P, getfacetset(grid, "left"), (x, t) -> [0.0,0.0,0.53], [1,2,3]))
#add!(boundaryconds, Dirichlet(:P, getfacetset(grid, "right"), (x, t) -> [0.0,0.0,-0.53], [1,2,3]))
close!(boundaryconds)
update!(boundaryconds, 0.0)

Expand Down
26 changes: 13 additions & 13 deletions docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,10 +86,10 @@ end
# to later assign Dirichlet boundary conditions
function importTestGrid()
grid = generate_grid(Tetrahedron, (5, 5, 5), zero(Vec{3}), ones(Vec{3}));
addfaceset!(grid, "myBottom", x -> norm(x[2]) 0.0);
addfaceset!(grid, "myBack", x -> norm(x[3]) 0.0);
addfaceset!(grid, "myRight", x -> norm(x[1]) 1.0);
addfaceset!(grid, "myLeft", x -> norm(x[1]) 0.0);
addfacetset!(grid, "myBottom", x -> norm(x[2]) 0.0);
addfacetset!(grid, "myBack", x -> norm(x[3]) 0.0);
addfacetset!(grid, "myRight", x -> norm(x[1]) 1.0);
addfacetset!(grid, "myLeft", x -> norm(x[1]) 0.0);
return grid
end;

Expand All @@ -98,16 +98,16 @@ end;
function create_values(interpolation_u, interpolation_p)
## quadrature rules
qr = QuadratureRule{RefTetrahedron}(4)
face_qr = FaceQuadratureRule{RefTetrahedron}(4)
facet_qr = FacetQuadratureRule{RefTetrahedron}(4)

## cell and facevalues for u
## cell and facetvalues for u
cellvalues_u = CellValues(qr, interpolation_u)
facevalues_u = FaceValues(face_qr, interpolation_u)
facetvalues_u = FacetValues(facet_qr, interpolation_u)

## cellvalues for p
cellvalues_p = CellValues(qr, interpolation_p)

return cellvalues_u, cellvalues_p, facevalues_u
return cellvalues_u, cellvalues_p, facetvalues_u
end;

# We now create the function for Ψ*
Expand Down Expand Up @@ -145,10 +145,10 @@ end;
# of the loading.
function create_bc(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myLeft"), (x,t) -> zero(Vec{1}), [1]))
add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBottom"), (x,t) -> zero(Vec{1}), [2]))
add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBack"), (x,t) -> zero(Vec{1}), [3]))
add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myRight"), (x,t) -> t*ones(Vec{1}), [1]))
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myLeft"), (x,t) -> zero(Vec{1}), [1]))
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBottom"), (x,t) -> zero(Vec{1}), [2]))
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBack"), (x,t) -> zero(Vec{1}), [3]))
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myRight"), (x,t) -> t*ones(Vec{1}), [1]))
close!(dbc)
Ferrite.update!(dbc, 0.0)
return dbc
Expand Down Expand Up @@ -284,7 +284,7 @@ function solve(interpolation_u, interpolation_p)

## Create the DofHandler and CellValues
dh = create_dofhandler(grid, interpolation_u, interpolation_p)
cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)
cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)

## Create the DirichletBCs
dbc = create_bc(dh)
Expand Down
Loading

0 comments on commit 2ca4f94

Please sign in to comment.