Boundary Conditions

Every PDE is accompanied with boundary conditions. There are different types of boundary conditions, and they need to be handled in different ways. Below we discuss how to handle the most common ones, Dirichlet and Neumann boundary conditions, and how to do it Ferrite.

Dirichlet Boundary Conditions

At a Dirichlet boundary the unknown field is prescribed to a given value. For the discrete FE-solution this means that there are some degrees of freedom that are fixed. To handle Dirichlet boundary conditions in Ferrite we use the ConstraintHandler. A constraint handler is created from a DoF handler:

ch = ConstraintHandler(dh)

We can now create Dirichlet constraints and add them to the constraint handler. To create a Dirichlet constraint we need to specify a field name, a part of the boundary, and a function for computing the prescribed value. Example:

dbc1 = Dirichlet(
    :u,                       # Name of the field
    getfaceset(grid, "left"), # Part of the boundary
    x -> 1.0,                 # Function mapping coordinate to a prescribed value
)

The field name is given as a symbol, just like when the field was added to the dof handler, the part of the boundary where this constraint is active is given as a face set, and the function computing the prescribed value should be of the form f(x) or f(x, t) (coordinate x and time t) and return the prescribed value(s).

Multiple sets

To apply a constraint on multiple face sets in the grid you can use union to join them, for example

left_right = union(getfaceset(grid, "left"), getfaceset(grid, "right"))

creates a new face set containing all faces in the "left" and "right" face sets, which can be passed to the Dirichlet constructor.

By default the constraint is added to all components of the given field. To add the constraint to selected components a fourth argument with the components should be passed to the constructor. Here is an example where a constraint is added to component 1 and 3 of a vector field :u:

dbc2 = Dirichlet(
    :u,                       # Name of the field
    getfaceset(grid, "left"), # Part of the boundary
    x -> [0.0, 0.0],          # Function mapping coordinate to prescribed values
    [1, 3],                   # Components
)

Note that the return value of the function must match with the components – in the example above we prescribe components 1 and 3 to 0 so we return a vector of length 2.

Adding the constraints to the constraint handler is done with add!:

add!(ch, dbc1)
add!(ch, dbc2)

Finally, just like for the dof handler, we need to use close! to finalize the constraint handler. Internally this will then compute the degrees-of-freedom that match the constraints we added.

If one or more of the constraints depend on time, i.e. they are specified as f(x, t), the prescribed values can be recomputed in each new time step by calling update! with the proper time, e.g.:

for t in 0.0:0.1:1.0
    update!(ch, t) # Compute prescribed values for this t
    # Solve for time t...
end
Examples

Most examples make use of Dirichlet boundary conditions, for example Heat Equation.

Neumann Boundary Conditions

At the Neumann part of the boundary we know something about the gradient of the solution.

As an example, the following code snippet can be included in the element routine, to evaluate the boundary integral:

for face in 1:nfaces(cell)
    if (cellid(cell), face) ∈ getfaceset(grid, "Neumann Boundary")
        reinit!(facevalues, cell, face)
        for q_point in 1:getnquadpoints(facevalues)
            dΓ = getdetJdV(facevalues, q_point)
            for i in 1:getnbasefunctions(facevalues)
                δu = shape_value(facevalues, q_point, i)
                fe[i] += δu * b * dΓ
            end
        end
    end
end

We start by looping over all the faces of the cell, next we check if this particular face is located on our faceset of interest called "Neumann Boundary". If we have determined that the current face is indeed on the boundary and in our faceset, then we reinitialize facevalues for this face, using reinit!. When reinit!ing facevalues we also need to give the face number in addition to the cell. Next we simply loop over the quadrature points of the face, and then loop over all the test functions and assemble the contribution to the force vector.

Examples

The following commented examples makes use of Neumann boundary conditions:

  • TODO

Periodic boundary conditions

Periodic boundary conditions ensure that the solution is periodic across two boundaries. To define the periodicity we first define the image boundary $\Gamma^+$ and the mirror boundary $\Gamma^-$. We also define a (unique) coordinate mapping between the image and the mirror: $\varphi:\ \Gamma^+\, \rightarrow\, \Gamma^-$. With the mapping we can, for every coordinate on the image, compute the corresponding coordinate on the mirror:

\[\boldsymbol{x}^- = \varphi(\boldsymbol{x}^+),\quad \boldsymbol{x}^- \in \Gamma^-,\, \boldsymbol{x}^+ \in \Gamma^+.\]

We now want to ensure that the solution on the image $\Gamma^+$ is mirrored on the mirror $\Gamma^-$. This periodicity constraint can thus be described by

\[u(\boldsymbol{x}^-) = u(\boldsymbol{x}^+).\]

Sometimes this is written as

\[\llbracket u \rrbracket = 0,\]

where $\llbracket \bullet \rrbracket := \bullet(\boldsymbol{x}^+) - \bullet(\boldsymbol{x}^-)$ is the "jump operator". Thus, this condition ensure that the jump, or difference, in the solution between the image and mirror boundary is the zero – the solution becomes periodic. For a vector valued problem the periodicity constraint can in general be written as

\[\boldsymbol{u}(\boldsymbol{x}^-) = \boldsymbol{R} \cdot \boldsymbol{u}(\boldsymbol{x}^+) \quad \Leftrightarrow \quad \llbracket \boldsymbol{u} \rrbracket = \boldsymbol{R} \cdot \boldsymbol{u}(\boldsymbol{x}^+) - \boldsymbol{u}(\boldsymbol{x}^-) = \boldsymbol{0}\]

