Assembling

Assembling functions for B-spline IgA.

This module contains functions to assemble stiffness matrices and load vectors for isogeometric discretizations. Furthermore, it provides support for computing and eliminating Dirichlet dofs from a linear system.

Assembling general problems

General variational forms can be assembled using the following function. See the section Assembling custom forms for further details.

pyiga.assemble.assemble(problem, kvs, args=None, bfuns=None, boundary=None, symmetric=False, format='csr', layout='blocked', **kwargs)

Assemble a matrix or vector in a function space.

Parameters:
  • problem –

    the description of the variational form to assemble. It can be passed in a number of formats (see Assembling custom forms for details):

    • string: a textual description of the variational form
    • VForm: an abstract description of the variational form
    • assembler class (result of compiling a VForm using pyiga.compile.compile_vform())
    • assembler object (result of instantiating an assembler class with a concrete space and input functions) – in this case kvs and args are ignored
  • kvs –

    the space or spaces in which to assemble the problem. Either:

    • a tuple of KnotVector instances, describing a tensor product spline space
    • if the variational form requires more than one space (e.g., for a Petrov-Galerkin discretization), then a tuple of such tuples, each describing one tensor product spline space (usually one for the trial space and one for the test space)
    • an HSpace instance for problems in hierarchical spline spaces
  • args (dict) –

    a dictionary which provides named inputs for the assembler. Most problems will require at least a geometry map; this can be given in the form {'geo': geo}, where geo is a geometry function defined using the pyiga.geometry module. Further values used in the problem description must be passed here.

    For convenience, any additional keyword arguments to this function are added to the args dict automatically.

  • bfuns –

    a list of used basis functions. By default, scalar basis functions ‘u’ and ‘v’ are assumed, and the arity of the variational form is determined automatically based on whether one or both of these functions are used. Otherwise, bfuns should be a list of tuples (name, components, space), where name is a string, components is an integer describing the number of components the basis function has, and space is an integer referring to which input space the function lives in. Shorter tuples are valid, in which case the components default to components=1 (scalar basis function) and space=0 (all functions living in the same, first, input space).

    This argument is only used if problem is given as a string.

For the meaning of the remaining arguments and the format of the output, refer to assemble_entries().

pyiga.assemble.assemble_vf(vf, kvs, symmetric=False, format='csr', layout='blocked', args=None, **kwargs)

Compile the given variational form (VForm) into a matrix or vector.

Any named inputs defined in the vform must be given in the args dict or as keyword arguments. For the meaning of the remaining arguments, refer to assemble_entries().

pyiga.assemble.assemble_entries(asm, symmetric=False, format='csr', layout='blocked')

Given an instance asm of an assembler class, assemble all entries and return the resulting matrix or vector.

Parameters:
  • asm – an instance of an assembler class, e.g. one compiled using pyiga.compile.compile_vform()
  • symmetric (bool) – (matrices only) exploit symmetry of the matrix to speed up the assembly
  • format (str) – (matrices only) the sparse matrix format to use; default ‘csr’
  • layout (str) –

    (vector-valued problems only): the layout of the generated matrix or vector. Valid options are:

    • ’blocked’: the matrix is laid out as a k_1 x k_2 block matrix, where k_1 and k_2 are the number of components of the test and trial functions, respectively
    • ’packed’: the interactions of the components are packed together, i.e., each entry of the matrix is a small k_1 x k_2 block
Returns:

ndarray or sparse matrix –

  • if the assembler has arity=1: an ndarray of vector entries whose shape is given by the number of degrees of freedom per coordinate direction. For vector-valued problem, an additional final axis is added which has the number of components as its length.
  • if the assembler has arity=2: a sparse matrix in the given format
class pyiga.assemble.Assembler(problem, kvs, args=None, bfuns=None, boundary=None, symmetric=False, updatable=[], **kwargs)

A high-level interface to an assembler class.

