Taiga.jl

Tensor-product applications in isogeometric analysis.

Note

Examples listed on this page use following packages

using Taiga
using LinearAlgebra
using SortedSequences, CartesianProducts, KroneckerProducts, NURBS, SpecialSpaces, AbstractMappings, UnivariateSplines

Postprocessing

Some basic postprocessing methods are provided as Plots recipes in most packages in the Feather ecosystem. Taiga.jl implements Bézier extraction and can save mappings, fields and control nets to a VTK file, e.g

mapping = hole_in_square_plate()
mapping = refine(mapping, method=hRefinement(3))
space = ScalarSplineSpace(mapping.space)

euclidean_distance = Field((x,y) -> sqrt(x^2 + y^2))
d = Field(space)
project!(d, onto=euclidean_distance ∘ mapping; method=QuasiInterpolation)

vtk_save_bezier("geometry_and_field", mapping; fields = Dict("d" => d))
vtk_save_control_net("control_net", mapping)

Bezier extraction on square plate with a hole

More detailed examples are provided in the examples directory.

Aggregation of Kronecker products

More often then not Kronecker product approximations will be composed of a sum of Kronecker products. In most cases these sums will not have the classical Kronecker sum structure. Nonetheless, it is still possible to evaluate matrix-vector products with these operators with complexity $\mathcal O(M\cdot N^{3/2})$ and $\mathcal O(M\cdot N^{4/3})$ in 2 and 3 dimensions, where $M$ is the number of Kronecker products in the collection.

KroneckerProductAggregate is a collection of Kronecker product matrices of conforming size that acts as a linear map and supports reduction over the collection, e.g.

julia> A = rand(7, 3) ⊗ rand(5, 3) ⊗ rand(3, 5);

julia> B = rand(7, 3) ⊗ rand(5, 3) ⊗ rand(3, 5);

julia> C = A + B;

julia> K = KroneckerProductAggregate(A, B);

julia> v = rand(size(C, 2));

julia> b = K * v;

julia> b ≈ C * v
true

julia> mul!(b, K, v); # in-place evaluation

julia> b ≈ C * v
true

A KroneckerProductAggregate supports LinearAlgebra.mul!, LinearAlgebra.adjoint, Base.isempty, Base.push! and most other methods that abstract matrices implement.

`KroneckerFactory`

Let us consider that we need to assemble the bilinear form associated with the operator L(u) = ∇ ⋅ (C∇u) where C = [6.0 0.1875; 0.125 12.0] is a conductivity matrix. Clearly, even on Cartesian grids an exact FastDiagonalization of the resulting matrix is not possible due to non-zero off-diagonal entries in C. Furthermore, if C is a function of the spatial coordinates, the resulting stiffness matrix does not have a Kronecker product structure at all and the material pullback needs to be approximated using some separable approximation.

KroneckerFactory can be used to obtain approximations of a bilinear form given a separable approximation of the data, i.e. weighted spline approximation of the pullback operator on mapped geometries or Wachspress approximation at quadrature points.

The syntax is close to that of IgaFormation.Sumfactory.

# domain
Ω = Interval(0.0, 1.0) ⨱ Interval(0.0, 1.0)

# partition
Δ = Partition(Ω, (3,5))

# space
p = (2,3)
S = ScalarSplineSpace(p, Δ)

# mock data
Dim = length(S)
nqpts = map((p,n) -> (p+1)*n + 2, p, num_elements.(S))
f = [ map(rand, nqpts) for k ∈ 1:Dim, l ∈ 1:Dim]

# derivative indicators
∇u = k -> ι(k, dim=Dim)
∇v = k -> ι(k, dim=Dim)

# assemble stiffness
∫ = KroneckerFactory(S, S)
for β in 1:Dim
    for α in 1:Dim
        ∫(∇u(α), ∇v(β); data=f[α, β])
    end
end

# aggregated Kronecker products
K = ∫.data

In the above, f[α, β] is a Tuple of either univariate weighting splines in each parametric direction or plain vectors with weighting coefficients. This data is used to weight the test functions.

Similarly to Sumfactory, the data collection can be reset

∫(∇u(α), ∇v(β); data=f[α, β], reset=true)

or simply

reset!(∫)

Linear solvers

Taiga provides minimal reference implementations of linear solvers:

TaigaCG Conjugate Gradient method
TaigaPCG preconditioned Conjugate Gradient method
TaigaIPCG inexactly preconditioned Conjugate Gradient method

These implementations are reasonably fast. The solution times are at least as good as Krylov.jl, whereby due to the minimal implementation there are less memory allocations. Furthermore, some solvers like TaigaIPCG are nowhere to be found and are essential to Taiga's conceptual framework.

julia> solver = TaigaPCG(L, P; atol=10e-8, rtol=10e-8, itmax=100)
Linear solver of type TaigaPCG (atol=1.0e-7, rtol=1.0e-7, itmax=100)

julia> x, stats = linsolve!(solver, b);

julia> stats
Linear solver statistics for TaigaPCG:
┌──────────────────┬───────────────────────────────┐
│ Metric           │ Value                         │
├──────────────────┼───────────────────────────────┤
│ converged        │ true                          │
│ niter            │ 47                            │
│ residual_norm    │ 9.25797e-8                    │
│ residual_norm_x₀ │ 0.0732086                     │
│ status           │ coverged: atol ✔, rtol ✘      │
└──────────────────┴───────────────────────────────┘

Index

Taiga.BezierExtractionContext — Type

BezierExtractionContext{Dim, T}

Bezier extraction context.

Fields:

splinespace::TensorProduct{Dim, SplineSpace{T}}: spline space
partition::CartesianProduct{Dim}: partition
C::Array{KroneckerProduct, Dim}: Bezier extraction operators
bezier_basis_dimension::NTuple{Dim}: Bezier basis dimension