API

box(a::Number, b::Number)
box(a::NTuple, b::NTuple)

Creates a CartesianProduct from two points a and b, which define the corners of the box. The component intervals are defined as $[\min(a_i, bᵢ), \max(a_i, b_i)]$. It accepts both numbers (for 1D) and NTuple (for D-dimensions).

interval(x, y)

Constructs a 1-dimensional CartesianProduct representing the closed interval $[x, y]$. The inputs are converted to floating-point numbers. It also accepts an existing 1D CartesianProduct as a single argument.

×(X::CartesianProduct, Y::CartesianProduct)

Computes the CartesianProduct of two CartesianProducts X and Y. The new dimension will be the sum of the dimensions of X and Y. Can be used as

X × Y

dim(X::CartesianProduct)

Returns the dimension of the space where the CartesianProduct is embedded. Can be applied directly to the type of the CartesianProduct.

Examples

julia> X = cartesian_product(0, 1);
	   dim(X)
1

julia> Y = cartesian_product(((0, 1), (4, 5)));
	   dim(Y)
2

dim(Ω::Domain)

Returns the dimension of the ambient space where the Domain Ω is embedded. It can also be applied to the type of the domain.

Example

julia> I = interval(0.0, 1.0);
	   dim(domain(I × I))
2

dim(Ωₕ::AbstractMeshType)

Returns the dimension of the space where Ωₕ is embedded.

This function is a required part of the AbstractMeshType interface. Can be applied to the type of an AbstractMeshType.

dim(Wₕ::AbstractSpaceType)

Returns the spatial dimension of the mesh associated with the functionpace Wₕ.

topo_dim(X::CartesianProduct)

Returns the topological dimension of a CartesianProduct. The depends on the dimension of the CartesianProduct and the number of collapsed dimensions.

topo_dim(Ω::Domain)

Returns the topological dimension Domain Ω.

topo_dim(Ωₕ::AbstractMeshType)

Returns the topological dimension Ωₕ.

markers(space_set, [time_set], pairs...)

Constructs a DomainMarkers object from a series of label => identifier pairs.

The identifier can be a Symbol, a Tuple of Symbols, or a Function. The full list of predefined boundary symbols can be found via get_boundary_symbols.

Example

julia> I = cartesian_product(0.0, 1.0);
	   tuples = (:corners => (:top, :right), :all_boundary => (:top, :right, :left, :bottom));
	   ids = (:left_boundary => :left, tuples..., :internal => x -> 0.2 < x < 0.8);
	   m = markers(I, ids...);

julia> length(m.symbols)
1

julia> length(m.tuples)
2

julia> length(m.conditions)
1

markers(Ω::Domain)

Returns the DomainMarkers object associated with the Domain Ω.

markers(Ωₕ::AbstractMeshType)

Returns the DomainMarkers) associated with the mesh Ωₕ. This function is a required part of the AbstractMeshType interface. Any concrete subtype of AbstractMeshType must implement this method and have a field called markers of type DomainMarkers.

domain(X::CartesianProduct, [markers...])

Returns a Domain from a CartesianProduct, assuming a single Marker with the label :boundary that marks the whole boundary of X. Alternatively, a list of Marker can be passed as argument in the form of :symbol => key (see examples and markers).

labels(Ω::Domain)

Returns a generator that yields the labels of all markers associated with the Domain Ω.

mesh(Ω::Domain, npts::Int, unif::Bool)
mesh(Ω::Domain, npts::NTuple{D}, unif::NTuple{D})

Returns a Mesh1D or a MeshnD ($D=2,3$) defined on the Domain Ω. The number of points for each coordinate projection mesh are given in the tuple npts. The distribution of points on the submeshes are encoded in the tuple unif.

For future reference, the mesh points are denoted as

• 1D mesh, with npts = $N_x$

\[x_i, \, i=1,\dots,N.\]

• 2D mesh, with npts = ($N_x$, $N_y$)

\[(x_i,y_j), \, i=1,\dots,N_x, \, j=1,\dots,N_y\]

• 3D mesh, with npts = ($N_x$, $N_y$, $N_z$)

\[(x_i,y_j,z_l), \, i=1,\dots,N_x, \, j=1,\dots,N_y, \, l=1,\dots,N_z.\]

Examples

I = interval(0, 1);
Ωₕ = mesh(domain(I), 10, true);

X = domain(interval(0, 1) × interval(4, 5));
Ωₕ = mesh(X, (10, 15), (true, false))

mesh(Wₕ::AbstractSpaceType)

Returns the underlying mesh object associated with the function space Wₕ.

points(Ωₕ::AbstractMeshType)

