# assembleFEMatrices

Assemble finite element matrices

## Syntax

``FEM = assembleFEMatrices(model)``
``FEM = assembleFEMatrices(model,matrices)``
``FEM = assembleFEMatrices(model,bcmethod)``
``FEM = assembleFEMatrices(___,state)``

## Description

example

````FEM = assembleFEMatrices(model)` returns a structural array containing all finite element matrices for a PDE problem specified as a `model`.```

example

````FEM = assembleFEMatrices(model,matrices)` returns a structural array containing only the specified finite element matrices.```

example

````FEM = assembleFEMatrices(model,bcmethod)` assembles finite element matrices and imposes boundary conditions using the method specified by `bcmethod`.```

example

````FEM = assembleFEMatrices(___,state)` assembles finite element matrices using the input time or solution specified in the `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.```

## Examples

collapse all

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] ```

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] ```

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.

`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 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.

`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);`

`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 = 1.1680e-04 ```

## Input Arguments

collapse all

Model object, specified as a `PDEModel` object, `ThermalModel` object, `StructuralModel` object, or `ElectroMagneticModel` object.

`assembleFEMatrices` does not support assembling FE matrices for 3-D magnetostatic analysis models.

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

Example: ```thermalmodel = createpde("thermal","steadystate")```

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

Example: ```emagmodel = createpde("electromagnetic","electrostatic")```

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 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`

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)```

## Output Arguments

collapse all

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.

## Algorithms

collapse all

Partial Differential Equation Toolbox™ solves equations of the form

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

and eigenvalue equations of the form

`$\begin{array}{l}-\nabla ·\left(c\otimes \nabla u\right)+au=\lambda du\\ \text{or}\\ -\nabla ·\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}}·\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.

### Imposing Dirichlet Boundary Conditions

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.

Note

Internally, the toolbox uses the `"nullspace"` approach to impose Dirichlet boundary conditions while computing the solution using `solvepde` and `solve`.

### Degrees of Freedom (DoFs)

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.

### Thermal, Structural, and Electromagnetic Analysis

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. 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```. Hysteretic (structural) damping contributes to the stiffness matrix `K`, which becomes complex.

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```.

Note

Assembling FE matrices does not work for harmonic analysis and 3-D magnetostatic analysis.

## Version History

Introduced in R2016a

expand all