Usually, you will simply call assemble() to assemble a matrix or vector. When assembling a sequence of problems with changing parameters, it may be more efficient to instantiate the assembler class only once and inform it of the changing inputs via the updatable argument. You may then call update() to change these input fields and assemble the problem via assemble(). Instead of calling update() explicitly, you can also simply pass the fields to be updated as additional keyword arguments to assemble().

updatable is a list of names of assembler arguments which are to be considered updatable. The other arguments have the same meaning as for assemble().

assemble(format='csr', layout='blocked', **upd_fields)

Assemble the problem.

Any additional keyword arguments are passed to update().

update(**kwargs)

Update all input fields given as name=func keyword arguments.

All fields updated in this manner must have been specified in the list of updatable arguments when creating the Assembler object.

Tensor product Gauss quadrature assemblers

Standard Gauss quadrature assemblers for mass and stiffness matrices. They take one or two arguments:

  • kvs (list of KnotVector): Describes the tensor product B-spline basis for which to assemble the matrix. One KnotVector per coordinate direction. In the 1D case, a single KnotVector may be passed directly.
  • geo (BSplineFunc or NurbsFunc; optional): Geometry transform, mapping from the parameter domain to the physical domain. If omitted, assume the identity map; a fast Kronecker product implementation is used in this case.
pyiga.assemble.mass(kvs, geo=None, format='csr')

Assemble a mass matrix for the given basis (B-spline basis or tensor product B-spline basis) with an optional geometry transform.

pyiga.assemble.stiffness(kvs, geo=None, format='csr')

Assemble a stiffness matrix for the given basis (B-spline basis or tensor product B-spline basis) with an optional geometry transform.

Fast low-rank assemblers

Fast low-rank assemblers based on the paper “A Black-Box Algorithm for Fast Matrix Assembly in Isogeometric Analysis”. They may achieve significant speedups over the classical Gauss assemblers, in particular for fine discretizations and higher spline degrees. They only work well if the geometry transform is rather smooth so that the resulting matrix has relatively low (numerical) Kronecker rank.

They take the following additional arguments:

  • tol: the stopping accuracy for the Adaptive Cross Approximation (ACA) algorithm
  • maxiter: the maximum number of ACA iterations
  • skipcount: terminate after finding this many successive near-zero pivots
  • tolcount: terminate after finding this many successive pivots below the desired accuracy
  • verbose: the amount of output to display. 0 is silent, 1 prints basic information, 2 prints detailed information
pyiga.assemble.mass_fast(kvs, geo=None, tol=1e-10, maxiter=100, skipcount=3, tolcount=3, verbose=2)

Assemble a mass matrix for the given tensor product B-spline basis with an optional geometry transform, using the fast low-rank assembling algorithm.

pyiga.assemble.stiffness_fast(kvs, geo=None, tol=1e-10, maxiter=100, skipcount=3, tolcount=3, verbose=2)

Assemble a stiffness matrix for the given tensor product B-spline basis with an optional geometry transform, using the fast low-rank assembling algorithm.

Right-hand sides

pyiga.assemble.inner_products(kvs, f, f_physical=False, geo=None)

Compute the \(L_2\) inner products between each basis function in a tensor product B-spline basis and the function f (i.e., the load vector).

Parameters:
  • kvs (seq) – a sequence of KnotVector, representing a tensor product basis
  • f – a function or BSplineFunc object
  • f_physical (bool) – whether f is given in physical coordinates. If True, geo must be passed as well.
  • geo – a BSplineFunc or NurbsFunc which describes the integration domain; if not given, the integrals are computed in the parameter domain
Returns:

ndarray – the inner products as an array of size kvs[0].ndofs × kvs[1].ndofs × … × kvs[-1].ndofs. Each entry corresponds to the inner product of the corresponding basis function with f. If f is not scalar, then each of its components is treated separately and the corresponding dimensions are appended to the end of the return value.

Boundary and initial conditions

