Main Content

Assemble finite element matrices

assembles finite element matrices using the input
time or solution specified in the
`FEM`

= assembleFEMatrices(___,`state`

)`state`

structure array. The
function uses the `time`

field of
the structure for time-dependent models and the
solution field `u`

for nonlinear
models. You can use this argument with any of the
previous syntaxes.

Create a PDE model for the Poisson equation on an L-shaped membrane with zero Dirichlet boundary conditions.

model = createpde(1); geometryFromEdges(model,@lshapeg); specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1); applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);

Generate a mesh and obtain the default finite element matrices for the problem and mesh.

```
generateMesh(model,'Hmax',0.2);
FEM = assembleFEMatrices(model)
```

`FEM = `*struct with fields:*
K: [401x401 double]
A: [401x401 double]
F: [401x1 double]
Q: [401x401 double]
G: [401x1 double]
H: [80x401 double]
R: [80x1 double]
M: [401x401 double]
T: [401x401 double]

Make computations faster by specifying which finite element matrices to assemble.

Create a transient thermal model and include the geometry of the built-in function `squareg`

.

thermalmodel = createpde('thermal','steadystate'); geometryFromEdges(thermalmodel,@squareg);

Plot the geometry with the edge labels.

pdegplot(thermalmodel,'EdgeLabels','on') xlim([-1.1 1.1]) ylim([-1.1 1.1])

Specify the thermal conductivity of the material and the internal heat source.

```
thermalProperties(thermalmodel,'ThermalConductivity',0.2);
internalHeatSource(thermalmodel,10);
```

Set the boundary conditions.

thermalBC(thermalmodel,'Edge',[1,3],'Temperature',100);

Generate a mesh.

generateMesh(thermalmodel);

Assemble the stiffness and mass matrices.

`FEM_KM = assembleFEMatrices(thermalmodel,'KM')`

`FEM_KM = `*struct with fields:*
K: [1541x1541 double]
M: [1541x1541 double]

Now, assemble the finite element matrices M, K, A, and F.

`FEM_MKAF = assembleFEMatrices(thermalmodel,'MKAF')`

`FEM_MKAF = `*struct with fields:*
M: [1541x1541 double]
K: [1541x1541 double]
A: [1541x1541 double]
F: [1541x1 double]

The four matrices M, K, A, and F correspond to discretized versions of the PDE coefficients m, c, a, and f. These four matrices also represent the domain of the finite-element model of the PDE. Instead of specifying them explicitly, you can use the `domain`

argument.

`FEMd = assembleFEMatrices(thermalmodel,'domain')`

`FEMd = `*struct with fields:*
M: [1541x1541 double]
K: [1541x1541 double]
A: [1541x1541 double]
F: [1541x1 double]

The four matrices Q, G, H, and R, correspond to discretized versions of q, g, h, and r in the Neumann and Dirichlet boundary condition specification. These four matrices also represent the boundary of the finite-element model of the PDE. Use the `boundary`

argument to assemble only these matrices.

`FEMb = assembleFEMatrices(thermalmodel,'boundary')`

`FEMb = `*struct with fields:*
H: [74x1541 double]
R: [74x1 double]
G: [1541x1 double]
Q: [1541x1541 double]

`nullspace`

and `stiff-spring`

MethodsCreate a PDE model for the Poisson equation on an L-shaped membrane with zero Dirichlet boundary conditions.

model = createpde(1); geometryFromEdges(model,@lshapeg); specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1); applyBoundaryCondition(model,'Edge',1:model.Geometry.NumEdges,'u',0);

Generate a mesh.

`generateMesh(model,'Hmax',0.2);`

Obtain the finite element matrices after imposing the boundary condition using the null-space approach. This approach eliminates the Dirichlet degrees of freedom and provides a reduced system of equations.

`FEMn = assembleFEMatrices(model,'nullspace')`

