mlinterface_pdf.book
V E R S I O N 3 . 2
COM SO L Mu l t i p h y s i c s
S c r i p t i n g G u i d e
™
COMSOL
How to contact COMSOL:
Benelux
COMSOL BV
Röntgenlaan 19
2719 DX Zoetermeer
The Netherlands
Phone: +31 (0) 79 363 4230
Fax: +31 (0) 79 361 4212
info@femlab.nl
www.femlab.nl
Denmark
COMSOL A/S
Diplomvej 376
2800 Kgs. Lyngby
Phone: +45 88 70 82 00
Fax: +45 88 70 80 90
info@comsol.dk
www.comsol.dk
Finland
COMSOL OY
Lauttasaarentie 52
FIN-00200 Helsinki
Phone: +358 9 2510 400
Fax: +358 9 2510 4010
info@comsol.fi
www.comsol.fi
France
COMSOL France
19 rue des bergers
F-38000 Grenoble
Phone: +33 (0)4 76 46 49 01
Fax: +33 (0)4 76 46 07 42
info@comsol.fr
www.comsol.fr
Germany
FEMLAB GmbH
Berliner Str. 4
D-37073 Göttingen
Phone: +49-551-99721-0
Fax: +49-551-99721-29
info@femlab.de
www.femlab.de
Norway
COMSOL AS
Verftsgata 4
NO-7485 Trondheim
Phone: +47 73 84 24 00
Fax: +47 73 84 24 01
info@comsol.no
www.comsol.no
Sweden
COMSOL AB
Tegnérgatan 23
SE-111 40 Stockholm
Phone: +46 8 412 95 00
Fax: +46 8 412 95 10
info@comsol.se
www.comsol.se
Switzerland
FEMLAB GmbH
Technoparkstrasse 1
CH-8005 Zürich
Phone: +41 (0)44 445 2140
Fax: +41 (0)44 445 2141
info@femlab.ch
www.femlab.ch
United Kingdom
COMSOL Ltd.
Studio G8 Shepherds Building
Rockley Road
London W14 0DA
Phone:+44-(0)-20 7348 9000
Fax: +44-(0)-20 7348 9020
info.uk@comsol.com
www.uk.comsol.com
United States
COMSOL, Inc.
1 New England Executive Park
Suite 350
Burlington, MA 01803
Phone: +1-781-273-3322
Fax: +1-781-273-6603
COMSOL, Inc.
1100 Glendon Avenue, 17th Floor
Los Angeles, CA 90024
Phone: +1-310-689-7250
Fax: +1-310-689-7527
COMSOL, Inc.
744 Cowper Street
Palo Alto, CA 94301
Tel: +1-650-324-9935
Fax: +1-650-324-9936
info@comsol.com
www.comsol.com
For a complete list of international
representatives, visit
www.comsol.com/contact
Company home page
www.comsol.com
COMSOL user forums
www.comsol.com/support/forums
COMSOL Multiphysics Scripting Guide
© COPYRIGHT 1994–2005 by COMSOL AB. All rights reserved
Patent pending
The software described in this document is furnished under a license agreement. The software may be used
or copied only under the terms of the license agreement. No part of this manual may be photocopied or
reproduced in any form without prior written consent from COMSOL AB.
COMSOL, COMSOL Multiphysics, and COMSOL Script are trademarks of COMSOL AB.
Other product or brand names are trademarks or registered trademarks of their respective holders.
Version: September 2005 COMSOL 3.2
C O N T E N T S
C h a p t e r 1 : C O M S O L M u l t i p h y s i c s S c r i p t i n g
Typographical Conventions . . . . . . . . . . . . . . . . . . . 1
Introduction to Programming in COMSOL Multiphysics 3
Running COMSOL Multiphysics from COMSOL Script or MATLAB. . . . 3
A First Example—Poisson’s Equation on the Unit Disk. . . . . . . . . 4
Exporting and Importing the FEM Structure. . . . . . . . . . . . . 5
Saving and Loading a FEM Structure . . . . . . . . . . . . . . . . 6
Creating the FEM Structure Using a Model M-file. . . . . . . . . . . 7
Data Structure and Function Overview 8
FEM Structure Overview . . . . . . . . . . . . . . . . . . . . 8
Data Structures and Functions . . . . . . . . . . . . . . . . . . 9
Geometry Modeling 11
Working with Geometry Objects . . . . . . . . . . . . . . . . 11
Creating a 1D Geometry . . . . . . . . . . . . . . . . . . . . 12
Creating a 2D Geometry Using Solid Modeling . . . . . . . . . . . 12
Creating a 2D Geometry Using Boundary Modeling . . . . . . . . . 14
Creating 3D Geometries Using Solid Modeling . . . . . . . . . . . 15
Working with a Geometry Model . . . . . . . . . . . . . . . . 16
Importing and Exporting Geometry and CAD Models from File. . . . . . 16
Working with the Analyzed Geometry . . . . . . . . . . . . . . 17
The Geometry Object Hierarchy. . . . . . . . . . . . . . . . . 18
Command-Line Geometry Modeling . . . . . . . . . . . . . . . 21
Working with Image Files and MRI Data . . . . . . . . . . . . . . 24
Spline Interpolation for Creating Geometry Objects . . . . . . . . . 30
Surface Interpolation . . . . . . . . . . . . . . . . . . . . . 32
Meshing 36
Specifying a Model 38
The Geometry and the FEM Problem . . . . . . . . . . . . . . . 40
Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
C O N T E N T S | i
ii | C O N T E N T S
Shape Functions . . . . . . . . . . . . . . . . . . . . . . . 48
Variables . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Materials/Coefficients Libraries . . . . . . . . . . . . . . . . . 55
Equations and Constraints . . . . . . . . . . . . . . . . . . . 56
Discretization . . . . . . . . . . . . . . . . . . . . . . . . 68
Initial Values. . . . . . . . . . . . . . . . . . . . . . . . . 70
The Extended Mesh . . . . . . . . . . . . . . . . . . . . . . 71
Multiple Geometries . . . . . . . . . . . . . . . . . . . . . 72
The Application Structure . . . . . . . . . . . . . . . . . . . 73
Overview of Postprocessing Functions. . . . . . . . . . . . . . . 80
Overview of Common Solver Properties and Their Values . . . . . . . 82
Associative Geometry and Parametrization . . . . . . . . . . . . . 88
The Structure of the Model M-file 90
Model M-files for Models With Multiple Geometries . . . . . . . . . 97
C h a p t e r 2 : S i m u l i n k a n d S t a t e – S p a c e E x p o r t
Simulink Export 100
Dynamic or Static Export? . . . . . . . . . . . . . . . . . . 100
General or Linearized Export? . . . . . . . . . . . . . . . . . 101
The Simulink Export Dialog Box . . . . . . . . . . . . . . . . 101
The Simulink Export Types . . . . . . . . . . . . . . . . . . 102
Relation to Solver Parameters and Solver Manager Settings . . . . . . 103
The Input/Output Variables . . . . . . . . . . . . . . . . . . 104
Using the COMSOL Multiphysics Model in Simulink. . . . . . . . . 105
State-Space Export 107
The State-Space Form . . . . . . . . . . . . . . . . . . . . 107
Model Reduction . . . . . . . . . . . . . . . . . . . . . . 108
The State-Space Export Dialog Box . . . . . . . . . . . . . . . 109
The Input/Output Variables . . . . . . . . . . . . . . . . . . 110
INDEX 113
1
C O M S O L M u l t i p h y s i c s S c r i p t i n g
This chapter explains how to use the COMSOL Multiphysics programming
language and scripting capabilities to build models.
Typographical Conventions
All COMSOL manuals use a set of consistent typographical conventions that
should make it easy for you to follow the discussion, realize what you can expect to
see on the screen, and know which data you must enter into various data-entry
fields. In particular, you should be aware of these particular conventions:
• A boldface font of the shown size and style indicates that the given word(s)
appear exactly that way on the COMSOL graphical user interface (for toolbar
buttons in the corresponding tooltip). For instance, we often refer to the Model
Navigator, which is the window that appears when you start a new modeling
session in COMSOL; the corresponding window on the screen has the title
Model Navigator. As another example, the instructions might say to click the
Multiphysics button, and the boldface font indicates that you can expect to see a
button with that exact label on the COMSOL user interface.
• The names of other items on the graphical user interface that do not have direct
labels contain a leading uppercase letter. For instance, we often refer to the Draw
1
2 | C H A P T E R
toolbar; this vertical bar containing many icons appears on the left side of the user
interface during geometry modeling. However, nowhere on the screen will you see
the term “Draw” referring to this toolbar (if it were on the screen, we would print
it in this manual as the Draw menu).
• The symbol > indicates a menu item or an item in a folder in the Model Navigator.
For example, Physics>Equation System>Subdomain Settings is equivalent to: On the
Physics menu, point to Equations System and then click Subdomain Settings.
COMSOL>Heat Transfer>Conduction means: Open the COMSOL folder, open the Heat
Transfer folder, and select Conduction.
• A Code (monospace) font indicates keyboard entries in the user interface. You might
see an instruction such as “Type 1.25 in the Current density edit field.” The
monospace font also indicates COMSOL Script codes.
• An italic font indicates the introduction of important terminology. Expect to find
an explanation in the same paragraph or in the Glossary. The names of books in the
COMSOL documentation set also appear using an italic font.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
I n t r o du c t i o n t o P r o g r amm ing i n
COMSOL Mu l t i p h y s i c s
This manual describes how to model using the COMSOL Multiphysics™
programming language in the COSMOL Script™ environment or at the MATLAB®
command line. The sections cover the steps of geometry modeling, meshing the
geometry, defining the PDE problem definition, and finally computing the solution
and performing postprocessing.
The COMSOL Multiphysics Command Reference serves as reference documentation
for all scripting-related functionality in COMSOL Multiphysics. The COMSOL Script
documentation and MATLAB documentation, respectively, serves as documentation
for the programming languages.
Running COMSOL Multiphysics from COMSOL Script or MATLAB
You can always use the COMSOL Script environment, which you start from the
COMSOL Multiphysics user interface by choosing COMSOL Script from the File menu.
If you want to use the features described in this manual in MATLAB, you must run
COMSOL with MATLAB. The easiest way to do so in Windows is to double-click on
the icon COMSOL with MATLAB, which should be visible on the desktop. Under UNIX
and Linux, start COMSOL with MATLAB by using the command comsol matlab.
See the COMSOL Multiphysics Installation Guide for more information about how
to run COMSOL Multiphysics.
Note: All the following procedures refer to COMSOL Script, unless otherwise
indicated. If you run COMSOL Multiphysics with MATLAB, the commands and
procedure are identical, but the scripting environment is then MATLAB instead of
COMSOL Script.
When working on the command line, you will encounter the FEM structure. It is a
data structure (a scalar structure array in COMSOL Script and MATLAB) containing
the definition of the PDE problem. This chapter reviews how to set up the FEM
structure to solve a problem. Step-by-step instructions show how to specify a model
by assigning values to the fields in the FEM structure.
I N T R O D U C T I O N T O P R O G R A M M I N G I N C O M S O L M U L T I P H Y S I C S | 3
4 | C H A P T E R
A First Example—Poisson’s Equation on the Unit Disk
A classic PDE with a well-known behavior is Poisson’s equation
on the unit disk Ω with f = 1. The exact solution is
which makes it possible to compare the numerical solution with the exact solution at
the node points of the mesh.
To set up this example on the command line, follow these steps:
1 Take as habit to clear the FEM structure when starting a new problem:
clear fem
2 The geom field in the FEM structure contains the problem’s geometry. You create a
unit circle with the command circ2:
fem.geom=circ2(0,0,1);
3 The meshinit function creates a triangular mesh on the geometry defined in the
fem.geom field.
fem.mesh=meshinit(fem);
4 To specify the PDE coefficients, use the coefficient form of the basic equations and
set both c and f to 1. All boundaries should have u=0 as boundary conditions, which
means setting h to 1. All coefficients you do not specify are zero.
fem.equ.f=1;
fem.equ.c=1;
fem.bnd.h=1;
The default name of the dependent variables i u, you can choose another name, for
example T, by entering fem.dim=’T’;.
5 Choose quadratic Lagrange elements. The default order for the lagrange elements
is 1 on the command line. This differs from the graphical user interface where the
default order is 2.
fem.shape=2;
6 The meshextend function creates the extended mesh object, which is necessary to
assemble the problem.
∇ ∇u( )⋅– f= Ω
u 0= ∂Ω⎩
⎨
⎧
u x y,( )
1 x
2
– y
2
–
4
—————————,=
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
fem.xmesh=meshextend(fem);
7 Solve the PDE and plot the solution:
fem.sol=femlin(fem);
postplot(fem,’tridata’,’u’,’triz’,’u’);
8 Compute the maximum error by evaluating the difference to the exact solution:
pd=posteval(fem,’u-(1-x^2-y^2)/4′);
er=max(max(pd.d))
9 If the error is not sufficiently small, you can refine the mesh:
fem.mesh=meshrefine(fem);
fem.xmesh=meshextend(fem);
10 Solve the problem on the new mesh and plot its solution:
fem.sol=femlin(fem);
postplot(fem,’tridata’,’u’,’triz’,’u’);
You can recompute the error with the steps just given.
Find this example for the COMSOL Multiphysics graphical user interface in the
COMSOL Multiphysics Modeling Guide in the section “Poisson’s Equation on the
Unit Disk” on page 249.
Exporting and Importing the FEM Structure
A good way to get to know the FEM structure is by exporting it from COMSOL
Multiphysics to the command line. You can do this anytime during modeling in
COMSOL Multiphysics. Also when working at the command line, you can look at the
model in COMSOL Multiphysics by importing the FEM structure from the command
line.
E X P O R T I N G T H E F E M S T R U C T U R E
Export the current FEM structure to the command line by selecting Export>FEM
Structure as ‘fem’ from the File menu. If you also want to specify the name of the
exported FEM structure you can select Export>FEM Structure from the File menu.
Including Default Values
By default, the export does not include fields with their default values in the exported
FEM structure. To include also the default values:
1 From the Options menu, open the Preferences dialog box.
2 On the Modeling tab, select the Include default values check box in the Model M-file/
FEM export area.
I N T R O D U C T I O N T O P R O G R A M M I N G I N C O M S O L M U L T I P H Y S I C S | 5
6 | C H A P T E R
3 Click OK.
Shrinking Coefficients
By default, the export “shrinks” the coefficients in the physics settings to a more
compact notation than the full syntax (see “Equations and Constraints” on page 56 for
more information). To get the full syntax in the exported FEM structure:
1 From the Options menu, open the Preferences dialog box.
2 On the Modeling tab, clear the Shrink coefficients check box in the Model M-file/FEM
export area.
3 Click OK.
I M P O R T I N G T H E F E M S T R U C T U R E
You can set up a complete PDE problem using a FEM structure on the command line
and then import it into COMSOL Multiphysics. Select Import>FEM Structure on the
File menu. Type the name of the FEM structure in the dialog box. COMSOL
Multiphysics tries to import the FEM structure, attempting not to terminate the
import due to an syntactic error. Instead it indicates any problem in the message log
and sets up the parts of the FEM structure it can handle.
Saving and Loading a FEM Structure
COMSOL Multiphysics provides special commands for saving and loading FEM
structures to and from a COMSOL Multiphysics Model MPH-file.
S A V I N G A F E M S T R U C T U R E
Use the command flsave to save a FEM structure to a Model MPH-file. For example,
to save the FEM structure fem to the model file model.mph, type
flsave model fem
You can then open the model file in the COMSOL Multiphysics graphical user
interface or load the FEM structure using the command flload.
To save a FEM structure or any part of a FEM structure as a MAT-file, use the save
command.
L O A D I N G A F E M S T R U C T U R E
Use the command flload to load a FEM structure from a Model MPH-file or a
MAT-file.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
If filename is a Model MPH-file, the command flload filename assigns a variable
fem to be the FEM structure in the model file.
If filename is a MAT-file created with an older version of COMSOL Multiphysics
(FEMLAB), the flload command loads the FEM structures, the Simulink structures,
or the geometry objects saved under the file name filename. In this case the
difference from the standard load command in MATLAB is that flload reobjectifies
the structures for future compatibility; that is, the original objects are recreated from
the saved structures.
Creating the FEM Structure Using a Model M-file
Another way to learn command-line modeling is to save a model from the COMSOL
Multiphysics as a Model M-file. It is a complete description of a PDE problem that you
develop in the COMSOL Multiphysics user interface. To create a FEM structure in the
command prompt, run the Model M-file from the command line. Find additional
information in the section “The Structure of the Model M-file” on page 90.
I N T R O D U C T I O N T O P R O G R A M M I N G I N C O M S O L M U L T I P H Y S I C S | 7
8 | C H A P T E R
Da t a S t r u c t u r e and Fun c t i o n
Ov e r v i ew
FEM Structure Overview
The data structures that define a PDE problem are stored in a single data structure—
the FEM structure. A structure array is an array of containers for a set of fields, which
you can access with a dot notation; the FEM Structure is a scalar structure array. It is
similar to the structure in the C language. The data structures in the fields of the FEM
structure define different aspects of the PDE problem, and the various processing
stages produce new data structures from old ones. The most important fields for
command-line modeling are listed below.
TABLE 1-1: FEM STRUCTURE FIELDS
FIELD DETAILS IN INTERPRETATION
fem.sdim page 41 Names of space coordinates (independent
variables)
fem.frame page 41 Frame names
fem.draw page 16 The solid objects, face objects (3D only),
curve objects (2D and 3D only), and point
objects used for geometry modeling
fem.geom page 11 Geometry object
fem.mesh page 47 Mesh object or structure. Contains the
finite element mesh
fem.appl page 74 Application modes
fem.sshape page 41 Geometry shape
fem.dim page 56 Names or number of dependent variables
fem.shape page 50 Shape functions (finite elements)
fem.const page 53 Definition of constants
fem.cporder Constraint point orders
fem.gporder Numerical integration orders
fem.expr page 53 Definition of variables as expressions
fem.var page 77 Application scalar variables
fem.lib page 55 Material library
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Data Structures and Functions
The following table indicates how the major commands in COMSOL Multiphysics
affect the FEM structure. Most functions in COMSOL Multiphysics are designed to
accept the full FEM structure as an input, and by default, they return only a part of the
FEM structure. Most functions can optionally return the full FEM structure.
fem.equ page 38 Variables, equations, constraints, and initial
values on subdomains
fem.bnd page 38 Variables, equations, constraints, and initial
values on boundaries (not present in 0D)
fem.edg page 38 Variables, equations, constraints, and initial
values on edges (only in 3D)
fem.pnt page 38 Variables, equations, constraints, and initial
values on vertices (only in 2D and 3D)
fem.ode page 38 ODEs and other global scalar variables,
equations and initial conditions
fem.form page 56 Form of equations. Can be coefficient
(default), general, or weak
fem.solform page 71 The solution form. Can be coefficient,
general, or weak (default)
fem.border page 65 Assembly on interior boundaries
fem.xmesh page 71,
meshextend
Extended mesh structure
fem.sol femsol Solution object.
Subfields: u, ut, lambda, tlist, plist
TABLE 1-2: MAJOR COMSOL MULTIPHYSICS FUNCTIONS
FUNCTION INPUT FIELDS OUTPUT FIELDS
geomcsg draw geom
geomanalyze many fields, especially draw, and
geom
draw, geom, appl,
equ, bnd, edg, pnt
meshinit geom mesh
meshrefine geom, mesh mesh
TABLE 1-1: FEM STRUCTURE FIELDS
FIELD DETAILS IN INTERPRETATION
D A T A S T R U C T U R E A N D F U N C T I O N O V E R V I E W | 9
10 | C H A P T E R
multiphysics many fields, especially fem.appl dim, form, equ, bnd,
edg, pnt, var, elemmph,
eleminitmph, shape,
sshape, border
femdiff sdim, dim, equ, bnd equ, bnd
meshextend most fields xmesh
assemble xmesh, const Assembled matrices
asseminit xmesh, const Solution object
femlin,
femnlin,
femtime, femeig
xmesh, const, init sol
adaption geom, mesh, xmesh, const, dim,
init
mesh, xmesh, sol
posteval,
postinterp
geom, mesh, xmesh, const, sol,
equ.ind, bnd.ind, edg.ind,
pnt.ind
Post data (values of
expressions)
postplot,
postmovie
geom, mesh, xmesh, const, sol,
equ.ind, bnd.ind, edg.ind,
pnt.ind
Graphics output
TABLE 1-2: MAJOR COMSOL MULTIPHYSICS FUNCTIONS
FUNCTION INPUT FIELDS OUTPUT FIELDS
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Geome t r y Mode l i n g
The basics of geometry modeling appear in the chapter “Geometry Modeling and
CAD Tools” on page 29 of the COMSOL Multiphysics User’s Guide. This section
adds information about geometry modeling on the command line.
Working with Geometry Objects
Use Geometry objects to form the geometry of the PDE problem. For each space
dimension, there is a base class of geometry objects containing necessary information
for the geometry modeling: geom0, geom1, geom2, and geom3.
In 3D, there are solid objects (solid3), face objects (face3), curve objects (curve3),
and point objects (point3), all of which are subclasses to the geom3 class. These have
subclasses representing primitive geometry objects such as blocks, spheres, and
cylinders.
In 2D, there is one subclass to geom2 for solid objects (solid2), one subclass for
rational Bézier curve objects (curve2), and one subclass for points (point2). You
can also create geometry objects representing the primitive objects: squares,
rectangles, circles, and ellipses.
A set of similar geometry classes exists for 1D geometries. The base class geom1 has the
subclasses solid1 and point1 for intervals and points.
I M P O R T I N G G E O M E T R Y O B J E C T S
You can import geometry objects from the COMSOL Script workspace into
COMSOL Multiphysics by choosing Import>Geometry Object from the File menu.
Type the variables names of the geometry objects in the dialog box.
E X P O R T I N G G E O M E T R Y O B J E C T S
You can export geometry objects from COMSOL Multiphysics to the COMSOL
Script workspace by choosing the Export>Geometry Object from the File menu. The
names of the geometry objects in COMSOL Multiphysics also serve as the names in
the scripting workspace.
G E O M E T R Y M O D E L I N G | 11
12 | C H A P T E R
Creating a 1D Geometry
For more information about 1D geometry modeling in COMSOL Multiphysics, see
“Creating a 1D Geometry Model” on page 42 of the COMSOL Multiphysics User’s
Guide. From the COMSOL Script prompt, you create the same 1D geometry model
using the constructors solid1, point1, and the function geomcsg.
The command
s=solid1([0 1 2])
creates a 1D solid object with two subdomains.
Entering
p=point1(0.5)
creates a 1D point object, and
g=geomcsg({s},{p})
merges the two geometry objects to a 1D analyzed geometry.
Creating a 2D Geometry Using Solid Modeling
For more information about 2D geometry modeling in COMSOL Multiphysics, see
the example “2D Solid Modeling Techniques” on page 47 of the COMSOL
Multiphysics User’s Guide. You can set up the same geometry using Boolean
operations from the COMSOL Script command line.
Creating Holes
To create the holes, follow these steps:
1 Create a solid rectangle with corners at (-1, -1) and (1, 1):
s1 = rect2(-1,1,-1,1)
geomplot(s1)
The function rect2 creates a solid rectangular object s1. The function geomplot
plots the rectangle (in this case a square).
2 Create a circular hole with a radius of 0.5 and center it at (0, 0):
s2 = circ2(0,0,0.5)
The function circ2 creates a circular object s2. Use it to drill a hole in the first
object by forming the difference:
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
s3 = s1-s2
or by entering
s3 = geomcomp({s1,s2},’ns’,{‘s1′,’s2′},’sf’,’s1-s2′)
3 Visualize the result:
geomplot(s3)
When entering s3 = s1-s2, you take advantage of the object-orientated capabilities
in COMSOL Multiphysics. Both s1 and s2 are solid objects, and the operation is
equivalent to calling the function geomcomp with the parameters shown earlier. The
intersection operator * and the union operator + work similarly. The function
geomcomp analyzes the specified geometric objects and returns the analyzed
geometry.
Trimming Solids
To remove the corners from the solid rectangle follow the steps below:
1 Rotate s1 around (0, 0):
s2 = rotate(s1,45*(pi/180),0,0)
2 Form the intersection:
s = s3*s2
3 Visualize the result:
geomplot(s)
The reason for using (pi/180) as an the argument for rotate is that the function
expects the angle of rotation in radians rather than degrees.
Adding Domains
To add a domain to a solid, use the union operator (+). Take the object, the solid s,
and attach an elliptical domain on its right side.
1 Create a solid ellipse centered at (0, 0) and with semiaxes (0.5, 0.25):
e = ellip2(0,0,0.5,0.25,45*pi/180)
The final argument, 45*pi/180, expresses the angle of rotation in radians with
respect to the semiaxes and coordinate axes.
2 Rotate the solid ellipse by -45 degrees to get an ellipse with semiaxes parallel to the
coordinate axes:
e = rotate(e,-45*pi/180)
3 Move the center of the ellipse from (0, 0) to (1, 0):
G E O M E T R Y M O D E L I N G | 13
14 | C H A P T E R
e = move(e,1,0)
The second and third parameters in the function call correspond to the x and y
displacements, respectively.
4 Use the union operation to add the solid ellipse to the solid s:
s = s+e
5 Visualize the result:
geomplot(s)
Creating a 2D Geometry Using Boundary Modeling
Use command-line boundary modeling to obtain result from the example “2D
Boundary Modeling” on page 51 of the COMSOL Multiphysics User’s Guide.
1 Create the six boundary curve objects:
w = 1/sqrt(2)
c1 = curve2([-0.5 -1 -1],[-0.5 -0.5 0],[1 w 1])
c2 = curve2([-1 -1 -0.5],[0 0.5 0.5],[1 w 1])
c3 = curve2([-0.5 0.5],[0.5 0.5])
c4 = curve2([0.5 1 1],[0.5 0.5 0],[1 w 1])
c5 = curve2([1 1 0.5],[0 -0.5 -0.5],[1 w 1])
c6 = curve2([0.5 -0.5],[-0.5 -0.5])
The objects c1, c2, c3, c4, c5, and c6 are all curve2 objects. The vector [1 w 1]
specifies the weights for a rational Bézier curve that is equivalent to a quarter-circle
arc. The weights are adjusted to create elliptical or circular arcs. For more
information on the available geometry objects, see “The Geometry Object
Hierarchy” on page 18 and the individual geometry object entries in the Command
Reference.
2 Create a solid object s of the object type solid2:
s = geomcoerce(‘solid’,{c1,c2,c3,c4,c5,c6})
3 Visualize the result:
geomplot(s)
Create a single-curve object by converting a solid2 object to a curve2 object:
c = curve2(s)
You can also accomplish the same thing with a call to geomcomp:
c = geomcomp({c1,c2,c3,c4,c5,c6})
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Creating 3D Geometries Using Solid Modeling
See the example “Solid Modeling Using Solid Objects and Boolean Operations” on
page 72 of the COMSOL Multiphysics User’s Guide. This section demonstrates the
use of Boolean operations on solid objects from the COMSOL Script command line.
The examples show how to create composite solids starting from primitive solids.
Using the programming language for creating geometry objects can be useful. For
instance, when you need a sequence of commands to compute exact coordinates of
geometry objects. To create the geometry perform the following steps:
1 Create a solid rectangle with corners at (0,0) and (1,2).
r = rect2(0,1,0,2)
2 Use the fillet utility to take off its edges.
rf = fillet(r,’radii’,0.125)
3 Create a circle object with radius 0.25.
c1 = circ2(0.5,1,0.25)
4 Form a triangle using three curve objects.
t1 = curve2([0.15 0.2],[0.75 0.75])
t2 = curve2([0.2 0.2],[0.75 0.8])
t3 = curve2([0.2 0.15],[0.8 0.75])
5 Create a solid triangle by a call to geomcoerce.
t = geomcoerce(‘solid’,{t1,t2,t3});
The function geomcoerce analyzes the specified geometric objects and coerces the
resulting geometric object into the specified class.
6 Extrude the rectangle along the z-axis from 0 to 0.2 and extrude circle c1 from 0
to 0.2.
blk = extrude(rf,’distance’,0.2);
hole = extrude(c1,’distance’,0.2);
7 Create the chamfer by revolving the triangle object t 360 degrees about the center
of the hole. To do so, define a work plane going through the center of the hole, then
define an axis of rotation.
wrkpln = [[0;1;0],[0;1;1],[1;1;0]];
ax = [[0;0.5],[1;0]];
chamf = revolve(t,’angles’,[0,2*pi],…
‘revaxis’,ax,’wrkpln’,wrkpln);
chamf = solid3(chamf);
G E O M E T R Y M O D E L I N G | 15
16 | C H A P T E R
8 Use the resulting 3D objects to create the composite solid object by subtracting the
chamfer and hole from the filleted block using the Boolean operations + and -:
plate = blk-(chamf+hole);
9 Plot the final geometry using geomplot.
geomplot(plate)
10 Use use the Headlight or Scenelight toolbar buttons to turn on lighting in COMSOL
Script.
Working with a Geometry Model
A geometry model consists of a number of geometry objects. Store the geometry
model in the draw structure fem.draw in the fields s, f, c, and p for solids, faces,
curves, and points. See the Command Reference entry on geomcsg for details.
I M P O R T I N G A G E O M E T R Y M O D E L
Import a geometry model, draw, from scripting workspace to COMSOL Multiphysics
by typing
clear fem;
fem.draw=draw;
fem.sdim=2;
and then choose Import>FEM Structure on the File menu. Type the name fem in the
dialog box.
E X P O R T I N G A G E O M E T R Y M O D E L
Export a geometry model from COMSOL Multiphysics to the command line by
selecting Export>FEM Structure on the File menu. The geometry model becomes
available in fem.draw in the scripting environment.
Importing and Exporting Geometry and CAD Models from File.
With COMSOL Multiphysics, you can import and export geometries from a variety of
file formats. See the COMSOL Multiphysics Command Reference entries on
geomimport and geomexport for details. Below is a short summary of the various file
formats.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
C O M S O L M U L T I P H Y S I C S F I L E S
A natural choice for storing geometries in 1D, 2D, and 3D is the native file format of
COMSOL’s geometry kernel (.mphtxt, .mphbin). Note that this file format is only
used for geometry and mesh objects. It is not the same as a Model MPH-file (.mph).
2 D C A D F O R M A T S
COMSOL Multiphysics supports import and export for the DXF® file format, the
native format of the CAD system AutoCAD®. Furthermore, you can import files on
the neutral GDS format.
3 D C A D F O R M A T S
It is possible to import surface meshes on the STL and VRML formats. Moreover,
COMSOL provides a CAD Import Module, which provides import of most 3D CAD
file formats such as: Parasolid®, SAT®, STEP, IGES, CATIA® V4, CATIA® V5, Pro/
ENGINEER®, Autodesk Inventor®, and VDA-FS. The CAD Import Module User’s
Guide covers these file formats in detail.
Working with the Analyzed Geometry
When a geometry is used for the actual PDE model it is called the analyzed geometry.
The analyzed geometry is stored in the FEM structure field fem.geom.
In addition to a geometry object, you can also store the name of a Geometry M-file,
mesh object, or a PDE Toolbox decomposed geometry matrix in fem.geom, but we
no longer recommend. Convert Geometry M-file names, meshes, and PDE Toolbox
decomposed geometry matrices to geometry objects by using the function
geomobject before storing them in fem.geom. This ensures efficient operation.
The Geometry M-file and mesh are described in the Command Reference entries on
geomfile and femmesh, respectively.
The function geominfo provides an information interface to the analyzed geometry.
You can plot the analyzed geometry with the geomplot function.
I M P O R T I N G T H E A N A L Y Z E D G E O M E T R Y
An analyzed geometry, geom, can be imported into COMSOL Multiphysics by typing
clear fem; fem.geom=geom; and then choosing Import>FEM structure from the File
menu. Type the name fem in the dialog box.
G E O M E T R Y M O D E L I N G | 17
18 | C H A P T E R
E X P O R T I N G T H E A N A L Y Z E D G E O M E T R Y
The analyzed geometry can be exported to the command line using the Export>FEM
Structure option from the File menu. The analyzed geometry is available in fem.geom.
The Geometry Object Hierarchy
The geometric objects store information necessary for modeling. A geometric entity
refers to a vertex, edge segment, face segment, or subdomain. For an n-dimensional
space, a subdomain is an n-dimensional geometric entity bounded by
(n−1)-dimensional boundaries. The crucial information of the geometry objects is the
description of these boundaries together with information on the up and down, or left
and right, subdomain of each boundary. The main difference between 1D, 2D, and
3D geometry objects is the amount of information necessary to properly describe the
boundaries. The higher the space dimension, the more information is necessary.
• In 1D, the boundaries are vertices (points) defined by one coordinate value. To each
vertex there is an associated up and down subdomain.
• In 2D, the subdomains are bounded by rational Bézier curves of degree one, two,
and three. These curves, or edge segments, are defined by control points and
weights, as described above, and for each curve the up and down subdomain is
given. There can also be isolated vertices and curves, unrelated to the subdomains.
• In 3D, the subdomains are bounded by faces or trimmed surfaces. A face is a 2D
geometric entity bounded by curves. To properly describe a face, it is necessary to
first of all have information on the underlying surface, that is, control points and
weights of the rational Bézier patch. Furthermore, the bounding curves must
defined. These can either be rational Bézier curves with descriptions, as in the 2D
case, or curves that are intersections between several surfaces. These can only be
described in terms of the associated rational Bézier patches of the intersection. The
up and down subdomain numbering of faces works analogously to the 1D and 2D
cases. If a face has the same subdomain on the up and down side, the face is referred
to as isolated. There can also be isolated vertices and curves. An isolated vertex can
either lie on a face or inside a subdomain.
I N H E R I T A N C E S T R U C T U R E
The base class of each dimension contains the properties described above. For each
base class there are a number of subclasses with specific properties describing common
geometric objects such as rectangles, cylinders, and right-angled blocks.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
3D
The following is the hierarchical inheritance structure for the 3D geometry classes:
The four subclasses to the geom3 class describe geometric entities of dimension 0, 1,
2, and 3. The solid3 and face3 classes have similar branches of subclasses used to
represent primitive objects of surface or solid types. The figure below shows the
inheritance structure for the 3D solid object branch.
An identical structure exists for the face object branch, where the class names are
tetrahedron2 and so on. For more information on the primitive objects, see the
corresponding entries on the constructor methods in the Command Reference.
3D objects can also be created from 2D objects such as by extrusion or revolution.
solid3 face3 curve3 point3
geom3
solid3
tetrahedron3 hexahedron3
block3 gencyl3 torus3 ellipsoid3
pyramid3 econe3 sphere3
cone3
cylinder3
G E O M E T R Y M O D E L I N G | 19
20 | C H A P T E R
2D
Below is the hierarchical inheritance structure of the 2D geometry classes
The base class geom2 has three subclasses corresponding to solid objects, curve objects,
and point objects. Curve and solid subclasses describe the primitive objects rectangle
and ellipse.
An object of the class solid2 consists of a boundary and an interior part as well as
optional internal borders and subdomains. An object of the class curve2 consists only
of a boundary, which can either be open or closed. An object of the class curve2 can
be coerced to a nonempty solid object of the class solid2 if there is a closed boundary
present in the object. A solid object or a curve object can always be coerced into a point
object.
geom2
curve2 point2solid2
rect2 ellip2
square2 circ2
rect1 ellip1
square1 circ1
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
1D
Below is the hierarchical inheritance structure of the 1D geometry classes
The 1D geometry objects are either intervals (solid objects) or points. A point object
can be coerced to a nonempty solid object if it contains more than one point.
Command-Line Geometry Modeling
The 3D geometry classes are listed below.
TABLE 1-3: GEOMETRY CLASSES AND THEIR SUPERCLASSES, 3D.
FUNCTION DESCRIPTION SUPERCLASS
block2 Rectangular block face object face3
block3 Rectangular block solid object solid3
cone2 Cone face object econe2
cone3 Cone solid object econe3
curve3 3D curve object geom3
cylinder2 Cylinder face object cone2
cylinder3 Cylinder solid object cone3
econe2 Eccentric cone face object gencyl2
econe3 Eccentric cone solid object gencyl3
ellipsoid2 Ellipsoid face object face3
ellipsoid3 Ellipsoid solid object solid3
face3 3D face object geom3
geom Geometry object
geom3 3D geometry object geom
hexahedron2 Hexahedron face object face3
hexahedron3 Hexahedron solid object solid3
geom1
point1solid1
G E O M E T R Y M O D E L I N G | 21
22 | C H A P T E R
The 2D geometry classes are listed below.
point3 3D point object geom3
pyramid2 Rectangular pyramid face object gencyl2
pyramid3 Rectangular pyramid solid object gencyl3
gencyl3 Straight homogeneous generalized
cylinder solid object with linearly varying
section
solid3
gencyl2 Straight homogeneous generalized
cylinder face object with linearly varying
section
face3
solid3 3D solid object geom3
sphere2 Sphere face object face3
sphere3 Sphere solid object solid3
tetrahedron2 Tetrahedron face object face3
tetrahedron3 Tetrahedron solid object solid3
torus2 Torus face object face3
torus3 Torus solid object solid3
TABLE 1-4: GEOMETRY CLASSES AND THEIR SUPERCLASSES, 2D.
FUNCTION DESCRIPTION SUPERCLASS
circ1 Circle curve object ellip1
circ2 Circle solid object ellip2
curve2 2D rational Bézier curve object geom2
ellip1 Ellipse curve object curve2
ellip2 Ellipse solid object solid2
geom Geometry object
geom2 2D geometry object geom
point2 2D point object geom2
rect1 Rectangle curve object curve2
rect2 Rectangle solid object solid2
solid2 2D solid object geom2
square1 Square curve object rect1
square2 Square solid object rect2
TABLE 1-3: GEOMETRY CLASSES AND THEIR SUPERCLASSES, 3D.
FUNCTION DESCRIPTION SUPERCLASS
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
The 1D geometry classes are listed below.
G E O M E T R Y F U N C T I O N S
The most important geometry methods and functions are shown in the following table
TABLE 1-5: GEOMETRY CLASSES AND THEIR SUPERCLASSES, 1D.
FUNCTION DESCRIPTION SUPERCLASS
geom Geometry object
geom1 1D geometry object geom
point1 1D point object geom1
solid1 1D solid object geom1
TABLE 1-6: GEOMETRY METHODS AND FUNCTIONS
METHOD 1D 2D 3D DESCRIPTION
arc1 ♦ Create elliptical or circular arc
arc2 ♦ Create elliptical or circular solid sector
chamfer ♦ Create flattened corners in geometry
object
elevate ♦ Elevate degrees of geometry object
Bézier curves
embed ♦ Embed a 2D object in a plane
extrude ♦ Create an extruded 3D object
fillet ♦ Create circular rounded corners in
geometry object
flim2curve ♦ Create curve object from image data
geomanalyze ♦ ♦ ♦ Analyze geometry objects in draw
structure
geomarrayr ♦ ♦ ♦ Create rectangular array of geometry
objects
geomcoerce ♦ ♦ ♦ Coerce geometry objects
geomcomp ♦ ♦ ♦ Analyze (compose) geometry objects
geomcsg ♦ ♦ ♦ Combines any number of solid objects,
face objects, curve objects, and point
objects and return an analyzed
geometry
geomdel ♦ ♦ ♦ Delete interior boundaries
geomexport ♦ ♦ ♦ Export geometry object to file
geomfile ♦ ♦ Geometry M-file
G E O M E T R Y M O D E L I N G | 23
24 | C H A P T E R
Working with Image Files and MRI Data
Magnetic resonance imaging (MRI) is an imaging technique used primarily in medical
settings to produce high quality images of the inside of the human body. The more
general name for the analytical technique is NMR (Nuclear Magnetic Resonance
spectroscopy). MRI data is typically represented as a sequence of 2D images.
Note: The functionality described in this section is only available when you run
COMSOL Multiphysics with MATLAB.
MRI import consists of two steps. First the images need to be imported, and the
corresponding 2D geometry objects need to be created. After that, 3D geometry
objects can be created using a lofting technique.
geomgetwrkpln ♦ Retrieve work plane information
geomimport ♦ ♦ ♦ Import geometry object from file
geominfo ♦ ♦ ♦ Retrieve geometry information
geomplot ♦ ♦ ♦ Plot a geometry object
geomspline ♦ Spline interpolation
geomsurf ♦ Surface interpolation
line1 ♦ Create open curve polygon
line2 ♦ Create solid polygon
loft ♦ Loft 2D geometry sections to a 3D
geometry
mirror ♦ ♦ Reflect geometry
move ♦ ♦ ♦ Move geometry object
poly1 ♦ Create curve polygon
poly2 ♦ Create solid polygon
revolve ♦ Create a revolved object
rotate ♦ ♦ Rotate geometry object
scale ♦ Scale geometry object
split ♦ ♦ Split solid, curve, or point object
TABLE 1-6: GEOMETRY METHODS AND FUNCTIONS
METHOD 1D 2D 3D DESCRIPTION
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
I M P O R T I N G I M A G E S
Images of the formats JPEG, TIFF, GIF, BMP, PNG, HDF, PCX, XWD, ICO, and
CUR can be converted into 2D or 3D geometry objects, working from the MATLAB
command prompt. To read graphic file format images use the MATLAB function
imread. This function converts images on graphic file format to a matrix format in the
MATLAB work space. Other useful functions are imwrite for saving images, image
and imagesc for displaying images on matrix format, and imfinfo for obtaining
information about graphics file format images. See the MATLAB documentation for
more information on the MATLAB functions for image analysis.
Note: The Image Processing Toolbox for MATLAB contains additional functionality
for advanced image processing.
In COMSOL Multiphysics, a 2D curve object can be created from an image using the
function flim2curve, and a 3D solid object from a sequence of images using the
functions flim2curve and loft. For more information, see below in this section and
the Command Reference entries flim2curve and loft.
A 2D curve describing the contour curves of an image is created by calling
flim2curve with appropriate property values. The function flim2curve utilizes the
MATLAB function imcontour for detecting contours in images. The following
example shows an image on matrix format, generated from MATLAB’s sample
function peaks, that is converted to a 2D curve object. Instead of creating the image
in MATLAB as in this example, you can import the image from a file using the imread
command.
p = (peaks+7)*5;
figure
image(p)
v = axis;
c = flim2curve(p,{[],[5:5:75]});
g = geomcsg({rect2(5,45,0,50)},{c});
s = solid2(g);
Visualize the geometry object with the following commands.
figure
geomplot(s,’pointmode’,’off’,’sublabels’,’on’);
axis(v)
G E O M E T R Y M O D E L I N G | 25
26 | C H A P T E R
axis ij
You can continue working on this model either from the MATLAB prompt or from
the COMSOL Multiphysics graphical user interface. To insert this geometry model
into the GUI, select Import and then Geometry Objects from the File menu. Once you
have done this, you can use the Split and Union buttons for creating the geometry to
be used in the simulation.
If the original image is noisy, the resulting geometry object may contain too many edge
segments to be practical to work with. The function flim2curve provides basic
functionality for filtering data by supporting the same input argument properties as
flmesh2spline (see the Command Reference entries for these functions). Also,
MATLAB provides functionality for filtering images before converting them to
geometry objects. See the MATLAB documentation on the functions filter2 and
conv2 for more information on how to do this.
For more information on reducing noise in point data, see “Spline Interpolation for
Creating Geometry Objects” on page 30.
If the processed image is large, the computation time for the MRI import functions
can be long. To reduce the computation time it is often possible to subsample the
image before processing it, that is, reduce the number of pixels in the image. The most
simple subsampling method is to include every second pixel in the x and y direction,
respectively. More advanced subsampling methods preprocess the image by low-pass
filtering (also referred to as softening or blurring) before reducing the number of
pixels. This often gives a higher quality of the subsampled image due to the fact that
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
low-pass filtering reduces the amount of high frequency components in the
subsampled image (due to rasterization effects).
I M P O R T I N G M R I D A T A
When creating a 3D geometry object from a sequence of images, for example, MRI
data, the work flows as follows.
1 Load MRI data to MATLAB. The 3D MRI data is typically represented as a
sequence of 2D images.
2 Create a set of 2D geometry sections from the images using the flim2curve
command.
3 Use the loft command to create an interpolating 3D surface between the 2D
geometry sections.
4 Define your FEM problem on the resulting geometry.
The example below shows how to create a 3D solid object from MRI-data.
Start by loading the MRI data from MATLAB.
load mri
This loads a sequence of images stored in multidimensional array called mri. It also
loads a colormap named map, for illustrating the images.
Create an index vector indicating which sections to include.
i = [1 6 12 17 22 27];
Visualize the sections.
figure
for k = 1:6
subplot(2,3,k)
image(D(:,:,1,i(k)))
title(sprintf(‘Image %d’,k))
axis off
end
G E O M E T R Y M O D E L I N G | 27
28 | C H A P T E R
colormap(map)
Before creating curve objects from the images, the threshold and keepfrac data
need to be defined. Use the threshold value to create a binary image where all parts of
the image with a lower value than the threshold value is turned white. This binary
image is used for creating contour curve as in the function imcontour. When the
threshold value is set to 1, the contour curves are created on the boundaries of the
regions of the images that are completely black. Use the keepfrac property in
flim2curve for determining how many points in the contours created from the
images should be used for creating a curve object. You must tune this value for all
images so that every created curve object consists of the same number of segments.
This is crucial for creating a 3D geometry model of the cross-section images.
th = [1 1 1 1 1 1];
kf = [0.11 0.10 0.112 0.115 0.129 0.165];
Loop over sections and save curve-data in the cell-array c.
for k = 1:6
[c{k},r] = …
flim2curve(D(:,:,1,i(k)),{th(k),[]},’KeepFrac’,kf(k));
end
Convert 2D curve objects to 2D solid objects.
for k = 1:6
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
c{k} = solid2(c{k});
end
Now, each of these solid objects could be used for creating a 2D mesh, and for solving
2D problems. Alternatively, the objects could be given specific names and then
imported into COMSOL Multiphysics. For example, create
s1=c{1}
and import the solid object s1 into COMSOL Multiphysics using Import>Geometry
Objects from the File menu. You can then continue modeling from COMSOL
Multiphysics, creating a mesh, specifying the boundary conditions and material
parameters, solving, and visualizing the results.
Instead, in this example, the 2D solids are seen as 2D cross sections of a 3D geometry,
as follows.
Create a mapping table that maps edges between sections. See loft help text or
Command Reference entry.
el = {18 18 18 19 19 18};
This means that edge 18 of section 1 is mapped on edge 18 of section 2. The rest of
the edge mappings between sections 1 and 2 are then automatically derived from this.
The mappings between sections 2,3; 3,4; 4,5; and 5,6; are: 18-18, 18-19, 19-19 and
19-18, respectively.
Create 3D positioning data for the sections.
dvr = {repmat(12.5,1,5),repmat(0,2,6),repmat(0,1,6)};
Loft between the sections using cubic lofting (the default lofting method), with the
lofting weights explicitly given.
lg = loft(c,’loftedge’,el,’loftsecpos’,dvr,…
‘loftweights’,repmat(0.1,2,5));
Visualize the geometry.
geomplot(lg)
G E O M E T R Y M O D E L I N G | 29
30 | C H A P T E R
The geometry object is imported into COMSOL Multiphysics by choosing
Import>Geometry Objects from the File menu. You can then continue and use this
geometry as part of a simulation: create a mesh, solve, and postprocess the solution.
Note that to be able to mesh this on a computer with limited memory resources, you
need to select a coarse mesh setting.
Spline Interpolation for Creating Geometry Objects
To create a curve object from a set of data points, you can use the function
geomspline. It creates C1 or C2 continuous splines for the interpolation between the
points.
Note: The functionality described in this section is only available when you run
COMSOL Multiphysics with MATLAB.
The following example illustrates the use of geomspline by creating a geometry object
with splines created from irregularly distributed points on a circle.
First create the circle data, sorted by the angle, and remove some of the points.
phi = 0:0.2:2*pi;
phi([1 3 6 7 10 20 21 25 28 32]) = [];
p = [cos(phi);sin(phi)];
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Add some noise to the data points.
randn(‘state’,17)
p = p+0.01*randn(size(p));
Different global parameterization techniques can be used in geomspline. The global
parameterization is a parameter that varies from 0 in the first interpolated point to 1 in
the last interpolated point. For a closed curve these two points coincide. For details on
the different parameterization techniques, see the reference entry on geomspline in
the Command Reference.
First use a uniform parameterization technique. This is the most straightforward
technique, but it does not always work well on irregularly distributed data.
plot(p(1,:),p(2,:),’ro’)
c = geomspline(p,’splinemethod’,’uniform’,’closed’,’on’);
hold on
geomplot(c,’pointmode’,’off’)
Then interpolate using centripetal parameterization and note the difference.
c = geomspline(p,’splinemethod’,’centripetal’,’closed’,’on’);
hold on
geomplot(c,’pointmode’,’off’,’edgecolor’,’b’)
For creating a solid geometry from this, on which you can base an analysis, enter
s = solid2(c);
To use this geometry for analysis in COMSOL Multiphysics, select Import>Geometry
Objects from the File menu. Do this when the draw mode is active, either when setting
−1 −0.5 0 0.5 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
G E O M E T R Y M O D E L I N G | 31
32 | C H A P T E R
up a 2D problem or when a work plane has been defined in a 3D model. Once this is
done, you can use any of the available geometry operations to continue working on
the model.
Occasionally, point data can be too dense or contain too much noise to be practical to
work with. The functions flcontour2curve and flmesh2spline contain basic
functionality to reduce the number of edge segments in a point set. To understand the
MATLAB contour data format, see the MATLAB functions imcontour, contour, and
contourc.
Note that you can also use geomspline for creating 3D curve objects. Using the data
created in the previous example, define the z coordinates to create the following 3D
curve object.
p = [p;linspace(0,3,22)];
c = geomspline(p,’splinemethod’,’centripetal’,’closed’,’off’)
Surface Interpolation
You can create a 3D face object or a 2D solid object from a 3D point set or a 2D point
set, respectively, using the function geomsurf. See the Command Reference entry on
geomsurf for more information.
2 D S O L I D O B J E C T S
using geomsurf, you can create a solid object with rectangular subdomains from a
rectangular point set. In the following example, a solid object with non uniformly
distributed subdomains is generated from a nonuniform grid.
[x,y] = meshgrid(0:0.1:1,0:0.2:1);
rand(‘state’,1);
x = x+repmat((0.1*rand(1,size(x,2))),size(x,1),1);
rand(‘state’,1);
y= y+repmat((0.2*rand(size(y,1),1)),1,size(y,2));
s = geomsurf(x,y);
figure
geomplot(s,’pointmode’,’off’)
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
This can be used for defining a domain with a nonuniform distribution of some
material parameter.
You can import this geometry object into COMSOL Multiphysics by choosing
Import>Geometry Objects from the File menu.
3 D F A C E O B J E C T S
By using the function geomsurf, you can create a surface from height data defined on
a grid by bilinear interpolation. The surface created is identical in shape to the surface
created with the command surf. The following example illustrates how to generate a
surface from simulated or measured data and how to insert the generated surface as
the top surface of a solid block.
1 The height data is defined as deviations, with a maximum deviation of 0.03 from a
reference plane situated at z equal to 0.7.
dev_max = 0.03;
ref_level = 0.7;
[x,y] = meshgrid(0:0.1:1,0:0.1:1);
randn(‘state’,1);
z = ref_level+dev_max*randn(size(x));
2 Create a surface from the height data:
G E O M E T R Y M O D E L I N G | 33
34 | C H A P T E R
f = geomsurf(x,y,z);
3 Create a solid block whose top surface lies above the generated surface.
b = block3(1,1,ref_level+3*dev_max);
4 Coerce the face and the solid block to a composite solid object.
g = geomcoerce(‘solid’,{f,b});
5 Split the object into its subdomains.
ss = split(g);
6 Visualize the solid object with the interpolated surface as top surface.
figure
geomplot(ss{1})
7 Use use the Headlight or Scenelight toolbar buttons to turn on lighting in COMSOL
Script.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
8 You can import this geometry object into COMSOL Multiphysics by choosing
Import>Geometry Objects from the File menu.
If the original data set is noisy, the resulting geometry object may contain too many
faces to be practical to work with. MATLAB provides functionality for filtering data
sets before converting them to geometry objects. See the MATLAB documentation
on the functions filter2 and conv2 for more information on how to do this.
G E O M E T R Y M O D E L I N G | 35
36 | C H A P T E R
Me s h i n g
The mesh of a PDE problem is stored in the mesh object or (if mesh cases are used) the
mesh structure fem.mesh. In the mesh, subdomains are divided into elements, and
boundaries are broken up into boundary elements.
The mesh is created from the analyzed geometry by the functions meshinit, meshmap,
meshextrude, and meshrevolve and can be altered by the functions meshrefine and
meshsmooth.
The adaption function creates mesh data as part of the solution process. You can plot
the mesh with the meshplot function. The Command Reference entry femmesh
contains details on the mesh data representation.
I M P O R T I N G A M E S H O B J E C T
You can choose between importing a mesh object into COMSOL Multiphysics along
with a compatible analyzed geometry or without a geometry. To import a mesh along
with an associated analyzed geometry, geom, type
clear fem
fem.geom=geom;
fem.mesh=mesh;
otherwise type
clear fem
fem.mesh=mesh;
and choose Import>FEM Structure from the File menu.
Using a Mesh to Represent the Geometry
If you do not have a geometry, you can use a mesh to represent the geometry. To use
a mesh as a geometry, type the following commands at the command line:
clear fem
fem.geom=geomobject(mesh);
fem.mesh=mesh;
and choose Import>FEM Structure from the File menu. An alternative is to use the mesh
as a geometry without importing the actual mesh. To do so, enter
clear fem
fem.geom=geomobject(mesh);
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
and choose Import>FEM Structure from the File menu. You can then create a new mesh
using the geometry represented by this mesh.
E X P O R T I N G A M E S H O B J E C T
The mesh object can be exported to the command line using the Export>FEM Structure
option from the File menu. The mesh object is available in fem.mesh.
I M P O R T I N G M E S H D A T A C R E A T E D B Y A N E X T E R N A L P R O G R A M
You can create a mesh object from mesh data generated by an external program by
using the commands femmesh and meshenrich. Use the command femmesh to create
an initial mesh object from the external mesh data and use meshenrich to enrich the
initial mesh object with necessary information such that a valid mesh object is created.
Example
Consider the text files coord.txt and tet.txt (found in the multiphysics/
directory) containing the mesh vertex coordinates and element vertex indices,
respectively, of a tetrahedral mesh created with an external program. To import this
data into a mesh object, follow these steps:
1 Load the external mesh data to the command line.
load coord.txt
load tet.txt
2 Create an initial mesh object from the external mesh data.
el = cell(1,0);
tet = tet+1; % Lowest mesh point index is zero in tet.txt
el{1} = struct(‘type’,’tet’,’elem’,tet’);
m = femmesh(coord,el);
3 Enrich the initialized mesh object with necessary information such that a valid mesh
object is formed.
m = meshenrich(m);
meshplot(m);
M E S H I N G | 37
38 | C H A P T E R
S p e c i f y i n g a Mode l
In this section we describe how to specify a model using the scripting environment,
either interactively or by editing an M-file. All data describing the model is contained
in a structure which is called the FEM structure and is denoted fem in our examples.
Specify a model by assigning values to the fields of the FEM structure. We start with
the case of a single geometry. If there are several geometries, you must use an extended
FEM structure (see page 72).
V A R I A B L E S A N D F I N I T E E L E M E N T S
The names of the dependent variables are specified in the field fem.dim. The finite
elements used can be given in the field fem.shape. Global scalar dependent variables,
available everywhere in the model, and corresponding equations can be added in the
field fem.ode. Additional variables can be defined in terms of other variables in
fem.expr. Constants are specified in fem.const.
E Q U A T I O N S , C O N S T R A I N T S , A N D I N I T I A L V A L U E S
The equations, constraints (boundary conditions), and initial conditions of the PDE
problem are stored in the following fields:
• fem.equ—contains information about equations on subdomains
• fem.bnd—contains information belonging to boundaries
• fem.edg—in 3D, this field is related to edges (curves)
• fem.pnt—in 2D and 3D, this field corresponds to vertices (points)
• fem.ode—specifies ODEs and other scalar equations independent of position
Equations can be given both in strong form (as a PDE) or in weak form (as an equality
of integrals). There are two formulations of the strong form, the coefficient form,
which is suitable for linear problems, and the general form, which is suitable for both
linear and nonlinear problems. The form is specified in the field fem.form.
T H E E X T E N D E D M E S H
After the PDE problem is specified, the COMSOL Multiphysics function meshextend
creates an extended mesh structure fem.xmesh. The extended mesh contains a
low-level description of the PDE problem. Remember to run meshextend when you
have made changes to the model, because the solution and postprocessing stages are
based on the extended mesh.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
S O L U T I O N
The solution of the PDE problem is computed by one of the following solvers:
• femlin
• femnlin
• femtime
• femeig
• adaption
The solution is stored in the solution object fem.sol. The solution object is described
in the reference entry femsol. The matrix fem.sol.u contains values for the degrees
of freedom. For femtime, the time derivative of the solution is available in the matrix
fem.sol.ut, but the accuracy of the time derivative will be good only if the time
derivative have been stored using the property Outcomp (see femsolver). For
femtime and femeig, fem.sol has the additional fields fem.sol.tlist and
fem.sol.lambda, respectively. In these cases each column in fem.sol.u is a solution
object corresponding to a time value in fem.sol.tlist or an eigenvalue in
fem.sol.lambda. When femlin or femnlin has been used as a parametric solver, the
solution object has the additional field fem.sol.plist. In this case fem.sol.u has
one column for each parameter value fem.sol.plist. Also, the parameter name is in
fem.sol.pname.
PO S T P R O C E S S I N G A N D P R E S E N T A T I O N
Given a solution/mesh pair, a variety of tools is provided for the visualization and
processing of the data:
• posteval and postinterp—evaluate expressions in selected points
• postint—integrate expressions
• postplot and postcrossplot—visualize expressions and variables
• postarrow, postarrowbnd, postcont, postflow, postiso, postlin,
postslice, postsurf, and posttet—shorthand commands for different
postplot plot types
• postmovie—animate a plot
• postanim—shorthand function for standard animations.
A P P L I C A T I O N M O D E S
Using the application modes is an alternative method for entering equations in the
FEM structure. The application modes contain predefined equations for a variety of
physics applications. The data is stored in the field fem.appl and the function
S P E C I F Y I N G A M O D E L | 39
40 | C H A P T E R
multiphysics converts the data to equations in the FEM structure. The application
modes are described in the COMSOL Multiphysics Modeling Guide.
The Geometry and the FEM Problem
The geometry is given in the field fem.geom. See the Command Reference entry for
geomcsg and “Geometry Modeling” on page 11 for details. For instance,
fem.geom = rect2(-2,2,-2,2) + circ2;
gives a plane geometry consisting of the unit circle inside a square with side length 4.
The domains make up the bounded region of space in which you solve the equations.
The subdomains are numbered 1, 2, 3,…, and so are the boundaries, edges, and
vertices (if they exist). In the regions of space that are outside the domains, no physics
is allowed. These regions are sometimes given the number 0.
Data pertaining to subdomains, boundaries, edges, or vertices are stored in separate
fields of the FEM structure, as shown in the following table.
(The field name equ has historic reasons: originally it contained only equations.)
D O M A I N G R O U P S
It is often convenient to treat several domains as one group using the same equations,
material properties, or boundary conditions. You can do this by forming domain
groups. The subdomain groups are defined in the optional field fem.equ.ind. This
can be a cell array of numeric vectors or a numeric vector. Two examples illustrate these
syntaxes:
fem.equ.ind = { [1 2] [6 3 4] };
This means that subdomain group 1 consists of subdomains 1 and 2, while subdomain
group 2 consists of subdomains 6, 3, and 4. Equivalently:
fem.equ.ind = [1 1 2 2 0 2];
TABLE 1-7: FEM STRUCTURE FIELDS
DOMAIN NAME FIELD IN FEM STRUCTURE
subdomain fem.equ
boundary fem.bnd
edge fem.edg
vertex fem.pnt
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
where for each subdomain, its group number has been specified. 0 means that the
subdomain (number 5 in this case) does not belong to any group.
If the field fem.equ.ind is not given, then each subdomain forms its own subdomain
group, that is, the subdomain groups are the same as the subdomains.
Similarly, groupings of boundaries, edges, and vertices can be defined in
fem.bnd.ind, fem.edg.ind, and fem.pnt.ind, respectively.
F E M . S D I M — S P A C E C O O R D I N A T E S
The default names of the (global) space coordinates are x, y, z in 3D, x and y in 2D,
and x in 1D. You can override this by specifying your own names in the optional cell
array fem.sdim, for example
fem.sdim = { ‘r’ ‘phi’ ‘z’ };
in 3D. The space coordinates are sometimes called the independent variables.
Problems with a moving mesh require a second set of coordinates to describe the
motion of the mesh. The coordinate sets are called frames. Specify two frames by
letting fem.sdim be a cell array of two cell arrays of coordinates, for example,
fem.sdim = {{‘X’ ‘Y’} {‘x’ ‘y’}};
in 2D. One the frames is the reference frame, the fixed frame where you draw the
geometry, and the other frame describes the mesh motion. The field fem.sshape (see
below) specifies which of the two is the reference frame. For postprocessing to use the
deformed mesh in its evaluations, make sure the reference frame becomes the first of
the two frames in fem.sdim.
F E M . F R A M E — F R A M E N A M E S
When there are more than one frame you must give each of them a name. Do this in
the field fem.frame, for example,
fem.frame = {‘ref’ ‘ale’};
If you do not specify fem.frame the frames get a default name which is the
concatenation of the coordinates.
F E M . S S H A P E — G E O M E T R Y S H A P E
The global coordinates are polynomials in the local element coordinates of a certain
degree k. This degree k can be specified in the field fem.sshape. For instance,
fem.sshape = 2;
S P E C I F Y I N G A M O D E L | 41
42 | C H A P T E R
uses quadratic shape functions for the global space coordinates. This makes it possible
for the mesh elements at the boundary to be curved, and thus come closer to the true
geometric boundary. The default k is equal to the maximum order of the shape
function objects in fem.shape. k is called the geometry shape order.
Moving Mesh
When solving problems with a moving mesh, it is not sufficient to only specify the
geometry shape order. fem.sshape should then contain definitions of the geometry
shape due to mesh motion. In this case fem.sshape is a cell array of structures each
defining a geometry shape. fem.sshape should contain all such shape definitions
which appear anywhere in the model. The field fem.equ.sshape will specify which
geometry shape definitions that should be used in each domain. The structures in
fem.sshape have the following fields:
• fem.sshape.frame. The frame which this structure defines. If absent the first
frame in fem.sdim will be used.
• fem.sshape.sorder. The geometry shape order.
• fem.sshape.dvolname. The name of the mesh element scale factor.
• fem.sshape.type = ‘fixed’ | ‘moving_abs’ | ‘moving_rel’ |
‘moving_expr’ The default is ‘fixed’ which means that the frame does not move.
moving_abs is used when the mesh equation describes the absolute mesh
movement. moving_rel is used when the mesh equation describes the relative mesh
movement. moving_expr is used when the mesh movement is given as an
expression. For the reference frame the type should be ‘fixed’.
• fem.sshape.sdimdofs. A cell array of strings with the variables defining the spatial
coordinates for relative mesh movement. Required and only used when
fem.sshape.type = ‘moving_rel’. The spatial coordinates in sdimdofs have to
be defined directly by a shape function.
• fem.sshape.refframe. The name of the reference frame. Required and only used
when fem.sshape.type = ‘moving_rel’ and indicates which frame the motion
is relative to.
• fem.sshape.sdimexprs. A cell array of strings with the expressions giving the
mesh displacement. Used when fem.sshape.type = ‘moving_expr’ and gives
the expressions for the mesh deformation. The default is a cell array with zeros.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
D O M A I N S O F U S A G E F O R G E O M E T R Y S H A P E S
In the fields fem.equ.sshape, fem.bnd.sshape, fem.edg.sshape, and
fem.pnt.sshape, you can specify where the geometry shape objects in fem.sshape
are to be used. For example,
sshape1.frame = ‘ref’;
sshape1.type = ‘fixed’;
sshape1.sorder = 2;
sshape1.dvolname=’dvol’;
sshape2.frame = ‘ale’;
sshape2.type = ‘moving_abs’;
sshape2.sorder = 2;
sshape2.dvolname = ‘dvol_ale’;
sshape3.frame = ‘ale’;
sshape3.type = ‘moving_rel’;
sshape3.sorder = 2;
sshape3.sdimdofs = {‘dx’ ‘dy’};
sshape3.refframe = ‘ref’;
sshape3.dvolname = ‘dvol_ale’;
fem.sshape = {sshape1 sshape2 sshape3};
fem.equ.sshape = { [1 3] [1 2] };
means that on the first subdomain group, only the first and third geometry shape
objects (sshape1 and sshape2) are active. On the second subdomain group, the first
and third geometry shape objects are active.
To make sure that the spatial coordinates are well-defined everywhere all frames have
to be defined by one sshape object in each domain. In this example sshape1 defines
the coordinates in the reference frame while sshape2 and sshape3 are two definitions
of the coordinates in the frame ‘ale’. Therefore sshape1 should be active everywhere
and one of sshape2 and sshape3 should be defined everywhere.
If the field fem.equ.sshape is not given, then all geometry shape objects in
fem.sshape apply in all subdomain groups.
Similarly, the field fem.bnd.sshape is a cell array, which specifies for each boundary
group which shape function objects are active. If fem.bnd.sshape is not given, then
it is inherited from fem.equ.sshape. This means that a geometry shape object which
is active in a subdomain, is also active on the boundary of that subdomain (as well as
boundaries lying within the subdomain). If more than one geometry shape object
defining the coordinates of the same frame on the subdomains adjacent to a boundary
then only is selected. Anyone with type = ‘moving_abs’ is chosen first.
S P E C I F Y I N G A M O D E L | 43
44 | C H A P T E R
In 3D, the field fem.edg.sshape similarly specifies the usage of geometry shape
objects on edge groups. If it is not given, it is inherited from the usage on subdomains
and boundaries.
In 2D and 3D, the field fem.pnt.sshape similarly specifies usage on vertices. It is
defaulted by inheritance from subdomains, boundaries, and edges.
D O M A I N S O F U S A G E F O R G E O M E T R Y S H A P E S B Y E A C H E Q U A T I O N
When there are more than one frame in the problem you must specify on which frame
each equation is formulated. Do this using the field fem.sshapedim. Say you have a
problem with four equations. Then you can, for example, have
fem.sshape = {sshape1 sshape2 sshape3};
fem.equ.sshapedim = {[1 1 3 3] [1 1 2 2]};
The objects in fem.sshape are same ones as in the previous section. The first two
equations live on the reference frame using sshape1 in both subdomain groups. The
third and fourth equations live on the frame ‘ale’ using sshape3 in the first
subdomain group and sshape2 in the second group.
G E O M E T R I C V A R I A B L E S
Besides the space coordinates, there are a number of geometric variables that can be
used in expressions. Which variables that are available depends on the dimension of the
domain on which the expression is to be evaluated. In the following tables, x, y, and z
denote the names of the space coordinates (that is, they should be replaced with the
real names optionally given in fem.sdim).
TABLE 1-8: GENERAL GEOMETRIC VARIABLES
VARIABLES MEANING
x y z The space coordinates. Of course only the first n of these are
available in n dimensions.
dom The domain number
sd Alias for dom on subdomains
dvol The scale factor from local to global coordinates. For domains
of dimension 3, this is the factor that volumes are multiplied by
when going from local to global coordinates (the modulus of
the Jacobi determinant). For domains of dimension 2, it is the
area scaling factor. For domains of dimension 1, it is the length
scaling factor. For domains of dimension 0, it equals 1.
h The local mesh size, that is, the diameter of the mesh element
(in the domain) that the evaluation point resides in. This is not
available in 0D.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
On boundaries, there is a notion of up and down normal directions. In 3D, the up
normal is defined as the cross product of the tangent vectors t1 and t2 (which point
in the directions of the surface parameters s1 and s2, respectively). In 2D, the up
normal points to the left as you traverse the curve in the direction of increasing curve
parameter s. In 1D, the up normal points to the right, that is, in the direction of
increasing space coordinate. The subdomains which lay in the up and down directions
are called the upper and lower subdomains, respectively
Note: Note that the normal vector un defined below is the outward normal seen
from the upper subdomain. That is, un points in the down direction.
TABLE 1-9: ADDITIONAL GEOMETRIC VARIABLES ON BOUNDARIES IN 3D
VARIABLES MEANING
xg yg zg The space coordinates computed using the geometry’s
representation of the surface (instead of the approximation
with polynomial shape functions)
s1 s2 Surface parameters in the Bézier surface representation
t1x t1y t1z
t2x t2y t2z
Components of the unit tangent vectors to the surface in the
directions of the parameters s1 and s2, respectively
unx uny unz Components of the outward unit normal vector seen from the
upper subdomain. That is, the downward pointing normal.
dnx dny dnz Components of the outward unit normal vector seen from the
lower subdomain. That is, the upward pointing normal.
nx ny nz For boundaries that have exactly one adjacent subdomain,
these are the components of the outward unit normal vector,
that is pointing out from this subdomain. For other
boundaries, these components are the same as dnx dny dnz.
TABLE 1-10: ADDITIONAL GEOMETRIC VARIABLES ON EDGES IN 3D
VARIABLES MEANING
xg yg zg The space coordinates computed using the geometry’s
representation of the curve (instead of the approximation with
polynomial shape functions)
S P E C I F Y I N G A M O D E L | 45
46 | C H A P T E R
s1 The curve arc length parameter. This runs from 0 to the length
of the curve.
t1x t1y t1z Components of the unit tangent vector to the curve in the
direction of increasing parameter s1
TABLE 1-11: ADDITIONAL GEOMETRIC VARIABLES ON BOUNDARIES IN 2D
VARIABLES MEANING
xg yg The space coordinates computed using the geometry’s
representation of the curve (instead of the approximation with
polynomial shape functions)
s The curve parameter. For geometry objects, this parameter
varies from 0 to 1 in the direction of the curve. For mesh as
geometry and geometry M-files, this is not guaranteed.
s1 Alias for s
tx ty Components of the unit tangent vector to the curve in the
direction of increasing parameter s
t1x t1y Alias for tx, ty
unx uny Components of the outward unit normal vector seen from the
upper subdomain. That is, the downward pointing normal.
dnx dny Components of the outward unit normal vector seen from the
lower subdomain. That is, the upward pointing normal.
nx ny For boundaries that have exactly one adjacent subdomain,
these are the components of the outward unit normal vector,
that is pointing out from this subdomain. For other
boundaries, these components are the same as dnx dny.
TABLE 1-12: ADDITIONAL GEOMETRIC VARIABLES ON BOUNDARIES IN 1D
VARIABLE MEANING
unx Outward normal seen from the upper subdomain, that is, -1
dnx Outward normal seen from the lower subdomain, that is, 1
nx Outward unit normal. For boundaries that have exactly one
adjacent subdomain, this is the outward unit normal, that is
pointing out from this subdomain. For other boundaries, nx is
the same as dnx.
TABLE 1-10: ADDITIONAL GEOMETRIC VARIABLES ON EDGES IN 3D
VARIABLES MEANING
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Mesh
For mesh terminology, see “Meshing” on page 261 in the COMSOL Multiphysics
User’s Guide. Given a FEM structure with a geometry fem.geom, you can create an
initial unstructured mesh by
fem.mesh = meshinit(fem);
You can fine tune the meshing algorithm by specifying a lot of options. The most
important is the overall maximum mesh size hmax. For example,
fem.mesh = meshinit(fem, ‘hmax’, 0.01);
ensures that all mesh element diameters are less than 0.01. You can also specify local
values of Hmax if you want the mesh to be finer on certain domains.
Given an unstructured mesh, you can refine it using meshrefine:
fem.mesh = meshrefine(fem);
Here you can also choose among different algorithms.
A 1D or 2D mesh can also be used as geometry data, for example you can create a
coarser mesh starting from a finer mesh:
fem.geom = fem.mesh;
fem.mesh = meshinit(fem, ‘hmax’, 0.2);
In 1D and 2D, the geometry that is used by the meshinit function can also be a
Geometry M-file, for example the cardioid function:
fem.geom = ‘cardg’;
fem.mesh = meshinit(fem);
A mesh can be visualized using the meshplot function:
meshplot(fem);
F E M . M E S H — T H E M E S H O B J E C T O R S T R U C T U R E
The mesh object is described in detail in the Command Reference entry for femmesh.
If mesh cases are introduced, fem.mesh is a structure containing the following fields:
fem.mesh.default is the mesh corresponding to mesh case 0. This mesh is also used
in all mesh cases that does not occur in fem.mesh.mind.
fem.mesh.case is a cell array of meshes corresponding to the mesh cases in
fem.mesh.mind. The default an empty cell array.
S P E C I F Y I N G A M O D E L | 47
48 | C H A P T E R
fem.mesh.mind is a cell array containing vectors of positive mesh case indices. That is,
the mesh fem.mesh.case{i} is used in the mesh cases fem.mesh.mind{i}. The
default fem.mesh.mind is {1 2 … n}, where n is the length of fem.mesh.case.
Shape Functions
Finite elements are defined using shape function objects. Their task is to express field
variables as linear combinations of certain basis functions (shape functions). The
coefficients are called degrees of freedom (DOF), and they form the discrete
representation of the field variables. The shape function objects are objects from
certain shape function classes. COMSOL Multiphysics has the following shape
functions classes:
• shlag—Lagrange element of arbitrary order.
• shherm—Hermite element of arbitrary order (>2).
• sharg_2_5—Fifth order Argyris element on triangular meshes in 2D.
• shcurl—Vector (edge) element on triangular or tetrahedral meshes.
• shdiv—Divergence element on triangular or tetrahedral meshes.
• shbub—Bubble element on 1D, triangular or tetrahedral meshes.
• shdisc—Discontinuous element.
For more information about these elements, see “Finite Elements” on page 543 in
User’s Guide and the entries shlag etc. in the Command Reference. There are also a
number of shape function classes in the Structural Mechanics Module for beams,
plates, and shells.
S H L A G — T H E L A G R A N G E E L E M E N T
To construct a Lagrange shape function object, type
s = shlag(k, ‘u’);
This creates a Lagrange element of order k for the variable u. For instance,
s = shlag(1, ‘phi’);
means that the variable phi will have linear elements.
S H H E R M — T H E H E R M I T E E L E M E N T
To construct a Hermite shape function object, type
s = shherm(sdim, k, ‘u’);
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
where sdim is the space dimension. This creates a Hermite element of order k for the
variable u. The order k must be at least 3. For instance,
s = shherm(1, 3, ‘v’);
means that the variable v will represented as piecewise third-degree polynomials on a
1D domain, with continuous derivative.
S H A R G _ 2 _ 5 — T H E A R G Y R I S E L E M E N T
To construct an Argyris shape function object, type
s = sharg_2_5(‘u’);
This creates an Argyris element of order 5 for the variable u in 2D.
S H C U R L — T H E VE C T O R E L E M E N T ( F I R S T S Y N T A X )
To construct a vector shape function object, type
s = shcurl(order,’u’);
This means that the vector variable u will have vector (edge) elements.
S H C U R L — T H E VE C T O R E L E M E N T ( S E C O N D S Y N T A X )
An alternative syntax of the constructor is
s = shcurl(order,cell array with component names);
For instance
s = shcurl(order, {‘a’ ‘b’ ‘c’} );
in 3D, or
s = shcurl(order, {‘Hx’ ‘Hy’} );
in 2D.
S H D I V — T H E D I V E R G E N C E E L E M E N T
To construct a divergence element shape function object, type
s = shdiv(order,’u’);
S H B U B — T H E B U B B L E E L E M E N T
To construct an bubble element shape function object, type
s = shbub(mdim,’u’);
S P E C I F Y I N G A M O D E L | 49
50 | C H A P T E R
S H D I S C — T H E D I S C O N T I N U O U S E L E M E N T
To construct a discontinues element shape function object, type
s = shdisc(mdim,order,’u’);
F E M . S H A P E — D E F I N I T I O N O F S H A P E F U N C T I O N S
The field fem.shape is a cell array with shape function objects. For example,
fem.shape = {shlag(1,’u’) shlag(2,’u’) shcurl(2,’A’)};
defines three shape function objects. You can choose on which domains these objects
will be active, by using the fields fem.equ.shape, fem.bnd.shape, fem.edg.shape,
fem.pnt.shape, see below.
Alternatively, fem.shape can be a numeric vector. In this case its length has to match
the dim variables (see page 56). Lagrange shape function objects are then constructed
for the dim variables. For example,
fem.dim = {‘u’ ‘v’ ‘p’};
fem.shape = [2 2 1];
is equivalent to
fem.dim = {‘u’ ‘v’ ‘p’};
fem.shape = {shlag(2,’u’) shlag(2,’v’) shlag(1,’p’)};
A further alternative is to give a single number k in fem.shape, which amounts to
having Lagrange elements of order k for all dim variables.
If mesh cases are used, the field fem.shape can be a structure containing the following
fields:
• fem.shape.default, which is the shape functions for mesh case 0. This has the
same syntax as described above. These shape functions are also used in all mesh cases
that does not occur in fem.shape.mind.
• fem.shape.case, which is a cell array of shape function definitions. One entry
fem.shape.case{i} has the same syntax as described above, and it applies to the
mesh cases fem.shape.mind{i}. The default is an empty cell array
fem.shape.case.
• fem.shape.mind{i}, which is a cell array of vectors of positive mesh case numbers.
The default is fem.shape.mind = {1 2 … n}, where n is the length of
fem.shape.case.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Finally, if fem.shape is not given, it is taken to be 1, that is, all dim variables will have
linear elements. However, in the COMSOL Multiphysics user interface the default
Lagrange element order is 2 in most application modes.
D O M A I N S O F U S A G E F O R S H A P E F U N C T I O N S
In the fields fem.equ.shape, fem.bnd.shape, fem.edg.shape, and
fem.pnt.shape, you can specify where the shape function objects in fem.shape are
to be used. For example,
fem.shape = {shlag(1,’u’) shlag(2,’u’) sharg_2_5(‘v’)};
fem.equ.shape = { [1 3] [2 3] [] 3 };
means that on the first subdomain group, only the first and third shape function
objects (shlag(1,’u’) and sharg_2_5(‘v’)) are active. On the second subdomain
group, the second and third shape function objects are active. On the third subdomain
group, no shape function objects are defined, while on the fourth subdomain group
only sharg_2_5(‘v’) is active. Thus, the variable u will be defined on subdomain
groups 1 and 2, having linear elements in subdomain group 1, and quadratic elements
in subdomain group 2. If these subdomain groups are adjacent, this will cause
problems because COMSOL Multiphysics does not take care of the “hanging nodes”
that appear. Thus you should not mix elements for the same variable in adjacent
subdomains.
If the field fem.equ.shape is not given, then all shape function objects in fem.shape
apply in all subdomain groups.
Similarly, the field fem.bnd.shape is a cell array, which specifies for each boundary
group which shape function objects are active. If fem.bnd.shape is not given, then it
is inherited from fem.equ.shape. This means that a shape function object which is
active in a subdomain, is also active on the boundary of that subdomain (as well as
boundaries lying within the subdomain).
In 3D, the field fem.edg.shape similarly specifies the usage of shape function objects
on edge groups. If it is not given, it is inherited from the usage on subdomains and
boundaries. That is, a shape function object which is active on some subdomain or
some boundary, is also active on all edges that touch this subdomain or boundary.
In 2D and 3D, the field fem.pnt.shape similarly specifies usage on vertices. It is
defaulted by inheritance from subdomains, boundaries, and edges.
S P E C I F Y I N G A M O D E L | 51
52 | C H A P T E R
F R A M E F O R S H A P E F U N C T I O N S
For problems with moving meshes and more than one frame, you have to indicate on
which frame each shape function lives. To do this specify an additional property frame
for the shape function. When adding this property you have to use the shape function’s
property value syntax. For example, if the shape function shlag(2,’u’) should live
on the frame ale, then the shape function becomes
shlag(‘order’,2,’basename’,’u’,’frame’,’ale’).
Variables
The variables that can be used in expressions are divided into the following categories:
• Geometric variables. See page 44.
• Field variables. These are the dependent variables that describe the physical
quantities that you model. The field variables can be functions of the space
coordinates and time. In the FEM discretization, they are expressed in terms of the
degrees of freedom.
• Constants. Besides the predefined constants (for instance pi, i, and j), you can use
your own. See below.
• Special variables. These are the phase angle phase, the time t for a time-dependent
problem, and the eigenvalue lambda for an eigenvalue problem. The variable phase
is only used in postprocessing; it is put equal to 0 in the solvers and the assembly
(unless you override it with fem.const). The eigenvalue lambda cannot be used in
equation coefficients; it can only be used in postprocessing.
• Weak form meta variables. These are formed by appending _test to a field variable
name. They can only be used in weak terms. See “The Weak Form” on page 65.
The field variables are further divided into the following groups:
• Shape function variables. These are the variables that are defined by the shape
function objects; see “Shape Functions” on page 48. The dim variables are a subset
of these; see page 56.
• Expression variables. These are variables that are defined as expressions in terms of
other variables; see below.
• Application mode variables and application scalar variables. These are like the
expression variables, but they are defined by an application mode. See “The
Application Structure” on page 73 and the application mode descriptions in the
Modeling Guide.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
• Material variables. These are like the expression variables, but they are defined by a
material library (see below).
• Boundary coupled shape variables. These are certain derivatives that a priori have no
well-defined values on the boundaries. The values of some of these variables are
formed as the average of the values in the adjacent subdomains. See below.
• Coupling variables. These are variables that are defined in one domain in terms of
other variables in possibly distant domains. This makes it possible to have nonlocal
dependencies in your model.
• Equation variables. These are shorthands for certain terms in equations in
coefficient and general form. Normally they should only be used in postprocessing.
See below.
• Boundary coupled equation variables. See below.
F E M . C O N S T — C O N S T A N T S
In the field fem.const you can define your own constants. You do this using a cell
array with alternating constant names (strings) and numeric values. For example,
fem.const = { ‘c0’ 3e8 ‘mu0’ 4*pi*1e-7 };
This gives c0 the value and mu0 the value .
An alternative syntax is to let fem.const be a structure with the fields corresponding
to the names of the constants.
fem.const.c0 = 3e8;
fem.const.mu0 = 4*pi*1e-7;
is equivalent to the cell array above.
E X P R E S S I O N V A R I A B L E S
You can define new field variables in terms of others in the fields fem.expr,
fem.equ.expr, fem.bnd.expr, fem.edg.expr, and fem.pnt.expr. This can be
convenient if your equations contain the same expression several times.
In the field fem.expr you put variable definitions that apply on all domain groups (of
all dimensions). fem.expr is a cell array of alternating variable names and expressions
or a structure where the fields are the variable names.
In fem.equ.expr you put variable definitions that are in force on subdomains. For
example,
fem.equ.expr = { ‘W’ ‘B*H/2’ ‘div’ ‘ux+vy’ };
3 10
8
⋅ 4π 10
7–
⋅
S P E C I F Y I N G A M O D E L | 53
54 | C H A P T E R
defines W=B*H/2, and div=ux+vy on all subdomain groups. The defining expressions
in the fem.equ.expr can be 1-dimensional cell arrays, for instance
fem.equ.expr = { ‘v’ {‘a+b’ ‘a-b’} ‘w’ ‘pi*x*b’ };
This defines v to be a+b on subdomain group 1, and a-b on subdomain group 2. The
use of an empty vector [] instead of an expression means that the variable is not
defined on the corresponding subdomain group. For example,
fem.equ.expr = { ‘v’ {‘a+b’ []} };
means that v = a+b on subdomain group 1, and that v is undefined on subdomain
group 2.
An alternative syntax is to let fem.equ.expr be a structure where the fields are the
variable names. For example,
fem.equ.expr.v = {‘a+b’ ‘a-b’};
fem.equ.expr.w = ‘pi*x*b’;
is equivalent to
fem.equ.expr = { ‘v’ {‘a+b’ ‘a-b’} ‘w’ ‘pi*x*b’ };
Similarly, variable definitions on boundaries, edges, and vertices are put in
fem.bnd.expr, fem.edg.expr, and fem.pnt.expr, respectively.
B O U N D A R Y C O U P L E D S H A P E V A R I A B L E S
COMSOL Multiphysics can extrapolate the gradients of the dependent variables to the
boundary of the domain making the boundary coupled shape variables available. You
can turn on the generation of these variables by using the property Cplbndsh to
meshextend:
fem.xmesh = meshextend(fem, ‘cplbndsh’,’on’);
This set of variables is documented in the COMSOL Multiphysics User’s Guide in
section “Boundary Coupled Shape Variables” on page 170.
E Q U A T I O N V A R I A B L E S
When you use solution form coefficient or solution form general, COMSOL
Multiphysics can make variables representing certain terms in the coefficient form
equation available — the equation variables. The field xfem.solform in the extended
FEM structure or the property Solform to the meshextend function controls the
solution form. In addition, the meshextend property Eqvars controls the generation
of the equation variables. It is set to on by default, but the variables will still not be
available unless the solution form is coefficient or general. This set of variables is
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
documented in the COMSOL Multiphysics User’s Guide in section “Equation
Variables” on page 172.
B O U N D A R Y C O U P L E D E Q U A T I O N V A R I A B L E S
When you use solution form coefficient or solution form general, COMSOL
Multiphysics can extrapolate the equation variables to the boundary — the boundary
coupled equation variables. Control the solution form by the field xfem.solform in
the extended FEM structure or the property Solform in the meshextend function. In
addition, the meshextend property Cplbndeq controls the generation of the equation
variables. It is set to on by default, but the variables will still not be available unless the
solution form is coefficient or general. This set of variables is documented in the
COMSOL Multiphysics User’s Guide in section “Boundary Coupled Equation
Variables” on page 173. The boundary coupled equation variables should be used with
caution in equations, in order to avoid circular definitions.
Materials/Coefficients Libraries
F E M . L I B — L I B R A R Y D A T A
The field fem.lib is a structure that you can use to define a materials/coefficients
library, which defines material parameters, cross sections used in structural mechanics,
and coordinate systems. The fields in fem.lib correspond either to materials, cross
sections or a coordinate system.
For example,
fem.lib.mat{1}.name = ‘Water’;
fem.lib.mat{1}.varname = ‘mat1’;
fem.lib.mat{1}.variables.C = ‘4200’;
fem.lib.mat{1}.variables.k = ‘k(T)’;
fem.lib.mat{1}.functions{1}.type = ‘inline’;
fem.lib.mat{1}.functions{1}.name = ‘k(T)’;
fem.lib.mat{1}.functions{1}.expr = ‘0.0015*T+0.1689’;
fem.lib.mat{1}.functions{1}.dexpr = {‘diff(0.0015*T+0.1689,T)’};
defines the material mat1. This creates the material variables mat1_C, which is
defined by the corresponding expression, and the local material function mat1_k(T).
In the graphical user interface you can read and write the material library from and to
a file (see “Using the Materials/Coefficients Library” on page 206 in the COMSOL
Multiphysics User’s Guide).
Another example is
S P E C I F Y I N G A M O D E L | 55
56 | C H A P T E R
fem.lib.sec{1}.name = ‘HEA 100’;
fem.lib.sec{1}.varname = ‘sec1’;
fem.lib.sec{1}.variables.A = ‘21.2*1E-4’;
fem.lib.sec{1}.variables.Iyy = ‘349.0*1E-8’;
fem.lib.sec{1}.variables.Izz = ‘134.0*1E-8’;
fem.lib.sec{1}.variables.J = ‘5.2647e-008’;
fem.lib.sec{1}.variables.heighty = ‘100*1E-3′;
fem.lib.sec{1}.variables.heightz = ’96*1E-3’;
which defines a cross section sec1.
Coordinate systems are defined as
fem.lib.coord1.type = ‘coordsys’;
fem.lib.coord1.name = ‘Coordnate system 1’;
fem.lib.coord1.T = {‘cos(pi/90)’ ‘-sin(pi/90)’;
‘sin(pi/90)’ ‘cos(pi/90)’};
The field T defines the transformation matrix of the coordinate system. This is an
n-by-n matrix, where n is the number of space dimensions, transforming the
coordinates in the local coordinate system to the global coordinate system. The field
name coord1 is arbitrary, which kind of object is defined is indicated by the type field.
Equations and Constraints
This section describes the syntax for equations and constraints.
F E M . F O R M — F O R M O F E Q U A T I O N S
The field fem.form can be one of the strings ‘coefficient’, ‘general’, or ‘weak’.
The default is ‘coefficient’. Weak equations can be used in all three forms; they are
added to the weak reformulation of the coefficient or general formulation. If you use
fem.form = ‘weak’, contributions from the coefficient and general syntaxes are
ignored.
F E M . D I M — N A M E S O F V A R I A B L E S I N E Q U A T I O N S
The names for the dependent variables u1, u2,…, uN that occur in the coefficient and
general forms can be given in the field fem.dim. For example,
fem.dim = {‘u’ ‘v’ ‘p’};
for a system with N=3. These are called the dependent variables, and they are subset
of the shape function variables (see “Variables” on page 52). In the case of a single
variable, you can skip the braces:
fem.dim = ‘v’;
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Alternatively, fem.dim can be a number:
fem.dim = 3;
which gives the variable names u1, u2, u3. However,
fem.dim = 1;
gives the variable name u, which is the default.
You can override the variables in fem.dim on subdomains, boundaries, edges, and
points by using the fields fem.equ.dim, fem.bnd.dim, fem.egd.dim, and
fem.pnt.dim, respectively.
In a FEM structure that you export from the COMSOL Multiphysics user interface,
you find the dependent variables in fem.equ.dim instead of fem.dim.
T H E C O E F F I C I E N T F O R M
Recall that a time-dependent problem in coefficient form can be written
where k and l range from 1 to N, and m ranges from 1 to M, where M <= N. The PDE
and boundary coefficients in these equations are stored in the subfields da, c, al, ga,
be, a, and f of fem.equ, and the subfields q, g, h, and r of fem.bnd. In these fields
the coefficients are represented as nested cell arrays. Coefficients that are not specified
are assumed to be 0.
Let us give some examples. For a system with N=2,
fem.equ.f = { {'x' 'y+1'} {'x*y' 3} };
means that on subdomain group 1, you have f1=x and f2=y+1, while on subdomain
group 2, you have f1=xy and f2=3. Thus the outer pair of braces corresponds to the
different subdomain groups, and the inner pairs of braces correspond to the equation
index l.
The coefficient a is a matrix with components alk, so it can be entered as
fem.equ.a = { {'x' 1; 'y' 0} };
da lk t∂
∂uk ∇+ clk uk∇– αlkuk– γl+( )⋅ βlk ∇uk alkuk+⋅+ fl= in Ω
n clk uk αlkuk γl–+∇( )⋅ qlkuk = gl hmlµm–+ on Ω∂
hmkuk rm= on Ω∂⎩
⎪
⎪
⎨
⎪
⎪
⎧
S P E C I F Y I N G A M O D E L | 57
58 | C H A P T E R
(if N=2) using the syntax for matrices. In this case the outermost pair of curly braces
only contains one entry. In such cases, that entry applies to all subdomain groups.
Thus,
on all subdomain groups.
Because the coefficient βlk is a vector, you need a third level of braces to represent β.
For example (in case N=2 and the space dimension is 3):
fem.equ.be = { { {1 0 0} {2 4 5}; {'x' 3 4} {1 2 3} } };
means that
on all subdomain groups, or, put in matrix form,
After these examples, let us describe the general syntax of these nested cell arrays. The
syntax has some shortcuts, but except for these the levels of braces correspond to the
following:
• Level 1: The different domains or domain groups.
• Level 2: The indices k and l, which enumerate equations or variables.
• Level 3: The indices i and j, which enumerate space coordinates.
Level 1 Coefficients
Let us call the content of any of the coefficient fields da, c, al, ga, be, a, f, q, g, h, r
a level 1 coefficient. A level 1 coefficient can be a cell array or something else. The
latter case is reduced to the former case by including the coefficient in braces. Thus,
fem.equ.c = 3;
is equivalent to
a x 1
y 0
=
β11
1
0
0
β12
2
4
5
β21
x
3
4
β22
1
2
3
,=,=,=,=
β
1 0 0 2 4 5
x 3 4 1 2 3
=
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
fem.equ.c = { 3 };
The entries in a level 1 coefficient cell array are called level 2 coefficients. If the cell
array has length 1, then its single entry is valid in all subdomain/boundary groups. If
the cell array has length >1, then its nth entry applies in the nth subdomain/boundary
group.
Level 2 Coefficients
A level 2 coefficient can be a cell array or something else. The second case is interpreted
as a cell array with one entry. Thus you need to consider only the case of a cell array.
The entries in this cell array are called level 3 coefficients. The interpretation of the cell
array depends on its dimensions, according to the following table:.
The case of a c coefficient given as a vector with N(N+1)/2 components is interpreted
as symmetric matrix in the sense that ckl = clk
T, where the T denotes transpose.
COEFFICIENT NUMBER OF CELL ARRAY
ELEMENTS AT LEVEL 2
MATRIX/VECTOR INTERPRETATION
da c α β a q h 1 diagonal N-by-N matrix with
the same element throughout
the diagonal
da c α β a q N diagonal N-by-N matrix with
the given elements on the
diagonal
da c α β a q N(N + 1)/2 symmetric (see below) N-by-N
matrix, given by listing the
columns in the upper triangular
part in a vector
da c α β a q N
2 full N-by-N matrix given as cell
array
γ f g 1 vector of length N with the
same element in all components
γ f g N vector with the given
components
h MN full M-by-N matrix given as cell
array
r M vector with the given
components
S P E C I F Y I N G A M O D E L | 59
60 | C H A P T E R
Level 3 Coefficients
The meaning of level 3 coefficients is shown in the following table:.
At level 3, the coefficients a, γ, and β are vectors of length n (the space dimension).
Hence they are represented as cell arrays. The entries in these cell arrays are called level
4 coefficients (see below).
The coefficients da, a, f, q, g, h, and r are scalars at level 3. Hence they are represented
directly as level 4 coefficients. (You can optionally surround these level 4 coefficients
with a pair of braces.)
Finally, the c coefficient is a scalar or matrix at level 3. If it is a scalar it can be
represented directly as a level 4 coefficient. If it is a matrix it is represented as a cell
array of level 4 coefficients, according to the following table (where n is the space
dimension).
FIELD NAME INTERPRETATION OF LEVEL 3 COEFFICIENT
fem.equ.da dalk (scalar)
fem.equ.c clk (scalar or matrix)
fem.equ.al αlk (vector)
fem.equ.ga γl (vector)
fem.equ.be βlk (vector)
fem.equ.a alk (scalar)
fem.equ.f fl (scalar)
fem.bnd.q qlk (scalar)
fem.bnd.g gl (scalar)
fem.bnd.h hmk (scalar)
fem.bnd.r rm (scalar)
NUMBER OF CELL ARRAY ELEMENTS MATRIX INTERPRETATION OF LEVEL 3 COEFFICIENT
FOR FEM.EQU.C
1 diagonal n-by-n matrix with the same element
throughout diagonal
n diagonal n-by-n matrix with the given elements
on the diagonal
n(n + 1)/2 symmetric n-by-n matrix, given by listing the
columns of the upper triangular part as a vector
n2 full n-by-n matrix, given as cell array
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Level 4 Coefficients
A level 4 coefficient represents the scalar quantity dalk, clkij, βlkj, γlj, αlkj, alk, fl, qlk,
gl, hmk, or rm. It can be given in one of the following forms:
• A numeric scalar (possibly complex), like 3.14.
• A string containing a COMSOL Multiphysics expression. You can use variables (see
“Variables” on page 52) and functions, including your own M-files in this
expression. The functions must be vectorized, that is, take numeric arrays as inputs,
and output a numeric array of the same size as the inputs. Functions with string
arguments are not supported.
T H E G E N E R A L F O R M
The quantities da, Γ, F, G, and R in the general form are specified in the fields
fem.equ.da, fem.equ.ga, fem.equ.f, fem.bnd.g, and fem.bnd.r, respectively.
The same syntax is used as in coefficient form. You must also specify the partial
derivatives c, α, β, a, f, q, and h that occur in the linearized problem. See “The Linear
or Linearized Model” on page 321 in the COMSOL Multiphysics User’s Guide. These
are most conveniently calculated by using the function femdiff:
fem = femdiff(fem);
but you can also specify them by hand. If you use your own M-files in expressions, you
may need to specify differentiation rules for them; see the Command Reference entry
on femdiff.
E X A M P L E : N A V I E R ’ S E Q U A T I O N S
We illustrate the syntax for equations in coefficient and general form with the Navier’s
equations for structural mechanics. The plane stress case of the Navier equation for
structural mechanics can be written as
where
x∂
∂
2G+µ( )
x∂
∂u
⎝ ⎠
⎛ ⎞–
y∂
∂
G
y∂
∂u
⎝ ⎠
⎛ ⎞–
x∂
∂ µ
y∂
∂v
⎝ ⎠
⎛ ⎞–
y∂
∂ G
x∂
∂v
⎝ ⎠
⎛ ⎞– = Fx in Ω
x∂
∂ G
y∂
∂u
⎝ ⎠
⎛ ⎞–
y∂
∂ µ
x∂
∂u
⎝ ⎠
⎛ ⎞–
x∂
∂
G
x∂
∂v
⎝ ⎠
⎛ ⎞–
y∂
∂ 2G+µ( )
y∂
∂v
⎝ ⎠
⎛ ⎞– = Fy in Ω
⎩
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎧
S P E C I F Y I N G A M O D E L | 61
62 | C H A P T E R
E is Young’s modulus and ν is the Poisson ratio. Using the coefficient form you can
identify the c and f coefficients.
The following commands set up the data structures.
clear fem
E = 1e2; nu = 0.3; G = E/2/(1+nu); mu = 2*G*nu/(1-nu);
Fx = 0; Fy = -1;
fem.dim = {‘u’ ‘v’};
fem.equ.c = { { {2*G+mu 0; 0 G} {0 mu; G 0}; …
{0 G; mu 0} {G 0; 0 2*G+mu} } };
fem.equ.f = { {Fx Fy} };
There is just one subdomain, so the length of the outermost cell array is one. The level
4 coefficients are just numeric values. For the c coefficient, alternatively use the
diagonal syntax at the space-coordinate level, and the symmetry syntax at level 2:
fem.equ.c = { { {2*G+mu G} {0 mu; G 0} {G 2*G+mu} } };
You can solve a full structural mechanics problem with geometry, mesh, and boundary
conditions by typing the lines below. The boundary conditions fix the left side of the
rectangle.
fem.bnd.h = { 0 0 0 1 };
fem.geom = rect2(0, 1, 0, 0.1);
fem.mesh = meshinit(fem);
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postplot(fem, ‘tridata’, ‘v’)
The postplot command plots the displacement v in the y-direction. Because there is
no specification of fem.shape, linear elements are used.
G =
E
2 1 ν+( )
———————
µ = 2G
ν
1 ν–
————
⎩
⎪
⎪
⎨
⎪
⎪
⎧
c
2G µ+ 0
0 G
0 µ
G 0
0 G
µ 0
G 0
0 2G µ+
= f
Fx
Fy
=
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
You can also use the general form when identifying coefficients:
The following commands set up the data structures.
clear fem
E = 1e2; nu = 0.3; G = E/2/(1+nu);
Fx = 0; Fy = -1;
fem.const = {‘G’ G ‘mu’ 2*G*nu/(1-nu)};
fem.dim = {‘u’ ‘v’};
fem.form = ‘general’;
fem.equ.ga = { { {‘-(2*G+mu)*ux-mu*vy’ ‘-G*uy-G*vx’} …
{‘-G*uy-G*vx’ ‘-mu*ux-(2*G+mu)*vy’} } };
fem.equ.f={{Fx Fy}};
The symbols ux, uy, vx, and vy refer to the x and y derivatives of u and v.
There is just one subdomain, so the length of the outermost cell array is one. Because
this is a 2D problem for two variables, you have two cell array entries at both level 2
and 3. The above example uses the string expression syntax for the level 4 coefficients.
You can now solve the same problem but using quadratic elements:
fem.bnd.r = { 0 0 0 {‘u’ ‘v’} };
fem.shape = 2;
fem.geom = rect2(0, 1, 0, 0.1);
fem = femdiff(fem);
fem.mesh = meshinit(fem);
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem,’v’)
For linear problems, you can normally choose the form which suits you best. In some
cases it is easier to write down a problem using coefficient form and in other cases the
general form is more convenient.
Γ
2G µ+( )–
x∂
∂u
µ
y∂
∂v
–
G
y∂
∂u
– G
x∂
∂v
–
G
y∂
∂u
– G
x∂
∂v
–
µ
x∂
∂u
– 2G+µ( )
y∂
∂v
–
= F
Fx
Fy
=
S P E C I F Y I N G A M O D E L | 63
64 | C H A P T E R
D E A C T I V A T E D E Q U A T I O N S A N D B O U N D A R Y C O N D I T I O N S
If a dim variable is not defined on the whole geometry (see “Domains of Usage for
Shape Functions” on page 51), then the corresponding equations and boundary
conditions (in general and coefficient form) are deactivated on the subdomains and
boundaries where the variable does not exist. More precisely, assume that dim variable
number k is not defined on a domain. Then the kth equation or Neumann boundary
condition is ignored on that domain. Also, the coefficients qlk and hmk are put equal
to 0 on the domain.
Example with a Slit
This example shows how to implement a slit by defining different variables on different
domains. Assume you want to solve Laplace’s equation on a unit disk with a slit on the
negative x-axis. For boundary conditions, use u = atan2(y,x) on the circular boundary
and a linearly varying value inside the slit.
To model the slit, use two equations and two dependent variables that are active only
for y≥0 and y≤0, respectively.
clear fem
fem.draw.s.objs = {circ2(0,0,1)};
fem.draw.c.objs = {curve2([-1 1],[0 0])};
fem.draw.p.objs = {point2(0,0)};
fem.geom = geomcsg(fem);
geomplot(fem, ‘edgelabels’, ‘on’, ‘sublabels’, ‘on’)
fem.dim = 2;
fem.equ.shape = {1 2};
fem.equ.c = 1;
This means that the variable u1 is defined on subdomain 1 (y≥0), and the variable u2
is defined on subdomain 2 (y≤0).
The boundary conditions are more tricky. Only the boundary conditions on
boundaries adjacent to active subdomains will have any effect. Assume you would like
to use u1 = atan2(y,x) on boundaries 4 and 6, and u2 = atan2(y-ε,x) on boundaries 3
and 5 to ensure a negative value also in the point (-1,0). On boundary 1, set u1 = −πx
and u2 = πx, and on boundary 2 set u1 = u2:
fem.bnd.h = {1 1 {1 -1} 1};
fem.bnd.r = { {‘atan2(y,x)’ 0} {0 ‘atan2(y-eps,x)’} 0 …
{‘-pi*x’ ‘pi*x’} };
fem.bnd.ind = [4 3 2 1 2 1];
Introduce an auxiliary variable U to be u1 on subdomain 1, and u2 on subdomain 2.
Then activate the assembly on interior boundaries, create the mesh and extended
mesh, solve, and plot the solution:
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
fem.equ.expr = {‘U’ {‘u1’ ‘u2’}};
fem.border = ‘on’;
fem.mesh = meshinit(fem,’hmax’,0.05);
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem,’U’)
F E M . B O R D E R — A C T I V A T E B O U N D A R Y C O N D I T I O N S O N I N T E R I O R
B O U N D A R I E S
The field fem.border determines whether boundary conditions will be used on
interior boundaries (in coefficient or general form). fem.border can be one of the
strings ‘on’ or ‘off’, or a cell array of such strings (of length matching the dim
variables). Alternatively, fem.border can be 1 or 0, or a numeric vector of ones and
zeros (of length matching the dim variables). In this case 1 is interpreted as ‘on’ and
0 as ‘off’. If fem.border{k}=’off’, then the kth Neumann boundary condition
will be ignored on interior boundaries for the kth dim variable. Furthermore, the
coefficients qlk and hmk will be put equal to 0 on these interior boundaries. If
fem.border is not a cell array or vector it applies to all dim variables. The default is
fem.border = ‘off’.
T H E W E A K F O R M
For a description of the weak problem formulation, see the chapter “The Weak Form”
in the COMSOL Multiphysics Modeling Guide. Equations in weak form are stored in
the fields fem.equ.weak, fem.bnd.weak, fem.edg.weak, and fem.pnt.weak.
The field fem.equ.weak contains the integrand in the integral over subdomains. The
field fem.equ.weak can be a string expression or a cell array. In the first case the
expression applies to all subdomain groups. In the second case, fem.equ.weak{k}
applies to the kth subdomain group. The entry fem.equ.weak{k} can be an
expression or a cell array of expressions. In the latter case, the expressions in the cell
array are added.
In the expressions representing the integrand, the test function corresponding to a
variable v is denoted v_test. Then v_test has the same shape functions as v.
Similarly, the fields fem.bnd.weak, fem.edg.weak, and fem.pnt.weak contain
integrands that are integrated over boundaries, edges, and points, respectively. All
these integrals are put on the right-hand side of the weak equation.
Time-Dependent Problems
For a time-dependent problem, the terms containing time derivatives are stored in
fem.equ.dweak, fem.bnd.dweak, fem.edg.dweak, and fem.pnt.dweak. These have
S P E C I F Y I N G A M O D E L | 65
66 | C H A P T E R
the same syntax as the weak fields, except that the time derivatives must enter linearly.
The time derivative of a variable v is denoted v_time, but you can also use the syntax
vt. The integrals defined by the dweak fields are put on the left-hand side of the weak
equation. You can also use vtt for the second time derivative.
Eigenvalue Problems
Eigenvalue problems are specified like time-dependent problems. COMSOL
Multiphysics then interprets vt as −λv, where λ is the eigenvalue.
Constraints in the Weak Problem Formulation
The constraints are stored in the fields fem.equ.constr, fem.bnd.constr,
fem.edg.constr, and fem.pnt.constr. These constraints are implemented point
wise, as described in “Discretization of the Equations” on page 559 in the COMSOL
Multiphysics User’s Guide. The Jacobians of these constraints will be calculated
correctly. In contrast, when using the general form, the Jacobian of R only accounts
for derivatives with respect to the dim variables (and not with respect to their
derivatives, for example).
The constraints on subdomains are given in the field fem.equ.constr. This field can
be an expression or a cell array. In the first case, the expression fem.equ.constr is
constrained to be zero on all subdomain groups. In the second case,
fem.equ.constr{k} applies to the kth subdomain group. The entry
fem.equ.constr{k} can be an expression or a (possibly empty) cell array of
expressions. These expressions are constrained to be zero on the kth subdomain
group.
Similarly, constraints on boundaries, edges, and vertices are defined in
fem.bnd.constr, fem.edg.constr, and fem.pnt.constr, respectively.
Example: Poisson’s Equation with a Point Source
As an example of a problem that needs some weak contribution, consider
in the unit disk Ω, where δ is Dirac’s δ-function at the origin. Multiplying with a test
function v and integrating by parts give:
x
2
2
∂
∂ u
–
y
2
2
∂
∂ u
– δ=
0 v 0 0,( )
x∂
∂v
x∂
∂u
y∂
∂v
y∂
∂u
+⎝ ⎠
⎛ ⎞ Ad
Ω
∫–=
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
using the boundary conditions u=0 and v=0. (By introducing a Lagrange multiplier
you can remove the boundary condition on v.) To solve the problem with COMSOL
Multiphysics, type the following commands:
clear fem
fem.geom = geomcoerce(‘solid’, {circ2 point2(0,0)});
geomplot(fem, ‘pointlabels’, ‘on’)
The second command above introduces a vertex at the origin. The plot shows that the
vertex number of the origin is 3.
fem.mesh = meshinit(fem);
fem.form = ‘weak’;
fem.pnt.weak = ‘u_test’;
fem.pnt.ind = {3};
fem.equ.weak = ‘-ux_test*ux-uy_test*uy’;
fem.bnd.constr = ‘u’;
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem,’u’,’triz’,’u’)
The exact solution
is infinite at the origin. Of course, the finite element approximation cannot represent
this singularity.
Alternatively, you can use the coefficient form, with a weak contribution only for the
point source:
clear fem
fem.geom = geomcoerce(‘solid’, {circ2 point2(0,0)});
fem.mesh = meshinit(fem);
fem.pnt.weak = ‘u_test’;
fem.pnt.ind = {3};
fem.equ.c = 1;
fem.bnd.h = 1;
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem,’u’,’triz’,’u’)
F E M . O D E — G L O B A L V A R I A B L E S A N D E Q U A T I O N S
Use the field fem.ode to introduce single named degrees of freedom which are
available as variables everywhere in the model and to set equations for these variables.
Despite the suggestive name, fem.ode does not have to specify an ordinary differential
equation. The variables introduced are available everywhere to be used in domain
u
1
4π
—— x
2
y
2
+( )log–=
S P E C I F Y I N G A M O D E L | 67
68 | C H A P T E R
equations, algebraic equations or ODEs. It is also possible to set equations for variables
introduced somewhere else.
Dependent variables are listed in the fem.ode.dim field, which is a cell array of strings.
Initial values can be given both for the variables—in the fem.ode.init field—and
their time derivatives—in the fem.ode.dinit field.
Scalar equations are specified either in fem.ode.f or in fem.ode.weak. The former,
if present, is a cell array of strings with the same number of entries as fem.ode.dim.
Each entry is set equal to zero and added to the global set of equations, using the test
function of the variable at the corresponding position in the dim field. The weak field
can contain an arbitrary number of entries on weak form. Expressions used in the
equations may reference ODE variables, constants and integration coupling variables
with global destination.
Note that the introduction of degrees of freedom is very much independent from the
specification of equations. The equation fields can be left blank if the degrees of
freedom are only needed in a domain equation. Conversely, the weak field can be used
on its own to set equations for integration coupling variables with global destination.
For example, to solve Volterra’s prey-predator, you can enter the following:
clear fem;
fem.geom = geom0(zeros(0,1));
fem.mesh = meshinit(fem);
fem.const = {‘a’,1,’b’,0.1,’c’,0.1,’d’,0.001};
fem.ode.dim = {‘r’,’f’};
fem.ode.f = {‘rt-r*(a-b*f)’,’ft-f*(d*r-c)’};
fem.ode.init = {‘120′,’8′};
fem.xmesh = meshextend(fem);
fem.sol = femtime(fem,’tlist’,0:0.1:100);
postcrossplot(fem,0,1,’pointxdata’,’r’,’pointdata’,’f’)
Discretization
The matrices K, L, M, N, and D occurring in the discretized version of the problem
can be computed by
[K, L, M, N, D] = assemble(fem, ‘u’, u0);
where u0 is the solution object for which the linearization is done. See “The Linear or
Linearized Model” on page 321 in the COMSOL Multiphysics User’s Guide. See the
Command Reference entry femsol for information about the solution object.
The following data structures control the discretization.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
G P O R D E R — O R D E R O F Q U A D R A T U R E F O R M U L A
The integrals occurring in the assembly of the matrices are computed numerically
using a quadrature formula. The order of this quadrature formula is specified in the
fields fem.equ.gporder, fem.bnd.gporder, fem.edg.gporder, and
fem.pnt.gporder.
The field fem.equ.gporder gives the order for integrals over subdomains. The field
fem.equ.gporder can be a number or a cell array. In the first case fem.equ.gporder
applies to all subdomain groups. In the second case, fem.equ.gporder{i} applies to
the ith subdomain group. The entry fem.equ.gporder{i} can be a number or a cell
array of numbers. In the latter case, fem.equ.gporder{i}{k} applies to the kth
equation in coefficient or general form, and the kth integrand fem.equ.weak{i}{k}
in the weak formulation.
If the field fem.gporder is present, another interpretation of the syntax is used. Then
fem.gporder can be a cell array with quadrature formula orders, and the numbers in
fem.equ.gporder are indices into fem.gporder. When mesh cases are used, the field
fem.gporder is instead a structure with the fields
• fem.gporder.default, which is a cell array of quadrature formula orders used in
mesh case 0. The numbers in fem.equ.gporder are indices into this cell array.
These quadrature formula orders are also used for all mesh cases that do not occur
in fem.gporder.mind.
• fem.gporder.case, which is a cell array of cell arrays of quadrature formula orders.
fem.gporder.case{i} applies to the mesh cases fem.gporder.mind{i}. The
numbers in fem.equ.gporder are indices into this cell array. The default is an
empty cell array fem.gporder.case.
• fem.gporder.mind, which is a cell array of vectors of positive mesh case numbers.
The default is fem.gporder.mind = {1 2 … n}, where n is the length of
fem.gporder.case.
Similarly, fem.bnd.gporder, fem.edg.gporder, and fem.pnt.gporder give the
order for integrals over boundaries, edges, and vertices, respectively.
The default value of the gporder fields is twice the maximum order of the shape
functions you are using.
C P O R D E R — D E N S I T Y O F P O I N T W I S E C O N S T R A I N T S
The pointwise constraints are enforced in the Lagrange points of a certain order. This
order is given in the fields fem.equ.cporder, fem.bnd.cporder, fem.edg.cporder,
and fem.pnt.cporder. These fields have the same syntax as the gporder fields.
S P E C I F Y I N G A M O D E L | 69
70 | C H A P T E R
The field fem.equ.cporder gives the order for constraints on subdomains. The field
fem.equ.cporder can be a number or a cell array. In the first case fem.equ.cporder
applies to all subdomain groups. In the second case, fem.equ.cporder{i} applies to
the ith subdomain group. The entry fem.equ.cporder{i} can be a number or a cell
array of numbers. In the latter case, fem.equ.cporder{i}{k} applies to the kth
constraint on subdomain group i.
Like for gporder, another interpretation of the syntax is used if the field fem.cporder
is present; see the previous section on gporder.
Similarly, fem.bnd.cporder, fem.edg.cporder, and fem.pnt.cporder give the
order for constraints on boundaries, edges, and vertices, respectively. The default value
of the cporder fields is equal to the maximum of the orders of the shape functions you
are using.
If you use Hermite or Argyris shape functions and the boundary is curved, a locking
phenomenon can occur due to the way constraints are implemented in COMSOL
Multiphysics. That is, if you choose cporder to be the same as the element order, then
the derivative degrees of freedom will be over-constrained at the mesh vertices on the
boundary. To prevent this, you can use a cporder that is the element order minus 1.
Initial Values
If you have a time-dependent problem, you need to specify initial values for your
variables. Also, for a nonlinear problem it helps to provide a good starting guess to the
solver. Initial values are specified in the fields fem.equ.init, fem.bnd.init,
fem.edg.init, and fem.pnt.init. There are also other ways to specify initial values;
see asseminit in the Command Reference. If initial values are not specified, they are
taken to be 0.
Initial values on subdomains are given in the field fem.equ.init. This field has the
same syntax as fem.equ.f. That is, fem.equ.init is a numeric scalar, an expression,
or a cell array. In the first two cases, this defines the initial value for all dim variables on
all subdomain groups. In the third case fem.equ.init{i} applies to the ith
subdomain group (but if the cell array fem.equ.init has length 1, then its single
entry applies to all subdomain groups). The entry fem.equ.init{i} can be a numeric
scalar, an expression, or a cell array of scalars or expressions. In the first two cases,
fem.equ.init{i} is the initial value for all dim variables. In the third case, the
components of the cell array fem.equ.init{i} correspond to the different dim
variables.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
In a similar way, fem.equ.dinit defines the initial value for the time derivative when
solving wave-equation problems.
The definitions given in fem.equ.init are applied only where they have meaning,
that is, only where the corresponding dim variable is defined. Similarly, initial
conditions on boundaries, edges, and vertices can be given in fem.bnd.init,
fem.edg.init, and fem.pnt.init, although it is normally not needed. If a variable
is defined only on a boundary (for instance, a Lagrange multiplier in a weak
constraint), then it is useful to define initial values for that variable using
fem.bnd.init.
The solution object corresponding to an initial value specification can be computed by
the function asseminit (assuming that the extended mesh has been generated):
sol = asseminit(fem);
See the Command Reference entry femsol for information about the solution object.
The Extended Mesh
Once you have completed your model specification, and before you use one of
COMSOL Multiphysics’ solvers, you must issue the command
fem.xmesh = meshextend(fem);
The function meshextend takes your model specification and puts it in a suitable form
for the solvers. The resulting object fem.xmesh is called the extended mesh. Thus, if
you change your model in any way, you must run meshextend before you solve again.
The function meshextend does the following:
• Converts the standard syntax to element syntax. The result is taken together with
the element syntax you have specified.
• Determines a partition of the domains into domain groups. Within each domain
group there are the same variables, the same equations, and the same constraints.
Similarly, the mesh elements are partitioned into mesh element groups. Within each
mesh element group, the mesh elements are of the same type and the same shape
functions are used for the variables.
• Collects all the node points that are introduced when using higher-order element.
• Collects all the degree of freedom names.
S P E C I F Y I N G A M O D E L | 71
72 | C H A P T E R
• Collects all degrees of freedom (DOF). A DOF is a pair consisting of a DOF name
and a node point. Each DOF is given a number and corresponds to an entry in the
solution object.
• Collects all the variables that are defined and checks for name collisions.
If you prefer to only use the element syntax, you can turn off the conversion from
standard syntax by
fem.xmesh = meshextend(fem, ‘standard’,’off’);
Also, you can turn off generation of boundary coupled variables by
fem.xmesh = meshextend(fem, ‘cplbndsh’,’off’, ‘cplbndeq’,’off’);
S O L U T I O N F O R M
You can generate the extended mesh using a form of the equation which is different
from fem.form. This is called the solution form as is specified in the field
fem.solform. The function meshextend reads this field and transforms the equations
to this form before generating the extended mesh. The possible values for
fem.solform are coefficient, general, and weak (the default).
Multiple Geometries
Up to now the discussion covers how to set up a FEM structure for a model with a
single geometry. If you have several geometries in a model, you must use an extended
FEM structure called xfem. The extended FEM structure has a field xfem.fem, which
is a cell array of ordinary FEM structures, one for each geometry. These FEM
structures can contain all the fields described above, except const, elem, eleminit,
xmesh, init, and sol. These fields are instead added to xfem. (There can be some
more fields in the extended FEM structure; see the Command Reference entry on
femstruct.)
Find a simple example of the use of two geometries below. First, define the single
geometry FEM structures fem1 and fem2:
fem1.geom = solid1([0 1]);
fem1.mesh = meshinit(fem1);
fem1.equ.c = 1;
fem1.equ.f = 1;
fem1.bnd.h = 1;
fem2.geom = rect2;
fem2.mesh = meshinit(fem2);
fem2.equ.c = 1;
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
fem2.equ.f = 1;
fem2.bnd.h = 1;
Thus you have uncoupled Poisson’s equations with homogeneous Dirichlet conditions
on the unit interval and the unit square. The dependent variable u has linear elements
by default. Now define the extended FEM structure, extend the mesh, solve, and plot:
xfem.fem = {fem1 fem2};
xfem.xmesh = meshextend(xfem);
xfem.sol = femlin(xfem);
postplot(xfem, ‘geomnum’,1, ‘liny’,’u’);
figure
postplot(xfem, ‘geomnum’,2, ‘tridata’,’u’, ‘tribar’,’on’);
Now the solution object xfem.sol.u contains values of the DOFs for both
geometries.
The Application Structure
If you want to use predefined equations for a certain application, you can use
application modes. Instead of specifying the equations explicitly, just enter values for
the relevant physical parameters. For setting up a problem using application modes at
the command line, application structures are used. You can have several application
structures at the same time, describing different physics phenomena in a multiphysics
problem. Once the application structures have been set up, the function
multiphysics can translate them to equations and boundary conditions in the FEM
structure.
The application structures are stored in the fem.appl field in the FEM structure. This
field can be a cell array of application structures or just a single application structure.
These structures have many fields in common with the FEM structure.
The table below lists the fields in the appl structure.
TABLE 1-13: APPL STRUCTURE FIELDS
FIELD INTERPRETATION
appl.assignsuffix Application mode variable name suffix
appl.assign Application mode variable name assignments
appl.bnd Boundary coefficients/application data
appl.border Assembly on interior boundaries
appl.dim Cell array with names of the dependent variable names
appl.edg Edge coefficients/application data
S P E C I F Y I N G A M O D E L | 73
74 | C H A P T E R
The fields in the FEM structure that are affected by the application structure are the
dim, form, equ, bnd, edg, pnt, var, elemmph, eleminitmph, shape, sshape, and
border fields.
D E F I N I N G T H E A P P L I C A T I O N M O D E
Which application mode the application structure refers to is defined in the mode field.
The mode field is a structure with two fields, the class field specifying the application
mode class and a type field specifying if the model is Cartesian or axisymmetric. For
example
appl.mode.class = ‘Electrostatics’;
specifies the Electrostatics application mode. To make the model axisymmetric use
appl.mode.type = ‘axi’;
The default is that the model is Cartesian, and in that case the type field can be
omitted. The class name can then be given directly to the mode field. For example,
appl.mode = ‘Electrostatics’;
specifies a Cartesian Electrostatics application mode.
The application modes are divided into the following categories:
• Acoustics
• Convection and Diffusion
• Electromagnetics
• Fluid Dynamics
• Heat Transfer
appl.equ Subdomain coefficients/application data
appl.mode Application mode
appl.name Application mode name
appl.pnt Point coefficients/application data
appl.prop Application mode properties
appl.shape Shape functions. A cell array of shape function objects or
strings.
appl.sshape Preferred element geometry order. An integer.
appl.var Application mode specific variables
TABLE 1-13: APPL STRUCTURE FIELDS
FIELD INTERPRETATION
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
• Structural Mechanics
• PDE modes
In the additional modules for chemical engineering, earth sciences, electromagnetics,
heat transfer, MEMS, and structural mechanics there are tailored application modes
from these application areas.
For an overview of the different application modes, see the COMSOL Multiphysics
Modeling Guide.
E X P A N D I N G A P P L I C A T I O N M O D E D A T A T O T H E F E M S T R U C T U R E
Once you have defined one or several application structures, the data in them can be
used to generate fields in the FEM structure. This is done by using the multiphysics
function.
The names of the dependent variables, the fem.equ.dim field, are given directly from
the corresponding field in the application structure. If there is no such field there, a
default value for the application mode is used. If several application modes are used to
generate the FEM structure fields, the dim field in the FEM structure is a cell array
with the dependent variables of all application modes.
For the field border there are important differences between the application structure
and the FEM structure. The border field in the application structure is used for
disabling boundary settings on interior boundaries when generating the boundary
conditions in the FEM structure. The field fem.border is set to 1 by multiphysics.
A P P L . E Q U — S U B D O M A I N S E T T I N G S
When using the PDE modes, the Classical PDEs or the weak modes, the subfields in
appl.equ have direct correspondence to the same fields in the FEM structure. You can
therefore define the equ field in one of these application modes in the same way as
when working directly with the FEM structure, as was described in “Equations and
Constraints” on page 56. One advantage with working with application structures is
that you can generate the weak form from a problem defined in general form, and if
you have defined an application in coefficient form you can generate a FEM structure
using the general or weak forms. It is possible to do a similar thing using the FEM
structure and the function flform, but then you lose the original problem
formulation. For more information about flform, see the Command Reference.
For physics modes there is a much clearer difference between the fields in the
application structure and the fields in the FEM structure. The quantities available for
each physics mode correspond to some physical property appearing in the equation
S P E C I F Y I N G A M O D E L | 75
76 | C H A P T E R
defined by the application mode. For example, the default values of the fields in
appl.equ defined by the Heat Transfer by Conduction application mode are
k: {‘400’}
ktensor: {{2×2 cell}}
ktype: {‘iso’}
Dts: {‘1’}
rho: {‘8700’}
C: {‘385’}
Q: {‘0’}
htrans: {‘0’}
Text: {‘0’}
Ctrans: {‘0’}
Tambtrans: {‘0’}
For instance, appl.equ.k is the thermal conductivity, appl.equ.rho is the density,
and appl.equ.C is the heat capacity. Now consider an example with two materials. To
use other values than the default for the physical properties, enter them as follows.
appl.equ.k = {‘300’ ‘500’};
This means that k will have the values 300 and 500 in subdomain groups 1 and 2,
respectively.
In addition to the fields discussed above, the appl.equ structure can have the fields
shape, usage, gporder, cporder, and init. See multiphysics for details. The fields
shape, gporder, cporder, and init have the same syntax as the corresponding FEM
structure fields.
A P P L . B N D — B O U N D A R Y S E T T I N G S
One difference between setting up boundary conditions in the application structures
and the FEM structure is the use of different boundary condition types. This is true
also for the PDE coefficient and weak application modes. This means that, for example,
if the boundary condition on a certain boundary is defined to be of Neumann type,
the h and r coefficients in the application structure are not used when generating fields
in the FEM structure.
For PDE modes there are two boundary condition types available. For Dirichlet
conditions the type is labeled ‘dir’ and for Neumann conditions it is ‘neu’. The
following example sets a Neumann condition on boundary 1 and Dirichlet conditions
on all other.
appl.mode = ‘FlPDEC’;
appl.bnd.r = {1 1};
appl.bnd.h = {1 1};
appl.bnd.type = {‘neu’ ‘dir’};
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
appl.bnd.ind = [1 2 2 2];
For physics modes it is of the same importance to define the boundary condition type.
For each type there can be a number of additional quantities that also need to be set
in order to specify the boundary condition completely.
The default values in the appl.bnd structure for the Heat Transfer by Conduction
application mode are
q0: {‘0’}
h: {‘0’}
Tinf: {‘0’}
Const: {‘0’}
Tamb: {‘0’}
T0: {‘0’}
type: {‘q0’}
To change boundary conditions, simply enter,
appl.bnd.type = {‘q0’ ‘T’};
appl.bnd.T0 = {‘0’ ‘2’};
The type q0 corresponds to the thermal insulation boundary condition and the type T
sets the temperature to the value given by appl.bnd.T0. Thus, in this case we have
thermal insulation on boundary group number 1, and the temperature T = 2 on
boundary group number 2. The ‘0’ in appl.bnd.T0 can be replaced by anything,
since when appl.bnd.type = ‘q0′ appl.bnd.T0 is not used.
In addition to the fields discussed above, the appl.bnd structure can have the fields
shape, gporder, and cporder. For certain application modes which do not have an
appl.equ field, the appl.bnd structure can also contain the fields shape and usage.
The fields shape, gporder, cporder, and init have the same syntax as the
corresponding FEM structure fields. See multiphysics for details.
A P P L . E D G A N D A P P L . P N T — E D G E A N D P O I N T S E T T I N G S
In the same way as appl.equ and appl.bnd are used to defined to define the equation
and boundary conditions, appl.edg can be used to define equations and conditions
on the edges (3D) and appl.pnt on points (2D and 3D).
A P P L I C A T I O N – S P E C I F I C V A R I A B L E S
For many of the application modes, additional variables are generated. These are
similar to the expression variables described in “Expression Variables” on page 53, but
the names and definitions of these application specific variables are predefined.
S P E C I F Y I N G A M O D E L | 77
78 | C H A P T E R
The application mode variables defined on subdomains, boundaries, edges, and
points are automatically generated and stored in var fields in the equ, bnd, edg, and
pnt fields in the FEM structure. There are no corresponding var fields in the
application structure though, because you should not change the expression that will
be used when generating the variables. You can however modify the fields in the FEM
structure, but that is not recommended. Use the expression variables in the expr fields
instead.
For the application scalar variables, which are defined everywhere and stored in
appl.var, you can change the value in the application structure. For example, this
feature is used in the AC Power Electromagnetics application mode, which defines a
frequency nu in the appl.var field. It can be set as
appl.var.nu = ’50’;
or as a cell array, that is,
appl.var = {‘nu’ ’50’};
A S S I G N E D V A R I A B L E N A M E S
Although the internal names used in the application structure are predefined, you can
assign another name for the variable to use in the FEM structure. This is important so
that variable name conflicts can be avoided. To change a name of an application mode
variable or an application scalar variable, use the assign field. If you set up the
assigned field name as
appl.assign = {‘nu’ ‘freq’}
before calling the multiphysics function, the generated FEM structure will have the
var field
fem.var = {‘freq’ ’50’}
There is also a field assignsuffix, which can be used to add a suffix to the name of
all application mode variables. Variables appearing in appl.assign will use the given
assigned name, while the other gets the suffix defined by appl.assignsuffix added
to their names.
A P P L I C A T I O N M O D E P R O P E R T I E S
Some application modes define properties to, for example, specify which type of
analysis to perform. They are specified in the prop field. As an example both
Magnetostatics and AC Power Electromagnetics are implemented in the same
application mode, FlPerpendicularCurrents. This is because they define the same
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
type of physics, one for static fields and one for time-harmonic fields. The application
mode has an analysis property to specify which type of analysis to perform.
appl.mode = ‘FlPerpendicularCurrents’;
appl.prop.analysis = ‘static’;
gives magnetostatics, while
appl.mode = ‘FlPerpendicularCurrents’;
appl.prop.analysis = ‘harmonic’;
gives AC power electromagnetics.
Default Element Property
In the application structure you can directly specify the contents of the fields shape,
sshape, ***.shape, ***.gporder, and ***.cporder, which are then moved into the
appropriate parts of the FEM structure. There is, however, a simpler way of defining
these properties in the application structure. You can use the elemdefault property.
The elemdefault property is used to generate default values of the shape, sshape,
gporder, and cporder fields when these are missing in the application structure.
For the PDE coefficient application mode, the default shape function used is the
quadratic Lagrange element, with suitable values for sshape, gporder, and cporder.
To solve a problem with linear Lagrange elements instead, simply enter
appl.prop.elemdefault = ‘Lag1’;
and the shape functions and other fields in the FEM structure will be generated in
correspondence with the default values for the linear Lagrange element. Because this
field is only available in the application structure, this is another advantage with using
application modes for equation modeling.
E X A M P L E : N A V I E R ’ S E Q U A T I O N R E V I S I T E D
For examples using physics modes, see the Command Reference entry on
multiphysics. The example below uses the PDE Coefficient Form application mode.
The example of solving Navier’s equation on a rectangle is demonstrated on page 61
using the equ and bnd fields of the FEM structure directly. Coefficient form as well as
general form are used, which requires that you redefine the PDE coefficients and
boundary conditions.
If the application structure is used for setting up the problem, you can use
multiphysics repeatedly for redefining the problem. First set the geometry and mesh
in the FEM structure as in the previous example.
S P E C I F Y I N G A M O D E L | 79
80 | C H A P T E R
clear fem
fem.geom = rect2(0, 1, 0, 0.1);
fem.mesh = meshinit(fem);
The constants can be defined in the fem.expr field (although fem.const is more
efficient).
fem.expr = {‘E’ 1e2 ‘nu’ 0.3 ‘G’ ‘E/2/(1+nu)’ …
‘mu’ ‘2*G*nu/(1-nu)’ ‘Fx’ ‘0’ ‘Fy’ ‘-1’};
To set up the application structure, the name of the application mode and the nature
of the subdomain and boundary settings must be specified. In this case, the PDE
coefficient application mode for two variables is used.
clear appl
appl.mode = ‘FlPDEC’;
appl.dim = {‘u’ ‘v’ ‘u_t’ ‘v_t’};
appl.equ.c = { { {‘2*G+mu’ ‘G’} {0 ‘mu’; ‘G’ 0}…
{‘G’ ‘2*G+mu’} } };
appl.equ.f = { {‘Fx’ ‘Fy’} };
appl.bnd.h = { 0 1 };
appl.bnd.type = {‘neu’ ‘dir’};
appl.bnd.ind = [1 1 1 2];
fem.appl = appl;
To generate the fields in the FEM structure enter
fem = multiphysics(fem);
The problem is solved and visualized with
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem, ‘v’)
The above solution uses quadratic Lagrange elements by default. To solve the problem
using linear Lagrange elements, make the following changes:
fem.appl.prop.elemdefault = ‘Lag1’;
fem = multiphysics(fem);
fem.xmesh = meshextend(fem);
fem.sol = femlin(fem);
postsurf(fem, ‘v’)
Overview of Postprocessing Functions
COMSOL Multiphysics has a set of functions specifically designed for postprocessing.
The functions convert solution data from the solvers. The term expression in the
following table refers to an expression containing any of the variables described in
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
“Using Variables and Expressions” on page 142 in the COMSOL Multiphysics User’s
Guide.
The term domain below refers to a geometric domain of a certain space dimension. For
instance, in 3D, there are subdomains (dimension 3), boundaries (dimension 2), edges
(dimension 1), and vertices (dimension 0).
There is also a set of visualization functions. The general plot command is postplot.
In order to conveniently plot an expression in some standard ways, COMSOL
Multiphysics provides a set of shorthand commands for postplot:
The postmovie function animates any postplot plots. The function postanim is a
shorthand function for standard animations.
The function postcrossplot plots expressions in cross sections of the geometry of
any space dimension. It can also plot expressions on any domain. For time-dependent
problems, postcrossplot also provides time extrusion plots, and, letting the
dimension of the cross section be zero, node plots.
FUNCTION DESCRIPTION
posteval Evaluates an expression on any domain and returns the result in
postdata format
postint Integrates an expression on any domain
postinterp Evaluates an expression on any domain of dimension >0 in points
specified by global coordinates (or parameter values on a geometry
face or edge)
FUNCTION DESCRIPTION WORKS IN
postarrow Arrow plot on subdomains 2D, 3D
postarrowbnd Arrow plot on boundaries 2D, 3D
postcont Contour plot 2D
postflow Flow line plot 2D, 3D
postiso Isosurface plot 3D
postlin Line plot 1D, 2D, 3D
postslice Slice plot 3D
postsurf Surface plot 2D, 3D
posttet Subdomain plot 3D
S P E C I F Y I N G A M O D E L | 81
82 | C H A P T E R
Overview of Common Solver Properties and Their Values
The following table shows the property/value pairs of some important solver
functions: adaption, assemble, femeig, femlin, femnlin, and femtime. From the
table you can determine where to look for an explanation of a property value pair. The
entry “r” in the table indicates that information about the property/value pair is
available in the Command Reference entry on that function. The entry “u” indicates
that the function uses that property/value pair. No entry indicates that the function
does not accept that property/value pair.
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
Amgauto Quality of multigrid
hierarchy (AMG)
integer between
1 and 10
u u u u u
Atol Absolute tolerance numeric scalar or
initial value
r
Blocksize Assembly block size positive integer u r u u u u
Complex Complex numbers on | off r
Complexfun Use complex-valued
functions with real
input
on | off u u u u u u
Conjugate Use complex
conjugate, Hermitian
transpose
on | off u u u u u
Consistent Consistent
initialization of DAE
system
off | on |
bweuler
r
Const Definition of
constants
cell array u r u u u u
Csolver Coarse solver (AMG,
GMG)
string u u u u u
Csolverpar Coarse solver
property/value list
cell array u u u u u
Droptol Drop tolerance
(LUINC,
TAUCS_LLT,
UMFPACK, SPOOLES)
scalar between 0
and 1
u u u u u
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Eigselect Eigenvalue adaption
parameter
numeric array r
Estrat Error estimation
strategy
0 | 1 r
Etol Eigenvalue tolerance non-negative
scalar
r
Fillratio Fill ratio (LUINC) non-negative
scalar
u u u u u
Geomnum Geometry number positive integer r
Hnlin Highly nonlinear
problem
off | on u r
In Input matrices cell array of
names and
matrices
r r r
Init Initial value u u u u u
Initialstep Initial time step positive scalar r
Initstep Initial damping
parameter
positive scalar u r
Iter Number of iterations
(preconditioners and
smoothers)
positive integer u u u u u
Iterative Use iterative solver off | on r
Itol Iterative tolerance positive scalar u r r
Itrestart Number of iterations
before restart
(GMRES)
positive integer u u u u u
Keep Manual control of
reassembly
string u u u u
Krylovdim Dimension of Krylov
space
positive integer u r
L2scale Scale factors for the
L2 error norm
vector of positive
numbers
r
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
S P E C I F Y I N G A M O D E L | 83
84 | C H A P T E R
L2staborder Orders in the stability
estimate for the L2
error estimate
vector of positive
integers
r
Linsolver Direct linear system
solver
string u u u u u
Masssingular Singular mass matrix yes | maybe r
Maxcoarsedof Maximum number of
DOFs on coarsest
level (AMG)
positive integer u u u u u
Maxdepth Maximum recursion
depth
(TAUCS_LLT_MF)
positive integer u u u u u
Maxiter Maximum number of
iterations
positive integer u r
Maxlinit Maximum number of
linear iterations
positive integer u u u u u
Maxorder Maximum BDF order integer between
1 and 5
r
Maxstep Maximum time step positive scalar r
Maxt Maximum number of
mesh elements
positive scalar r
Mcase Mesh case(s) vector of integers u u u u u u
Mcasekeep Keep generated mesh
cases (GMG)
on | off u u u u u
Meshscale Mesh scaling factors
(GMG)
vector of positive
numbers
u u u u u
Method Constraint handling
method
eliminate |
elimlagr |
lagrange |
spring
u u u u
Mgassem Assemble on coarse
levels (GMG)
on | off | numeric
array
u u u u u
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Mgauto Method for mesh case
generation (GMG)
off | explicit |
meshscale | shape
| both
u u u u u
Mgcycle Multigrid cycle type
(AMG, GMG)
V | W u u u u u
Mggeom Geometries to use in
multigrid (GMG)
vector of
integers
u u u u u
Mglevels Maximum number of
levels (AMG, GMG)
positive integer u u u u u
Minstep Minimum step size positive scalar u r
Modified Modified incomplete
Cholesky
(TAUCS_LLT)
on | off u u u u u
Neigs Number of
eigenvalues sought
positive integer
Ngen Maximum number of
meshes
positive integer r
Nonlin Nonlinear problem off | on r
Ntol Nonlinear tolerance positive scalar u r
Nullfun Null space function auto | flnullorth |
flspnull
u u u u u
Out Output arguments cell array of
strings
r r r r r r
Outcomp Solution components
to store in output
cell arrray of
strings
u u u u u
Pinitstep Initial stepsize for
continuation
positive real u r
Plist List of parameter
values for
continuation
real vector u r
Pmaxstep Maximum stepsize for
continuation
positive real u r
Pminstep Minimum stepsize for
continuation
positive real u r
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
S P E C I F Y I N G A M O D E L | 85
86 | C H A P T E R
Pname Name of parameter string u r
Porder Predictor order for
continuation
0 | 1 u r
Postsmooth Postsmoother (AMG,
GMG)
string u u u u u
Postsmoothpar Postsmoother
property/value list
(AMG, GMG)
cell array u u u u u
Prefun Preconditioner string u u u u u
Preorder Preordering algorithm
(SPOOLES)
mmd | nd | ms |
both
u u u u u
Prepar Preconditioner
property/value list
cell array u u u u u
Presmooth Presmoother (AMG,
GMG)
string u u u u u
Presmoothpar Presmoother
property/value list
(AMG, GMG)
cell array u u u u u
Relax Relaxation factor
(SSOR, SOR, SORU,
SSORVEC, SORVEC,
SORUVEC, JAC)
positive scalar u u u u u
Report Show progress
window
off | on r r r r r r
Resorder Order of decrease of
equation residuals
scalar | vector r
Respectpattern Do not drop original
nonzeros (LUINC)
on | off u u u u u
Rhob Factor in linear error
estimate (iterative
solvers)
real scalar > 1 u u u u u
Rmethod Refinement method regular |
longest |
meshinit
r
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
Rowscale Row equilibration on | off u u u u
Rstep Restriction for
stepsize update
real scalar > 1 u u
Rtol Relative tolerance positive scalar r
Seconditer Number of secondary
iterations (SSORVEC,
SORVEC, SORUVEC)
positive integer
Shapechg Change in shape
function orders
(GMG)
vector of
integers
u u u u u
Shift Eigenvalue search
location
scalar u r
Solcomp Degree of freedom
names to solve for
cell array of
strings
u r u u u u
Solver Solver type stationary |
eigenvalue
r
Stop Deliver partial
solution when failing
off | on r r r r r
T Time value real scalar r
Thresh Pivot threshold
(UMFPACK,
SPOOLES, LUINC)
scalar between 0
and 1
u u u u u
Tlist Time list numeric vector r
Tout Output times tout | tsteps r
Tpfun Element pick function fltpft |
fltpworst |
fltpqty
r
Tppar Parameter passed to
element pick
function
nonnegative
scalar
r
Tsteps Time-stepping mode free |
intermediate |
strict
r
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
S P E C I F Y I N G A M O D E L | 87
88 | C H A P T E R
Associative Geometry and Parametrization
For repeated simulations with modifications in a geometry between each simulation
iteration, you can carry out the geometry modeling such that there is an association
between the building blocks in the geometry modeling and the final geometry used
for the numerical simulation. This concept is called associative geometry.
To set up a model for geometric parameterization, use the function geomanalyze. It
creates one decomposed and analyzed geometry from all the geometry objects
contained in the Draw structure, and it also sets up the associative geometry
information.
You can, for example, create the initial geometry for a model using the command
fem = geomanalyze(fem, {rect2 circ2},’ns’,{‘Rectangle’,’Circle’});
It stores one rectangle object and one circle object in the draw structure; it stores one
decomposed geometry object in the fem.geom field; and it computes the mappings
between the objects in fem.draw and the object in fem.geom.
You can now create a mesh, enter the physics settings, and solve the model. Once you
have studied the results for this first simulation, you might want to move the original
circle and carry out the same simulation again and compare the results. The concept
of associative geometry minimizes the risk that you must reenter physics settings.
You could, for example, do so with the following commands:
circObj = drawgetobj(fem,’Circle’);
fem = drawsetobj(fem,’Circle’,move(circObj,.1,0));
fem = geomanalyze(fem);
U Value of variables not
solved for and
linearization point
solution object |
numeric vector |
scalar
u r r r r r
Umfalloc Memory allocation
factor (UMFPACK)
positive scalar u u u u u
Uscale Scaling of variables auto | init |
none | cell array |
numeric vector
u u u u u
TABLE 1-14: PROPERTY/VALUE PAIRS SHARED BETWEEN SOLVER ALGORITHMS
PROPERTY DESCRIPTION VALUE TYPE
a
d
a
p
t
i
o
n
a
s
s
e
m
b
l
e
f
e
m
e
i
g
f
e
m
l
i
n
f
e
m
n
l
i
n
f
e
m
t
i
m
e
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
This code extracts the circle object, modifies it with the move operation, and sets it in
the fem.draw field. The geomanalyze command updates all fields that are geometry
dependent in the FEM structure, so that you can solve the same problem on the
updated geometry. Of course, there might be cases where complicated geometry
operations result in another geometry mapping instead of the one you expect.
S P E C I F Y I N G A M O D E L | 89
90 | C H A P T E R
T h e S t r u c t u r e o f t h e Mode l M – f i l e
You can save an entire modeling from a GUI session as a sequence of COMSOL
Multiphysics commands. This command sequence is called a Model M-file. You can:
• Run these commands directly from COMSOL Script or MATLAB or load it into
COMSOL Multiphysics
• Modify the Model M-file using COMSOL or MATLAB commands
If you want to open a modified Model M-file in COMSOL Multiphysics, it is
important to follow the syntactic rules outlined the following sections.
The Model M-files consist of several parts describing the COMSOL Multiphysics
model. The semantic and syntactic rules are outlined in the following sections. These
rules do not need to be considered if you modify a Model M-file from the command
line. This section describes the way COMSOL Multiphysics generates a Model M-files.
The order of commands that are recorded in the Model M-file must be the same as in
this document. The order of the commands is very important. Starting with the
geometry section and ending with the multiphysics, meshextend, initial value, solver,
and postprocessing sections (in that order). All of the commands in between here can
occur in any reasonable order. The detailed rules for this order are indicated in the
sections below.
VE R S I O N S E C T I O N
The version section assigns a version to the FEM structure.
fem.version=version;
The version must be a legal COMSOL version structure, that is, it must contain the
fields
version.name=’COMSOL 3.2′;
version.major=0;
version.build=numeric;
where numeric must be an integer less than or equal to the number visible in the About
COMSOL window that you open from the Help menu.
If the version field is missing, it assumes the current COMSOL version.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
S P A C E D I M E N S I O N S E C T I O N
The space dimension section determines the names of the space coordinates
(independent variables).
fem.sdim=sdim;
This section can be omitted.
D R A W S E C T I O N
The draw section determines the geometry objects in the drawing. It is generated each
time you leave draw mode and change the geometry.
fem.draw=draw;
The geometry model, draw, is documented in the Command Reference entry
geomcsg.
This section can be omitted.
Separate File for Geometry Objects in Model M-Files
Sometimes the geometry objects in the draw section in Model M-files, which defines
fem.draw, are not generated from commands but saved in a separate file. This happens
in the following cases:
• When you have used Reset Model in the File menu
• When you have imported the geometry from COMSOL Script or MATLAB
This is because in these cases the information on how to generate the geometry from
commands is lost.
When you save a Model M-file, for example, mymodel.m, COMSOL Multiphysics also
saves a file mymodel.mphm with the geometry objects. In mymodel.m there is a variable
flbinaryfile that specifies the geometry file.
flbinaryfile = ‘mymodel.mphm’;
The geometry objects are defined by statements like
g1 = flbinary(‘g1′,’draw’,flbinaryfile);
The function flbinary loads the object, which is identified by the tags g1 and draw,
from mymodel.mphm.
G E O M E T R Y S E C T I O N
The geometry section defines the analyzed geometry of the model. It is generated each
time that you leave draw mode and change the geometry.
T H E S T R U C T U R E O F T H E M O D E L M – F I L E | 91
92 | C H A P T E R
fem.geom=geomcsg(fem);
You can assign any valid analyzed geometry such as a geometry object, a mesh (in 1D
and 2D), a Decomposed Geometry matrix list from the PDE Toolbox (2D), or a
Geometry M-file to the fem.geom field. This section cannot be omitted. You can
assign any valid analyzed geometry such as a geometry object, a mesh (in 1D and 2D),
a Decomposed Geometry matrix list from the PDE Toolbox (2D), or a Geometry
M-file to the fem.geom field.
This section cannot be omitted.
A P P L I C A T I O N M O D E S E C T I O N
The application mode section determines the application modes.
fem.appl{i}=appl;
When a Model M-file has been loaded, a number of GUI settings are determined by
this section. The appl{i} cell array item contains the fields:
• mode—create the application mode object
• dim—set the dependent variable names
• border—set the flags for assembly on interior boundaries
• name—set the name
• usage—set the subdomain usage flags
• var—determines the Application Mode Variables dialog box settings
• assign—the assigned names
• assignsuffix—application mode variable name suffix
• prop—application mode properties
• shape—set up shape functions and preferred shape order
• sshape—set up shape functions and preferred shape order
• equ—determines the Subdomain Settings dialog box settings
• bnd—determines the Boundary Settings dialog box settings
• edg—determines the Edge Settings dialog box settings (3D only)
• pnt—determines the Point Settings dialog box settings (2D and 3D).
M E S H S E C T I O N
The mesh section contains a command for meshing the geometry.
fem.mesh=mesh
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
The mesh can be computed by a mesh command or directly created by the Model
M-file as a sequence of commands for generating the mesh.
The property/value pairs in the meshinit command are used to restore the settings
in the Mesh Parameters dialog box when the Model M-file is loaded.
The first mesh section in a Model M-file must contain a meshinit statement
fem.mesh=meshinit(fem,…);
This section can be omitted if there is no solver section. At least one mesh section must
occur before a solver section.
See the Command Reference entries on meshinit, meshrefine, and meshsmooth for
more information.
I N T E R I O R B O U N D A R I E S S E C T I O N
The interior boundaries section determines the interior boundaries settings.
fem.border=border;
This section can be omitted. In that case, fem.border=’off’ is assumed.
F O R M S E C T I O N
The form section determines the PDE form. If the Model M-file calls multiphysics,
the outform field determines the PDE form returned from the multiphysics call.
fem.outform=outform;
The multiphysics function returns a FEM structure on the desired PDE form, and
the resulting form is stored in fem.form.
This section can be omitted. In that case, fem.form=’coefficient’ is assumed. The
form section must be compatible with the application mode section. If any application
mode form is in general form, you must have fem.form=’general’.
The form which is used when solving the equation can be different from the form
which is used to enter them. The function meshextend uses the field fem.solform to
determine the solution form of the equations. fem.solform can be omitted. In that
case fem.solform=’weak’ is assumed.
See the section “Using the PDE Modes” on page 224 in the COMSOL Multiphysics
Modeling Guide for more information about the form of a PDE.
T H E S T R U C T U R E O F T H E M O D E L M – F I L E | 93
94 | C H A P T E R
S I M P L I F Y S E C T I O N
The simplify section determines if expressions should be simplified or not. This is used
by the multiphysics call when combining the composite system.
fem.simplify=simplify;
This section can be omitted. In that case, fem.simplify=’on’ is assumed.
See the Command Reference entry multiphysics for more information about
simplify.
F U N C T I O N D E F I N I T I O N S E C T I O N
The function section defines functions used when deriving the composite system.
fem.functions=functions;
This section can be omitted if no functions are needed.
M A T E R I A L L I B R A R Y S E C T I O N
The material library section defines material parameters, coordinate systems, and
cross sections used in the model.
fem.lib=lib
This section can be omitted if no material library is used in the model.
G E O M E T R Y S H A P E O R D E R S E C T I O N
The geometry shape order section defines the order of geometry approximation,
sshape.The sshape is determined by the outsshape field.
fem.outsshape=outsshape;
If omitted, a default value is computed based on the shape functions used in the model.
The resulting geometry shape order is stored in the sshape field.
C O N S T A N T S S E C T I O N
The constants section determines names of constants and their values as defined by the
user.
fem.const=const;
This section can be omitted.
A P P L I C A T I O N S C A L A R V A R I A B L E S S E C T I O N
The application scalar variables section determines variables defined by the
application modes. This is defined for each application mode separately.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
fem.appl{i}.var=var;
This section can be omitted.
E X P R E S S I O N S E C T I O N
The expression section contains expression variables and their defining expressions as
specified by the user.
fem.expr=expr;
These expressions are defined everywhere on the geometry. Domain-based expressions
are defined in the point settings, edge settings, boundary conditions, and PDE
coefficient sections, in the sub-field expr.
This section can be omitted.
S H A P E F U N C T I O N S E C T I O N
The shape function section determines the shape functions used when solving the PDE
problem. This is defined for each application mode separately:
fem.appl{i}.shape=shape;
PO I N T S E T T I N G S S E C T I O N
The point settings section determines the point settings used by the solvers to solve
the PDE problem. These settings are made for each application mode separately:
fem.appl{i}.pnt=pnt;
For 1D models, there is no point section.
E D G E S E T T I N G S S E C T I O N
The edge settings section determines the edge settings used by the solvers to solve the
PDE problem. These settings are made for each application mode separately:
fem.appl{i}.edg=edg;
The edge section only exists for 3D models.
B O U N D A R Y C O N D I T I O N S E C T I O N
The boundary condition section determines the boundary conditions used by the
solvers to solve the PDE problem.
These settings are made for each application mode separately:
fem.appl{i}.bnd=bnd;
T H E S T R U C T U R E O F T H E M O D E L M – F I L E | 95
96 | C H A P T E R
P D E C O E F F I C I E N T S E C T I O N
The PDE coefficient section determines the PDE coefficients used by the solvers to
solve the PDE problem. These settings are made for each application mode separately:
fem.appl{i}.equ=equ;
When a Model M-file has been loaded, the settings made in this section determines the
contents of the Subdomain Settings dialog box for the equation system.
C O U P L I N G V A R I A B L E E L E M E N T S S E C T I O N
The coupling variable elements section defines coupling variable elements. Such
variables can be defined in the dialog boxes for coupling variables.
fem.elemcpl=elemcpl;
M U L T I P H Y S I C S S E C T I O N
The multiphysics section combines the applications in fem.appl to the composite
system.
fem=multiphysics(fem,…);
M E S H E X T E N D S E C T I O N
In the meshextend section the finite element mesh is extended to the desired element
types.
fem.xmesh=meshextend(fem,…);
This section cannot be omitted.
See the Command Reference entry on meshextend for details.
I N I T I A L V A L U E S E C T I O N
The initial value section computes a solution object corresponding to the initial
conditions that is used to solve the PDE problem.
init=asseminit(fem,…);
In some cases when, for example, a mesh change, a geometry change, or a change in
the included dependent variables has occurred, the solution object need to be mapped
onto the new extended mesh before the initial solution object can be computed.
This section can be omitted.
See the Command Reference entry asseminit for more information on initial values.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
S O L V E R S E C T I O N
In the solver section of the Model M-file, one of the solvers is called to produce a
solution structure.
fem.sol=sol;
For some solvers, the full FEM structure must be returned, for example, for the
adaptive solver.
fem=adaption(fem,…);
The property/value pairs in the solver commands are used to restore the settings in
the user interface when the Model M-file is loaded.
If there is a postprocessing section, the solver section must be present.
Immediately following the solver call there is an assignment of the type
fem0=fem;
This copy of the FEM structure is used when using an old solution when computing
the initial value vector, for example, when clicking the Restart button. This line can be
omitted unless following commands refers to the variable fem0.
PO S T P R O C E S S I N G S E C T I O N
The postprocessing section reproduces all the plot types available in COMSOL
Multiphysics.
Standard Plots
The postplot command corresponds to plots available from the Plot Parameters
dialog box (except animations).
postplot(fem,…
‘tridata’,expr,…
‘contdata’,expr,…
…)
The property/value pairs in the postplot command are used to restore the settings
in COMSOL Multiphysics when opening the Model M-file.
This section can be omitted. At least one solver section must precede the
postprocessing section.
Model M-files for Models With Multiple Geometries
A Model M-file for an extended multiphysics model using multiple geometries appear
slightly different from a single geometry model. In the multiple geometry case, there
T H E S T R U C T U R E O F T H E M O D E L M – F I L E | 97
98 | C H A P T E R
is an additional level in the FEM structure. This structure is called the extended FEM
structure, xfem, to differentiate it from the single geometry FEM structure. The top
level of the structure contains a fem field that is a cell array of single geometry
structures. The xfem structure also contains other fields, such as the extended mesh
and other global properties of the model.
For a multiple geometry model, the xfem structure is passed to meshextend,
asseminit, the solvers, and all postprocessing functions. Only one geometry is active
at a time. Prior to a section with geometry specific commands, you have to activate a
geometry. You can do this by the following assignment:
fem = xfem.fem{g}
where g is the geometry number. After such a line, the following geometry-specific
commands are identical to those of a single geometry model. When another geometry
is activated or before the xfem structure is used, the changes to the active geometry
must be stored in the xfem structure again.
xfem.fem{g} = fem
The extended FEM structure can have the fields xfem.fem, xfem.const, xfem.elem,
xfem.eleminit, xfem.xmesh, xmesh.init, and xfem.sol (and some more fields;
see the Command Reference entry on femstruct). The field xfem.fem is a cell array
containing an ordinary FEM structure xfem.fem{g} for each geometry. These
ordinary FEM structures do not contain the fields const, elem, eleminit, xmesh,
init, and sol.
1 : C O M S O L M U L T I P H Y S I C S S C R I P T I N G
2
S i m u l i n k a n d S t a t e – S p a c e E x p o r t
This chapter explains how to use the Simulink and state-space export functions in
COMSOL Multiphysics to include finite element models in dynamic system
simulations. For model examples of the Simulink export, see the chapter
“Multidisciplinary Models” in the COMSOL Multiphysics Model Library.
99
100 | C H A P T E
S imu l i n k Expo r t
Simulink is a general software package for simulation of dynamic systems. It provides
a graphical user interface for intuitive building of models, a block library with a large
number of components, and the possibility to create new blocks. Furthermore, as
Simulink is built on MATLAB, it has access to advanced mathematical and numerical
routines as well as all other products in the MATLAB family.
It is of great interest to integrate static and time-dependent COMSOL Multiphysics
models using Simulink export. This enhances COMSOL Multiphysics with
functionality such as control loop validation and design, signal processing, and
optimization, all easily accessible through Simulink’s user interface. Many complex
problems in these areas involve physical systems, and the use of COMSOL
Multiphysics for the modeling of systems gives much more accurate results than
empirical models.
For examples using Simulink export, see the models “Magnetic Brake” on page 266
and “Controlling Temperature” on page 290 in the COMSOL Multiphysics Model
Library.
Dynamic or Static Export?
By default, COMSOL Multiphysics exports a dynamic model. This means that the
output of the COMSOL Multiphysics Subsystem block can depend on both the input
variables and the Simulink state vector. The COMSOL Multiphysics solution vector is
a part of the Simulink state vector.
If the time scales in a time-dependent COMSOL Multiphysics model are very fast
compared to the time scales in the Simulink model, a static model can be exported.
This means that the terms of the equations that contain time derivatives are neglected,
so that the output of the COMSOL Multiphysics Subsystem block can be computed
from the input variables by solving a stationary PDE problem. If, for example, you use
COMSOL Multiphysics to model electromagnetic phenomena as a part of a
mechanical system, the dynamics of the COMSOL Multiphysics model can be
neglected. The static model can be either linear or nonlinear.
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
General or Linearized Export?
By default, COMSOL Multiphysics uses general export. General export means that
the COMSOL Multiphysics solvers will be called in every simulation step from
Simulink. This is not necessary for linear models. To perform a simulation of a linear
model or to perform a nonlinear simulation close to an equilibrium point, you can use
linearized export. The linearized model is exported as a set of matrices, which gives
faster simulations because the COMSOL Multiphysics solver does not need to be
called during the simulation. It is also possible to simplify the exported model even
further by using model reduction.
The Simulink Export Dialog Box
Open the Export Simulink Model dialog box from the File menu under Export.
Use this dialog box to create a Simulink structure in the MATLAB workspace.
Simulink uses the Simulink structure to view your COMSOL Multiphysics model as
an S-function. Define the variable name for the Simulink structure in the Structure
name edit field. You need to enter this variable name in the COMSOL Multiphysics
Subsystem block in Simulink.
When you have defined your Simulink model parameters, click Apply or OK to create
and save the model as a Simulink structure.
S I M U L I N K E X P O R T | 101
102 | C H A P T E
The Simulink Export Types
There are four different types of export corresponding to combinations of dynamic/
static and general/linearized export. Select the export type in the Simulink export type
list. There is just one COMSOL Multiphysics Subsystem block in Simulink, but it
functions differently depending on your choice of export type.
G E N E R A L D Y N A M I C E X P O R T
Perform export of a general dynamic model, where the COMSOL Multiphysics
degrees of freedom are part of the Simulink state vector. The COMSOL Multiphysics
solver is called several times for each time step to compute the time derivative of the
state vector, with inputs from Simulink. Only linear, time-independent constraints and
the elimination constraint handling method are supported.
G E N E R A L S T A T I C E X P O R T
Export a general static model, where the COMSOL Multiphysics degrees of freedom
are not part of the Simulink state vector. To compute the outputs of the COMSOL
Multiphysics Subsystem block, the COMSOL Multiphysics linear or nonlinear
stationary solver is called for each time step, with inputs from Simulink. Select the Use
nonlinear stationary solver check box to use the nonlinear solver instead of the default
linear solver.
L I N E A R I Z E D D Y N A M I C E X P O R T
Export a linearized dynamic model. The COMSOL Multiphysics model is linearized
about an equilibrium solution, and the matrices in the state-space form are computed.
The COMSOL Multiphysics degrees of freedom are part of the Simulink state vector.
At each time step, the matrices in the state-space form are used to compute the time
derivative of the state vector instead of calling the COMSOL Multiphysics solver. Only
linear, time-independent constraints and the elimination constraint handling method
are supported.
The linearization point should be an equilibrium point (stationary solution). You
control it by using the settings in the Value of variables not solved for and linearization
point area on the Initial Value tab of the Solver Manager dialog box. It can be useful to
first compute a stationary solution and store that solution as the linearization point.
Note that the inputs and the outputs are deviations from the equilibrium values.
You can use model reduction to approximate the linearized model with a model that
has fewer degrees of freedom; select the Model reduction check box to enable it. The
model reduction uses eigenmodes and static modes of the original PDE problem and
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
projects the linearized model onto the subspace of the eigenmodes. You control the
number of eigenmodes using the Number of eigenmodes edit fields and whether or not
to use the static modes using the Include static modes check box.
The model reduction uses all settings in the Solver Parameters dialog box for the
eigenvalue solver. As an exception, the number of eigenvectors is taken from the Export
Simulink Model dialog box.
When using model reduction on wave equation models that have been rewritten as a
system of first-order equations (wave extension), the algorithm needs to know the
names of the original (non time derivative) solution components. For structural
mechanics application modes, these names are determined automatically. If you are
using wave extension in another context, you need to perform the export from the
MATLAB command prompt with the function femsim (see the Command
Reference).
L I N E A R I Z E D S T A T I C E X P O R T
Export a linearized static model. The COMSOL Multiphysics model is linearized
about an equilibrium solution and a transfer matrix is computed. The COMSOL
Multiphysics degrees of freedom are not part of the Simulink state vector. To compute
the outputs of the COMSOL Multiphysics Subsystem block, the transfer matrix is
used.
Relation to Solver Parameters and Solver Manager Settings
The Simulink export considers the settings from the Solver Parameters and Solver
Manager dialog boxes.
For the dynamic block types, the export only considers certain solver parameters: the
setting Matrix symmetry on the General tab and the settings on the Advanced tab in the
Solver parameters dialog box. The time-dependent simulation in Simulink considers
the settings in the Scaling of variables and Manual control of assembly areas during the
simulation. Simulink itself provides the settings for the time stepping.
For the static block types, the solver type setting in COMSOL Multiphysics controls if
the static solver is linear or nonlinear. All settings related to the linear or nonlinear
solvers in Solver Parameters are considered for the static export blocks.
For the dynamic block types, COMSOL Multiphysics computes the initial value of the
Simulink simulation using the settings in the Initial value area on the Initial Value tab
of the Solver Manager dialog box. For choices that include initial value expressions,
S I M U L I N K E X P O R T | 103
104 | C H A P T E
these software takes these expressions from the Init tab in the Subdomain Settings,
Boundary Settings, Edge Settings, and Point Settings dialog boxes.
The linearized block types take the linearization points according to the settings in the
Value of variables not solved for and linearization point area on the Initial Value tab of
Solver Manager dialog box.
Also, the exported models take the settings on the Solve for tab into account. The
variables not solved for are taken from the Value of variables not solved for and
linearization point area on the Initial Value tab of the Solver Manager dialog box.
The Input/Output Variables
To make use of a PDE model within Simulink, it needs inputs and outputs. Use the
Variables tab to define the inputs and outputs. The input and output variable names
appear in the COMSOL Multiphysics Subsystem block as inputs and outputs,
respectively. For a multigeometry model, you control what geometry the inputs and
outputs belong to by opening the Export Simulink Model dialog box from that
geometry.
C R E A T I N G I N P U T V A R I A B L E S
Enter input variables in the Variables edit field in the Input area. Use a space-separated
list of variable names. Use the variables in expressions, coupling variables, or in the
physics settings.
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
Note: For dynamic block types, the input variables must not occur in the Dirichlet
boundary conditions (constraints). Avoid this problem by approximating Dirichlet
boundary conditions with a Neumann boundary condition using the “stiff spring”
method. See the model “Controlling Temperature” on page 290 of the COMSOL
Multiphysics Model Library, a heat flux condition with a heat transfer coefficient is
used instead of a Dirichlet condition for the temperature.
C R E A T I N G O U T P U T V A R I A B L E S
Define output variables in the Output area. Initially, there are no output variables. Add
an output variable by clicking Add. Enter the output variable name in the Variable name
edit field. Define the value of the variable name by an expression. Click the Subdomain
option button to evaluate the expression in any location by providing the coordinates,
or click the Point option button to evaluate the expression in a geometry vertex
(point). Then select a quantity from the Predefined quantities list or type in an
expression in the Expression edit field.
M O D I F Y I N G O R D E L E T I N G A V A R I A B L E
You can modify a variable by selecting it in Variable number list and then modifying it.
Similarly, to remove a variable, select it in the Variable number list and then click Delete.
Using the COMSOL Multiphysics Model in Simulink
The result of the COMSOL Multiphysics Simulink export is a structure in the
MATLAB main workspace. To use the model in Simulink, open the Blocksets &
Toolboxes library in Simulink, double-click the COMSOL Multiphysics icon, and drag
the COMSOL Multiphysics Subsystem block to your Simulink model. Double-click your
S I M U L I N K E X P O R T | 105
106 | C H A P T E
copy of the block and enter the name of the Simulink structure. This sets up the input
and the output ports of the block as in the example below.
A COMSOL Multiphysics Subsystem block in Simulink with two inputs and one output.
Because COMSOL Multiphysics models usually are stiff, we recommend using an
implicit stiff ODE solver like ode15s in Simulink’s Simulation Parameters dialog box.
Note: Simulink does not support complex-valued inputs or outputs to the COMSOL
Multiphysics Subsystem block.
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
S t a t e – S pa c e E xpo r t
Use state-space export to create a linearized state-space model corresponding to your
COMSOL Multiphysics model. The state-space export can create a state-space object
suitable for use in the Control System Toolbox for MATLAB. You can also export the
matrices of the state-space form directly to the MATLAB workspace.
The State-Space Form
The vector x is called the state vector. The state-space form of the model is
where A, B, C, and D are constant matrices. The constraints (Dirichlet boundary
conditions) are eliminated by writing the solution vector U as
where U0 is the linearization point, and Null is the null-space matrix for the
constraints.
You can also export the system on the equivalent form
This form is more suitable for large systems, because the matrices M and MA usually
become much more sparse than A.
If the mass matrix M is small, it is possible to approximate the dynamic state-space
model with a static model, where M = 0:
Equivalently, y = Hu, where H = D – CA-1B is the transfer matrix.
td
dx Ax Bu+=
y Cx Du+=⎩
⎪
⎨
⎪
⎧
U Null x U0+=
M
td
dx
MAx MBu+=
y Cx Du+=⎩
⎪
⎨
⎪
⎧
0 Ax Bu+=
y Cx Du+=⎩
⎨
⎧
S T A T E – S P A C E E X P O R T | 107
108 | C H A P T E
Model Reduction
Finite element models often contain a large number of degrees of freedom, and the
exported state-space models then become very large. Model reduction makes it
possible to reduce the size of the system while keeping some of the important model
properties.
First, the mode-reduction algorithm creates a matrix T whose columns consist of
• A given number of the eigenmodes with lowest eigenvalues
• The static solutions for unit values on each of the input variables
Next, it introduces state variables as that projects the state-space equations
onto the subspace spanned by T:
or, equivalently,
.
This is a new state-space form with the number of degrees of freedom equal to the
number of modes in the basis T. You can then compute the solution vector by
,
where U0 is the linearization point.
The quality of the reduced system depends heavily on the properties of the system.
Including eigenmodes with small eigenvalues ensures correct modeling of the low
frequencies of the system, and including static solutions ensures correct stationary
values. If the high frequencies of the system have a large influence on the outputs,
however, the results may be unacceptable.
x̃ x Tx̃=
T
T
MT
td
d x̃ T
T
MATx̃ T
T
MBu+=
y CTx̃ Du+=⎩
⎪
⎨
⎪
⎧
M̃
td
d
x̃ MA˜ x̃ MB˜ u+=
y C̃x̃ Du+=⎩
⎪
⎨
⎪
⎧
U Null x̃˜ U0+=
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
The State-Space Export Dialog Box
Open the Export State-Space Model dialog box from the File menu under Export.
Use this dialog box to create a structure containing the state-space form matrices or a
state-space object for the MATLAB Control System toolbox in the MATLAB
workspace. Define the variable name for the state-space object in the Structure/object
name edit field. Use this variable name to refer to the result of the state-space export
in MATLAB.
When you have defined the export parameters, click Apply or OK to perform the export
to the MATLAB workspace.
Select the Use sparse matrices check box if you want to export the state-space model as
sparse matrices. This keeps the sparsity of the matrices but is not recommended if you
want to perform further analysis, for example, in the Control System Toolbox, because
sparse matrices may not be supported.
The linearization point should be an equilibrium point (stationary solution). You
control it by using the settings in the Value of variables not solved for and linearization
point area on the Initial Value tab of the Solver Manager dialog box. The inputs and the
outputs are deviations from the equilibrium values.
The export only considers certain solver parameters: the Matrix symmetry setting on
the General tab and the settings on the Advanced tab in the Solver parameters dialog
box.
S T A T E – S P A C E E X P O R T | 109
110 | C H A P T E
You can use model reduction to approximate the linearized model with a model that
has fewer degrees of freedom; select the Model reduction check box to enable it. The
model reduction uses eigenmodes and a static mode of the original PDE problem, and
the linearized model is projected onto the subspace of the eigenmodes. Use the
Number of eigenmodes edit field to control the number of eigenmodes and the Include
static modes check box to determine whether or not to use the static modes in the
state-space model.
The model reduction uses all settings in the Solver Parameters dialog box for the
eigenvalue solver. As an exception, the number of eigenmodes are taken from the
Export State-Space Model dialog box.
The Input/Output Variables
The state-space model needs input and output. Use the Variables tab to define the
input and output. For a multigeometry model, you control what geometry the inputs
and outputs belong to by opening the Export State-Space Model dialog box from that
geometry.
C R E A T I N G I N P U T V A R I A B L E S
Enter input variables in Input. Use a space-separated list of variable names. Use the
variables in expressions, coupling variables, or in the physics settings.
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
C R E A T I N G O U T P U T V A R I A B L E S
Define output variable in Output. Initially, there are no output variables. Click Add to
add an output variable. Type the output variable name in the Variable name edit field.
Click the Subdomain option button to evaluate the expression in any location by
providing the coordinates, or click the Point option button to evaluate the expression
in a geometry vertex (point). Then define the value of the variable name by select a
quantity from the Predefined quantities list or by typing an expression in the Expression
edit field
M O D I F Y I N G A N D D E L E T I N G V A R I A B L E S
You can modify a variable by selecting it in the Variable number list and then modifying
it. Similarly, to delete a variable, select it in the Variable number list and then click
Delete.
S T A T E – S P A C E E X P O R T | 111
112 | C H A P T E
R 2 : S I M U L I N K A N D S T A T E – S P A C E E X P O R T
I N D E X
_test 65
_time 66
1D geometry object 21
2D geometry object 20
3D geometry object 19
A a PDE coefficient 57
adaption 82
al PDE coefficient 57
analyzed geometry 17
exporting 18
importing 17
appl 73
application mode 74
in Model M-file 92
structure 73
application mode variable 78
application scalar variable 78
application scalar variables section
in Model M-file 94
application structure 73
Argyris element 49
assemble 68, 82
assigned variable names 78
associative geometry 88
B be PDE coefficient 57
bilinear interpolation 33
block library
in Simulink 100
BMP 25
bnd 40
border 65, 75
boundary coefficients 57, 77
boundary condition section
in Model M-file 95
boundary conditions
types of 76
boundary coupled equation variable 55
boundary coupled shape variable 54
C coefficient form 57
coerce geometry object 20
const 53
constants 53
constants section
in Model M-file 94
constr 66
constraints 56, 66, 102
contour 32
contour curves in images 25
contourc 32
Control System Toolbox 107
conv2 26, 35
coordinate system 55
coupling variable elements section
in Model M-file 96
cporder 69, 79
cross section 55
CUR 25
curve object 11
curve parameter 46
curve2 11
curve3 11
curved mesh elements 42
D da PDE coefficient 57
data structures 8
deactivated equations and boundary con-
ditions 64
Decomposed Geometry matrix 92
degree of freedom 72
delta function 66
dependent variables 56
I N D E X | 113
114 | I N D E X
difference 12, 15
dim 56
dnx geometric variable on boundary in
3D 45
dny geometric variable on boundary in
3D 45
dnz geometric variable on boundary in
3D 45
DOF 72
dom general geometric variable 44
domain groups 40
down direction 45
draw section
in Model M-file 91
dvol general geometric variable 44
dweak 65
dynamic model 100, 102
dynamic systems 100
E edg 40
edge elements 49
edge segment 18
edge settings section 95
in Model M-file 95
eigenvalue 52
elemdefault 79
element
default settings 79
equ 40
equation variable 54
equations 56
equilibrium solution 102
exporting
FEM structure 5
expr 53
expression 61
expression section
in Model M-file 95
expression variables 53, 78
extended FEM structure 72, 98
extended mesh 71
external mesh
importing 37
extrude 15
F f PDE coefficient 57
face 18
face segment 18
face3 11
FEM structure 3, 8, 38
exporting 5
importing 6
loading 6
saving 6
femeig 82
femlin 82
femnlin 82
femtime 82
field variables 52
fillet 15
filter2 26, 35
finite elements 48
flcontour2curve 32
flform 75
flim2curve 25, 26
flmesh2spline 26, 32
form 56
form section
in Model M-file 93
frame 41
for equation 44
for shape function 52
reference 41
frame 41
function definition section
in Model M-file 94
G g boundary coefficient 57
ga PDE coefficient 57
general dynamic model 102
general form 61
general model 101
general static model 102
geom 40
geomanalyze command 88
geomcoerce 14, 15
geomcomp 13
geometric entity 18
geometric variables 44
geometry 40
geometry function 23
geometry methods 23
Geometry M-file 17, 92
geometry model 16, 91
exporting 16
importing 16
geometry objects 11, 18, 21, 92
exporting 11
importing 11
in Model M-files 91
inheritance structure 18
geometry section
in Model M-file 91
geometry shape order 42
in Model M-file 94
geometry, mesh-based 36
geomspline 30
geomsurf 32
GIF 25
gporder 69, 79
H h boundary coefficient 57
h general geometric variable 44
HDF 25
Hermite element 49
I ICO 25
image
filtering 26
importing 25
reducing noise 26
image 25
image contour curves 25
image processing 25
Image Processing Toolbox 25
imagesc 25
imcontour 25, 32
imfinfo 25
importing
external mesh 37
FEM structure 6
MRI data 27
imread 25
imwrite 25
including default values
in FEM structure 5
ind 40
inheritance structure
for geometry classes 18, 20, 21
init 70
initial value 70
initial value section
in Model M-file 96
input variables 104
interior boundaries section
in Model M-file 93
intersection 13, 15
J JPEG 25
L Lagrange element 48
lambda 52
level 1 coefficient 58
level 2 coefficient 59
level 3 coefficient 60
level 4 coefficient 61
lib 55
linearization point 102, 109
I N D E X | 115
116 | I N D E X
linearized dynamic simulation 102
linearized model 101
linearized state-space model 107
linearized static model 103
loading a FEM structure 6
loft 25
lower subdomain 45
M material library section
in Model M-file 94
material variables 55
materials/coefficients library 55
MATLAB 100, 107
mesh 36, 47, 92
importing 36
mesh 47
mesh cases 47
mesh element scale factor 42
mesh object 36
importing as geometry 17
mesh section
in Model M-file 92
mesh, using as geometry 36
meshextend 71
meshextend section
Model M-file 96
meta variables 52
Model M-file 90, 95
application mode section 92
application scalar variables section 94
boundary condition section 95
constants section 94
coupling variable elements section 96
draw section 91
expression section 95
form section 93
function definition section 94
geometry section 91
geometry shape order section 94
initial value section 96
interior boundaries section 93
material library section 94
mesh section 92
meshextend section 96
multiphysics section 96
multiple geometry model 97
order of commands 90
PDE coefficient section 96
point settings section 95
postprocessing section 97
shape function section 95
simplify section 94
solver section 97
space dimension section 91
version section 90
Model M-files
geometry objects 91
model reduction 101, 102, 108, 110
moving mesh 41
MRI data
import of 27
importing 27
multiphysics 75
multiphysics section
Model M-file 96
multiple geometries 72
multiple geometry model
in Model M-file 97
N Navier’s equations 61, 79
node point 71
noise
reducing, in point data 32
normal vector 45
numerical quadrature 69
nx geometric variable on boundary in 3D
45
ny geometric variable on boundary in 3D
45
nz geometric variable on boundary in 3D
45
O output variables 105
P PCX 25
PDE coefficient
c PDE coefficient 57
PDE coefficient section
in Model M-file 96
PDE coefficients 57, 96
PDE Toolbox 92
PDE Toolbox decomposed geometry
matrix
importing 17
phase 52
PNG 25
pnt 40
point data
for creating geometry models 32
point settings section
in Model M-file 95
point source 66
point1 11
point2 11
point3 11
pointwise constraint 69
postprocessing function 80
postprocessing section
in Model M-file 97
primitive geometry objects 11
property/value pairs 82
Q q boundary coefficient 57
quadrature formula 69
R r boundary coefficient 57
rational Bézier curve
object 11
reference frame 41
revolve 15
S s geometric variable on boundary in 2D
46
s1 geometric variable on boundary in 3D
45
s2 geometric variable on boundary in 3D
45
saving a FEM structure 6
Sd general geometric variable 44
sdim 41
several geometries 72
S-function 101
shape 50, 79
shape function classes 48
shape function objects 48, 50
shape function section
in Model M-file 95
shape function variables 52
sharg_2_5 49
shbub 49
shcurl 49
shdisc 50
shdiv 49
shherm 48
shlag 48
shrinking coefficients
in FEM structure 6
simplify section
in Model M-file 94
Simulink 100
Simulink export 100
dynamic model 100, 102
examples 100
general dynamic model 102
general model 101
general static model 102
input variables 104
linearized dynamic simulation 102
I N D E X | 117
118 | I N D E X
linearized model 101
linearized static model 103
model reduction 101, 102
output variables 105
S-function 101
solver manager 103
solver parameters 103
static model 100, 102
slit 64
solid
ellipse 13
rectangle 12
solid object 11
solid1 11
solid2 11
solid3 11
solution form
coefficient 54, 55
general 54, 55
solution object 39
solver manager 103
solver parameters 103
solver section
in Model M-file 97
space coordinates 41
space dimension section
in Model M-file 91
special variables 52
sshape 41, 79
sshapedim 44
state-space export 107
model reduction 110
state-space form 102
static model 100, 102
subdomain 18
settings 75
surface parameters 45
T t 52
t1x geometric variable on boundary in
3D 45
t1y geometric variable on boundary in
3D 45
t1z geometric variable on boundary in
3D 45
t2x geometric variable on boundary in
3D 45
t2y geometric variable on boundary in
3D 45
t2z geometric variable on boundary in
3D 45
tangent vector 45
test function 65
TIFF 25
time 52
transfer matrix
in state-space model 107
trimmed surfaces 18
trimming solids 13
tx geometric variable on boundary in 2D
46
ty geometric variable on boundary in 2D
46
typographical conventions 1
U union 13, 14, 15
unx geometric variable on boundary in
3D 45
uny geometric variable on boundary in
3D 45
unz geometric variable on boundary in
3D 45
up direction 45
upper subdomain 45
V variable 52
application mode 77, 78
application scalar 77, 78
assigned name 78
expression 53, 77, 78
field 52
geometric 44
material 55
shape function 52
special 52
weak form meta 52
variables 53
vector element 49
version section
in Model M-file 90
vertex 18
W weak 65
weak form 65
weak form, generating 75
weights
of control polygon 14
X x space coordinate 41
xfem 72
xg geometric variable on boundary in 3D
45
xmesh 71
XWD 25
Y y space coordinate 41
yg geometric variable on boundary in 3D
45
Z z space coordinate 41
zg geometric variable on boundary 45
I N D E X | 119
120 | I N D E X
CONTENTS
Chapter 1: COMSOL Multiphysics Scripting
Introduction to Programming in COMSOL Multiphysics 3
Data Structure and Function Overview 8
Geometry Modeling 11
Meshing 36
Specifying a Model 38
The Structure of the Model M-file 90
Chapter 2: Simulink and State-Space Export
Simulink Export 100
State-Space Export 107
COMSOL Multiphysics Scripting
Typographical Conventions
Introduction to Programming in COMSOL Multiphysics
Running COMSOL Multiphysics from COMSOL Script or MATLAB
A First Example-Poisson’s Equation on the Unit Disk
Exporting and Importing the FEM Structure
Saving and Loading a FEM Structure
Creating the FEM Structure Using a Model M-file
Data Structure and Function Overview
FEM Structure Overview
Data Structures and Functions
Geometry Modeling
Working with Geometry Objects
Creating a 1D Geometry
Creating a 2D Geometry Using Solid Modeling
Creating a 2D Geometry Using Boundary Modeling
Creating 3D Geometries Using Solid Modeling
Working with a Geometry Model
Importing and Exporting Geometry and CAD Models from File.
Working with the Analyzed Geometry
The Geometry Object Hierarchy
Command-Line Geometry Modeling
Working with Image Files and MRI Data
Spline Interpolation for Creating Geometry Objects
Surface Interpolation
Meshing
Specifying a Model
The Geometry and the FEM Problem
Mesh
Shape Functions
Variables
Materials/Coefficients Libraries
Equations and Constraints
Discretization
Initial Values
The Extended Mesh
Multiple Geometries
The Application Structure
Overview of Postprocessing Functions
Overview of Common Solver Properties and Their Values
Associative Geometry and Parametrization
The Structure of the Model M-file
Model M-files for Models With Multiple Geometries
Simulink and State-Space Export
Simulink Export
Dynamic or Static Export?
General or Linearized Export?
The Simulink Export Dialog Box
The Simulink Export Types
Relation to Solver Parameters and Solver Manager Settings
The Input/Output Variables
Using the COMSOL Multiphysics Model in Simulink
State-Space Export
The State-Space Form
Model Reduction
The State-Space Export Dialog Box
The Input/Output Variables
INDEX
A
B
C
D
E
F
G
H
I
J
L
M
N
O
P
Q
R
S
T
U
V
W
X
Y
Z