GridapROMs.TProduct

module TProduct

Infrastructure for tensor product finite element spaces and operators on Cartesian meshes. The key idea is that a D-dimensional Cartesian domain $[a_1,b_1] \times \cdots \times [a_D,b_D]$ can be discretised as the tensor product of D 1D meshes, allowing bilinear forms to be assembled as rank tensors of 1D matrices rather than full D-dimensional sparse matrices.

Geometry

• TProductDiscreteModel: D-dimensional CartesianDiscreteModel together with its D 1D factor models.
• TProductTriangulation: triangulation on a TProductDiscreteModel.
• TProductMeasure: quadrature measure composed of D 1D measures.

FE spaces

• TensorProductReferenceFE: a ReferenceFE that carries the 1D factor reffes and triggers TProductFESpace construction via the standard FESpace(model,reffe;kwargs...) interface.
• TProductFESpace: wraps an OrderedFESpace and D 1D OrderedFESpaces.

Cell data

• TProductFEBasis: FE basis in separated form; supports gradient and PartialDerivative.
• GenericTProductCellField, GenericTProductDiffCellField: cell fields in separated form.

Rank tensors

• Rank1Tensor: $a_1 \otimes \cdots \otimes a_D$.
• GenericRankTensor: $\sum_{k=1}^K a_1^k \otimes \cdots \otimes a_D^k$.
• BlockRankTensor: multi-field variant.

Assembly

• TProductSparseMatrixAssembler: assembles 1D matrices and packs them into an AbstractRankTensor.
• TProductBlockSparseMatrixAssembler: multi-field variant.

const MatrixOrTensor = Union{AbstractMatrix,AbstractRankTensor}

abstract type AbstractRankTensor{D,K} end

Type representing a tensor a of dimension D and rank K, i.e. assuming the form

$a = \sum\limits_{k=1}^K a_1^k \otimes \cdots \otimes a_D^k$

Subtypes:

• Rank1Tensor
• GenericRankTensor

struct BlockRankTensor{A<:AbstractRankTensor,N} <: AbstractArray{A,N}
  array::Array{A,N}
end

Multi-field version of a AbstractRankTensor

struct GenericRankTensor{D,K,A<:AbstractArray} <: AbstractRankTensor{D,K}
  decompositions::Vector{Rank1Tensor{D,A}}
end

Structure representing a generic rank-K tensor, i.e. assuming the form

$a = \sum\limits_{k=1}^K a_1^k \otimes \cdots \otimes a_D^k$

struct GenericTProductCellField{DS,A,B} <: TProductCellField
  single_fields::A
  trian::B
  domain_style::DS
end

A TProductCellField in separated form: single_fields is a vector of D 1D CellFields and trian is the underlying TProductTriangulation. All 1D fields must share the same DomainStyle.

struct GenericTProductDiffCellField{O,A,B,C} <: TProductDiffCellField
  op::O
  cell_data::A
  diff_cell_data::B
  summation::C
end

A differentiated TProductCellField. Stores the original field cell_data, its differentiated counterpart diff_cell_data (e.g. the gradient of each 1D factor), the differentiation operator op, and an optional summation cache used during multi-field or mixed block assembly.

The type alias GradientTProductCellField is provided for the op = gradient case.

struct Rank1Tensor{D,A<:AbstractArray} <: AbstractRankTensor{D,1}
  factors::Vector{A}
end

Structure representing rank-1 tensors, i.e. assuming the form

$a = a_1 \otimes \cdots \otimes a_D$

abstract type TProductCellField <: CellField end

Abstract supertype for cell fields defined on a TProductTriangulation. Concrete subtypes store D 1D CellFields in separated form.

See also: GenericTProductCellField, GenericTProductDiffCellField, TProductFEBasis.

struct TProductCellPoint{DS,A,B} <: CellDatum
  point::A
  single_points::B
end

Cell quadrature points on a TProductTriangulation, storing the D-dimensional CellPoint point and a vector of D 1D CellPoints single_points. Used internally when evaluating TProductFEBasis or TProductCellField on a TProductMeasure.

struct TProductDiscreteModel{D,A,B} <: DiscreteModel{D,D}
  model::A
  models_1d::B
end

A D-dimensional CartesianDiscreteModel together with D 1D CartesianDiscreteModels whose Cartesian product reproduces it.

Use TProductTriangulation and TProductMeasure to build the corresponding integration objects, and TProductFESpace (or the TensorProductReferenceFE interface) to build the FE space.

Construction

TProductDiscreteModel(args...; kwargs...)

Accepts the same arguments as CartesianDiscreteModel: a domain tuple and a partition tuple. The 1D components are split automatically from the D-dimensional CartesianDescriptor.

Example

model = TProductDiscreteModel((0,1,0,1),(10,10))  # 10×10 Cartesian mesh on [0,1]²

struct TProductFEBasis{DS,BS,A,B} <: FEBasis
  basis::A
  trian::B
  domain_style::DS
  basis_style::BS
end

FE basis in separated tensor product form. basis is a vector of D 1D FE bases (or MultiFieldCellFields in the multi-field case), one per spatial direction. trian is the TProductTriangulation.

Construct via get_tp_fe_basis or get_tp_trial_fe_basis. Supports gradient and PartialDerivative differentiation, and integration against a TProductMeasure.

struct TProductFESpace{S} <: SingleFieldFESpace
  space::S
  spaces_1d::Vector{<:SingleFieldFESpace}
  trian::TProductTriangulation
end