where $\boldsymbol{R}$ is a rotation matrix. If the mapping between mirror and image is simply a translation (e.g. sides of a cube) this matrix will be the identity matrix.

In Ferrite this type of periodic Dirichlet boundary conditions can be added to the ConstraintHandler by constructing an instance of PeriodicDirichlet. This is usually done it two steps. First we compute the mapping between mirror and image faces using collect_periodic_faces. Here we specify the mirror set and image sets (the sets are usually known or can be constructed easily ) and the mapping $\varphi$. Second we construct the constraint using the PeriodicDirichlet constructor. Here we specify which components of the function that should be constrained, and the rotation matrix $\boldsymbol{R}$ (when needed). When adding the constraint to the ConstraintHandler the resulting dof-mapping is computed.

Here is a simple example where periodicity is enforced for components 1 and 2 of the field :u between the mirror boundary set "left" and the image boundary set "right". Note that no rotation matrix is needed here since the mirror and image are parallel, just shifted in the $x$-direction (as seen by the mapping φ):

# Create a constraint handler from the dof handler
ch = ConstraintHandler(dofhandler)

# Compute the face mapping
φ(x) = x - Vec{2}((1.0, 0.0))
face_mapping = collect_periodic_faces(grid, "left", "right", φ)

# Construct the periodic constraint for field :u
pdbc = PeriodicDirichlet(:u, face_mapping, [1, 2])

# Add the constraint to the constraint handler
add!(ch, pdbc)

# If no more constraints should be added we can close
close!(ch)
Note

PeriodicDirichlet constraints are imposed in a strong sense, so note that this requires a periodic mesh such that it is possible to compute the face mapping between faces on the mirror and boundary.

Examples

Periodic boundary conditions are used in the following examples Computational homogenization, Stokes flow.

Heterogeneous "periodic" constraint

It is also possible to define constraints of the form

\[\llbracket u \rrbracket = \llbracket f \rrbracket \quad \Leftrightarrow \quad u(\boldsymbol{x}^+) - u(\boldsymbol{x}^-) = f(\boldsymbol{x}^+) - f(\boldsymbol{x}^-),\]

where $f$ is a prescribed function. Although the constraint in this case is not technically periodic, PeriodicDirichlet can be used for this too. This is done by passing a function to PeriodicDirichlet, similar to Dirichlet, which, given the coordinate $\boldsymbol{x}$ and time t, computes the prescribed values of $f$ on the boundary.

Here is an example of how to implement this type of boundary condition, for a known function f:

pdbc = PeriodicDirichlet(
    :u,
    face_mapping,
    (x, t) -> f(x),
    [1, 2],
)
Note

One application for this type of boundary conditions is multiscale modeling and computational homogenization when solving the finite element problem for the subscale. In this case the unknown $u$ is split into a macroscopic part $u^{\mathrm{M}}$ and a microscopic/fluctuation part $u^\mu$, i.e. $u = u^{\mathrm{M}} + u^{\mu}$. Periodicity is then usually enforced for the fluctuation part, i.e. $\llbracket u^\mu \rrbracket = 0$. The equivalent constraint for $u$ then becomes $\llbracket u \rrbracket = \llbracket u^{\mathrm{M}} \rrbracket$.

As an example, consider first order homogenization where the macroscopic part is constructed as $u^{\mathrm{M}} = \bar{u} + \boldsymbol{\nabla} \bar{u} \cdot [\boldsymbol{x} - \bar{\boldsymbol{x}}]$ for known $\bar{u}$ and $\boldsymbol{\nabla} \bar{u}$. This could be implemented as

pdbc = PeriodicDirichlet(
    :u,
    face_mapping,
    (x, t) -> ū + ∇ū  ⋅ (x - x̄)
)

Initial Conditions

When solving time-dependent problems, initial conditions, different from zero, may be required. For finite element formulations of ODE-type, i.e. $\boldsymbol{u}'(t) = \boldsymbol{f}(\boldsymbol{u}(t),t)$, where $\boldsymbol{u}(t)$ are the degrees of freedom, initial conditions can be specified by the apply_analytical! function. For example, specify the initial pressure as a function of the y-coordinate

ρ = 1000; g = 9.81    # density [kg/m³] and gravity [N/kg]
grid = generate_grid(Quadrilateral, (10,10))
dh = DofHandler(grid); add!(dh, :u, 2); add!(dh, :p, 1); close!(dh)
u = zeros(ndofs(dh))
apply_analytical!(u, dh, :p, x -> ρ * g * x[2])

See also Time Dependent Problems for one example.

Note about solving DAE: A Differential Algebraic Equations (DAE) is an equation of the form $\boldsymbol{r}(\boldsymbol{u}(t),\boldsymbol{u}'(t),t)=\boldsymbol{0}$, which usually cannot be expressed as a true ODE. They occur often, but not always, in forms where some time derivatives are missing

\[u_1'(t) = f(\boldsymbol{u}(t),t) 0 = g(\boldsymbol{u}(t),t)`\]

In for such equations, it is usually necessary to specify initial conditions for both $\boldsymbol{u}(0)$ and $\boldsymbol{u}'(0)$, and these must be consistent, i.e. $\boldsymbol{r}(\boldsymbol{u}(0),\boldsymbol{u}'(0),0)=\boldsymbol{0}$.