# assempde

Assemble finite element matrices and solve elliptic PDE

## Syntax

• ``u = assempde(model,c,a,f)``
example
• ``u = assempde(b,p,e,t,c,a,f)``
example
• ``````[Kc,Fc,B,ud] = assempde(___)``````
example
• ``````[Ks,Fs] = assempde(___)``````
example
• ``````[K,M,F,Q,G,H,R] = assempde(___)``````
example
• ``````[K,M,F,Q,G,H,R] = assempde(___,[],sdl)``````
• ``u = assempde(K,M,F,Q,G,H,R)``
example
• ``````[Ks,Fs] = assempde(K,M,F,Q,G,H,R)``````
• ``````[Kc,Fc,B,ud] = assempde(K,M,F,Q,G,H,R)``````

## Description

example

````u = assempde(model,c,a,f)` solves the PDE$-\nabla \cdot \left(c\nabla u\right)+au=f,$with geometry, boundary conditions, and finite element mesh in `model`, and coefficients `c`, `a`, and `f`. If the PDE is a system of equations (`model.PDESystemSize` > 1), then `assempde` solves the system of equations$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$```

example

````u = assempde(b,p,e,t,c,a,f)` solves the PDE with boundary conditions `b`, and finite element mesh (`p`,`e`,`t`).```

example

``````[Kc,Fc,B,ud] = assempde(___)```, for any of the previous input syntaxes, assembles finite element matrices using the reduced linear system form, which eliminates any Dirichlet boundary conditions from the system of linear equations. You can calculate the solution `u` at node points by the command `u = B*(Kc\Fc) + ud`. See Definitions.```

example

``````[Ks,Fs] = assempde(___)``` assembles finite element matrices that represent any Dirichlet boundary conditions using a stiff-spring approximation. You can calculate the solution `u` at node points by the command `u = Ks\Fs`. See Definitions.```

example

``````[K,M,F,Q,G,H,R] = assempde(___)``` assembles finite element matrices that represent the PDE problem. This syntax returns all the matrices involved in converting the problem to finite element form. See Algorithms.```
``````[K,M,F,Q,G,H,R] = assempde(___,[],sdl)``` restricts the finite element matrices to those that include the subdomain specified by the subdomain labels in `sdl`. The empty argument is required in this syntax for historic and compatibility reasons.```

example

````u = assempde(K,M,F,Q,G,H,R)` returns the solution `u` based on the full collection of finite element matrices.```
``````[Ks,Fs] = assempde(K,M,F,Q,G,H,R)``` returns finite element matrices that approximate Dirichlet boundary conditions using the stiff-spring approximation. See Algorithms.```
``````[Kc,Fc,B,ud] = assempde(K,M,F,Q,G,H,R)``` returns finite element matrices that eliminate any Dirichlet boundary conditions from the system of linear equations. See Algorithms.```

## Examples

collapse all

### Solve a Scalar PDE

Solve an elliptic PDE on an L-shaped region.

Create a scalar PDE model. Incorporate the geometry of an L-shaped region.

```model = createpde; geometryFromEdges(model,@lshapeg);```

Apply zero Dirichlet boundary conditions to all edges.

`applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);`

Generate a finite element mesh.

`generateMesh(model);`

Solve the PDE

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

with parameters c = 1, a = 0, and f = 5.

```c = 1; a = 0; f = 5; u = assempde(model,c,a,f);```

Plot the solution.

`pdeplot(model,'xydata',u)`

### 3-D Elliptic Problem

Solve a 3-D elliptic PDE using a PDE model.

Create a PDE model container, import a 3-D geometry description, and view the geometry.

```model = createpde; importGeometry(model,'Block.stl'); h = pdegplot(model,'FaceLabels','on'); h(1).FaceAlpha = 0.5;```

Set zero Dirichlet conditions on faces 1 through 4 (the edges). Set Neumann conditions with g = –1 on face 6 and g = 1 on face 5.

```applyBoundaryCondition(model,'Face',1:4,'u',0); applyBoundaryCondition(model,'Face',6,'g',-1); applyBoundaryCondition(model,'Face',5,'g',1);```

Set coefficients c = 1, a = 0, and f = 0.1.

```c = 1; a = 0; f = 0.1;```

Create a mesh and solve the problem.

```generateMesh(model); u = assempde(model,c,a,f);```

Plot the solution on the surface.

`pdeplot3D(model,'colormapdata',u);`

### 2-D PDE Using [p,e,t] Mesh

Solve a 2-D PDE using the older syntax for mesh.

Create a circle geometry.

`g = @circleg;`

Set zero Dirichlet boundary conditions.

`b = @circleb1;`

Create a mesh for the geometry.

`[p,e,t] = initmesh(g);`

Solve the PDE

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

with parameters c = 1, a = 0, and f = sin(x).

