Assemble finite element matrices and solve elliptic PDE

** assempde IS
NOT RECOMMENDED.** Use

`solvepde`

instead.

solves
the PDE`u`

= assempde(`model`

,`c`

,`a`

,`f`

)

$$-\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.$$

`[`

, for any of the previous input
syntaxes, assembles finite element matrices using the `Kc`

,`Fc`

,`B`

,`ud`

]
= assempde(___)*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 Reduced Linear System.

`[`

assembles finite element matrices
that represent any Dirichlet boundary conditions using a `Ks`

,`Fs`

]
= assempde(___)*stiff-spring* approximation.
You can calculate the solution `u`

at node points
by the command `u = Ks\Fs`

.
See Stiff-Spring Approximation.

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

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

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

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.

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.

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

Was this topic helpful?