pyiga.assemble.compute_dirichlet_bcs(kvs, geo, bdconds)

Compute indices and values for Dirichlet boundary conditions on several boundaries at once.

Parameters:
  • kvs – a tensor product B-spline basis
  • geo (BSplineFunc or NurbsFunc) – the geometry transform
  • bdconds – a list of (bdspec, dir_func) pairs, where bdspec specifies the boundary to apply a Dirichlet boundary condition to and dir_func is the function providing the Dirichlet values. For the exact meaning, refer to compute_dirichlet_bc(). As a shorthand, it is possible to pass a single pair ("all", dir_func) which applies Dirichlet boundary conditions to all boundaries.
Returns:

A pair (indices, values) suitable for passing to RestrictedLinearSystem.
pyiga.assemble.compute_dirichlet_bc(kvs, geo, bdspec, dir_func)

Compute indices and values for a Dirichlet boundary condition using interpolation.

Parameters:
  • kvs – a tensor product B-spline basis
  • geo (BSplineFunc or NurbsFunc) – the geometry transform
  • bdspec –

    a pair (axis, side). axis denotes the axis along which the boundary condition lies, and side is either 0 for the “lower” boundary or 1 for the “upper” boundary. Alternatively, one of the following six strings can be used for bdspec:

    value Meaning
    "left" x low
    "right" x high
    "bottom" y low
    "top" y high
    "front" z low
    "back" z high
  • dir_func – a function which will be interpolated to obtain the Dirichlet boundary values. Assumed to be given in physical coordinates. If it is vector-valued, one Dirichlet dof is computed per component, and they are numbered according to the “blocked” matrix layout. If dir_func is a scalar value, a constant function with that value is assumed.
Returns:

A pair of arrays (indices, values) which denote the indices of the dofs within the tensor product basis which lie along the Dirichlet boundary and their computed values, respectively.
pyiga.assemble.compute_initial_condition_01(kvs, geo, bdspec, g0, g1, physical=True)

Compute indices and values for an initial condition including function value and derivative for a space-time discretization using interpolation. This only works for a space-time cylinder with constant (in time) geometry. To be precise, the space-time geometry transform geo should have the form

\[G(\vec x, t) = (\widetilde G(\vec x), t).\]
Parameters:
  • kvs – a tensor product B-spline basis
  • geo (BSplineFunc or NurbsFunc) – the geometry transform of the space-time cylinder
  • bdspec – a pair (axis, side). axis denotes the time axis of geo, and side is either 0 for the “lower” boundary or 1 for the “upper” boundary.
  • g0 – a function which will be interpolated to obtain the initial function values
  • g1 – a function which will be interpolated to obtain the initial derivatives.
  • physical (bool) – whether the functions g0 and g1 are given in physical (True) or parametric (False) coordinates. Physical coordinates are assumed by default.
Returns:

A pair of arrays (indices, values) which denote the indices of the dofs within the tensor product basis which lie along the initial face of the space-time cylinder and their computed values, respectively.
pyiga.assemble.combine_bcs(bcs)

Given a sequence of (indices, values) pairs such as returned by compute_dirichlet_bc(), combine them into a single pair (indices, values).

Dofs which occur in more than one indices array take their value from an arbitrary corresponding values array.

class pyiga.assemble.RestrictedLinearSystem(A, b, bcs, elim_rows=None)

Represents a linear system with some of its dofs eliminated.

Parameters:
  • A – the full matrix
  • b – the right-hand side (may be 0)
  • bcs – a pair of arrays (indices, values) which contain the indices and values, respectively, of dofs to be eliminated from the system
  • elim_rows – for Petrov-Galerkin discretizations, the equations to be eliminated from the linear system may not match the dofs to be eliminated. In this case, an array of indices of rows to be eliminated may be passed in this argument.

Once constructed, the restricted linear system can be accessed through the following attributes:

A

the restricted matrix

b

the restricted and updated right-hand side