Returns the points of Ωₕ either as a vector (1D case) or a tuple of vectors (nD case).

• 1D mesh, with npts = $N_x$

\[x_i, \, i=1,\dots,N_x\]

• 2D mesh, with npts = ($N_x$, $N_y$)

\[([x_i]_{i=1}^{N_x}, [y_j]_{j=1}^{N_y})\]

• 3D mesh, with npts = ($N_x$, $N_y$, $N_z$)

\[([x_i]_{i=1}^{N_x}, [y_j]_{j=1}^{N_y}, [z_l]_{l=1}^{N_z}).\]

hₘₐₓ(Ωₕ::AbstractMeshType)

Returns the maximum diagonal of mesh Ωₕ.

• 1D mesh

\[h_{max} \vcentcolon = \max_{i=1,\dots,N} x_i - x_{i-1}.\]

• 2D mesh

\[\max_{i,j} \Vert (h_{x,i}, h_{y,j}) \Vert_2\]

• 3D mesh

\[\max_{i,j,l} \Vert (h_{x,i}, h_{y,j}, h_{z,l}) \Vert_2\]

npoints(Ωₕ::AbstractMeshType, [::Type{Tuple}])

Returns the number of points of mesh Ωₕ. If Tuple is passed as the second argument, it returns a tuple with the number of points of each submesh composing Ωₕ.

change_points!(Ωₕ::AbstractMeshType, [Ω::Domain], pts)

Changes the coordinates of the internal points of the mesh Ωₕ to the new coordinates specified in pts. This function assumes the points in pts are ordered and that the first and last of them coincide with the bounds of the Ω. if the domain Ω is passed as an argument, the markers of the mesh are also recalculated after this change.

iterative_refinement!(Ωₕ::AbstractMeshType, [domain_markers::DomainMarkers])

Refines the given mesh Ωₕ by halving each existing cell (in every direction). If an object of type DomainMarkers is passed as an argument, it also updates the markers according to accordingly after the refinement.

gridspace(Ωₕ::AbstractMeshType{D}; cache_avg = false, cache_bwd = true)

Constructor for a ScalarGridSpace defined on the mesh Ωₕ. This builds the weights for the inner products mentioned in ScalarGridSpace as well as the differentiation matrices associated with the grid points of mesh Ωₕ. The keyword arguments cache_avg and cache_bwd can be used to indicate if the average and backward difference matrices should be precomputed and stored in the space (default is true for cache_bwd and false for cache_avg).

element(Wₕ::AbstractSpaceType, [α::Number])

Returns a VectorElement for grid space Wₕ with uninitialized components. if $\alpha$ is provided, the components are initialized to $\alpha$.

element(Wₕ::AbstractSpaceType, v::AbstractVector)

Returns a VectorElement for a grid space Wₕ with the same coefficients of v.

mesh(Wₕ::AbstractSpaceType)

Returns the underlying mesh object associated with the function space Wₕ.

ndofs(Wₕ::AbstractSpaceType, [::Type{Tuple}])

Returns the total number of degrees of freedom (DOFs) in the function space Wₕ. If Tuple is passed, it returns a tuple with the number of DOFs in each dimension.

avgₕ(Wₕ::AbstractSpaceType, f)

Returns a VectorElement with the average of function f with respect to the cell_measure of mesh(Wₕ) around each grid point. It is defined as follows

• 1D case

\[\textrm{avg}ₕ f(x_i) = \frac{1}{|\square_i|} \int_{\square_i} f(x) dx, \, i = 1,\dots,N\]

• 2D case

\[\textrm{avg}ₕ f(x_i, y_j) = \frac{1}{|\square_{i,j}|} \iint_{\square_{i,j}} f(x,y) dA, \, i = 1,\dots,N_x, j = 1,\dots,N_y\]

• 3D case

\[\textrm{avg}ₕ f(x_i, y_j, z_l) = \frac{1}{|\square_{i,j,l}|} \iiint_{\square_{i,j,l}} f(x,y,z) dV, \, i = 1,\dots,N_x, j = 1,\dots,N_y, l = 1,\dots,N_z\]

Please check the implementations of functions cell_measure (for the 1-dimensional case) and cell_measure (for the n-dimensional cases).

avgₕ!(uₕ::VectorElement, f)

In-place version of averaging operator avgₕ.

Rₕ(Wₕ::AbstractSpaceType, f)

Standard nodal restriction operator. It returns a VectorElement with the result of evaluating the function f at the points of mesh(Wₕ). It is defined as follows

• 1D case

\[\textrm{R}ₕ f(x_i) = f(x_i), \, i = 1,\dots,N\]

• 2D case

