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

- Main Page
- Related Pages
- Files

- File List
- File Members

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

    00001 function validStiff2DP1()
    00002 % function validStiff2DP1()
    00003 %  Validation function for the assembly of the stiffness matrix for
    00004 %  `P_1`-Lagrange finite element method
    00005 %
    00006 %   The Stiffness Matrix `\Stiff` is given by
    00007 %   ``\Stiff_{i,j}=\int_\DOMH \DOT{\GRAD\FoncBase_i(\q)}{\GRAD\FoncBase_j(\q)}d\q,\ \forall (i,j)\in\ENS{1}{\nq}^2``
    00008 %   where `\FoncBase_i` are `P_1`-Lagrange basis functions.
    00009 %   This Matrix is computed by functions StiffAssembling2DP1{Version} where {Version} is one of
    00010 %     'base', 'OptV0', 'OptV1' and 'OptV2'.
    00011 %     - Test 1: Computation of the Stiffness Matrix using all the versions giving errors and computation times
    00012 %     - Test 2: Computation of the integral ``\int_\DOM \DOT{\GRAD u(\q)}{\GRAD v(\q)}d\q \approx \DOT{\Stiff \vecb{U}}{\vecb{V}}``
    00013 %       where `\vecb{U}_i=u(\q^i)` and `\vecb{V}_i=v(\q^i)`.
    00014 %       Functions `u` and `v` are those defined in #valid_FEMmatrices.
    00015 %     - Test 3: One retrieves the order 2 of `P_1`-Lagrange integration
    00016 %       ``|\int_\DOM \DOT{\GRAD u}{\GRAD v} -\DOT{\GRAD \Pi_h(u)}{\GRAD \Pi_h(v)}d\DOM| \leq C h^2``
    00017 %
    00018 % See also:
    00019 %   #StiffAssembling2DP1base, #StiffAssembling2DP1OptV0,
    00020 %   #StiffAssembling2DP1OptV1, #StiffAssembling2DP1OptV2,
    00021 %   #valid_FEMmatrices, #SquareMesh, #GetMaxLengthEdges
    00022 %
    00023 % Results:
    00024 %  \image html images/validStiff2DP1.png "figure : validStiff2DP1 Test 3 result"
    00025 %
    00026 % @author François Cuvelier @date 2012-11-26
    00027 % Copyright:
    00028 %   See \ref license
    00029
    00030   disp('*******************************************')
    00031   disp('*     Stiff Assembling P1 validations     *')
    00032   disp('*******************************************')
    00033
    00034   Th=SquareMesh(50);
    00035
    00036 % TEST 1
    00037   disp('-----------------------------------------')
    00038   disp('  Test 1: Matrices errors and CPU times  ')
    00039   disp('-----------------------------------------')
    00040   tic();
    00041   Mbase=StiffAssembling2DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.areas);
    00042   T(1)=toc();
    00043   tic();
    00044   MOptV0=StiffAssembling2DP1OptV0(Th.nq,Th.nme,Th.q,Th.me,Th.areas);
    00045   T(2)=toc();
    00046   Test1.error(1)=norm(Mbase-MOptV0,Inf);
    00047   Test1.name{1}='StiffAssembling2DP1OptV0';
    00048   fprintf('    Error P1base vs OptV0 : %e\n',Test1.error(1))
    00049   tic();
    00050   MOptV1=StiffAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas);
    00051   T(3)=toc();
    00052   Test1.error(2)=norm(Mbase-MOptV1,Inf);
    00053   Test1.name{2}='StiffAssembling2DP1OptV1';
    00054   fprintf('    Error P1base vs OptV1 : %e\n',Test1.error(2))
    00055   tic();
    00056   MOptV2=StiffAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas);
    00057   T(4)=toc();
    00058   Test1.error(3)=norm(Mbase-MOptV2,Inf);
    00059   Test1.name{3}='StiffAssembling2DP1OptV2';
    00060   fprintf('    Error P1base vs OptV2 : %e\n',Test1.error(3))
    00061
    00062   fprintf('    CPU times base (ref) : %3.4f (s)\n',T(1))
    00063   fprintf('    CPU times OptV0       : %3.4f (s) - Speed Up X%3.3f\n',T(2),T(1)/T(2))
    00064   fprintf('    CPU times OptV1       : %3.4f (s) - Speed Up X%3.3f\n',T(3),T(1)/T(3))
    00065   fprintf('    CPU times OptV2       : %3.4f (s) - Speed Up X%3.3f\n',T(4),T(1)/T(4))
    00066   checkTest1(Test1)
    00067
    00068   M=Mbase;
    00069
    00070 % TEST 2
    00071   disp('-----------------------------------------------------')
    00072   disp('  Test 2: Validations by integration on [0,1]x[0,1]  ')
    00073   disp('-----------------------------------------------------')
    00074   Test=valid_FEMmatrices();
    00075   for kk=1:length(Test)
    00076     U=Test(kk).u(Th.q(1,:),Th.q(2,:));
    00077     V=Test(kk).v(Th.q(1,:),Th.q(2,:));
    00078     Test(kk).error=abs(Test(kk).Stiff-U*M*V');
    00079     fprintf('    functions %d : u(x,y)=%s, v(x,y)=%s,\n           -> Stiff error=%e\n',kk,Test(kk).cu,Test(kk).cv,Test(kk).error);
    00080   end
    00081   checkTest2(Test)
    00082
    00083 % TEST 3
    00084   disp('--------------------------------')
    00085   disp('  Test 3: Validations by order  ')
    00086   disp('--------------------------------')
    00087   n=length(Test);
    00088   u=Test(n).u;
    00089   v=Test(n).v;
    00090   ExSol=Test(n).Stiff;
    00091
    00092   for k=1:10
    00093     Th=SquareMesh(50*k+50);
    00094     fprintf('    Matrix size : %d\n',Th.nq);
    00095     h(k)=GetMaxLengthEdges(Th.q,Th.me);
    00096     tic();
    00097     M=StiffAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas);
    00098     TT(k)=toc();
    00099     U=u(Th.q(1,:),Th.q(2,:));
    00100     V=v(Th.q(1,:),Th.q(2,:));
    00101     Error(k)=abs(ExSol-U*M*V');
    00102     fprintf('      StiffAssembling2DP1OptV2 CPU times : %3.3f(s)\n',TT(k));
    00103     fprintf('      Error                            : %e\n',Error(k));
    00104   end
    00105
    00106   loglog(h,Error,'+-r',h,h*1.1*Error(1)/h(1),'-sm',h,1.1*Error(1)*(h/h(1)).^2,'-db')
    00107   legend('Error','O(h)','O(h^2)')
    00108   xlabel('h')
    00109   title('Test 3 : Stiffness Matrix')
    00110   checkTest3(h,Error)
    00111 end
    00112
    00113 function checkTest1(Test)
    00114   I=find(Test.error>1e-14);
    00115   if isempty(I)
    00116     disp('------------------------')
    00117     disp('  Test 1 (results): OK')
    00118     disp('------------------------')
    00119   else
    00120     disp('----------------------------')
    00121     disp('  Test 1 (results): FAILED')
    00122     disp('----------------------------')
    00123   end
    00124 end
    00125
    00126 function checkTest2(Test)
    00127   N=length(Test);
    00128   cntFalse=0;
    00129   for k=1:N
    00130     if ( ismember(Test(k).degree,[0 1]) )
    00131       if (Test(k).error>1e-12)
    00132         cntFalse=cntFalse+1;
    00133       end
    00134     end
    00135   end
    00136   if (cntFalse==0)
    00137     disp('------------------------')
    00138     disp('  Test 2 (results): OK')
    00139     disp('------------------------')
    00140   else
    00141     disp('----------------------------')
    00142     disp('  Test 2 (results): FAILED')
    00143     disp('----------------------------')
    00144   end
    00145 end
    00146
    00147 function checkTest3(h,error)
    00148   % order 2
    00149   P=polyfit(log(h),log(error),1);
    00150   if abs(P(1)-2)<1e-2
    00151     disp('------------------------')
    00152     disp('  Test 3 (results): OK')
    00153     fprintf('    -> found numerical order %f. Must be 2\n',P(1))
    00154     disp('------------------------')
    00155   else
    00156     disp('----------------------------')
    00157     disp('  Test 3 (results): FAILED')
    00158     fprintf('    -> found numerical order %f. Must be 2\n',P(1))
    00159     disp('----------------------------')
    00160   end
    00161 end

All Files Functions