Input

Scenario

Scenarios to simulate: instances(Scenario)

• F_CAVI = 1 –> lid-driven cavity
• F_COUE = 2 –> Couette flow
• F_POIS = 3 –> Hagen–Poiseuille flow
• F_PULL = 4 –> tether pulling from a flat patch
• F_BEND = 5 –> bending of a flat patch

The first letter indicates mesh Topology: F for a flat patch and C for a cylinder. At present, no cylindrical scenarios are implemented.

To create a new scenario:

• add and export a new entry in the Scenario enum
• edit get_scenario_bc_info in src/input/Bc.jl
• write the corresponding apply_<SCENARIO>_bcs function in src/input/Mesh.jl
• add the new scenario to the implemented list in src/input/Params.jl

Topology

Mesh topologies: instances(Topology)

• FLAT = 1
• CYLINDER = 2

Motion

Types of mesh motion: instances(Motion)

• STATIC = 1 –> mesh does not move
• EUL = 2 –> Eulerian mesh motion
• LAG = 3 –> Lagrangian mesh motion
• ALEV = 4 –> ALE-viscous mesh motion
• ALEVB = 5 –> ALE-viscous-bending mesh motion

The choice of motion is stored in Params.motion variable in Params.jl, and is relevant both to the specification of boundary conditions in Bc.jl and the finite element analysis.

Boundary

Enumerate four mesh boundaries: instances(Boundary), never to be used as an index.

• BOTTOM = 1
• RIGHT = 2
• TOP = 3
• LEFT = 4

Corner

Enumerate four mesh corners: instances(Corner), never to be used as an index.

• BOTTOM_LEFT = 1
• BOTTOM_RIGHT = 2
• TOP_LEFT = 3
• TOP_RIGHT = 4

Neumann

Enumerate types of Neumann boundary conditions: instances(Neumann)

• SHEAR = 1 –> apply in-plane force/length in the τ direction
• STRETCH = 2 –> apply in-plane force/length in the ν direction
• MOMENT = 3 –> apply boundary moment

Associate finite element calculations are found in calc_bdry_element_residual

Curve

Enumerate types of curves: instances(Curve)

• CLAMPED
• CLOSED

These curves are used to generate B-spline basis functions—see Spline.jl

Params

Parameters that must be specified for any simulation.

For Scenario-specific data that must be included via keyword arguments, see check_params. The following information is contained:

• motion::Motion → type of motion
• scenario::Scenario → type of scenario
• num1el::Int64 → number of elements in $ζ^1$
• num2el::Int64 → number of elements in $ζ^2$
• output::Bool → write output (yes by default)
• length::Float64 → length, in real space
• kb::Float64 → mean bending modulus
• kg::Float64 → Gaussian bending modulus
• ζv::Float64 → membrane viscosity
• pn::Float64 → normal pressure (body force)
• poly::Int64 → polynomial order of basis functions
• gp1d::Int64 → number of Gauss points, in 1-D
• nders::Int64 → number of 1-D B-spline derivatives
• nen::Int64 → number of element nodes
• αdb::Float64 → Dohrmann–Bochev parameter
• αm::Float64 → mesh parameter
• εk::Float64 → tolerance numerical calculation of matrix K
• εnr::Float64 → Netwon–Raphson tolerance

check_params(p::Params; args...)

Check relevant parameters, including scenario-specific keyword arguments.

For any simulation, the fields t0, t0_id, and Δts are required. If data is output (as is the default, for p.output), the output path out_path and file name out_file are necessary. Finally, the information needed to solve the various Scenarios is checked.

Dof.Unknown

Possible degrees of freedom (or fundamental unknowns) at each node.

For different mesh motions, different fundamental unknowns are required. A Lagrangian simulation requires only vx, vy, vz, and λ, while an ALE simulation requires all eight fundamental unknowns. Currently, the mapping from the choice of Motion to the corresponding fundamental unknowns is specified in the function get_dofs and stored in Mesh.dofs.

This enum is encapsulated in a module; access elements with e.g. Dof.vx.

Types of fundamental unknowns:

• vx –> $x$-velocity, $v_x$
• vy –> $y$-velocity, $v_y$
• vz –> $z$-velocity, $v_z$
• vmx –> $x$-mesh velocity, $v^{\text{m}}_x$
• vmy –> $y$-mesh velocity, $v^{\text{m}}_y$
• vmz –> $z$-mesh velocity, $v^{\text{m}}_z$
• λ –> surface tension, $λ$
• pm –> mesh pressure, $p^{\text{m}}$