```c = 1; a = 0; f = 'sin(x)'; u = assempde(b,p,e,t,c,a,f);```

Plot the solution.

`pdeplot(p,e,t,'xydata',u)`

### Finite Element Matrices

Obtain the finite-element matrices that represent the problem using a reduced linear algebra representation of Dirichlet boundary conditions.

Create a scalar PDE model. Import a simple 3-D geometry.

```model = createpde; importGeometry(model,'Block.stl');```

Set zero Dirichlet boundary conditions on all the geometry faces.

`applyBoundaryCondition(model,'Face',1:model.Geometry.NumFaces,'u',0);`

Generate a mesh for the geometry.

`generateMesh(model);`

Obtain finite element matrices `K`, `F`, `B`, and `ud` that represent the equation

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

with parameters c = 1, a = 0, and f = log(1 + x + y/(1+z)).

```c = 1; a = 0; f = 'log(1+x+y./(1+z))'; [K,F,B,ud] = assempde(model,c,a,f);```

You can obtain the solution `u` of the PDE at mesh nodes by executing the command

`u = B*(K\F) + ud;`

Generally, this solution is slightly more accurate than the stiff-spring solution, as calculated in Stiff-Spring Finite Element Solution.

### Stiff-Spring Finite Element Solution

Obtain the stiff-spring approximation of finite element matrices.

Create a scalar PDE model. Import a simple 3-D geometry.

```model = createpde; importGeometry(model,'Block.stl');```

Set zero Dirichlet boundary conditions on all the geometry faces.

`applyBoundaryCondition(model,'Face',1:model.Geometry.NumFaces,'u',0);`

Generate a mesh for the geometry.

`generateMesh(model);`

Obtain finite element matrices `Ks` and `Fs` that represent the equation

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

with parameters c = 1, a = 0, and f = log(1 + x + y/(1+z)).

```c = 1; a = 0; f = 'log(1+x+y./(1+z))'; [Ks,Fs] = assempde(model,c,a,f);```

You can obtain the solution `u` of the PDE at mesh nodes by executing the command

`u = Ks\Fs;`

Generally, this solution is slightly less accurate than the reduced linear algebra solution, as calculated in Finite Element Matrices.

### Full Collection of Finite Element Matrices

Obtain the full collection of finite element matrices for an elliptic problem.

Import geometry and set up an elliptic problem with Dirichlet boundary conditions. The `'Torus.stl'` geometry has only one face, so you need set only one boundary condition.

```model = createpde(); importGeometry(model,'Torus.stl'); applyBoundaryCondition(model,'Face',1,'u',0); c = 1; a = 0; f = 1; generateMesh(model);```

Create the finite element matrices that represent this problem.

`[K,M,F,Q,G,H,R] = assempde(model,c,a,f);`

Most of the resulting matrices are quite sparse. `G`, `M`, `Q`, and `R` are all zero sparse matrices.

```howsparse = @(x)nnz(x)/numel(x); disp(['Maximum fraction of nonzero entries in K or H is ',num2str(max(howsparse(K),howsparse(H)))])```
`Maximum fraction of nonzero entries in K or H is 0.00086641`

To find the solution to the PDE, call `assempde` again.

`u = assempde(K,M,F,Q,G,H,R);`

## Input Arguments

collapse all

### `model` — PDE model`PDEModel` object

PDE model, specified as a `PDEModel` object.

Example: `model = createpde(1)`

### `c` — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. `c` represents the c coefficient in the scalar PDE

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specify`c` in various ways, detailed in c Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: `'cosh(x+y.^2)'`

Data Types: `double` | `char` | `function_handle`
Complex Number Support: Yes

### `a` — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. `a` represents the a coefficient in the scalar PDE

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specify`a` in various ways, detailed in a or d Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: `2*eye(3)`

Data Types: `double` | `char` | `function_handle`
Complex Number Support: Yes

### `f` — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. `f` represents the f coefficient in the scalar PDE

$-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specify`f` in various ways, detailed in f Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: `char('sin(x)';'cos(y)';'tan(z)')`

Data Types: `double` | `char` | `function_handle`
Complex Number Support: Yes

### `b` — Boundary conditionsboundary matrix | boundary file

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a string naming the file.

For more information on boundary conditions, see Forms of Boundary Condition Specification.

Example: `b = 'circleb1'` or equivalently ```b = @circleb1```

Data Types: `double` | `char` | `function_handle`

### `p` — Mesh nodesoutput of `initmesh` | output of `meshToPet`

Mesh nodes, specified as the output of `initmesh` or `meshToPet`. For the structure of a `p` matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: `[p,e,t] = initmesh(g)`

Data Types: `double`

### `e` — Mesh edgesoutput of `initmesh` | output of `meshToPet`

Mesh edges, specified as the output of `initmesh` or `meshToPet`. For the structure of `e`, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: `[p,e,t] = initmesh(g)`