A SingleFieldFESpace on a tensor product mesh, storing the D-dimensional OrderedFESpace space and a vector of D 1D OrderedFESpaces spaces_1d. The TProductTriangulation trian is stored explicitly to ensure compatibility with MultiField scenarios.

All standard SingleFieldFESpace interface methods are delegated to space. The 1D spaces are used exclusively for tensor product assembly via TProductSparseMatrixAssembler and basis extraction via get_tp_fe_basis / get_tp_trial_fe_basis.

The preferred construction path is via TensorProductReferenceFE:

model = TProductDiscreteModel((0,1,0,1),(10,10))
reffe = TensorProductReferenceFE(model,lagrangian,Float64,1)
V = FESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")

Alternatively, the reffe tuple form is still supported:

Ω  = Triangulation(model)
dΩ = Measure(Ω,2)
V  = FESpace(Ω,(lagrangian,Float64,1);conformity=:H1,dirichlet_tags="boundary")

struct TProductMeasure{A,B} <: Measure
  measure::A
  measures_1d::B
end

A Measure whose quadrature is the Cartesian product of D 1D quadratures stored in measures_1d. The full D-dimensional measure measure is also kept for standard Gridap integration.

Integrating a TProductCellField against a TProductMeasure returns a vector of D DomainContributions, one per spatial direction, which are later assembled into a AbstractRankTensor by TProductSparseMatrixAssembler.

Construct via Measure(trian::TProductTriangulation,degree;kwargs...).

struct TProductSparseMatrixAssembler{A<:SparseMatrixAssembler} <: SparseMatrixAssembler
  assems_1d::Vector{A}
end

A SparseMatrixAssembler for tensor product FE spaces. Wraps D 1D SparseMatrixAssemblers, one per spatial direction.

Assembly operates direction-by-direction: integrating a bilinear form against a TProductMeasure yields a vector of D DomainContributions, and the assembler collects and assembles each into a 1D sparse matrix. The results are combined into an AbstractRankTensor:

• A Rank1Tensor for forms without derivatives (e.g. mass matrix).
• A GenericRankTensor for forms involving gradient or PartialDerivative (e.g. stiffness matrix), where the rank equals the spatial dimension D.

Construction

TProductSparseMatrixAssembler(trial::TProductFESpace,test::TProductFESpace)
TProductSparseMatrixAssembler(mat,trial::TProductFESpace,test::TProductFESpace)
TProductSparseMatrixAssembler(mat,vec,trial,test[,strategy])

For multi-field scenarios use TProductBlockSparseMatrixAssembler.

struct TProductTriangulation{Dt,Dp,A,B,C} <: Triangulation{Dt,Dp}
  model::A
  trian::B
  trians_1d::C
end

A Triangulation whose cells are the Cartesian product of D 1D triangulations stored in trians_1d. The full D-dimensional triangulation trian and the background TProductDiscreteModel model are also stored for standard Gridap compatibility.

Construct via Triangulation(model::TProductDiscreteModel) or by wrapping an existing Triangulation with a vector of 1D triangulations.

TProductBlockSparseMatrixAssembler(trial::MultiFieldFESpace,test::MultiFieldFESpace
  ) -> TProductSparseMatrixAssembler

Returns a TProductSparseMatrixAssembler in a MultiField scenario

get_1d_tags(model::TProductDiscreteModel,tags) -> Vector{Vector{Int8}}

Fetches the tags of the tensor product 1D models corresponding to the tags of the D-dimensional model tags. The length of the output is D

get_arrays_1d(a::GenericRankTensor{D,D}) -> Vector

For a gradient-assembled GenericRankTensor built by tproduct_array(gradient,arrays,grads), recovers the original arrays_1d (mass-like factors) as a vector of length D.

get_decomposition(a::AbstractRankTensor,k::Integer) -> Vector{<:AbstractArray}

For a tensor a of dimension D and rank K assuming the form

$a = \sum\limits_{k=1}^K a_1^k \otimes \cdots \otimes a_D^k$

returns the decomposition relative to the kth rank:

$[a_1^k, \hdots , a_D^k]$

get_factor(a::AbstractRankTensor,d::Integer,k::Integer) -> AbstractArray

Returns the d-th factor array of the k-th rank-1 decomposition of a. For a Rank1Tensor k must equal 1.

get_factors(a::Rank1Tensor) -> Vector

Returns the vector of D factor arrays of the rank-1 tensor a.

get_gradients_1d(a::GenericRankTensor{D,D}) -> Vector

For a gradient-assembled GenericRankTensor built by tproduct_array(gradient,arrays,grads), recovers the original gradients_1d (stiffness-like factors) as a vector of length D.

get_tp_fe_basis(f::TProductFESpace) -> TProductFEBasis

get_tp_trial_fe_basis(f::TProductFESpace) -> TProductFEBasis

tproduct_array(arrays_1d::Vector{<:AbstractArray}) -> Rank1Tensor
tproduct_array(op,arrays_1d::Vector{<:AbstractArray},gradients_1d::Vector{<:AbstractArray},args...) -> GenericRankTensor

Returns a AbstractRankTensor storing the arrays arrays_1d (usually matrices) arising from an integration routine on D 1-d triangulations whose tensor product gives a D-dimensional triangulation. In the absence of the field gradients_1d, the output is a Rank1Tensor; when provided, the output is a GenericRankTensor

tproduct_array(arrays_1d::Vector{<:BlockArray}) -> BlockRankTensor
tproduct_array(op,arrays_1d::Vector{<:BlockArray},gradients_1d::Vector{<:BlockArray},args...) -> BlockRankTensor

Generalization of the previous functions to multi-field scenarios