Skip to content

Commit

Permalink
Ferrite facet definition and OrderedSets (#26)
Browse files Browse the repository at this point in the history
  • Loading branch information
KnutAM authored May 21, 2024
1 parent 647fdbc commit 00c5ae0
Show file tree
Hide file tree
Showing 40 changed files with 377 additions and 380 deletions.
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,8 @@ At its core, `FerriteAssembly.jl` lets you overload some of the following functi
```julia
element_routine!(Ke, re, state, ae, material, cellvalues, cellbuffer)
element_residual!(re, state, ae, material, cellvalues, cellbuffer)
face_routine!(Ke, re, ae, material, facevalues, facebuffer)
face_residual!(re, ae, material, facevalues, facebuffer)
facet_routine!(Ke, re, ae, material, facetvalues, facetbuffer)
facet_residual!(re, ae, material, facetvalues, facetbuffer)
```
define a `domainbuffer` with `setup_domainbuffer`,

Expand Down
8 changes: 4 additions & 4 deletions docs/figure_factory/grid_domains.jl
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,9 @@ function save_as_vtk()
end
end

vtk_grid(joinpath(@__DIR__, "..", "..", "DomainSpecs"), grid) do vtk
vtk_cell_data(vtk, cell_type_data, "cell_type_data")
vtk_cell_data(vtk, cell_physical_domains, "cell_physical_domains")
vtk_cell_data(vtk, cell_work_domains, "cell_work_domains")
VTKFile(joinpath(@__DIR__, "..", "..", "DomainSpecs"), grid) do vtk
write_cell_data(vtk, cell_type_data, "cell_type_data")
write_cell_data(vtk, cell_physical_domains, "cell_physical_domains")
write_cell_data(vtk, cell_work_domains, "cell_work_domains")
end
end
2 changes: 1 addition & 1 deletion docs/src/Customization.md
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ The interface for creating an assembler is that the assembler must support the `
as well as the methods,

* `FerriteAssembly.work_single_cell!(assembler, buffer::AbstractCellBuffer)`
* `FerriteAssembly.work_single_face!(assembler, buffer::FaceBuffer)`
* `FerriteAssembly.work_single_facet!(assembler, buffer::FacetBuffer)`

The internal methods for builtin assemblers can be used as examples.

Expand Down
12 changes: 6 additions & 6 deletions docs/src/DomainBuffers/ItemBuffer.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,9 @@ CurrentModule = FerriteAssembly

# AbstractItemBuffer
Depending on the domain that is worked over, different item buffers are available,
e.g. `CellBuffer` and `FaceBuffer`. These are set up during call to `setup_domainbuffer`.
e.g. `CellBuffer` and `FacetBuffer`. These are set up during call to `setup_domainbuffer`.

For each item (cell or face), the values in the buffer are updated to the current item,
For each item (cell or facet), the values in the buffer are updated to the current item,
and can be accessed with the following functions:

```@docs
Expand All @@ -29,9 +29,9 @@ To allocate the user cache, overload `allocate_cell_cache`:
allocate_cell_cache(::Any, ::Any)
```

## FaceBuffer
The `FaceBuffer` is similar to the `CellBuffer`, except that it has `FaceValues` instead of `CellValues`
and that state variables are not supported. To allocate user cache, overload `allocate_face_cache`:
## FacetBuffer
The `FacetBuffer` is similar to the `CellBuffer`, except that it has `FacetValues` instead of `CellValues`
and that state variables are not supported. To allocate user cache, overload `allocate_facet_cache`:
```@docs
allocate_face_cache
allocate_facet_cache
```
16 changes: 8 additions & 8 deletions docs/src/Workers/Assemblers.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
Assemblers are used to assemble the system matrices and residual vectors. To calculate the contribution to these matrices, one or more of the following functions should be overloaded for the specific `material`:
* [`element_routine!`](@ref FerriteAssembly.element_routine!)
* [`element_residual!`](@ref FerriteAssembly.element_residual!)
* [`face_residual!`](@ref FerriteAssembly.face_residual!)
* [`face_routine!`](@ref FerriteAssembly.face_routine!)
* [`facet_residual!`](@ref FerriteAssembly.facet_residual!)
* [`facet_routine!`](@ref FerriteAssembly.facet_routine!)


## Available assemblers
Expand All @@ -13,17 +13,17 @@ The following assemblers can be used to assemble the system matrix and vector:
- [`ReAssembler`](@ref)

### Ferrite assemblers
Ferrite assemblers are supported and are used to assemble both a matrix and a vector. They are created by calling Ferrite's `start_assemble` function. They will primarily call `element_routine!` or `face_routine!`.
However, if not defined for the given material, automatic differentiation will be used to get the matrix contribution by calling `element_residual!` or `face_residual!` instead.
Ferrite assemblers are supported and are used to assemble both a matrix and a vector. They are created by calling Ferrite's `start_assemble` function. They will primarily call `element_routine!` or `facet_routine!`.
However, if not defined for the given material, automatic differentiation will be used to get the matrix contribution by calling `element_residual!` or `facet_residual!` instead.

### ReAssembler
The `ReAssembler` assembles only the residual vector, and will thus call `element_residual!` or `face_residual!`.
The `ReAssembler` assembles only the residual vector, and will thus call `element_residual!` or `facet_residual!`.
```@docs
ReAssembler
```

### KeReAssembler
The `KeReAssembler` assembles both the residual vector and the system matrix, similar to Ferrite's assemblers. It will also primarily call `element_routine!` or `face_routine!` and fall back to automatic differentiation of `element_residual!` and `face_residual!` if the former methods are not defined.
The `KeReAssembler` assembles both the residual vector and the system matrix, similar to Ferrite's assemblers. It will also primarily call `element_routine!` or `facet_routine!` and fall back to automatic differentiation of `element_residual!` and `facet_residual!` if the former methods are not defined.
```@docs
KeReAssembler
```
Expand All @@ -32,8 +32,8 @@ KeReAssembler
```@docs
FerriteAssembly.element_routine!
FerriteAssembly.element_residual!
FerriteAssembly.face_routine!
FerriteAssembly.face_residual!
FerriteAssembly.facet_routine!
FerriteAssembly.facet_residual!
```

## Residual scaling
Expand Down
2 changes: 1 addition & 1 deletion docs/src/Workers/Integrators.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,5 +21,5 @@ Using the `Integrator` requires the following functions
to be overloaded, depending on the type of domain
```@docs
FerriteAssembly.integrate_cell!
FerriteAssembly.integrate_face!
FerriteAssembly.integrate_facet!
```
2 changes: 1 addition & 1 deletion docs/src/Workers/Workers.md
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,6 @@ work!

A worker needs to define what to do for different domains,
specifically it needs to define
`work_single_cell!` and `work_single_face!`. Further details
`work_single_cell!` and `work_single_facet!`. Further details
about writing custom workers can be found in the
[Customizations](@ref Assembler-interface) section.
8 changes: 4 additions & 4 deletions docs/src/design.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,10 @@ The worker is the simplest object, and determines what is done for each item in

## Domain buffer
Two domain buffer types are implemented, a regular and a threaded version. To understand the structure, it sufficies to to understand the regular buffer. The threaded version just contains additional values to have per-task buffers etc.
The key fields of the `DomainBuffer` type are (1) a `set` of items to iterate over (for example a vector of `Int` for a list of cells or a vector of `FaceIndex` for a list of faces) and an `itembuffer` (e.g., `CellBuffer` or `FaceBuffer`) containing the data required to do something for each entity.
The key fields of the `DomainBuffer` type are (1) a `set` of items to iterate over (for example a vector of `Int` for a list of cells or a vector of `FacetIndex` for a list of facets) and an `itembuffer` (e.g., `CellBuffer` or `FacetBuffer`) containing the data required to do something for each entity.

### ItemBuffer
The most important fields in an item buffer are the user-defined `material` and the `FEValues` (e.g. `CellValues` or `FaceValues`).
The most important fields in an item buffer are the user-defined `material` and the `FEValues` (e.g. `CellValues` or `FacetValues`).

### Requirements for being a domain
In addition to having the same `FEValues` and `material`,
Expand Down Expand Up @@ -48,10 +48,10 @@ several functions that happen at the innermost part of the iteration of items in
* `element_residual!`: Calculate `re`, and `state`
* `create_cell_state`: For the given material and initial cell dof values, return the initial value of the cell state. Typically one value per integration point.
* `allocate_cell_cache`: Define workspace for calculation on a single cell
* `allocate_face_cache`: Same as above, but for a face
* `allocate_facet_cache`: Same as above, but for a facet

### State variables
State variables (both current and old) are stored as a `Dict{Int}` in the domain buffer (for cells, not for faces, this may be added in the future).
State variables (both current and old) are stored as a `Dict{Int}` in the domain buffer (for cells, not for facets, this may be added in the future).
The key corresponds to the cell id, and the state can be gotten with the
`get_state` and `get_old_state` functions.

Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ CurrentModule = FerriteAssembly
The goal of [FerriteAssembly](https://github.com/KnutAM/FerriteAssembly.jl)
is make it convenient to iterate over a grid and perform some task for each item.
The primary use is to iterate over all cells and assemble into global matrices and
vectors. However, it makes many other tasks easy as well, for example iterating over faces for adding boundary conditions, or integrating a function over a given domain.
vectors. However, it makes many other tasks easy as well, for example iterating over facets for adding boundary conditions, or integrating a function over a given domain.

The design is built around two main types of objects
1. **Workers**: What to do for a given item
Expand Down
6 changes: 3 additions & 3 deletions docs/src/internals.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,15 @@ FerriteAssembly._copydofs!
FerriteAssembly._create_cell_state
FerriteAssembly.CellBuffer
FerriteAssembly.AutoDiffCellBuffer
FerriteAssembly.FaceBuffer
FerriteAssembly.FacetBuffer
FerriteAssembly.get_itembuffer
FerriteAssembly.skip_this_domain
FerriteAssembly.can_thread
FerriteAssembly.fast_getindex
FerriteAssembly.work_single_cell!
FerriteAssembly.work_single_face!
FerriteAssembly.work_single_facet!
FerriteAssembly.autogenerate_cellvalues
FerriteAssembly.autogenerate_facevalues
FerriteAssembly.autogenerate_facetvalues
```

## Threading model
Expand Down
4 changes: 2 additions & 2 deletions docs/src/literate_howto/local_constraints.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ add!(dh, :u, ip)
close!(dh)

ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfaceset(grid, "left"), Returns(zero(Vec{2}))))
add!(ch, Dirichlet(:u, getfaceset(grid, "right"), Returns(1.0), [1,]))
add!(ch, Dirichlet(:u, getfacetset(grid, "left"), Returns(zero(Vec{2}))))
add!(ch, Dirichlet(:u, getfacetset(grid, "right"), Returns(1.0), [1,]))
close!(ch)
update!(ch, 0.0)