Data Types: `double`

### `t` — Mesh elementsoutput of `initmesh` | output of `meshToPet`

Mesh elements, specified as the output of `initmesh` or `meshToPet`. Mesh elements are the triangles or tetrahedra that form the finite element mesh. For the structure of a `t` matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: `[p,e,t] = initmesh(g)`

Data Types: `double`

### `K` — Stiffness matrixsparse matrix | full matrix

Stiffness matrix, specified as a sparse matrix or full matrix. Generally, you obtain `K` from a previous call to `assema` or `assempde`. For the meaning of stiffness matrix, see Elliptic Equations.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `M` — Mass matrixsparse matrix | full matrix

Mass matrix, specified as a sparse matrix or full matrix. Generally, you obtain `M` from a previous call to `assema` or `assempde`. For the meaning of mass matrix, see Elliptic Equations.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `F` — Finite element `f` representationvector

Finite element `f` representation, specified as a vector. Generally, you obtain `F` from a previous call to `assema` or `assempde`. For the meaning of this representation, see Elliptic Equations.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `Q` — Neumann boundary condition matrixsparse matrix | full matrix

Neumann boundary condition matrix, specified as a sparse matrix or full matrix. Generally, you obtain `Q` from a previous call to `assemb` or `assempde`. For the meaning of this matrix, see Elliptic Equations.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `G` — Neumann boundary condition vectorsparse vector | full vector

Neumann boundary condition vector, specified as a sparse vector or full vector. Generally, you obtain `G` from a previous call to `assemb` or `assempde`. For the meaning of this vector, see Elliptic Equations.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `H` — Dirichlet boundary condition matrixsparse matrix | full matrix

Dirichlet boundary condition matrix, specified as a sparse matrix or full matrix. Generally, you obtain `H` from a previous call to `assemb` or `assempde`. For the meaning of this matrix, see Algorithms.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `R` — Dirichlet boundary condition vectorsparse vector | full vector

Dirichlet boundary condition vector, specified as a sparse vector or full vector. Generally, you obtain `R` from a previous call to `assemb` or `assempde`. For the meaning of this vector, see Algorithms.

Example: `[K,M,F,Q,G,H,R] = assempde(model,c,a,f)`

Data Types: `double`
Complex Number Support: Yes

### `sdl` — Subdomain labelsvector of positive integers

Subdomain labels, specified as a vector of positive integers. For 2-D geometry only. View the subdomain labels in your geometry using the command

`pdegplot(g,'SubdomainLabels','on')`

Example: `sdl = [1,3:5];`

Data Types: `double`

## Output Arguments

collapse all

### `u` — PDE solutionvector

PDE solution, returned as a vector.

• If the PDE is scalar, meaning only one equation, then `u` is a column vector representing the solution u at each node in the mesh. `u(i)` is the solution at the `i`th column of `model.Mesh.Nodes` or the `i`th column of `p`.

• If the PDE is a system of N > 1 equations, then `u` is a column vector with N*`Np` elements, where `Np` is the number of nodes in the mesh. The first `Np` elements of `u` represent the solution of equation 1, then next `Np` elements represent the solution of equation 2, etc.

To obtain the solution at an arbitrary point in the geometry, use `pdeInterpolant`.

To plot the solution, use `pdeplot` for 2-D geometry, or see Plot 3-D Solutions.

### `Kc` — Stiffness matrixsparse matrix

Stiffness matrix, returned as a sparse matrix. See Elliptic Equations.

`u1 = Kc\Fc` returns the solution on the non-Dirichlet points. To obtain the solution `u` at the nodes of the mesh,

`u = B*(Kc\Fc) + ud`

Generally, `Kc`, `Fc`, `B`, and `ud` make a slower but more accurate solution than `Ks` and `Fs`.

### `Fc` — Load vectorvector

Load vector, returned as a vector. See Elliptic Equations.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`, `Fc`, `B`, and `ud` make a slower but more accurate solution than `Ks` and `Fs`.

### `B` — Dirichlet nullspacesparse matrix

Dirichlet nullspace, returned as a sparse matrix. See Algorithms.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`, `Fc`, `B`, and `ud` make a slower but more accurate solution than `Ks` and `Fs`.

### `ud` — Dirichlet vectorvector

Dirichlet vector, returned as a vector. See Algorithms.

`u = B*(Kc\Fc) + ud`

Generally, `Kc`, `Fc`, `B`, and `ud` make a slower but more accurate solution than `Ks` and `Fs`.

### `Ks` — Stiffness matrix corresponding to the stiff-spring approximation for Dirichlet boundary conditionsparse matrix

Finite element matrix for stiff-spring approximation, returned as a sparse matrix. See Algorithms.

To obtain the solution `u` at the nodes of the mesh,

`u = Ks\Fs`.

