1 function inner_product_matrices = evol_TB_2D_weak_form_inner_product_matrices(model,model_data)
2 %
function inner_product_matrices = evol_TB_2D_weak_form_inner_product_matrices(model,model_data)
4 % This
function computes the L2- and H1-inner product matrices
for the
5 % (COMSOL-) evol_TB_2D_weak_form model. This is done by setting the weak expression
6 % to u*u (L2) or u*u + ux*ux + uy*uy (H1).
8 % Input: model, model_data
9 % output: inner_product_matrices:
struct containing the full and eliminated
10 % versions of the L2- and H1- inner product matrices.
12 % Oliver Zeeb, 2013/11/20
14 comsol_model=model_data.comsol_model;
16 %extract the expression in wfeq1 so thaht we can set it back to
this value at the end of
this function
17 w_wfeq1 = comsol_model.physics(
'w').feature(
'wfeq1').getString(
'weak');
19 %deactivate the unneeded weak expressions
20 comsol_model.physics(
'w').feature(
'wfeq2').active(
false);
21 comsol_model.physics(
'w').feature(
'wfeq3').active(
false);
22 comsol_model.physics(
'w').feature(
'wfeq4').active(
false);
23 comsol_model.physics(
'w').feature(
'wfeq5').active(
false);
24 comsol_model.physics(
'w').feature(
'weak1').active(
false);
25 comsol_model.physics(
'w').feature(
'init1').set(
'u',
'0');
27 %only the stiffness matrix is needed --> deactivate assembling of the other
29 comsol_model.sol(
'sol1').feature(
'v1').set(
'scaleval',
'1');
30 comsol_model.sol(
'sol1').feature(
'v1').set(
'scalemethod',
'manual');
31 comsol_model.sol(
'sol1').feature(
't1').active(
false);
32 comsol_model.sol(
'sol1').feature(
'a1').set(
'Nullf',
false);
33 comsol_model.sol(
'sol1').feature(
'a1').set(
'Null',
false);
34 comsol_model.sol(
'sol1').feature(
'a1').set(
'Lc',
false);
35 comsol_model.sol(
'sol1').feature(
'a1').set(
'ud',
false);
36 comsol_model.sol(
'sol1').feature(
'a1').set(
'Kc',
false);
37 comsol_model.sol(
'sol1').feature(
'a1').set(
'Ec',
false);
38 comsol_model.sol(
'sol1').feature(
'a1').set(
'uscale',
false);
39 comsol_model.sol(
'sol1').feature(
'a1').set(
'Dc',
false);
40 comsol_model.sol(
'sol1').feature(
'a1').set(
'K',
true);
42 nr_dofs_full = size(model_data.ind_vectors.Null,1);
45 %%%%%%%%%%%%%%%%%%%%%%%%
47 %%%%%%%%%%%%%%%%%%%%%%%%
48 comsol_model.physics(
'w').feature(
'wfeq1').set(
'weak',
'-u*test(u)');
49 %run study --> assemble L2 matrices
50 comsol_model.study(
'std1').run;
51 % extract L2 (full version)
52 str.KCol = 1+double(comsol_model.sol(
'sol1').feature(
'a1').getSparseMatrixCol(
'K'));
53 str.KRow = 1+double(comsol_model.sol(
'sol1').feature(
'a1').getSparseMatrixRow(
'K'));
54 str.KVal = comsol_model.sol(
'sol1').feature(
'a1').getSparseMatrixVal(
'K');
55 str.K = sparse(nr_dofs_full,nr_dofs_full);
56 for k=1:length(str.KRow)
57 str.K(str.KRow(k), str.KCol(k)) = str.KVal(k);
59 inner_product_matrices.L2 = str.K;
60 inner_product_matrices.L2_eliminated=inner_product_matrices.L2(model_data.ind_vectors.com_DOF_ind,model_data.ind_vectors.com_DOF_ind);
65 %%%%%%%%%%%%%%%%%%%%%%%%
67 %%%%%%%%%%%%%%%%%%%%%%%%
68 comsol_model.physics('w').feature('wfeq1').set('weak', '-u*test(u) - ux*test(ux) - uy*test(uy)');
69 %run study --> assemble H1 matrices
70 comsol_model.study('std1').run;
71 % extract H1 (full version)
72 str.KCol = 1+
double(comsol_model.sol('sol1').feature('a1').getSparseMatrixCol('K'));
73 str.KRow = 1+
double(comsol_model.sol('sol1').feature('a1').getSparseMatrixRow('K'));
74 str.KVal = comsol_model.sol('sol1').feature('a1').getSparseMatrixVal('K');
75 str.K = sparse(nr_dofs_full,nr_dofs_full);
76 for k=1:length(str.KRow)
77 str.K(str.KRow(k), str.KCol(k)) = str.KVal(k);
79 inner_product_matrices.H1 = str.K;
80 inner_product_matrices.H1_eliminated=inner_product_matrices.H1(model_data.ind_vectors.com_DOF_ind,model_data.ind_vectors.com_DOF_ind);
84 %TEST: Check whether the H1 and L2 Matrices are ok:
85 % comsol_model.physics('w').feature('init1').set('u', 'x+y');
86 % comsol_model.study('std1').run;
87 % U_init = mphgetu(comsol_model);
88 % L2_norm_sqr = U_init' * inner_product_matrices.L2 * U_init % L2-Norm zum quadrat --> should be 1.1667
89 % H1_norm_sqr = U_init' * inner_product_matrices.H1 * U_init % H1 norm zum quadrat --> should be 3.1667
91 % comsol_model.physics('w').feature('init1').set('u', '1');
92 % comsol_model.study('std1').run;
93 % U_init = mphgetu(comsol_model);
94 % L2_norm_sqr = U_init' * inner_product_matrices.L2 * U_init % L2-Norm zum quadrat --> should be 1
95 % H1_norm_sqr = U_init' * inner_product_matrices.H1 * U_init % H1 norm zum quadrat --> should be
99 % activate the weak forms again and set back the original values:
100 comsol_model.physics('w').feature('wfeq1').set('weak', w_wfeq1);
101 comsol_model.physics('w').feature('wfeq2').active(true);
102 comsol_model.physics('w').feature('wfeq3').active(true);
103 comsol_model.physics('w').feature('wfeq4').active(true);
104 comsol_model.physics('w').feature('wfeq5').active(true);
105 comsol_model.physics('w').feature('weak1').active(true);
107 % activate the solver again
108 comsol_model.sol('sol1').feature('a1').active(false);
109 comsol_model.sol('sol1').feature('t1').active(true);
1.8.8