Dof.Position

Cartesian components of the mesh position at each node.

This enum is encapsulated in a module; access elements with e.g. Dof.xm.

Options:

• xm –> $x$-position
• ym –> $y$-position
• zm –> $z$-position

KnotVector(ζs::Vector{Float64}, nel::Int64, poly::Int64, curve::Curve)

Vector of knots, with which all 1-D B-spline functions are calculated.

The knot vector struct contains four fields:

• ζs –> list of all knots, including repeated ones
• nel –> number of 1-D elements, or knot spans
• poly –> polynomial order
• curve –> type of Curve, which is either CLAMPED or CLOSED

All of these data are not required to generate a knot vector. Two constructors are available—see knot_vector

knot_vector(ζs::Vector{Float64}, poly::Int64, curve::Curve)

External constructor of the KnotVector struct, where knots need not be uniform.

The list of knots ζs, polynomial order poly, and Curve type curve are passed in. If they represent a valid knot vector, then the struct is created.

knot_vector(nel::Int64, poly::Int64, curve::Curve)

External constructor of the KnotVector struct with uniform knots.

The number of elements nel, polynomial order poly, and Curve type curve are passed in. If the curve is CLOSED, then no knots are repeated. If the curve is CLAMPED, then the first and last (poly+1) knots are repeated. In both cases, the generated list of knots ζs ranges from 0 to 1.

get_knot_span_index(kv::KnotVector, ζ::Float64)::Int64

Return knot span index of ζ in CLAMPED or CLOSED KnotVector kv.

Algorithm A2.1 in Chap. 2, §5 of Piegl and Tiller (1997) ensures kv.ζs[idx] ≤ ζ < kv.ζs[idx+1], where idx is the returned value. The algorithm, valid only for CLAMPED knots, is here extended to CLOSED knots as well.

get_bspline_vals(kv::KnotVector, ζ::Float64)::Vector{Float64}

Return kv.poly+1 nonzero B-spline functions at ζ for the given KnotVector kv.

Algorithm A2.2 in Chap. 2, §5 of Piegl and Tiller (1997), valid for CLAMPED knot vectors with (kv.poly+1) repeated knots at the start and end, is also used for CLOSED knot vectors. The local-to-global mapping returned by get_bspline_indices is necessary to determine which global basis functions are nonzero at the provided ζ.

get_bspline_ders(kv::KnotVector, ζ::Float64, num_ders::Int64)::Matrix{Float64}

Return B-spline functions and num_ders derivatives at ζ for the given KnotVector kv.

Algorithm A2.3 in Chap. 2, §5 of Piegl and Tiller (1997), valid for CLAMPED knot vectors with (kv.poly+1) repeated knots at the start and end, is also used for CLOSED knot vectors. Here num_ders is required to not be greater than kv.poly, and also greater than or equal to zero. The local-to-global mapping returned by get_bspline_indices is necessary to determine which global basis functions are nonzero at the provided ζ.

The returned matrix has kv.poly+1 rows and num_ders+1 columns. The first column contains the functions themselves, and is thus identical to the output of get_bspline_vals. The second column contains the first derivatives, and so on.

get_bspline_indices(kv::KnotVector, ζ::Float64)::Vector{Int64}

Return global indices of nonzero B-splines at the value ζ in the KnotVector kv.

By definition, there are only kv.poly+1 nonzero B-spline functions at any ζ. The functions get_bspline_vals and get_bspline_ders calculate these nonzero quantities, and the mapping returned here places them within the global ordering. The global indices are the same for all ζ within a single knot span, and thus we begin by calling get_knot_span_index.

Both CLAMPED knot vectors, with (kv.poly+1) repeats at the start and end, as well as CLOSED knot vectors, are handled.

get_bspline_indices(kv::KnotVector, ks_id::Int64)::Vector{Int64}

Return global indices of nonzero B-splines at ks_id$^{\textrm{th}}$ knot span in the KnotVector kv.

By definition, there are only kv.poly+1 nonzero B-spline functions over any knot span. The mapping returned here places these functions within the global ordering. Both CLAMPED knot vectors, with (kv.poly+1) repeats at the start and end, as well as CLOSED knot vectors, are handled. Note that ks_id is not the element ID eid; rather, ks_id = eid + kv.poly

get_1d_bspline_cps(kv::KnotVector, x::Function)::Vector{Float64}

Return the global control points $\{ x^{}_K \}$ for the function $x(ζ)$.