`FEMn = `*struct with fields:*
Kc: [321x321 double]
Fc: [321x1 double]
B: [401x321 double]
ud: [401x1 double]
M: [321x321 double]

Obtain the solution to the PDE using the `nullspace`

finite element matrices.

un = FEMn.B*(FEMn.Kc\FEMn.Fc) + FEMn.ud;

Compare this result to the solution given by `solvepde`

. The two solutions are identical.

u1 = solvepde(model); norm(un - u1.NodalSolution)

ans = 0

Obtain the finite element matrices after imposing the boundary condition using the stiff-spring approach. This approach retains the Dirichlet degrees of freedom, but imposes a large penalty on them.

`FEMs = assembleFEMatrices(model,'stiff-spring')`

`FEMs = `*struct with fields:*
Ks: [401x401 double]
Fs: [401x1 double]
M: [401x401 double]

Obtain the solution to the PDE using the stiff-spring finite element matrices. This technique gives a less accurate solution.

us = FEMs.Ks\FEMs.Fs; norm(us - u1.NodalSolution)

ans = 0.0098

Assemble finite element matrices for the first and last time steps of a transient structural problem.

Create a transient structural model for solving a solid (3-D) problem.

structuralmodel = createpde('structural','transient-solid');

Create the geometry and include it in the model. Plot the geometry.

gm = multicylinder(0.01,0.05); addVertex(gm,'Coordinates',[0,0,0.05]); structuralmodel.Geometry = gm; pdegplot(structuralmodel,'FaceLabels','on','FaceAlpha',0.5)

Specify the Young's modulus and Poisson's ratio.

structuralProperties(structuralmodel,'Cell',1,'YoungsModulus',201E9, ... 'PoissonsRatio',0.3, ... 'MassDensity',7800);

Specify that the bottom of the cylinder is a fixed boundary.

structuralBC(structuralmodel,'Face',1,'Constraint','fixed');

Specify the harmonic pressure on the top of the cylinder.

structuralBoundaryLoad(structuralmodel,'Face',2,'Pressure',5E7,'Frequency',50);

Specify the zero initial displacement and velocity.

structuralIC(structuralmodel,'Displacement',[0;0;0],'Velocity',[0;0;0]);

Generate a linear mesh.

generateMesh(structuralmodel,'GeometricOrder','linear'); tlist = linspace(0,1,300);

Assemble the finite element matrices for the initial time step.

state.time = tlist(1); FEM_domain = assembleFEMatrices(structuralmodel,state)

`FEM_domain = `*struct with fields:*
K: [6609x6609 double]
A: [6609x6609 double]
F: [6609x1 double]
Q: [6609x6609 double]
G: [6609x1 double]
H: [252x6609 double]
R: [252x1 double]
M: [6609x6609 double]

Pressure applied at the top of the cylinder is the only time-dependent quantity in the model. To model the dynamics of the system, assemble the boundary-load finite element matrix G for the initial, intermediate, and final time steps.

```
state.time = tlist(1);
FEM_boundary_init = assembleFEMatrices(structuralmodel,'G',state)
```

`FEM_boundary_init = `*struct with fields:*
G: [6609x1 double]

```
state.time = tlist(floor(length(tlist)/2));
FEM_boundary_half = assembleFEMatrices(structuralmodel,'G',state)
```

`FEM_boundary_half = `*struct with fields:*
G: [6609x1 double]

```
state.time = tlist(end);
FEM_boundary_final = assembleFEMatrices(structuralmodel,'G',state)
```

`FEM_boundary_final = `*struct with fields:*
G: [6609x1 double]

Assemble finite element matrices for a heat transfer problem with temperature-dependent thermal conductivity.

Create a steady-state thermal model.

thermalmodelS = createpde('thermal','steadystate');

Create a 2-D geometry by drawing one rectangle the size of the block and a second rectangle the size of the slot.

