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