The global control points $\{ x^{}_K \}$ are calculated for the provided KnotVector kv such that

\[x(ζ) \, \approx \, \sum_{K = 1}^{\texttt{nn}} N^{}_K (ζ) \, x^{}_K ~,\]

where $x(ζ)$ is the provided x::Function and here $\texttt{nn}$ is the number of global 1-D nodes. To determine the unknown $\{ x^{}_K \}$, a set of $\texttt{nn}$ collocation points $\{ ζ^{}_J \}$ is chosen via collocate_ζ. A matrix equation is then obtained as

\[[x (ζ^{}_J)] \, = \, [N^{}_K (ζ^{}_J)] \, [x^{}_K] ~,\]

where the $N^{}_K (ζ^{}_J)$ are determined with get_bspline_indices and get_bspline_vals. The $x^{}_K$ are then calculated and returned as a vector.

get_2d_bspline_cps(kv1::KnotVector, kv2::KnotVector, x::Function)::Matrix{Float64}

Return the global control points $\{ x^{}_K \}$ for the scalar-valued function $z(ζ^1, ζ^2)$.

The procedure of get_1d_bspline_cps is repeated, except in this case collocation points are chosen across the 2-D parametric domain $(ζ^1, ζ^2)$.

collocate_ζ(kv::KnotVector)::Vector{Float64}

Take a KnotVector kv and return a list of points ζ used for collocation.

The length of the returned list of ζ values is equal to the number of global basis functions associated with the provided knot vector. Both CLAMPED and CLOSED knot vectors are accounted for.

get_unique_1d_elements(kv::KnotVector)

Get all unique elements for a given knot vector kv.

This function returns several quantities, in the following order:

• uel_num::Int64 -> number of unique elements
• num_el::Int64 -> number of elements
• uel_ids::Vector{Int64} -> mapping from element id to unique element id
• uel_list::Vector{Tuple{Float64, Float64}} -> start and end ζ for each unique element

GaussPointsξ(ngp::Int64)

Weights $\hat{w}_k$ and points $ξ_k$ on the interval $[-1, +1]$ with which a 1-D integral is approximated.

The struct has three fields:

• ngp –> number of Gaussian quadrature points
• ξs –> set of points $ξ_k$
• ws –> set of weights $\hat{w}_k$

Calculations from Wikipedia: Gaussian quadrature.

GaussPointsζ(gpsξ::GaussPointsξ, ζlo::Float64, ζhi::Float64)

Weights $w_k$ and points $ζ_k$ with which a 1-D integral is approximated over ζlo ≤ ζ ≤ ζhi.

Calculated via the transform from ξ to ζ via the relations

ζ$_k$ = ξ$_k$ (ζhi - ζlo) / 2 + (ζhi + ζlo) / 2

$w_k$ = $\hat{w}_k$ (ζhi - ζlo) / 2

The struct has three fields:

• ngp –> number of Gaussian quadrature points
• ζs –> set of points $ζ_k$
• ws –> set of weights $w_k$

GpBasisFnsζ(gpwζ::Float64, ζ::Float64, kv::KnotVector)

1-D basis functions $N(ζ)$ at the provided quadrature point ζ with weight gpwζ.

The B-spline values and derivatives are determined with get_bspline_ders and stored here. The struct contains four fields:

• w –> Gauss point weight, passed in as gpw
• N –> (poly+1)×1 vector of nonzero basis functions $N(ζ)$
• dN –> (poly+1)×1 vector of nonzero first derivatives $N'(ζ)$
• ddN –> (poly+1)×1 vector of nonzero second derivatives $N''(ζ)$

GpBasisFnsζα(fns1::GpBasisFnsζ, fns2::GpBasisFnsζ)

2-D basis functions $N(ζ^α)$ at a single Gauss point $ζ^α$, as set by 1-D GpBasisFnsζ fns1 and fns2.

Here $ζ^1$ and $ζ^2$ respectively correspond to the ζ values passed to fns1 and fns2, even though this information is not stored. The 2-D basis functions are calculated as the tensor product of 1-D functions. The struct contains four fields:

• w –> Gauss point weight, given by fns1.w × fns2.w
• N –> NEN×1 vector of nonzero basis functions $N(ζ^α)$, ordered by the tensor product structure
• ∂Nα –> NEN×2 matrix of first derivatives: rows respectively contain $N_{, 1}$ and $N_{, 2}$
• ∂∂Nαβ –> NEN×3 matrix of second derivatives: rows respectively contain $N_{, 1 1}$, $N_{, 2 2}$, and $N_{, 1 2}$