complete(u)

Given a solution u of the restricted linear system, complete it with the values of the eliminated dofs to a solution of the original system.

extend(u)

Given a vector u containing only the free dofs, pad it with zeros to all dofs.

restrict(u)

Given a vector u containing all dofs, return its restriction to the free dofs.

restrict_matrix(B)

Given a matrix B which operates on all dofs, return its restriction to the free dofs.

restrict_rhs(f)

Given a right-hand side vector f, return its restriction to the non-eliminated rows.

If elim_rows was not passed, this is equivalent to restrict().

Multipatch problems

class pyiga.assemble.Multipatch(patches, automatch=False)

Represents a multipatch structure, consisting of a number of patches together with their discretizations and the information about shared dofs between patches. Currently, only conforming patches are supported.

Parameters:
  • patches – a list of (kvs, geo) pairs, each of which describes a patch. Here kvs is a tuple of KnotVector instances describing the tensor product spline space used for discretization and geo is the geometry map.
  • automatch (bool) – if True, attempt to automatically determine the matching interfaces between all patches and finalize then. If False, the user has to manually join the patches by calling join_boundaries() as often as needed, followed by finalize().
assemble_system(problem, rhs, args=None, bfuns=None, symmetric=False, format='csr', layout='blocked', **kwargs)

Assemble both the system matrix and the right-hand side vector for a variational problem over the multipatch geometry.

Here problem represents a bilinear form and rhs a linear functional. See assemble() for the precise meaning of the arguments.

Returns:A pair (A, b) consisting of the sparse system matrix and the right-hand side vector.
compute_dirichlet_bcs(bdconds)

Performs the same operation as the global function compute_dirichlet_bcs(), but for a multipatch problem.

The sequence bdconds should contain triples of the form (patch, bdspec, dir_func).

Returns:A pair (indices, values) suitable for passing to RestrictedLinearSystem.
finalize()

After all shared dofs have been declared using join_boundaries() or join_dofs(), call this function to set up the internal data structures.

global_to_patch(p)

Compute a sparse binary matrix which maps global dofs to local dofs in patch p. This is just the transpose of patch_to_global() and also its left-inverse.

Parameters:p (int) – the index of the patch
Returns:a sparse matrix
join_boundaries(p1, bdspec1, p2, bdspec2, flip=None)

Join the dofs lying along boundary bdspec1 of patch p1 with those lying along boundary bdspec2 of patch p2.

See compute_dirichlet_bc() for the format of the boundary specification.

If flip is given, it should be a sequence of booleans indicating for each coordinate axis of the boundary if the coordinates of p2 have to be flipped along that axis.

join_dofs(p1, I1, p2, I2)

Join the dofs I1 of patch p1 with the dofs I2 of patch p2.

numdofs

Number of dofs after eliminating shared dofs.

May only be called after finalize().

numpatches

Number of patches in the multipatch structure.

patch_to_global(p, j_global=False)

Compute a sparse binary matrix which maps dofs local to patch p to the corresponding global dofs.

Parameters:
  • p (int) – the index of the patch
  • j_global (bool) – if False, the matrix has only as many columns as p has dofs; if True, the number of columns is the sum of the number of local dofs over all patches
Returns:

a CSR sparse matrix
patch_to_global_idx(p)

Return an array which maps local tensor product indices for patch p to global indices.

Integration

pyiga.assemble.integrate(kvs, f, f_physical=False, geo=None)

Compute the integral of the function f over the geometry geo or a simple tensor product domain.

Parameters:
  • kvs (seq) – a sequence of KnotVector; determines the parameter domain and the quadrature rule
  • f – a function or BSplineFunc object
  • f_physical (bool) – whether f is given in physical coordinates. If True, geo must be passed as well.
  • geo – a BSplineFunc or NurbsFunc which describes the integration domain; if not given, the integral is computed in the parameter domain
Returns:

float – the integral of f over the specified domain