| OptFEM2DP1 Toolbox 1.2b4 Matlab/Octave Optimized P1-Lagrange Finite Element Method in 2D

- Main Page
- Related Pages
- Files

- File List
- File Members

valid/validStiffElas2DP1.m
Go to the documentation of this file.

    00001 function validStiffElas2DP1(Num)
    00002 % function validStiffElas2DP1(Num)
    00003 %  Validation function for the Assembly of the Stiffness Elasticity
    00004 %  Matrix for `P_1`-Lagrange finite elements.
    00005 %
    00006 %   The Stiffness Elasticity Matrix is given by
    00007 %   ``\StiffElas_{m,l}=\int_{\DOMH} \Odv^t(\BasisFuncTwoD_m) \Ocv(\BasisFuncTwoD_l)dT, \ \forall (m,l)\in\ENS{1}{2\,\nq}^2,``
    00008 %   where `\BasisFuncTwoD_m` are `P_1`-Lagrange vector basis functions.
    00009 %   Here `\Ocv=(\Occ_{xx},\Occ_{yy},\Occ_{xy})^t` and `\Odv=(\Odc_{xx},\Odc_{yy},2\Odc_{xy})^t`
    00010 %   are the elastic stress and strain tensors respectively.
    00011 %
    00012 %   This Matrix is computed by functions StiffElasAssembling2DP1{Version} where {Version} is one of
    00013 %     'base', 'OptV0', 'OptV1' and 'OptV2'.
    00014 %     - Test 1: computation of the Elastic Stiffness Matrix using previous functions giving errors and computation times
    00015 %     - Test 2: compute ``\int_\DOM \Odv^t(\vecb{u}) \Ocv(\vecb{v})d\q \approx \DOT{\StiffElas \vecb{U}}{\vecb{V}}``
    00016 %       with `\foncdefsmall{\vecb{u}=(u_1,u_2)}{\DOM}{\R^2}` and `\foncdefsmall{\vecb{v}=(v_1,v_2)}{\DOM}{\R^2}`.
    00017 %       We have `\vecb{U}_{2i-1}=u_1(\q^i),` `\vecb{U}_{2i}=u_2(\q^i),` `\vecb{V}_{2i-1}=v_1(\q^i),` `\vecb{V}_{2i}=v_2(\q^i)`
    00018 %       if Num in {0,2} and `\vecb{U}_{i}=u_1(\q^i),` `\vecb{U}_{i+\nq}=u_2(\q^i),` `\vecb{V}_{i}=v_1(\q^i),` `\vecb{V}_{i+\nq}=v_2(\q^i)`
    00019 %       if Num in {1,3}. Use functions `\vecb{u}` and `\vecb{u}` defined in #valid_StiffElas.
    00020 %     - Test 3: retrieve order 2 of `P_1`-Lagrange integration
    00021 %       ``|\int_\DOM \Odv^t(\vecb{u}) \Ocv(\vecb{v})d\q - \DOT{\StiffElas \vecb{U}}{\vecb{V}}| \leq C h^2``
    00022 %
    00023 %   Num :
    00024 %    - 0 global alternate numbering with local alternate numbering (classical method),
    00025 %    - 1 global block numbering with local alternate numbering,
    00026 %    - 2 global alternate numbering with local block numbering,
    00027 %    - 3 global block numbering with local block numbering.
    00028 % See also:
    00029 %   #StiffElasAssembling2DP1base, #StiffElasAssembling2DP1OptV0,
    00030 %   #StiffElasAssembling2DP1OptV1, #StiffElasAssembling2DP1OptV2, #SquareMesh
    00031 %
    00032 % @author François Cuvelier @date 2012-11-26
    00033 % Copyright:
    00034 %   See \ref license
    00035
    00036   disp('********************************************')
    00037   disp('*   Stiff Elas Assembling P1 validations   *')
    00038   disp('********************************************')
    00039   switch Num
    00040   case 0
    00041     s=sprintf('Global alternate numbering / local alternate numbering');
    00042   case 1
    00043     s=sprintf('Global block numbering / local alternate numbering');
    00044   case 2
    00045     s=sprintf('Global alternate numbering / local block numbering');
    00046   case 3
    00047     s=sprintf('Global block numbering / local block numbering');
    00048   otherwise
    00049     error('invalid Num value')
    00050   end
    00051   fprintf('  Numbering Choice : %s\n',s);
    00052
    00053
    00054
    00055   Th=SquareMesh(20);
    00056
    00057 % TEST 1
    00058   disp('-----------------------------------------')
    00059   disp('  Test 1: Matrices errors and CPU times  ')
    00060   disp('-----------------------------------------')
    00061   tic();
    00062   Mbase=StiffElasAssembling2DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.areas,1,1,Num);
    00063   T(1)=toc();
    00064   tic();
    00065   MOptV0=StiffElasAssembling2DP1OptV0(Th.nq,Th.nme,Th.q,Th.me,Th.areas,1,1,Num);
    00066   T(2)=toc();
    00067   Test1.error(1)=norm(Mbase-MOptV0,Inf);
    00068   Test1.name{1}='StiffElasAssembling2DP1OptV0';
    00069   fprintf('    Error P1base vs OptV0 : %e\n',Test1.error(1))
    00070   tic();
    00071   MOptV1=StiffElasAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas,1,1,Num);
    00072   T(3)=toc();
    00073   Test1.error(2)=norm(Mbase-MOptV1,Inf);
    00074   Test1.name{2}='StiffElasAssembling2DP1OptV1';
    00075   fprintf('    Error P1base vs OptV1 : %e\n',Test1.error(2))
    00076   tic();
    00077   MOptV2=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,1,1,Num);
    00078   T(4)=toc();
    00079   Test1.error(3)=norm(Mbase-MOptV2,Inf);
    00080   Test1.name{3}='StiffElasAssembling2DP1OptV2';
    00081   fprintf('    Error P1base vs OptV2 : %e\n',Test1.error(3))
    00082
    00083   fprintf('    CPU times base (ref) : %3.4f (s)\n',T(1))
    00084   fprintf('    CPU times OptV0       : %3.4f (s) - Speed Up X%3.3f\n',T(2),T(1)/T(2))
    00085   fprintf('    CPU times OptV1       : %3.4f (s) - Speed Up X%3.3f\n',T(3),T(1)/T(3))
    00086   fprintf('    CPU times OptV2       : %3.4f (s) - Speed Up X%3.3f\n',T(4),T(1)/T(4))
    00087   checkTest1(Test1)
    00088
    00089
    00090
    00091
    00092 % TEST 2
    00093   disp('-----------------------------------------------------')
    00094   disp('  Test 2: Validations by integration on [0,1]x[0,1]  ')
    00095   disp('-----------------------------------------------------')
    00096   Test=valid_StiffElas();
    00097   qf=Th.q;
    00098   for kk=1:length(Test)
    00099     KK=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,Test(kk).lambda,Test(kk).mu,Num);
    00100     u=Test(kk).u;
    00101     v=Test(kk).v;
    00102     switch Num
    00103     case {0,2}
    00104       U=[u{1}(qf(1,:),qf(2,:));u{2}(qf(1,:),qf(2,:))];
    00105       U=U(:);
    00106       V=[v{1}(qf(1,:),qf(2,:));v{2}(qf(1,:),qf(2,:))];
    00107       V=V(:);
    00108     case {1,3}
    00109       U=[u{1}(qf(1,:),qf(2,:)) u{2}(qf(1,:),qf(2,:))]';
    00110       V=[v{1}(qf(1,:),qf(2,:)) v{2}(qf(1,:),qf(2,:))]';
    00111     end
    00112     %whos
    00113     Test(kk).error=abs(Test(kk).Stiff-U'*KK*V);
    00114     fprintf('    functions %d : u(x,y)=(%s,%s), v(x,y)=(%s,%s),\n           -> Stiff error=%e\n', ...
    00115             kk,Test(kk).cu{1},Test(kk).cu{2},Test(kk).cv{1},Test(kk).cv{2},Test(kk).error);
    00116   end
    00117   checkTest2(Test)
    00118
    00119 % TEST 3
    00120   disp('--------------------------------')
    00121   disp('  Test 3: Validations by order  ')
    00122   disp('--------------------------------')
    00123   n=length(Test);
    00124   u=Test(n).u;
    00125   v=Test(n).v;
    00126   lambda=Test(n).lambda;
    00127   mu=Test(n).mu;
    00128   ExSol=Test(n).Stiff;
    00129
    00130   for k=1:10
    00131     Th=SquareMesh(20*k+50);
    00132     qf=Th.q;
    00133     fprintf('    Matrix size : %d\n',Th.nq);
    00134     h(k)=GetMaxLengthEdges(Th.q,Th.me);
    00135     tic();
    00136     M=StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,lambda,mu,Num);
    00137     TT(k)=toc();
    00138     switch Num
    00139     case {0,2}
    00140       U=[u{1}(qf(1,:),qf(2,:));u{2}(qf(1,:),qf(2,:))];
    00141       U=U(:);
    00142       V=[v{1}(qf(1,:),qf(2,:));v{2}(qf(1,:),qf(2,:))];
    00143       V=V(:);
    00144     case {1,3}
    00145       U=[u{1}(qf(1,:),qf(2,:)) u{2}(qf(1,:),qf(2,:))]';
    00146       V=[v{1}(qf(1,:),qf(2,:)) v{2}(qf(1,:),qf(2,:))]';
    00147     end
    00148     Error(k)=abs(ExSol-U'*M*V);
    00149     fprintf('      StiffElasAssembling2DP1OptV2 CPU times : %3.3f(s)\n',TT(k));
    00150     fprintf('      Error                            : %e\n',Error(k));
    00151   end
    00152
    00153   loglog(h,Error,'+-r',h,h*1.1*Error(1)/h(1),'-sm',h,1.1*Error(1)*(h/h(1)).^2,'-db')
    00154   legend('Error','O(h)','O(h^2)')
    00155   xlabel('h')
    00156   title(sprintf('Test 3 : Stiff Elas. Matrix (Num=%d)',Num))
    00157   checkTest3(h,Error)
    00158
    00159   %figure(3)
    00160   %spy(M)
    00161   %title(sprintf('Matrix sparsity for %s numbering',s))
    00162 end
    00163
    00164 function checkTest1(Test)
    00165   I=find(Test.error>1e-14);
    00166   if isempty(I)
    00167     disp('------------------------')
    00168     disp('  Test 1 (results): OK')
    00169     disp('------------------------')
    00170   else
    00171     disp('----------------------------')
    00172     disp('  Test 1 (results): FAILED')
    00173     disp('----------------------------')
    00174   end
    00175 end
    00176
    00177 function checkTest2(Test)
    00178   N=length(Test);
    00179   cntFalse=0;
    00180   for k=1:N
    00181     if (ismember(Test(k).degree,[0 1]))
    00182       if (Test(k).error>1e-12)
    00183         cntFalse=cntFalse+1;
    00184       end
    00185     end
    00186   end
    00187   if (cntFalse==0)
    00188     disp('------------------------')
    00189     disp('  Test 2 (results): OK')
    00190     disp('------------------------')
    00191   else
    00192     disp('----------------------------')
    00193     disp('  Test 2 (results): FAILED')
    00194     disp('----------------------------')
    00195   end
    00196 end
    00197
    00198 function checkTest3(h,error)
    00199   % order 2
    00200   P=polyfit(log(h),log(error),1);
    00201   if abs(P(1)-2)<5e-2
    00202     disp('------------------------')
    00203     disp('  Test 3 (results): OK')
    00204     fprintf('    -> found numerical order %f. Must be 2\n',P(1))
    00205     disp('------------------------')
    00206   else
    00207     disp('----------------------------')
    00208     disp('  Test 3 (results): FAILED')
    00209     fprintf('    -> found numerical order %f. Must be 2\n',P(1))
    00210     disp('----------------------------')
    00211   end
    00212 end

All Files Functions