LineGpBasisFns(kv::KnotVector, ngp::Int64)

1-D B-spline basis functions and derivatives at all quadrature points on the domain of the KnotVector kv.

Basis functions are determined at each of the ngp Gauss points, over each element. When using B-splines, scenarios often arise where the basis functions are identical across many elements. We accordingly store only the 1-D basis functions over unique elements, as well as a mapping from elements to unique elements. The basis functions at the edges of the 1-D parametric domain are stored separately, as they are required subsequently for application of boundary conditions. The struct has five fields:

• nel –> number of elements in the KnotVector kv
• uel_ids –> mapping from elem_id to unique_elem_id
• ufns –> num_unique_elems×ngp matrix of unique 1-D basis functions, of type GpBasisFnsζ
• ζmin_fns –> 1-D basis functions at the min value of the parametric domain
• ζmax_fns –> 1-D basis functions at the max value of the parametric domain

Note that this struct contains no information about how basis functions vary in the orthogonal parametric direction, and thus cannot be directly used to apply boundary conditions. See BdryGpBasisFns.

BdryGpBasisFns(line_gp_fns::LineGpBasisFns, perp_edge_fns::GpBasisFnsζ, bdry::Boundary)

2-D B-spline basis functions and derivatives at all quadrature points on a boundary.

The LineGpBasisFns struct line_gp_fns contains all 1-D basis functions and derivatives along the boundary, with ngp Gauss points for each element. Here perp_edge_fns are the basis functions and derivatives at the boundary in the orthogonal parametric direction, as contained in a GpBasisFnsζ struct. For example, on the RIGHT boundary, line_gp_fns captures 1-D derivatives in the $ζ^2$ direction, while perp_edge_fns contains 1-D derivatives in the $ζ^1$ direction on the right edge. With knowledge of the underlying tensor product structure (Piegl and Tiller, 1997), the 2-D basis functions along the specified Boundary bdry are generated. The struct has four fields:

• nel –> number of elements associated with LineGpBasisFns line_gp_fns
• uel_ids –> mapping from elem_id to unique_elem_id
• ufns –> num_unique_elems×ngp matrix of unique 2-D basis functions, of type GpBasisFnsζα
• bdry –> Boundary of the parametric domain

AreaGpBasisFns(line_gp_fns1::LineGpBasisFns, line_gp_fns2::LineGpBasisFns)

2-D B-spline basis functions and derivatives at all quadrature points in the mesh area.

The two LineGpBasisFns structs passed in, line_gp_fns1 and line_gp_fns2, respectively contain 1-D derivatives in the $ζ^1$ and $ζ^2$ directions. With the tensor product structure of B-splines over a surface (Piegl and Tiller, 1997), it is straightforward to calculate the 2-D basis functions and their derivatives. As is the case for LineGpBasisFns, a mapping from elements to unique elements is generated; only unique basis function calculations are stored. This struct has three fields:

• nel –> number of area elements
• uel_ids –> mapping from elem_id to unique_elem_id
• ufns –> num_unique_elems×ngp matrix of unique 2-D basis functions, of type GpBasisFnsζα

Mesh(p::Params; args...)

The generated mesh, which includes the degree-of-freedom numbering and all basis functions for the given Scenario.

This struct contains the following scenario-independent fields:

• num1el –> number of elements in $ζ^1$ direction
• num2el –> number of elements in $ζ^2$ direction
• numel –> number of area elements
• num1np –> number of nodal points in $ζ^1$ direction
• num2np –> number of nodal points in $ζ^2$ direction
• numnp –> number of nodal points on the 2-d mesh
• IX –> matrix containing nonzero node numbers for each area element
• bdry_elems –> mapping from Boundary to element numbers
• crnr_elems –> mapping from Corner to element number
• bdry_nodes –> mapping from Boundary to node numbers
• bdry_inner_nodes –> mapping from Boundary to inner node numbers
• crnr_nodes –> mapping from Corner to node number
• crnr_inner_nodes –> mapping from Corner to inner node number
• kv1 –> KnotVector along $ζ^1$ direction
• kv2 –> KnotVector along $ζ^2$ direction
• area_gp_fns –> AreaGpBasisFns: 2-D area basis functions
• bdry_gp_fns –> mapping from Boundary to BdryGpBasisFns
• crnr_gp_fns –> mapping from Corner to 2-D basis functions GpBasisFnsζα

