COMSOL 6.4 - Set System Matrices in the Model

Use the function mphinputmatrix to set a linear matrix system to a model:

mphinputmatrix(model,<str>,<soltag>,<soltypetag>)

This command sets the matrices of a linear system stored in the MATLAB® structure <str> into the model. The linear system is associated to the solver sequence <soltag> and is to be solved by the solver <soltypetag>.

mphinputmatrix only supports the solver types Stationary, Eigenvalue, and Time.

A valid structure <str> for a stationary solver includes the following fields:


Field	Description
K	Stiffness matrix
L	Load vector
M	Constraint vector
N	Constraint Jacobian

A valid structure <str> for a time-dependent or an eigenvalue solver includes the following fields:


Expression	Description
K	Stiffness matrix
L	Load vector
M	Constraint vector
N	Constraint Jacobian
D	Damping matrix
E	Mass matrix

You can also include the Constraint force Jacobian vector, defined in the field NF.

Once the linear system is loaded in the model, you can directly run the solver.


	The system matrices are not stored in the model when it is saved in the MPH-file format or loaded to the COMSOL Desktop.

Setting a Model with a Modified Matrix System

This example deals with heat transfer in solids physics. The geometry and physics settings are already set in the model and saved in the MPH-file format. This Model MPH file comes with the COMSOL installation.

At the MATLAB prompt you load the model and add an additional line heat source to the model directly in the system matrix by manually changing the load vector. Then compute the solution of the modified system in COMSOL.

Load the base Model MPH file and display the geometry:

model = mphopen('model_tutorial_llmatlab.mph');

mphgeom(model)

This results in the following MATLAB figure:

Draw the line to be used as a line heat source in the model and plot the modified geometry:

comp1 = model.component('comp1');

b1 = comp1.geom('geom1').feature.create('b1', 'BezierPolygon');

b1.set('p', {'1e-2' '5e-2'; '1e-2' '5e-2'; '1e-2' '1e-2'});

mphgeom(model,'geom1','edgelabels','on','facealpha',0.5);

In the figure below you can see that the added line as the index 21:

Generate a mesh with finer mesh settings:

mesh1 = comp1.mesh('mesh1');

mesh1.feature.create('ftet1', 'FreeTet');

mesh1.feature('size').set('hauto', 3);

mesh1.run;

mphmesh(model)

Set the solver sequence associated to a stationary study node:

std1 = model.study.create('std1');

std1.feature.create('stat', 'Stationary');

sol1 = model.sol.create('sol1');

sol1.study('std1');

st1 = sol1.feature.create('st1', 'StudyStep');

st1.set('studystep', 'stat');

v1 = sol1.feature.create('v1', 'Variables');

v1.set('control', 'stat');

sol1.feature.create('s1', 'Stationary');

Set the dependent variable discretization with linear shape function:

Shape = comp1.physics('ht').prop('ShapeProperty');

Shape.set('order_temperature', 1, 1);

The heat transfer interface automatically computes for internal DOFs in order to evaluate fluxes accurately at the boundaries. Deactivate the internal DOFs with this command:

Shape.set('boundaryFlux_temperature', false);

Now extract the matrices of the linear system associated to the solver sequence sol1:

ME = mphmatrix(model,'sol1','Out',{'K' 'L' 'M' 'N'},...

'initmethod','sol','initsol','zero');

To retrieve the degrees of freedom that belong to edge 21, you need to get the geometric mesh data:

[stats,data] = mphmeshstats(model);

With the mesh data structure data, you can get the element indices that belong to edge 2. Use the MATLAB find function to list all the indices:

elem_idx = find(data.elementity{1}==21)'

With the function mphxmeshinfo, retrieve the finite element mesh information associated to solver sequence sol1:

info = mphxmeshinfo(model,'soltag','sol1','studysteptag','v1');

In the info structure you can get the DOFs indices that belong to the edge element defined with the indices elem_idx:

dofs = info.elements.edg.dofs;

edgdofs_idx = [];

for i = 1:length(elem_idx)

edgdofs_idx = [edgdofs_idx; dofs(:,elem_idx(i))];

end

edgdofs_idx might contain duplicate DOFs indices. This is because the information is from the element level; the duplicate indices correspond to the connecting node between two adjacent elements.

First remove the duplicate entities:

unique_idx = unique(edgdofs_idx);

Edit the load vector for the DOF that belong to edge 21. The total applied power is 50 W:

ME.L(unique_idx+1) = 50/length(unique_idx);

Now that the linear system has been modified, set it back in the model:

mphinputmatrix(model,ME,'sol1','s1')

mphmatrix only assembles the matrix system for the DOFs solved in the specified solver configuration. mphinputmatrix inserts the matrix system as defined by the user. When inserting matrices in an existing model, the solution format may not be compatible with the inserted system matrices.

In order to have a compatible xmesh solution format compatible with the size of the inserted matrices, add a new equation-form physics interface, solving only for one variable.

gForm = comp1.physics.create('g', 'GeneralFormPDE', {'u'});

gForm.prop('ShapeProperty').set('order', 1);

gForm.prop('ShapeProperty').set('boundaryFlux', false);

Disable the Heat Transfer interface.

comp1.physics('ht').active(false);

Compute the solution of the added system:

model.sol('sol1').runAll;

Display the solution:

pg1 = model.result.create('pg1', 'PlotGroup3D');

pg1.feature.create('surf1', 'Surface');

mphplot(model,'pg1','rangenum',1)