Discussion Forum

I did not find any help on this topic on the web, so I thought I should publish my simple results:

When working with FEM-models it is necessary to check grid-convergence of the solution. The L2-norm (and the H1 too) is a good measure for this, and it is not difficult to make a matalb-script to find the norm(s). Thanks to Bertil Nistad @ Comsol support for invaluable help making my scripts/functions.

First define a function for calculating the norm:

function norm=femL2H1norm(femA,femB,l2expr,h1expr)
% femA and femB are the fem structures we want to compare, femB has the finest mesh (most elements), so we will compare the results om femB's mesh.
% l2expr and h1expr are the expressions evaluated; typically 'u^2' and 'ux^2+uy^2', respectively. Note the string-notation. Replace u with your solution-variable in this string.

% initializing:
norm = [0 0] % if the returned value of the function is zero then something is wrong...
femA.xmesh=meshextend(femA, ...
'mcase',[1]); % This might not be necessary, but it makes the method more general.
femB.xmesh=meshextend(femB, ...
'mcase',[1]);
% Start mapping fem to current mesh:
% New mesh:
femA.mesh = femB.mesh;
% Mapping solution A to extended mesh
initA = asseminit(femB,'init',femA.sol,'xmesh',femA.xmesh,'blocksize','auto');
% Update models
femA.sol = initA;
% Multiphysics
femA = multiphysics(femA);
% Update extended mesh
femA.xmesh=femB.xmesh;
% update femA.sol to contain the *difference*:
init = asseminit(femA, 'init', femA.sol.u - femB.sol.u);
femA.sol = init;

% Compute the norms:
norm(1) = sqrt(postint(femA, l2expr));
norm(2)= norm(1) + sqrt(postint(femA, h1expr));

end

This function can be expanded and modified a lot. Read more about postint in the Reference guide. You could also use Matlabs nargin- and nargout-functions to create default values.
However, we also might want to know the number of elements in our mesh, to make it possible to plot the norm versus number of elements in mesh. This is easily done in the script where we call femL2H2norm(). The script should be based on the model .m-file. After the mesh-initializing or refinement add the lines:

% Initialize mesh
fem.mesh=meshinit(fem, ...
'hauto',4);
% initialize counter:
numMeshElem = 0;
% Number of elements in mesh:
el = get(fem.mesh,'el');
M = el{3}.elem;
numMeshElem(end) = size(M,2);
clear M el


This satisfied my need for generating some convergence plots for my models.