\[\textrm{R}ₕ f(x_i, y_j)= f(x_i, y_j), \, i = 1,\dots,N_x, j = 1,\dots,N_y\]

• 3D case

\[\textrm{R}ₕ f(x_i, y_j, z_l)= f(x_i, y_j, z_l), \, i = 1,\dots,N_x, j = 1,\dots,N_y, l = 1,\dots,N_z\]

Rₕ!(uₕ::VectorElement, f; markers)

In-place version of the restriction operator Rₕ.

diff₋ₓ(arg)

Alias for backward_difference(arg, Val(1)). Computes the backward difference in the x-direction.

diff₋ᵧ(arg)

Alias for backward_difference(arg, Val(2)). Computes the backward difference in the y-direction.

diff₋₂(arg)

Alias for backward_difference(arg, Val(3)). Computes the backward difference in the z-direction.

diff₋ₕ(arg)

Computes the backward gradient of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, diff₋ₕ(uₕ) is equivalent to (backward_difference(uₕ, Val(1)), backward_difference(uₕ, Val(2))).

diff₊ₓ(arg)

Alias for forward_difference(arg, Val(1)). Computes the forward difference in the x-direction.

diff₊ᵧ(arg)

Alias for forward_difference(arg, Val(2)). Computes the forward difference in the y-direction.

diff₊₂(arg)

Alias for forward_difference(arg, Val(3)). Computes the forward difference in the z-direction.

diff₊ₕ(arg)

Computes the forward gradient of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, diff₊ₕ(uₕ) is equivalent to (forward_difference(uₕ, Val(1)), forward_difference(uₕ, Val(2))).

D₋ₓ(arg)

Alias for backward_finite_difference(arg, Val(1)). Computes the backward difference in the x-direction.

D₋ᵧ(arg)

Alias for backward_finite_difference(arg, Val(2)). Computes the backward difference in the y-direction.

D₋₂(arg)

Alias for backward_finite_difference(arg, Val(3)). Computes the backward difference in the z-direction.

∇₋ₕ(arg)

Computes the backward gradient of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, ∇₋ₕ(uₕ) is equivalent to (backward_finite_difference(uₕ, Val(1)), backward_finite_difference(uₕ, Val(2))).

D₊ₓ(arg)

Alias for forward_finite_difference(arg, Val(1)). Computes the forward difference in the x-direction.

D₊ᵧ(arg)

Alias for forward_finite_difference(arg, Val(2)). Computes the forward difference in the y-direction.

D₊₂(arg)

Alias for forward_finite_difference(arg, Val(3)). Computes the forward difference in the z-direction.

∇₊ₕ(arg)

Computes the forward gradient of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, ∇₊ₕ(uₕ) is equivalent to (forward_finite_difference(uₕ, Val(1)), forward_finite_difference(uₕ, Val(2))).

jump₋ₓ(arg)

Alias for backward_jump(arg, Val(1)). Computes the backward difference in the x-direction.

jump₋ᵧ(arg)

Alias for backward_jump(arg, Val(2)). Computes the backward difference in the y-direction.

jump₋₂(arg)

Alias for backward_jump(arg, Val(3)). Computes the backward difference in the z-direction.

jump₋ₕ(arg)

Computes the vectorial backward jump of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, jump₋ₕ(uₕ) is equivalent to (backward_jump(uₕ, Val(1)), backward_jump(uₕ, Val(2))).

jump₊ₓ(arg)

Alias for forward_jump(arg, Val(1)). Computes the forward difference in the x-direction.

jump₊ᵧ(arg)

Alias for forward_jump(arg, Val(2)). Computes the forward difference in the y-direction.

jump₊₂(arg)

Alias for forward_jump(arg, Val(3)). Computes the forward difference in the z-direction.

jump₊ₕ(arg)

Computes the vectorial forward jump of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, jump₊ₕ(uₕ) is equivalent to (forward_jump(uₕ, Val(1)), forward_jump(uₕ, Val(2))).

M₋ₓ(arg)

Alias for backward_average(arg, Val(1)). Computes the backward difference in the x-direction.

M₋ᵧ(arg)

Alias for backward_average(arg, Val(2)). Computes the backward difference in the y-direction.

M₋₂(arg)

Alias for backward_average(arg, Val(3)). Computes the backward difference in the z-direction.

M₋ₕ(arg)

Computes the vectorial backward average of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, M₋ₕ(uₕ) is equivalent to (backward_average(uₕ, Val(1)), backward_average(uₕ, Val(2))).

M₊ₓ(arg)

Alias for forward_average(arg, Val(1)). Computes the forward difference in the x-direction.

M₊ᵧ(arg)