Expand Down
24 changes: 12 additions & 12 deletions docs/src/literate_howto/robin_bc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,9 @@
# from the domain to the outside, where the latter has a constant temperature $u_\mathrm{b}$.

# ## Implementing Robin BC
# To implement this, we will overload the `face_routine!` function, which allows us to calculate
# both a residual and stiffness contribution from faces. Hence, we first define a "material",
# `RobinBC`, and implement the `face_routine!` that encodes the contribution to the weak form as
# To implement this, we will overload the `facet_routine!` function, which allows us to calculate
# both a residual and stiffness contribution from facets. Hence, we first define a "material",
# `RobinBC`, and implement the `facet_routine!` that encodes the contribution to the weak form as
# outlined above.
using Ferrite, FerriteAssembly

Expand All @@ -23,8 +23,8 @@ struct RobinBC{T}
ub::T
end

function FerriteAssembly.face_routine!(
Ke, re, ae, rbc::RobinBC, fv::FaceValues, facebuffer
function FerriteAssembly.facet_routine!(
Ke, re, ae, rbc::RobinBC, fv::FacetValues, facetbuffer
)
for q_point in 1:getnquadpoints(fv)
= getdetJdV(fv, q_point)
Expand Down Expand Up @@ -57,25 +57,25 @@ a = zeros(ndofs(dh));

# ## Using Robin BC
# For the Robin boundary condition, we also need to define the integration via
# FaceValues, in addition to an instance of `RobinBC`, for which we set the
# FacetValues, in addition to an instance of `RobinBC`, for which we set the
# "outside temperature" to -1.0.
face_qr = FaceQuadratureRule{RefQuadrilateral}(2);
fv = FaceValues(face_qr, ip);
facet_qr = FacetQuadratureRule{RefQuadrilateral}(2);
fv = FacetValues(facet_qr, ip);
rbc = RobinBC(1.0, -1.0);

# Finally, we set up a domain buffer with the actual materials and faceset,
# Finally, we set up a domain buffer with the actual materials and facetset,
# and create an assembler. Naturally, the boundary conditions are combined with
# other contributions, such as the assembly over cells. Typically, the same
# assembler can be used. If not, it is important to set the keyword arguments
# such that `K` and `r` are not zeroed between different calls to `work!`.
domainbuffer = setup_domainbuffer(DomainSpec(dh, rbc, fv; set=getfaceset(grid, "right")))
domainbuffer = setup_domainbuffer(DomainSpec(dh, rbc, fv; set=getfacetset(grid, "right")))
assembler = start_assemble(K, r)
work!(assembler, domainbuffer; a=a);

# ## Visualizing the BC
# To visualize the output, we can export the vtk the residual, showing how
# contributions have been added to the boundary.
vtk_grid("RobinBC", grid) do vtk
vtk_point_data(vtk, dh, r)
VTKFile("RobinBC", grid) do vtk
write_solution(vtk, dh, r)
end;

14 changes: 7 additions & 7 deletions docs/src/literate_howto/surface_integral.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,20 +3,20 @@ using Ferrite, FerriteAssembly
# This how-to shows integration of the normal flux
# on a surface. As usual, we need the basic Ferrite
# building blocks, in this case a grid, dofhandler,
# and facevalues
# and facetvalues
grid = generate_grid(Hexahedron, (5,5,5))
ip = Lagrange{RefHexahedron,1}()
dh = DofHandler(grid); add!(dh, :u, ip); close!(dh)
qr = FaceQuadratureRule{RefHexahedron}(2);
fv = FaceValues(qr, ip);
qr = FacetQuadratureRule{RefHexahedron}(2);
fv = FacetValues(qr, ip);

# We also need a solution vector to integrate,
# unless we only calculate geometric properties.
a = zeros(ndofs(dh))
apply_analytical!(a, dh, :u, norm); # f(x)=norm(x)

# And then we decide which faces to integrate over
domainbuffer = setup_domainbuffer(DomainSpec(dh, nothing, fv; set=getfaceset(grid, "right")));
# And then we decide which facets to integrate over
domainbuffer = setup_domainbuffer(DomainSpec(dh, nothing, fv; set=getfacetset(grid, "right")));

# ## Using `SimpleIntegrator`
# Using the simple interface, `SimpleIntegrator`, we simply
Expand All @@ -32,11 +32,11 @@ println("Flux: ∫qₙ dA = ", s_integrator.val)
# To demonstrate how the full-fledged interface can be used, we
# perform the same task by using the `Integrator`.
# To do that, we have to define an integrand, let's call it `NormalFlux`,
# and overload the face integration function
# and overload the facet integration function
mutable struct NormalFlux{T}
qn::T
end
function FerriteAssembly.integrate_face!(nf::NormalFlux, ae, material, fv, facebuffer)
function FerriteAssembly.integrate_facet!(nf::NormalFlux, ae, material, fv, facetbuffer)
for q_point in 1:getnquadpoints(fv)
dA = getdetJdV(fv, q_point)
∇u = function_gradient(fv, q_point, ae)
Expand Down
8 changes: 4 additions & 4 deletions docs/src/literate_tutorials/heat_equation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,8 @@ work!(assembler, buffer);
# ## Solve the problem.
# To actually solve the problem, we also need Dirichlet boundary conditions.
ch = ConstraintHandler(dh)
faceset = union((getfaceset(grid,k) for k in ("left", "right", "bottom", "top"))...)
add!(ch, Dirichlet(:u, faceset, Returns(0.0)))
facetset = union((getfacetset(grid,k) for k in ("left", "right", "bottom", "top"))...)
add!(ch, Dirichlet(:u, facetset, Returns(0.0)))
close!(ch);
apply_zero!(K, r, ch)
# where we use `apply_zero!` since we assembled assuming a zero temperature.
Expand All @@ -73,8 +73,8 @@ apply_zero!(K, r, ch)

# Finally, we can solve the problem and save the results
a = -K\r
vtk_grid("heat_equation", grid) do vtk
vtk_point_data(vtk, dh, a)
VTKFile("heat_equation", grid) do vtk
write_solution(vtk, dh, a)
end;

#md # ## [Plain program](@id heat_equation_plain_program)
Expand Down
26 changes: 13 additions & 13 deletions docs/src/literate_tutorials/iga.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,9 @@ function create_mesh(;nels=(20,10))
nurbsmesh = generate_nurbs_patch(:plate_with_hole, nels)
grid = BezierGrid(nurbsmesh)

addfaceset!(grid, "left", (x) -> x[1] -4.0)
addfaceset!(grid, "bot", (x) -> x[2] 0.0)
addfaceset!(grid, "right", (x) -> x[1] 0.0)
addfacetset!(grid, "left", (x) -> x[1] -4.0)
addfacetset!(grid, "bot", (x) -> x[2] 0.0)
addfacetset!(grid, "right", (x) -> x[1] 0.0)

return grid
end
Expand All @@ -34,10 +34,10 @@ grid = create_mesh();
orders=(2,2)
ip = BernsteinBasis{2,orders}()
qr_cell = QuadratureRule{2,RefCube}(4)
qr_face = QuadratureRule{1,RefCube}(3)
qr_facet = QuadratureRule{1,RefCube}(3)

cv = BezierCellValues( CellValues(qr_cell, ip^2) );
fv = BezierFaceValues( FaceValues(qr_face, ip^2) );
fv = BezierFacetValues( FacetValues(qr_facet, ip^2) );

# Distribute dofs as normal
dh = DofHandler(grid)
Expand All @@ -50,11 +50,11 @@ r = zeros(ndofs(dh))
K = create_sparsity_pattern(dh);

# Starting with Dirichlet conditions:
# 1) Bottom face should only be able to move in x-direction
# 1) Bottom facet should only be able to move in x-direction
# 2) Right boundary should only be able to move in y-direction
ch = ConstraintHandler(dh);
add!(ch, Dirichlet(:u, getfaceset(grid, "bot"), Returns(0.0), 2))
add!(ch, Dirichlet(:u, getfaceset(grid, "right"), Returns(0.0), 1))
add!(ch, Dirichlet(:u, getfacetset(grid, "bot"), Returns(0.0), 2))
add!(ch, Dirichlet(:u, getfacetset(grid, "right"), Returns(0.0), 1))
close!(ch)
update!(ch, 0.0);

Expand All @@ -63,7 +63,7 @@ update!(ch, 0.0);
# taking the negative value since r = fint - fext.
traction = Vec((-10.0, 0.0))
lh = LoadHandler(dh)
add!(lh, Neumann(:u, fv, getfaceset(grid, "left"), Returns(-traction)));
add!(lh, Neumann(:u, fv, getfacetset(grid, "left"), Returns(-traction)));

# FerriteAssembly setup
material = ElasticPlaneStrain(;E=100.0, ν=0.3)
Expand Down Expand Up @@ -109,10 +109,10 @@ projector = L2Projector(ip, grid)
σ_nodes = IGA.igaproject(projector, cellstresses.s, qr_cell; project_to_nodes=true);

# Output results to VTK
vtkgrid = vtk_grid("plate_with_hole.vtu", grid)
vtk_point_data(vtkgrid, dh, a)
vtk_point_data(vtkgrid, σ_nodes, "sigma", grid)
vtk_save(vtkgrid);
VTKFile("plate_with_hole.vtu", grid) do vtk
write_solution(vtk, dh, a)
write_node_data(vtk, σ_nodes, "sigma", grid)
end

using Test #src
@test sum(norm, σ_nodes) 3087.2447327126742 #src
Expand Down
8 changes: 4 additions & 4 deletions docs/src/literate_tutorials/incompressible_elasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -83,12 +83,12 @@ function solve(;ν, ipu, ipp)

## Boundary conditions
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfaceset(dh.grid, "left"), (x,t) -> zero(Vec{2})))
add!(ch, Dirichlet(:u, getfacetset(dh.grid, "left"), (x,t) -> zero(Vec{2})))
close!(ch)
update!(ch, 0.0)

lh = LoadHandler(dh)
add!(lh, Neumann(:u, 3, getfaceset(dh.grid, "right"), Returns(Vec{2}((0.0, 1/16)))))
add!(lh, Neumann(:u, 3, getfacetset(dh.grid, "right"), Returns(Vec{2}((0.0, 1/16)))))

## Cellvalues
qr = QuadratureRule{RefTriangle}(3)
Expand All @@ -114,8 +114,8 @@ function solve(;ν, ipu, ipp)
## Export the results
filename = "cook_" * (isa(ipu, Lagrange{2,RefTetrahedron,1}) ? "linear" : "quadratic") *
"_linear"
vtk_grid(filename, dh) do vtkfile
vtk_point_data(vtkfile, dh, u)
VTKFile(filename, dh) do vtk
write_solution(vtk, dh, u)
end
return u
end
Expand Down
Loading

0 comments on commit 00c5ae0

Please sign in to comment.