OptFEM3DP1 Toolbox  20130614_170440
Matlab/Octave Optimized P1-Lagrange Finite Element Method in 3D
 All Files Functions Variables Pages
validMassVF3DP1.m
Go to the documentation of this file.
1 function validMassVF3DP1(Num)
2 % Copyright:
3 % See \ref license
4 
5  disp('*******************************************')
6  disp('* MassVF 3D Assembling P1 validations *')
7  disp('*******************************************')
8  switch Num
9  case 0
10  s=sprintf('Global alternate numbering / local alternate numbering');
11  case 1
12  s=sprintf('Global block numbering / local alternate numbering');
13  case 2
14  s=sprintf('Global alternate numbering / local block numbering');
15  case 3
16  s=sprintf('Global block numbering / local block numbering');
17  otherwise
18  error('invalid Num value')
19  end
20  fprintf(' Numbering Choice : %s\n',s);
21 
22 
23 
24  Th=CubeMesh(10);
25 
26 % TEST 1
27  disp('----------------------------------StiffElasAssemblingP1OptV2-------')
28  disp(' Test 1: Matrices errors and CPU times ')
29  fprintf(' Matrix size : %d\n',3*Th.nq);
30  disp('-----------------------------------------')
31  tic();
32  Mbase=MassVFAssembling3DP1base(Th.nq,Th.nme,Th.me,Th.volumes,Num);
33  T(1)=toc();
34  tic();
35  MOptV0=MassVFAssembling3DP1OptV0(Th.nq,Th.nme,Th.me,Th.volumes,Num);
36  T(2)=toc();
37  Test1.error(1)=norm(Mbase-MOptV0,Inf);
38  Test1.name{1}='MassVFAssembling3DP1OptV0';
39  fprintf(' Error P1base vs OptV0 : %e\n',Test1.error(1))
40  tic();
41  MOptV1=MassVFAssembling3DP1OptV1(Th.nq,Th.nme,Th.me,Th.volumes,Num);
42  T(3)=toc();
43  Test1.error(2)=norm(Mbase-MOptV1,Inf);
44  Test1.name{2}='MassVFAssembling3DP1OptV1';
45  fprintf(' Error P1base vs OptV1 : %e\n',Test1.error(2))
46  tic();
47  MOptV2=MassVFAssembling3DP1OptV2(Th.nq,Th.nme,Th.me,Th.volumes,Num);
48  T(4)=toc();
49  Test1.error(3)=norm(Mbase-MOptV2,Inf);
50  Test1.name{3}='MassVFAssembling3DP1OptV2';
51  fprintf(' Error P1base vs OptV2 : %e\n',Test1.error(3))
52  MOptV3=MassVFAssembling3DP1OptV3(Th.nq,Th.nme,Th.me,Th.volumes,Num);
53  T(5)=toc();
54  Test1.error(4)=norm(Mbase-MOptV3,Inf);
55  Test1.name{4}='MassVFAssembling3DP1OptV3';
56  fprintf(' Error P1base vs OptV3 : %e\n',Test1.error(4))
57 
58  fprintf(' CPU times base (ref) : %3.4f (s)\n',T(1))
59  fprintf(' CPU times OptV0 : %3.4f (s) - Speed Up X%3.3f\n',T(2),T(1)/T(2))
60  fprintf(' CPU times OptV1 : %3.4f (s) - Speed Up X%3.3f\n',T(3),T(1)/T(3))
61  fprintf(' CPU times OptV2 : %3.4f (s) - Speed Up X%3.3f\n',T(4),T(1)/T(4))
62  fprintf(' CPU times OptV3 : %3.4f (s) - Speed Up X%3.3f\n',T(5),T(1)/T(5))
63  checkTest1(Test1)
64 
65 
66 
67 
68 % TEST 2
69  disp('-----------------------------------------------------------')
70  disp(' Test 2: Validations by integration on [0,1]x[0,1]x[0,1] ')
71  disp('-----------------------------------------------------------')
72  Test=valid_MECA();
73  qf=Th.q;
74  for kk=1:length(Test)
75  KK=MassVFAssembling3DP1OptV2(Th.nq,Th.nme,Th.me,Th.volumes,Num);
76  u=Test(kk).u;
77  v=Test(kk).v;
78  switch Num
79  case {0,2}
80  U=[u{1}(qf(1,:),qf(2,:),qf(3,:));u{2}(qf(1,:),qf(2,:),qf(3,:));u{3}(qf(1,:),qf(2,:),qf(3,:))];
81  U=U(:);
82  V=[v{1}(qf(1,:),qf(2,:),qf(3,:));v{2}(qf(1,:),qf(2,:),qf(3,:));v{3}(qf(1,:),qf(2,:),qf(3,:))];
83  V=V(:);
84  case {1,3}
85  U=[u{1}(qf(1,:),qf(2,:),qf(3,:)) u{2}(qf(1,:),qf(2,:)) u{3}(qf(1,:),qf(2,:),qf(3,:))]';
86  V=[v{1}(qf(1,:),qf(2,:),qf(3,:)) v{2}(qf(1,:),qf(2,:),qf(3,:)) v{3}(qf(1,:),qf(2,:),qf(3,:))]';
87  end
88  %whos
89  Test(kk).error=abs(Test(kk).Mass-U'*KK*V);
90  fprintf(' functions %d : u(x,y,z)=(%s,%s,%s), v(x,y,z)=(%s,%s,%s),\n -> MassVF error=%e\n', ...
91  kk,Test(kk).cu{1},Test(kk).cu{2},Test(kk).cu{3},Test(kk).cv{1},Test(kk).cv{2},Test(kk).cv{3},Test(kk).error);
92  end
93  checkTest2(Test)
94 
95 % TEST 3
96  disp('--------------------------------')
97  disp(' Test 3: Validations by order ')
98  disp('--------------------------------')
99  n=length(Test);
100  u=Test(n).u;
101  v=Test(n).v;
102  ExSol=Test(n).Mass;
103  N=10;
104  for k=1:N
105  Th=CubeMesh(5*k+5);
106  qf=Th.q;
107  fprintf('(%d/%d) Matrix size : %d\n',k,N,Th.nq);
108  h(k)=GetMaxLengthEdges(Th.q,Th.me);
109  tic();
110  M=MassVFAssembling3DP1OptV2(Th.nq,Th.nme,Th.me,Th.volumes,Num);
111  TT(k)=toc();
112  switch Num
113  case {0,2}
114  U=[u{1}(qf(1,:),qf(2,:),qf(3,:));u{2}(qf(1,:),qf(2,:),qf(3,:));u{3}(qf(1,:),qf(2,:),qf(3,:))];
115  U=U(:);
116  V=[v{1}(qf(1,:),qf(2,:),qf(3,:));v{2}(qf(1,:),qf(2,:),qf(3,:));v{3}(qf(1,:),qf(2,:),qf(3,:))];
117  V=V(:);
118  case {1,3}
119  U=[u{1}(qf(1,:),qf(2,:),qf(3,:)) u{2}(qf(1,:),qf(2,:)) u{3}(qf(1,:),qf(2,:),qf(3,:))]';
120  V=[v{1}(qf(1,:),qf(2,:),qf(3,:)) v{2}(qf(1,:),qf(2,:),qf(3,:)) v{3}(qf(1,:),qf(2,:),qf(3,:))]';
121  end
122  Error(k)=abs(ExSol-U'*M*V);
123  fprintf(' MassVFAssembling3DP1OptV2 CPU times : %3.3f(s)\n',TT(k));
124  fprintf(' Error : %e\n',Error(k));
125  end
126 
127  loglog(h,Error,'+-r',h,h*1.1*Error(1)/h(1),'-sm',h,1.1*Error(1)*(h/h(1)).^2,'-db')
128  legend('Error','O(h)','O(h^2)')
129  xlabel('h')
130  title(sprintf('Test 3 : 3D MassVF Matrix (Num=%d)',Num))
131  checkTest3(h,Error)
132 
133  %figure(3)
134  %spy(M)
135  %title(sprintf('Matrix sparsity for %s numbering',s))
136 end
137 
138 function checkTest1(Test)
139  I=find(Test.error>1e-14);
140  if isempty(I)
141  disp('------------------------')
142  disp(' Test 1 (results): OK')
143  disp('------------------------')
144  else
145  disp('----------------------------')
146  disp(' Test 1 (results): FAILED')
147  disp('----------------------------')
148  end
149 end
150 
151 function checkTest2(Test)
152  N=length(Test);
153  cntFalse=0;
154  for k=1:N
155  if (ismember(Test(k).degree,[0 1]))
156  if (Test(k).error>1e-12)
157  cntFalse=cntFalse+1;
158  end
159  end
160  end
161  if (cntFalse==0)
162  disp('------------------------')
163  disp(' Test 2 (results): OK')
164  disp('------------------------')
165  else
166  disp('----------------------------')
167  disp(' Test 2 (results): FAILED')
168  disp('----------------------------')
169  end
170 end
171 
172 function checkTest3(h,error)
173  % order 2
174  P=polyfit(log(h),log(error),1);
175  if abs(P(1)-2)<5e-2
176  disp('------------------------')
177  disp(' Test 3 (results): OK')
178  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
179  disp('------------------------')
180  else
181  disp('----------------------------')
182  disp(' Test 3 (results): FAILED')
183  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
184  disp('----------------------------')
185  end
186 end