r1 = [3 4 -.5 .5 .5 -.5 -.8 -.8 .8 .8]; r2 = [3 4 -.05 .05 .05 -.05 -.4 -.4 .4 .4]; gdm = [r1; r2]';

Subtract the second rectangle from the first to create the block with a slot.

g = decsg(gdm,'R1-R2',['R1'; 'R2']');

Convert the `decsg`

format into a geometry object. Include the geometry in the model and plot the geometry.

geometryFromEdges(thermalmodelS,g); figure pdegplot(thermalmodelS,'EdgeLabels','on'); axis([-.9 .9 -.9 .9]);

Set the temperature on the left edge to 100 degrees. Set the heat flux out of the block on the right edge to -10. The top and bottom edges and the edges inside the cavity are all insulated: there is no heat transfer across these edges.

thermalBC(thermalmodelS,'Edge',6,'Temperature',100); thermalBC(thermalmodelS,'Edge',1,'HeatFlux',-10);

Specify the thermal conductivity of the material as a simple linear function of temperature `u`

.

```
k = @(~,state) 0.7+0.003*state.u;
thermalProperties(thermalmodelS,'ThermalConductivity',k);
```

Generate a mesh.

generateMesh(thermalmodelS);

Calculate the steady-state solution.

Rnonlin = solve(thermalmodelS);

Because the thermal conductivity is nonlinear (depends on the temperature), compute the system matrices corresponding to the converged temperature. Assign the temperature distribution to the `u`

field of the `state`

structure array. Because the `u`

field must contain a row vector, transpose the temperature distribution.

state.u = Rnonlin.Temperature.';

Assemble finite element matrices using the temperature distribution at the nodal points.

`FEM = assembleFEMatrices(thermalmodelS,'nullspace',state)`

`FEM = `*struct with fields:*
Kc: [1277x1277 double]
Fc: [1277x1 double]
B: [1320x1277 double]
ud: [1320x1 double]
M: [1277x1277 double]

Compute the solution using the system matrices to verify that they yield the same temperature as `Rnonlin`

.

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

Compare this result to the solution given by `solve`

.

norm(u - Rnonlin.Temperature)

ans = 9.8994e-05

`model`

— Model object`PDEModel`

object | `ThermalModel`

object | `StructuralModel`

object | `ElectroMagneticModel`

objectModel object, specified as a
`PDEModel`

object,
`ThermalModel`

object,
`StructuralModel`

object, or
`ElectroMagneticModel`

object.

**Example: **```
model =
createpde(1)
```

**Example: **```
thermalmodel =
createpde('thermal','steadystate')
```

**Example: **```
structuralmodel =
createpde('structural','static-solid')
```

**Example: **```
emagmodel =
createpde('electromagnetic','electrostatic')
```

`bcmethod`

— Method for including boundary conditions`'none'`

(default) | `'nullspace'`

| `'stiff-spring'`

Method for including boundary conditions, specified as `'none'`

,
`'nullspace'`

, or
`'stiff-spring'`

. For more
information, see Algorithms.

**Example: **`FEM = assembleFEMatrices(model,'nullspace')`

**Data Types: **`char`

| `string`

`matrices`

— Matrices to assemblematrix identifiers |

`'boundary'`

| `'domain'`

Matrices to assemble, specified as:

Matrix identifiers, such as

`'F'`

,`'MKF'`

,`'K'`

, and so on — Assemble the corresponding matrices. Each uppercase letter represents one matrix:`K`

,`A`

,`F`

,`Q`

,`G`

,`H`

,`R`

,`M`

, and`T`

. You can combine several letters into one character vector or string, such as`'MKF'`

.`'boundary'`

— Assemble all matrices related to geometry boundaries.`'domain'`

— Assemble all domain-related matrices.

**Example: **```
FEM =
assembleFEMatrices(model,'KAF')
```

**Data Types: **`char`

| `string`

`state`

— Time for time-dependent models and solution for nonlinear modelsstructure array

Time for time-dependent models and solution for nonlinear models, specified in a structure array. The array fields represent the following values:

`state.time`

contains a nonnegative number specifying the time value for time-dependent models.`state.u`

contains a solution matrix of size*N*-by-*Np*that can be used to assemble matrices in a nonlinear problem setup, where coefficients are functions of`state.u`

. Here,*N*is the number of equations in the system, and*Np*is the number of nodes in the mesh.

**Example: **```
state.time = tlist(end);
FEM =
assembleFEMatrices(model,'boundary',state)
```

`FEM`

— Finite element matricesstructural array

Finite element matrices, returned as a structural array. Use the `bcmethod`

and `matrices`

arguments to
specify which finite element matrices you want to
assemble.

The fields in the structural array depend
on `bcmethod`

:

If the value is

`'none'`

, then the fields are`K`

,`A`

,`F`

,`Q`

,`G`

,`H`

,`R`

, and`M`

.If the value is

`'nullspace'`

, then the fields are`Kc`

,`Fc`

,`B`

,`ud`

, and`M`

.If the value is

`'stiff-spring'`

, then the fields are`Ks`

,`Fs`

, and`M`

.

The fields in the structural array also
depend on `matrices`

:

If the value is

`boundary`

, then the fields are all matrices related to geometry boundaries.If the value is

`domain`

, then the fields are all domain-related matrices.If the value is a matrix identifier or identifiers, such as

`'F'`

,`'MKF'`

,`'K'`

, and so on, then the fields are the corresponding matrices.

For more information, see Algorithms.

Partial Differential Equation Toolbox™ solves equations of the form

$$m\frac{{\partial}^{2}u}{\partial {t}^{2}}+d\frac{\partial u}{\partial t}-\nabla \xb7\left(c\otimes \nabla u\right)+au=f$$

and eigenvalue equations of the form

$$\begin{array}{l}-\nabla \xb7\left(c\otimes \nabla u\right)+au=\lambda du\\ \text{or}\\ -\nabla \xb7\left(c\otimes \nabla u\right)+au={\lambda}^{2}mu\end{array}$$

with the Dirichlet boundary conditions, **hu** = **r**, and
Neumann boundary conditions, $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g$$.

`assembleFEMatrices`

returns the following full finite element matrices and
vectors that represent the corresponding PDE problem:

`K`

is the stiffness matrix, the integral of the discretized version of the`c`

coefficient.`M`

is the mass matrix, the integral of the discretized version of the`m`

or`d`

coefficients.`M`

is nonzero for time-dependent and eigenvalue problems.`A`

is the integral of the discretized version of the`a`

coefficient.`F`

is the integral of the discretized version of the`f`

coefficient. For thermal, electromagnetic, and structural problems,`F`

is a source or body load vector.`Q`

is the integral of the discretized version of the`q`

term in a Neumann boundary condition.`G`

is the integral of the discretized version of the`g`

term in a Neumann boundary condition. For structural problems,`G`

is a boundary load vector.The

`H`

and`R`

matrices come directly from the Dirichlet conditions and the mesh.

The `'nullspace'`

technique eliminates
Dirichlet conditions from the problem using a linear algebra
approach. It generates the combined finite-element matrices
`Kc`

, `Fc`

, `B`

, and vector `ud`

corresponding to the reduced system
`Kc*u = Fc`

, where ```
Kc =
B'*(K + A + Q)*B
```

, and ```
Fc =
B'*(F + G)
```

. The `B`

matrix spans the null space of the columns of
`H`

(the Dirichlet condition
matrix representing `h*ud = r`

). The
`R`

vector represents the
Dirichlet conditions in `H*ud = R`

. The
`ud`

vector has the size of the
solution vector. Its elements are zeros everywhere except at
Dirichlet degrees-of-freedom (DoFs) locations where they
contain the prescribed values.

From the `'nullspace'`

matrices, you can
compute the solution `u`

as

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

.

If you assembled a particular set of matrices, for example
`G`

and `M`

, you
can impose the boundary conditions on `G`

and `M`

as follows. First, compute the
nullspace of columns of `H`

.

```
[B,Or] = pdenullorth(H);
ud = Or*((H*Or\R)); % Vector with known value of the constraint DoF.
```

Then use the `B`

matrix as follows. To
eliminate Dirichlet degrees of freedom from the load vector
`G`

, use:

GwithBC = B'*G

To eliminate Dirichlet degrees of freedom from mass matrix, use:

M = B'*M*B

You can eliminate Dirichlet degrees of freedom from other vectors and matrices using the same technique.

The `'stiff-spring'`

technique converts
Dirichlet boundary conditions to Neumann boundary conditions
using a stiff-spring approximation. It returns a matrix
`Ks`

and a vector
`Fs`

that together represent a
different type of combined finite element matrices. The
approximate solution is `u = Ks\Fs`

.
Compared to the `'nullspace'`

technique,
the `'stiff-spring'`

technique generates
matrices more quickly, but generally gives less accurate
solutions.

If the number of nodes in a model is
`NumNodes`

, and the number of equations is
`N`

, then the length of column
vectors `u`

and `ud`

is
`N*NumNodes`

. The toolbox assigns
the IDs to the degrees of freedom in `u`

and `ud`

:

Entries from 1 to

`NumNodes`

correspond to the first equation.Entries from

`NumNodes+1`

to`2*NumNodes`

correspond to the second equation.Entries from

`2*NumNodes+1`

to`3*NumNodes`

correspond to the third equation.

The same approach applies to all other entries, up to
`N*NumNodes`

.

For example, in a 3-D structural model, the length of a solution
vector `u`

is
`3*NumNodes`

. The first
`NumNodes`

entries correspond to
the `x`

-displacement at each node, the next
`NumNodes`

entries correspond to
the `y`

-displacement, and the next
`NumNodes`

entries correspond to
the `z`

-displacement.

In thermal analysis, the `m`

and
`a`

coefficients are zeros. The
thermal conductivity maps to the `c`

coefficient. The product of the mass density and the
specific heat maps to the `d`

coefficient.
The internal heat source maps to the `f`

coefficient. The temperature on a boundary corresponds to
the Dirichlet boundary condition term `r`

with `h = 1`

. Various forms of boundary
heat flux, such as the heat flux itself, emissivity, and
convection coefficient, map to the Neumann boundary
condition terms `q`

and
`g`

.

In structural analysis, the `a`

coefficient is
zero. The Young's modulus and Poisson's ratio map to the
`c`

coefficient. The mass density
maps to the `m`

coefficient. The body loads
map to the `f`

coefficient. Displacements,
constraints, and components of displacement along the axes,
map to the Dirichlet boundary condition terms
`h`

and `r`

.
Boundary loads, such as pressure, surface tractions, and
translational stiffnesses, correspond to the Neumann
boundary condition terms `q`

and
`g`

. When you specify the damping
model by using the Rayleigh damping parameters
`Alpha`

and
`Beta`

, the discretized damping
matrix `C`

is computed by using the mass
matrix `M`

and the stiffness matrix
`K`

as ```
C = Alpha*M +
Beta*K
```

.

In electrostatic and magnetostatic analyses, the
`m`

, `a`

, and
`d`

coefficients are zeros. The
relative permittivity and relative permeability map to the
`c`

coefficient. The charge
density and current density map to the `f`

coefficient. The voltage and magnetic potential on a
boundary correspond to the Dirichlet boundary condition term
`r`

with ```
h =
1
```

.

`ElectromagneticModel`

| `PDEModel`

| `solve`

| `solvepde`

| `StructuralModel`

| `ThermalModel`

You have a modified version of this example. Do you want to open this example with your edits?

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)