The generated mesh depends on the problem being solved. The function generate_scenario uses the scenario-specific results from Bc.jl to set the following scenario-dependent fields:

• topology –> surface Topology
• dofs –> list of active Unknowns with global ordering
• ndf –> number of active Unknowns
• ID –> matrix mapping node number to indices of unknowns at that node
• ID_inv –> inverse of ID: map unknown_id to (node_number, dof_number)
• nmdf –> total number of mesh degrees of freedom
• LM –> map from area element to indices of all unknowns of that element
• inh_dir_bcs –> list of inhomogeneous Dirichlet boundary conditions
• inh_neu_bcs –> list of inhomogeneous Neumann boundary conditions

generate_scenario(numel, numnp, IX, bdry_nodes, ..., p::Params; args...)

Apply scenario-specific boundary conditions towards mesh generation.

Checks to see whether a helper function is available in Bc.jl for the given Scenario, and then constructs the ID matrix that maps from nodes to global degrees of freedom. The inverse ID_inv is also generated, as is the LM matrix—which provides the local degrees of freedom for each element.

get_basis_fns(el_id::Int64, gp_id::Int64, mesh::Mesh)::GpBasisFnsζα

Return GpBasisFnsζα of the gp_id$^{\text{th}}$ 2-D Gauss point of the el_id$^{\text{th}}$ area element.

get_basis_fns(bdry::Boundary, el_id::Int64, gp_id::Int64, mesh::Mesh)::GpBasisFnsζα

Return GpBasisFnsζα of the gp_id$^{\text{th}}$ 1-D Gauss point of the el_id$^{\text{th}}$ element on the Boundary bdry.

get_gpw(el_id::Int64, gp_id::Int64, mesh::Mesh)::Float64

Return the Gauss point weight of the call to get_basis_fns.

get_gpw(bdry::Boundary, el_id::Int64, gp_id::Int64, mesh::Mesh)::Vector{Float64}

Return the Gauss point weight of the call to get_basis_fns.

get_N(el_id::Int64, gp_id::Int64, mesh::Mesh)::SVector{NEN,Float64}

Return local basis functions $[𝐍^e]$ of the call to get_basis_fns.

get_N(bdry::Boundary, el_id::Int64, gp_id::Int64, mesh::Mesh)

Return local basis functions $[𝐍^e]$ of the call to get_basis_fns.

get_∂Nα(el_id::Int64, gp_id::Int64, mesh::Mesh)::SMatrix{NEN,ζDIM,Float64}

Return local basis function derivatives $[𝐍^e]_{, α}$ of the call to get_basis_fns.

get_∂Nα(bdry::Boundary, el_id::Int64, gp_id::Int64, mesh::Mesh)

Return local basis function derivatives $[𝐍^e]_{, α}$ of the call to get_basis_fns.

get_∂∂Nαβ(el_id::Int64, gp_id::Int64, mesh::Mesh)::SMatrix{NEN,VOIGT,Float64}

Return local basis function derivatives $[𝐍^e]_{, α β}$ of the call to get_basis_fns.

get_∂∂Nαβ(bdry::Boundary, el_id::Int64, gp_id::Int64, mesh::Mesh)

Return local basis function derivatives $[𝐍^e]_{, α β}$ of the call to get_basis_fns.

get_scenario_bc_info(numnp, IX, bdry_nodes, ..., p::Params; args...)

Return Scenario-specific information for the provided Params.

Each of the functions called returns (dofs, ndf, ID, inh_dir_bcs, inh_neu_bcs), in that order.

get_f_cavi_bc_info(numnp, bdry_nodes, crnr_nodes)

Generate Mesh data for the flat cavity Scenario F_CAVI.

In this scenario, a STATIC mesh Motion is required; relevant Unknowns are vx, vy, and λ. The following boundary conditions are prescribed:

• $v_x = 1.0$ along the top edge (corners excluded)
• $v_x = 0.0$ along all other edges (corners included)
• $v_y = 0.0$ along all edges
• $λ = 0.0$ at (or near) the center of the domain

Note that $λ$ is pinned to remove the global indeterminacy of the surface tension, up to a constant.

get_f_coue_bc_info(numnp, bdry_nodes)

Generate Mesh data for the flat Couette Scenario F_COUE.

In this scenario, a STATIC mesh Motion is required; relevant Unknowns are vx, vy, and λ. The following boundary conditions are prescribed:

