Initial and 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).

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

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

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.

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)

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],
)

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.