Alias for forward_average(arg, Val(2)). Computes the forward difference in the y-direction.

M₊₂(arg)

Alias for forward_average(arg, Val(3)). Computes the forward difference in the z-direction.

M₊ₕ(arg)

Computes the vectorial forward average of arg, returning a tuple of operators/elements for each spatial dimension.

For a 2D space, M₊ₕ(uₕ) is equivalent to (forward_average(uₕ, Val(1)), forward_average(uₕ, Val(2))).

innerₕ(uₕ::VectorElement, vₕ::VectorElement)

Returns the discrete $L^2$ inner product of the grid functions uₕ and vₕ. Also accepts MatrixElement as any of the arguments.

• 1D case

\[(\textrm{u}_h, \textrm{v}_h)_h \vcentcolon = \sum_{i=1}^N |\square_{i}| \textrm{u}_h(x_i) \textrm{v}_h(x_i)\]

• 2D case

\[(\textrm{u}_h, \textrm{v}_h)_h \vcentcolon = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} |\square_{i,j}| \textrm{u}_h(x_i,y_j) \textrm{v}_h(x_i,y_j)\]

• 3D case

\[(\textrm{u}_h, \textrm{v}_h)_h \vcentcolon = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \sum_{l=1}^{N_z} |\square_{i,j,l}| \textrm{u}_h(x_i,y_j) \textrm{v}_h(x_i,y_j)\]

inner₊ₓ(uₕ::VecOrMatElem, vₕ::VecOrMatElem)

Returns the discrete modified $L^2$ inner product of the grid functions uₕ and vₕ associated with the first variable. It accepts arguments of type VectorElement or MatrixElement, in any order.

For VectorElements, it is defined as

• 1D case

\[(\textrm{u}_h, \textrm{v}_h)_+ \vcentcolon = \sum_{i=1}^{N_x} h_{i} \textrm{u}_h(x_i) \textrm{v}_h(x_i)\]

• 2D case

\[(\textrm{u}_h, \textrm{v}_h)_{+x} \vcentcolon = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y} h_{x,i} h_{y,j+1/2} \textrm{u}_h(x_i,y_j) \textrm{v}_h(x_i,y_j)\]

• 3D case

\[(\textrm{u}_h, \textrm{v}_h)_{+x} \vcentcolon = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}\sum_{l=1}^{N_z} h_{x,i} h_{y,j+1/2} h_{z,l+1/2} \textrm{u}_h(x_i,y_j,z_l) \textrm{v}_h(x_i,y_j,z_l).\]

inner₊ᵧ(uₕ::VecOrMatElem, vₕ::VecOrMatElem)

Returns the discrete modified $L^2$ inner product of the grid functions uₕ and vₕ associated with the second variable. It accepts

• 2D case

\[(\textrm{u}_h, \textrm{v}_h)_{+y} \vcentcolon = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y} h_{x,i} h_{y,j+1/2} \textrm{u}_h(x_i,y_j) \textrm{v}_h(x_i,y_j)\]

• 3D case

\[(\textrm{u}_h, \textrm{v}_h)_{+y} \vcentcolon = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}\sum_{l=1}^{N_z} h_{x,i+1/2} h_{y,j} h_{z,l+1/2} \textrm{u}_h(x_i,y_j,z_l) \textrm{v}_h(x_i,y_j,z_l).\]

inner₊₂(uₕ::VecOrMatElem, vₕ::VecOrMatElem)

Returns the discrete modified $L^2$ inner product of the grid functions uₕ and vₕ associated with the z variable

\[(\textrm{u}_h, \textrm{v}_h)_{+z} \vcentcolon = \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}\sum_{l=1}^{N_z} h_{x,i+1/2} h_{y,j+1/2} h_{z,l} \textrm{u}_h(x_i,y_j,z_l) \textrm{v}_h(x_i,y_j,z_l).\]

inner₊(uₕ::VecOrMatElem, vₕ::VecOrMatElem, [::Type{Tuple}])
inner₊(uₕ::NTuple{D}, vₕ::NTuple{D})

Returns the discrete modified $L^2$ inner product of the grid functions uₕ and vₕ. It accepts arguments of type VectorElement or MatrixElement, in any order.

If the Tuple argument is given, it returns D-tuple of all $\textrm{inner}_{x_i,+}$ applied to its input arguments, where D is the topological dimension of the mesh associated with the elements.

If NTuples of VectorElement or MatrixElement are passed as input arguments, it returns the sum of all inner products $(\textrm{u}_h[i],\textrm{v}_h[i])_{+x_i}$.

For VectorElements, the definition is given by

• 1D case