• $v_x = 0.0$ along the bottom edge (corners included)
• $v_x = 3.0$ along the top edge (corners included)
• $v_y = 0.0$ along all edges
• $\bm{\bar{F}} = 4\bm{\nu}$ on the left and right edges

get_f_pois_bc_info(numnp, bdry_nodes)

Generate Mesh data for the flat Hagen–Poiseuille Scenario F_POIS.

In this scenario, a STATIC mesh Motion is required; relevant Unknowns are vx, vy, and λ. The following boundary conditions are prescribed:

• $v_x = 0.0$ along the top and bottom edges (corners included)
• $v_y = 0.0$ on all edges
• $\bm{\bar{F}} = 4 \bm{\nu}$ on the left edge
• $\bm{\bar{F}} = 8 \bm{\nu}$ on the right edge

get_f_pull_bc_info(numnp, IX, bdry_nodes, bdry_inner_nodes, p::Params; args...)

Generate Mesh data for the flat tube-pulling Scenario F_PULL.

Relevant unknowns are determined via get_dofs, depending on the choice of mesh Motion. The following boundary conditions are prescribed for all mesh motions:

• $v_z = 0$ on all boundaries, and all inner boundaries
• $\bm{\bar{F}} = (k_{\text{b}} / 4)\bm{\nu}$ on all boundaries
• $v_x = 0$, $v_y = 0$ at center of each edge
• $v_x = 0$, $v_y = 0$, $v_z = \,$args[:pull_speed]$\bm{e}_z$ on nodes of the central element

The following are prescribed when the Motion is not Lagrangian:

• $v^{\text{m}}_x = 0$, $v^{\text{m}}_y = 0$ at center of each edge
• $v^{\text{m}}_x = 0$, $v^{\text{m}}_y = 0$, $v^{\text{m}}_z = \,$args[:pull_speed}$\bm{e}_z$ on nodes of the central element

The following are prescribed when the Motion is either ALE-viscous or ALE-viscous-bending:

• $\bm{v}^{\text{m}} = \bm{0}$ on all boundaries

The following are prescribed when the Motion is ALE-viscous-bending:

• $v^{\text{m}}_z = 0$ on all inner boundaries

get_f_bend_bc_info(numnp, bdry_nodes, p::Params; args...)

Generate Mesh data for the flat bending Scenario F_BEND.

Relevant unknowns are determined via get_dofs, depending on the choice of mesh Motion. The following boundary conditions are prescribed for all mesh motions:

• LEFT edge:
  • $\bm{v} = \bm{0}$
  • $M = \,$args[:bend_mf] (see Params.jl)
  • if ALE mesh motion: $\bm{v}^{\text{m}} = \bm{0}$
• RIGHT edge:
  • $f_x = 0$ and $f_y = 0$
  • $v_z = 0$
  • $M = \,$args[:bend_mf] (see Params.jl)
  • if ALE mesh motion: $v^{\text{m}}_z = 0$
• BOTTOM edge:
  • $f_x = 0$ and $f_z = 0$
  • $v_y = 0$
  • $M = 0$
  • if ALE mesh motion: $v^{\text{m}}_y = 0$
• TOP edge:
  • $f_x = 0$ and $f_z = 0$
  • $v_y = 0$
  • $M = 0$
  • if ALE mesh motion: $v^{\text{m}}_y = 0$

get_dofs(motion::Motion)::Dict{Dof.Unknown, Int64}

Return global ordering of Unknown degrees of freedom, depending on the mesh Motion.

The following is the mapping from motion to unknowns:

• Lagrangian –> $v_x$, $v_y$, $v_z$, $λ$
• Eulerian –> $v_x$, $v_y$, $v_z$, $v^{\text{m}}_x$, $v^{\text{m}}_y$, $v^{\text{m}}_z$, $λ$
• ALE –> $v_x$, $v_y$, $v_z$, $v^{\text{m}}_x$, $v^{\text{m}}_y$, $v^{\text{m}}_z$, $λ$, $p^{\text{m}}$

Note that the flat, 2-D Scenarios (F_CAVI, F_COUE, F_POIS) have a reduced number of degrees of freedom, because there are no unknowns in the $z$-direction. Thus, this function is not called in those scenarios.

prepare_input(p::Params; args...)

Generate and return the mesh, surface position xms, and Unknown control points cps.

The mesh is generated by instantiating a Mesh with the provided parameters and arguments, which include the Scenario and number of elements. The surface positions and initial unknowns are set with calls to get_1d_bspline_cps and get_2d_bspline_cps.