程序代写代做代考 assembly data structure gui algorithm matlab mlinterface_pdf.book

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


—— 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+=⎩




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