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

- Main Page
- Related Pages
- Files

- File List
- File Members

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

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

All Files Functions