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 objectPDEModel
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 assemble'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 modelsTime 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 matricesFinite 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
and eigenvalue equations of the form
with the Dirichlet boundary conditions, hu = r, and Neumann boundary conditions, .
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.