Generally, `Ks` and `Fs` make a quicker but less accurate solution than `Kc`, `Fc`, `B`, and `ud`.

### `Fs` — Load vector corresponding to the stiff-spring approximation for Dirichlet boundary conditionvector

Load vector corresponding to the stiff-spring approximation for Dirichlet boundary condition, returned as a vector. See Algorithms.

To obtain the solution `u` at the nodes of the mesh,

`u = Ks\Fs`.

Generally, `Ks` and `Fs` make a quicker but less accurate solution than `Kc`, `Fc`, `B`, and `ud`.

### `K` — Stiffness matrixsparse matrix

Stiffness matrix, returned as a sparse matrix. See Elliptic Equations.

`K` represents the stiffness matrix alone, unlike `Kc` or `Ks`, which are stiffness matrices combined with other terms to enable immediate solution of a PDE.

Typically, you use `K` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `M` — Mass matrixsparse matrix

Mass matrix. returned as a sparse matrix. See Elliptic Equations.

Typically, you use `M` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `F` — Load vectorvector

Load vector, returned as a vector. See Elliptic Equations.

`F` represents the load vector alone, unlike `Fc` or `Fs`, which are load vectors combined with other terms to enable immediate solution of a PDE.

Typically, you use `F` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `Q` — Neumann boundary condition matrixsparse matrix

Neumann boundary condition matrix, returned as a sparse matrix. See Elliptic Equations.

Typically, you use `Q` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `G` — Neumann boundary condition vectorsparse vector

Neumann boundary condition vector, returned as a sparse vector. See Elliptic Equations.

Typically, you use `G` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `H` — Dirichlet matrixsparse matrix

Dirichlet matrix, returned as a sparse matrix. See Algorithms.

Typically, you use `H` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

### `R` — Dirichlet vectorsparse vector

Dirichlet vector, returned as a sparse vector. See Algorithms.

Typically, you use `R` in a subsequent call to a solver such as `assempde` or `hyperbolic`.

expand all

### Reduced Linear System

This form of the finite element matrices eliminates Dirichlet conditions from the problem using a linear algebra approach detailed in Systems of PDEs. The finite element matrices reduce to the solution

`u = B*(Kc\Fc) + ud`

where `B` spans the null space of the columns of `H` (the Dirichlet condition matrix representing hu = r). `R` is the Dirichlet condition vector for `Hu = R`. `ud` is the vector of boundary condition solutions for the Dirichlet conditions. `u1` = `Kc\Fc` returns the solution on the non-Dirichlet points.

### Stiff-Spring Approximation

This form of the finite element matrices converts Dirichlet boundary conditions to Neumann boundary conditions using a stiff-spring approximation. Using this approximation, `assempde` returns a matrix `Ks` and a vector `Fs` that represent the combined finite element matrices. The approximate solution `u` is

`u = Ks\Fs`.

See Elliptic Equations. For details of the stiff-spring approximation, see Systems of PDEs.

### Algorithms

`assempde` performs the following steps to obtain a solution `u` to an elliptic PDE:

1. Generate the finite element matrices [`K`,`M`,`F`,`Q`,`G`,`H`,`R`]. This step is equivalent to calling `assema` to generate the matrices `K`, `M`, and `F`, and also calling `assemb` to generate the matrices `Q`, `G`, `H`, and `R`.

2. Generate the combined finite element matrices [`Kc`,`Fc`,`B`,`ud`]. The combined stiffness matrix is for the reduced linear system, ```Kc = K + M + Q```. The corresponding combined load vector is ```Fc = F + G```. The `B` matrix spans the null space of the columns of `H` (the Dirichlet condition matrix representing hu = r). The `R` vector represents the Dirichlet conditions in `Hu = R`. The `ud` vector represents boundary condition solutions for the Dirichlet conditions.

3. Calculate the solution `u` via

`u = B*(Kc\Fc) + ud`.

`assempde` uses one of two algorithms for assembling a problem into combined finite element matrix form. A reduced linear system form leads to immediate solution via linear algebra. You choose the algorithm by the number of outputs. For the reduced linear system form, request four outputs:

`[Kc,Fc,B,ud] = assempde(_)`

For the stiff-spring approximation, request two outputs:

`[Ks,Fs] = assempde(_)`

For details, see Definitions.

As explained in Elliptic Equations, the full finite element matrices and vectors are the following.

• `K` is the stiffness matrix, the integral of the `c` coefficient against the basis functions.

• `M` is the mass matrix, the integral of the `a` coefficient against the basis functions.

• `F` is the integral of the `f` coefficient against the basis functions.

• `Q` is the integral of the `q` boundary condition against the basis functions.

• `G` is the integral of the `g` boundary condition against the basis functions.

• The `H` and `R` matrices come directly from the Dirichlet conditions and the mesh. See Systems of PDEs.