\[(\textrm{u}_h, \textrm{v}_h)_+ \vcentcolon = \sum_{i=1}^{N_x} h_{i} \textrm{u}_h(x_i) \textrm{v}_h(x_i)\]

• 2D case

\[(\textrm{u}_h, \textrm{v}_h)_+ \vcentcolon = (\textrm{u}_h, \textrm{v}_h)_{+x} + (\textrm{u}_h, \textrm{v}_h)_{+y}\]

• 3D case

\[(\textrm{u}_h, \textrm{v}_h)_+ \vcentcolon = (\textrm{u}_h, \textrm{v}_h)_{+x} + (\textrm{u}_h, \textrm{v}_h)_{+y} + (\textrm{u}_h, \textrm{v}_h)_{+z}.\]

See the definitions of inner₊ₓ, inner₊ᵧ and inner₊₂ for more details.

normₕ(uₕ::VectorElement)

Returns the discrete $L^2$ norm of the grid function uₕ, defined as

\[\Vert \textrm{u}_h \Vert_h \vcentcolon = \sqrt{(\textrm{u}_h, \textrm{u}_h)_h}\]

norm₊(uₕ::VectorElement)
norm₊(uₕ::NTuple{D,VectorElement})

Returns the discrete modified $L^2$ norm of the grid function uₕ. It also accepts a NTuple of VectorElements.

For VectorElements uₕ, it is defined as

\[\Vert \textrm{u}_h \Vert_+ = \sqrt{(\textrm{u}_h,\textrm{u}_h)_+}.\]

and for NTuples of VectorElements it returns

\[\Vert \textrm{u}_h \Vert_+ \vcentcolon = \sqrt{ \sum_{i=1}^D(\textrm{u}_h[i],\textrm{u}_h[i])_{+,x_i}}.\]

snorm₁ₕ(uₕ::VectorElement)

Returns the discrete version of the standard $H^1$ seminorm of VectorElement uₕ.

\[|\textrm{u}_h|_{1h} \vcentcolon = \Vert \nabla_h \textrm{u}_h \Vert_h\]

norm₁ₕ(uₕ::VectorElement)

Returns the discrete version of the standard $H^1$ norm of VectorElement uₕ.

\[\Vert \textrm{u}_h \Vert_{1h} \vcentcolon = \sqrt{ \Vert \textrm{u}_h \Vert_h^2 + \Vert \nabla_h \textrm{u}_h \Vert_h^2 }\]

dirichlet_constraints(cartesian_product, [I::CartesianProduct{1}], pairs...)

Creates Dirichlet boundary constraints.

Each pair is of the form :label => func, where :label identifies the boundary region and func defines the Dirichlet values. If the optional time domain I is provided, func should be a time-dependent function func(x, t).

The cartesian_product can be a CartesianProduct mesh domain or an ScalarGridSpace from which the mesh can be extracted. The :label must match a label in the mesh definition.

dirichlet_constraints(X::CartesianProduct, f::Function)

Creates a single Dirichlet boundary constraint with function `f` with the label `:dirichlet`.

form(Wₕ::AbstractSpaceType, Vₕ::AbstractSpaceType, f)

Returns a bilinear form from a given expression and trial and test spaces.

form(Wₕ::AbstractSpaceType, f)

Returns a linear form from a given expression fand a test spaceWₕ`.

assemble(a::BilinearForm, [dirichlet_labels])

Returns the assembled matrix of a bilinear form with imposed constraints.

assemble(l::LinearForm)

Returns the assembled linear form as a vector.

assemble(l::LinearForm, dirichlet_conditions::DomainMarkers, [dirichlet_labels])

Returns the assembled linear form with imposed constraints as a vector of numbers.

assemble!(A::AbstractMatrix, a::BilinearFormType [dirichlet_labels])

Copies the assembled matrix of a bilinear form and imposes the Dirichlet constraints to a given matrix A.

assemble!(vec::AbstractVector, l::LinearForm; dirichlet_conditions::DomainMarkers, [dirichlet_labels])

In-place assemble of a linear form with imposed constraints into a given vector.

symmetrize!(A, F, Ωₕ, labels)

Modifies the linear system Ax = F to make A symmetric after applying Dirichlet conditions. It updates the vector F and zeros out the columns of A corresponding to Dirichlet nodes.

The algorithm goes as follows: for any given row i where Dirichlet boundary conditions have been applied

- calculate `dᵢ = cᵢ .* F`, where `cᵢ` is the `i`-th column of `A`;
- replace `F` by substracting `dᵢ` to `F` (except for the `i`-th component)
- replace all elements in the `i`-th column of `A` (except the `i`-th by zero).