OptFEM2DP1  1.0
Matlab/Octave Optimized P1-Lagrange Finite Element Method in 2D
 All Files Functions Variables Groups Pages
validMassWP1.m
Go to the documentation of this file.
1 function validMassWP1()
2 % function validMassWP1()
3 % Validation function for the assembly of the weighted mass matrix
4 % for `P_1`-Lagrange finite elements
5 %
6 % The Weighted Mass Matrix `\MasseF{w}` is given by
7 % ``\MasseF{w}_{i,j}=\int_\DOMH w(\q)\FoncBase_i(\q) \FoncBase_j(\q) d\q,\ \forall (i,j)\in\ENS{1}{\nq}^2``
8 % where `\FoncBase_i` are `P_1`-Lagrange basis functions.
9 % This Matrix is computed by functions MassWAssemblingP1{Version} where {Version} is one of
10 % 'base', 'OptV0', 'OptV1' and 'OptV2'.
11 % - Test 1: Computation of the MassW Matrix using all the versions giving errors and cputimes
12 % - Test 2: Computation of the integral ``\int_\DOM w(x,y) u(x,y) v(x,y) dxdy \approx \DOT{\MasseF{w} \vecb{U}}{\vecb{V}}``
13 % where `\vecb{U}_i=u(\q^i)` and `\vecb{V}_i=v(\q^i)`.
14 % Functions `u`, `v` and `w` are those defined in #valid_FEMmatrices.
15 % - Test 3: One retrieves the order 2 of `P_1`-Lagrange integration ``|\int_\DOM u\,v\,w -\Pi_h(u)\,\Pi_h(v)\,\Pi_h(w)d\DOM| \leq C h^2``
16 %
17 % See also:
18 % #MassWAssemblingP1base, #MassWAssemblingP1OptV0,
19 % #MassWAssemblingP1OptV1, #MassWAssemblingP1OptV2,
20 % #valid_FEMmatrices, #SquareMesh, #GetMaxLengthEdges
21 %
22 % @author François Cuvelier @date 2012-11-26
23 %
24 % OptFEM2DP1 [V1.0e] - Copyright (C) 2013 CJS (LAGA)
25 %
26 % This file is part of OptFEM2DP1.
27 % OptFEM2DP1 is free software: you can redistribute it and/or modify
28 % it under the terms of the GNU General Public License as published by
29 % the Free Software Foundation, either version 3 of the License, or
30 % (at your option) any later version.
31 %
32 % OptFEM2DP1 is distributed in the hope that it will be useful,
33 % but WITHOUT ANY WARRANTY; without even the implied warranty of
34 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
35 % GNU General Public License for more details.
36 %
37 % You should have received a copy of the GNU General Public License
38 % along with this program. If not, see <http://www.gnu.org/licenses/>.
39 
40  disp('*******************************************')
41  disp('* MassW Assembling P1 validations *')
42  disp('*******************************************')
43 
44  Th=SquareMesh(50);
45 
46 
47 % TEST 1
48  disp('-----------------------------------------')
49  disp(' Test 1: Matrices errors and CPU times ')
50  disp('-----------------------------------------')
51  w=@(x,y) cos(x+y);
52  Tw=w(Th.q(1,:),Th.q(2,:));
53  tic();
54  Mbase=MassWAssemblingP1base(Th.nq,Th.nme,Th.me,Th.areas,Tw);
55  T(1)=toc();
56  tic();
57  MOptV0=MassWAssemblingP1OptV0(Th.nq,Th.nme,Th.me,Th.areas,Tw);
58  T(2)=toc();
59  Test1.error(1)=norm(Mbase-MOptV0,Inf);
60  Test1.name{1}='MassWAssemblingP1OptV0';
61  fprintf(' Error P1base vs OptV0 : %e\n',Test1.error(1))
62  tic();
63  MOptV1=MassWAssemblingP1OptV1(Th.nq,Th.nme,Th.me,Th.areas,Tw);
64  T(3)=toc();
65  Test1.error(2)=norm(Mbase-MOptV1,Inf);
66  Test1.name{2}='MassWAssemblingP1OptV1';
67  fprintf(' Error P1base vs OptV1 : %e\n',Test1.error(2))
68  tic();
69  MOptV2=MassWAssemblingP1OptV2(Th.nq,Th.nme,Th.me,Th.areas,Tw);
70  T(4)=toc();
71  Test1.error(3)=norm(Mbase-MOptV2,Inf);
72  Test1.name{3}='MassWAssemblingP1OptV2';
73  fprintf(' Error P1base vs OptV2 : %e\n',Test1.error(3))
74 
75  fprintf(' CPU times base (ref) : %3.4f (s)\n',T(1))
76  fprintf(' CPU times OptV0 : %3.4f (s) - Speed Up X%3.3f\n',T(2),T(1)/T(2))
77  fprintf(' CPU times OptV1 : %3.4f (s) - Speed Up X%3.3f\n',T(3),T(1)/T(3))
78  fprintf(' CPU times OptV2 : %3.4f (s) - Speed Up X%3.3f\n',T(4),T(1)/T(4))
79  checkTest1(Test1)
80 
81 % TEST 2
82  disp('-----------------------------------------------------')
83  disp(' Test 2: Validations by integration on [0,1]x[0,1] ')
84  disp('-----------------------------------------------------')
85  Test=valid_FEMmatrices();
86  for kk=1:length(Test)
87  U=Test(kk).u(Th.q(1,:),Th.q(2,:));
88  V=Test(kk).v(Th.q(1,:),Th.q(2,:));
89  Tw=Test(kk).w(Th.q(1,:),Th.q(2,:));
90  M=MassWAssemblingP1OptV2(Th.nq,Th.nme,Th.me,Th.areas,Tw);
91  Test(kk).error=abs(Test(kk).MassW-U*M*V');
92  fprintf(' function %d : u(x,y)=%s, v(x,y)=%s, w(x,y)=%s\n -> MassW error=%e\n', ...
93  kk,Test(kk).cu,Test(kk).cv,Test(kk).cw,Test(kk).error);
94  end
95  checkTest2(Test)
96 
97 % TEST 3
98  disp('--------------------------------')
99  disp(' Test 3: Validations by order ')
100  disp('--------------------------------')
101  n=length(Test);
102  u=Test(n).u;
103  v=Test(n).v;
104  w=Test(n).w;
105  ExSol=Test(n).MassW;
106 
107  for k=1:10
108  Th=SquareMesh(50*k+50);
109  Tw=w(Th.q(1,:),Th.q(2,:));
110  fprintf(' Matrix size : %d\n',Th.nq);
111  h(k)=GetMaxLengthEdges(Th.q,Th.me);
112  tic();
113  M=MassWAssemblingP1OptV2(Th.nq,Th.nme,Th.me,Th.areas,Tw);
114  TT(k)=toc();
115  U=u(Th.q(1,:),Th.q(2,:));
116  V=v(Th.q(1,:),Th.q(2,:));
117  Error(k)=abs(ExSol-U*M*V');
118  fprintf(' MassWAssemblingP1OptV2 CPU times : %3.3f(s)\n',TT(k));
119  fprintf(' Error : %e\n',Error(k));
120  end
121 
122  loglog(h,Error,'+-r',h,h*1.1*Error(1)/h(1),'-sm',h,1.1*Error(1)*(h/h(1)).^2,'-db')
123  legend('Error','O(h)','O(h^2)')
124  xlabel('h')
125  title('Test 3 : MassW Matrix')
126  checkTest3(h,Error)
127 end
128 
129 function checkTest1(Test)
130  I=find(Test.error>1e-14);
131  if isempty(I)
132  disp('------------------------')
133  disp(' Test 1 (results): OK')
134  disp('------------------------')
135  else
136  disp('----------------------------')
137  disp(' Test 1 (results): FAILED')
138  disp('----------------------------')
139  end
140 end
141 
142 function checkTest2(Test)
143  N=length(Test);
144  cntFalse=0;
145  for k=1:N
146  if (Test(k).degree<=1)
147  if (Test(k).error>1e-14)
148  cntFalse=cntFalse+1;
149  end
150  end
151  end
152  if (cntFalse==0)
153  disp('------------------------')
154  disp(' Test 2 (results): OK')
155  disp('------------------------')
156  else
157  disp('----------------------------')
158  disp(' Test 2 (results): FAILED')
159  disp('----------------------------')
160  end
161 end
162 
163 function checkTest3(h,error)
164  % order 2
165  P=polyfit(log(h),log(error),1);
166  if abs(P(1)-2)<1e-2
167  disp('------------------------')
168  disp(' Test 3 (results): OK')
169  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
170  disp('------------------------')
171  else
172  disp('----------------------------')
173  disp(' Test 3 (results): FAILED')
174  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
175  disp('----------------------------')
176  end
177 end