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

- Main Page
- Related Pages
- Files

- File List
- File Members

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

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

All Files Functions