OptFEM3DP1 Toolbox  V1.0
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();
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();
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  fprintf(' functions %d : u(x,y,z)=(%s,%s,%s), v(x,y,z)=(%s,%s,%s),\n',n, ...
100  Test(n).cu{1},Test(n).cu{2},Test(n).cu{3}, ...
101  Test(n).cv{1},Test(n).cv{2},Test(n).cv{3});
102  fprintf(' Lame coefficients : lambda = %f, mu =%f\n',lambda,mu);
103  fprintf(' Integral of epsilon(u)^t C epsilon(v) over unit cube : %.16f\n',ExSol);
104  N=10;
105  for k=1:N
106  Th=CubeMesh(5*k+10);
107  qf=Th.q;
108  fprintf('(%d/%d) Matrix size : %d\n',k,N,Th.nq);
109  h(k)=GetMaxLengthEdges(Th.q,Th.me);
110  tic();
111  M=StiffElasAssembling3DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,lambda,mu,Num);
112  TT(k)=toc();
113  switch Num
114  case {0,2}
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  U=U(:);
117  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,:))];
118  V=V(:);
119  case {1,3}
120  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,:))]';
121  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,:))]';
122  end
123  Error(k)=abs(ExSol-U'*M*V);
124  fprintf(' StiffElasAssemblingP1OptV2 CPU times : %3.3f(s)\n',TT(k));
125  fprintf(' Error : %e\n',Error(k));
126  end
127 
128  loglog(h,Error,'+-r',h,h*1.1*Error(1)/h(1),'-sm',h,1.1*Error(1)*(h/h(1)).^2,'-db')
129  legend('Error','O(h)','O(h^2)')
130  xlabel('h')
131  title(sprintf('Test 3 : Stiff Elas. Matrix (Num=%d)',Num))
132  checkTest3(h,Error)
133 
134  %figure(3)
135  %spy(M)
136  %title(sprintf('Matrix sparsity for %s numbering',s))
137 end
138 
139 function checkTest1(Test)
140  I=find(Test.error>1e-14);
141  if isempty(I)
142  disp('------------------------')
143  disp(' Test 1 (results): OK')
144  disp('------------------------')
145  else
146  disp('----------------------------')
147  disp(' Test 1 (results): FAILED')
148  disp('----------------------------')
149  end
150 end
151 
152 function checkTest2(Test)
153  N=length(Test);
154  cntFalse=0;
155  for k=1:N
156  if (ismember(Test(k).degree,[0 1]))
157  if (Test(k).error>1e-12)
158  cntFalse=cntFalse+1;
159  end
160  end
161  end
162  if (cntFalse==0)
163  disp('------------------------')
164  disp(' Test 2 (results): OK')
165  disp('------------------------')
166  else
167  disp('----------------------------')
168  disp(' Test 2 (results): FAILED')
169  disp('----------------------------')
170  end
171 end
172 
173 function checkTest3(h,error)
174  % order 2
175  P=polyfit(log(h),log(error),1);
176  if abs(P(1)-2)<5e-2
177  disp('------------------------')
178  disp(' Test 3 (results): OK')
179  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
180  disp('------------------------')
181  else
182  disp('----------------------------')
183  disp(' Test 3 (results): FAILED')
184  fprintf(' -> found numerical order %f. Must be 2\n',P(1))
185  disp('----------------------------')
186  end
187 end