Documentation

assempde

Assemble finite element matrices and solve elliptic PDE

assempde IS NOT RECOMMENDED. Use solvepde instead.

Syntax

  • [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

(cu)+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

(cu)+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 Reduced Linear System.

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 Stiff-Spring Approximation.

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

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

(cu)+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

(cu)+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

(cu)+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

(cu)+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 modelPDEModel 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

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifyc 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

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifya 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

(cu)+au=f,

or in the system of PDEs

(cu)+au=f.

You can specifyf 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 ith column of model.Mesh.Nodes or the ith 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 and Their Gradients.

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.

More About

collapse all

Reduced Linear System

This form of the finite element matrices eliminates Dirichlet conditions from the problem using a linear algebra approach. 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.

See Systems of PDEs for details on the approach used to eliminate Dirichlet conditions.

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 Reduced Linear System and Stiff-Spring Approximation.

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.

Introduced before R2006a

Was this topic helpful?