fe_hmm2d_para

PURPOSE ^

FE_HMM2D_PARA two-dimensional heterogeneous multiscale finite element method method for the heat equation.

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

FE_HMM2D_PARA   two-dimensional heterogeneous multiscale finite element method method for the heat equation.
   FE_HMM2D_PARA solves the parabolic equation
              d/dt u =  div(a*grad(u))  +  f  in Omega
                   u = u_D                 on the Dirichlet boundary
              d/dn u = g                   on the Neumann boundary
   with a highly oscillatory tensor a(x)
   on a geometry described by triangles and parallelograms, respectively
   and presents the solution graphically.

   Time stepping is done using an implicit Euler method.


   The code is available at http://anmc.epfl.ch/ and described in 
   further detail in 

   A. Abdulle and A. Nonnenmacher
   "A short and versatile finite element multiscale code for
   homogenization problems"
   Computer Methods in Applied Mechanics and Engineering,
   http://dx.doi.org/10.1016/j.cma.2009.03.019

   Please cite this article in any publication describing research
   performed using the software.



   The structure of the directories is as follows:
   fe-hmm          
   +---examples
   | +---elliptic    contains main file and files for the problem specific setup
   | | +---data      contains geometry data files
   | +---parabolic   contains main file and files for the problem specific setup  
   | | +---data      contains geometry data files
   | +---tools       contains problem independent algorithm files
   
   To adapt the program to a given setup of the equation, the user has to 
   modify the geometry data files in ./data and the M-files
   <f_t.m>, <g_t.m>, <u_d_t.m> and <tensor_a.m>. They have to be in the same
   directory as <fe_hmm2d_para.m>.


   This code is based in part on
   [1] J. Alberty, C. Carstensen & S. Funken
       "Remarks around 50 lines of Matlab: short finite element implementation",
        Numerical Algorithms, Springer, 1999, 20, 117-137

   Email           : assyr.abdulle@epfl.ch and achim.nonnenmacher@epfl.ch
   Last updated    : 02/08/2010 with MATLAB 7.8

   FE_HMM2D is Copyright (C) 2009 A. Abdulle and A. Nonnenmacher. 
   The software is provided free for non-commercial use unter the terms of 
   the GNU General Public License. See "copyright.m" for full details.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %FE_HMM2D_PARA   two-dimensional heterogeneous multiscale finite element method method for the heat equation.
0002 %   FE_HMM2D_PARA solves the parabolic equation
0003 %              d/dt u =  div(a*grad(u))  +  f  in Omega
0004 %                   u = u_D                 on the Dirichlet boundary
0005 %              d/dn u = g                   on the Neumann boundary
0006 %   with a highly oscillatory tensor a(x)
0007 %   on a geometry described by triangles and parallelograms, respectively
0008 %   and presents the solution graphically.
0009 %
0010 %   Time stepping is done using an implicit Euler method.
0011 %
0012 %
0013 %   The code is available at http://anmc.epfl.ch/ and described in
0014 %   further detail in
0015 %
0016 %   A. Abdulle and A. Nonnenmacher
0017 %   "A short and versatile finite element multiscale code for
0018 %   homogenization problems"
0019 %   Computer Methods in Applied Mechanics and Engineering,
0020 %   http://dx.doi.org/10.1016/j.cma.2009.03.019
0021 %
0022 %   Please cite this article in any publication describing research
0023 %   performed using the software.
0024 %
0025 %
0026 %
0027 %   The structure of the directories is as follows:
0028 %   fe-hmm
0029 %   +---examples
0030 %   | +---elliptic    contains main file and files for the problem specific setup
0031 %   | | +---data      contains geometry data files
0032 %   | +---parabolic   contains main file and files for the problem specific setup
0033 %   | | +---data      contains geometry data files
0034 %   | +---tools       contains problem independent algorithm files
0035 %
0036 %   To adapt the program to a given setup of the equation, the user has to
0037 %   modify the geometry data files in ./data and the M-files
0038 %   <f_t.m>, <g_t.m>, <u_d_t.m> and <tensor_a.m>. They have to be in the same
0039 %   directory as <fe_hmm2d_para.m>.
0040 %
0041 %
0042 %   This code is based in part on
0043 %   [1] J. Alberty, C. Carstensen & S. Funken
0044 %       "Remarks around 50 lines of Matlab: short finite element implementation",
0045 %        Numerical Algorithms, Springer, 1999, 20, 117-137
0046 %
0047 %   Email           : assyr.abdulle@epfl.ch and achim.nonnenmacher@epfl.ch
0048 %   Last updated    : 11/23/2009 with MATLAB 7.4
0049 %
0050 %   FE_HMM2D is Copyright (C) 2009 A. Abdulle and A. Nonnenmacher.
0051 %   The software is provided free for non-commercial use unter the terms of
0052 %   the GNU General Public License. See "copyright.m" for full details.
0053 %
0054 
0055 %
0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 
0058 
0059 % check if tools directory was added to the path
0060 if(exist('shapefunction_quad_lin.m')~=2)
0061    error('HMMFEM:path', strcat('Please add "tools" directory to the Matlab path.\n',...
0062        'You can do that with the command:\n',...
0063        'addpath(''../../tools'');'));
0064 end
0065 
0066 
0067 % Variables
0068 epsilon =5e-2;
0069 NMicro  =4;
0070 
0071 delta= epsilon;
0072 
0073 
0074 % Initialisation
0075 Coordinates=load('./data/coordinates.dat');
0076 Elements3  =load('./data/elements3.dat');
0077 Elements4  =load('./data/elements4.dat');
0078 Neumann    =load('./data/neumann.dat');
0079 Dirichlet  =load('./data/dirichlet.dat');
0080 
0081 T = 1; dt = 0.1; N = T/dt;
0082 U = zeros(size(Coordinates,1),N+1);
0083 
0084 % Initial Condition
0085 U(:,1) = load('./data/u_0.dat');
0086 
0087 bctype='periodic';
0088 % bctype='dirichlet';
0089 
0090 FreeNodes=setdiff(1:size(Coordinates,1),unique(Dirichlet));
0091 A = sparse(size(Coordinates,1),size(Coordinates,1));
0092 B = sparse(size(Coordinates,1),size(Coordinates,1));
0093 b = sparse(size(Coordinates,1),1);
0094 
0095 % HMM initialization, microgrid with quadrilaterals, periodic constraints
0096 MicroElements=micromesh_elements(NMicro);
0097 Constraints=make_constraints(NMicro, bctype);
0098 
0099 % Assembly
0100 % For meshes with many DOF this can be very slow due to Matlab's way of
0101 % handling sparse matrices. An optimized code can be found at the very
0102 % bottom of this file.
0103 for j=1:size(Elements3,1)
0104     M=hmm_stima_tri(NMicro,epsilon,MicroElements,...
0105         Coordinates(Elements3(j,:),:),Constraints, delta, bctype);
0106     A(Elements3(j,:), Elements3(j,:)) = A(Elements3(j,:), Elements3(j,:))+M;
0107 end
0108 
0109 for j=1:size(Elements4,1)
0110     M=hmm_stima_quad(NMicro,epsilon,MicroElements,...
0111         Coordinates(Elements4(j,:),:),Constraints, delta, bctype);
0112     A(Elements4(j,:), Elements4(j,:)) = A(Elements4(j,:), Elements4(j,:))+M;
0113 end
0114 
0115 % Assembly mass matrix
0116 for j = 1:size(Elements3,1)
0117   B(Elements3(j,:),Elements3(j,:)) = B(Elements3(j,:),Elements3(j,:)) ...
0118       + det([1,1,1;Coordinates(Elements3(j,:),:)'])*[2,1,1;1,2,1;1,1,2]/24;
0119 end
0120 
0121 for j = 1:size(Elements4,1)
0122   B(Elements4(j,:),Elements4(j,:)) = B(Elements4(j,:),Elements4(j,:)) ...
0123       + det([1,1,1;Coordinates(Elements4(j,1:3),:)'])*...
0124       ([4 2 1 2; 2 4 2 1; 1 2 4 2; 2 1 2 4 ])/36;
0125 end
0126 
0127 
0128 % time steps
0129 for n = 2:N+1
0130     b = sparse(size(Coordinates,1),1);
0131 
0132     % Volume Forces
0133     for j=1:size(Elements3,1)
0134         b(Elements3(j,:)) = b(Elements3(j,:))+ ...
0135             dt*det([1,1,1; Coordinates(Elements3(j,:),:)'])* ...
0136             f_t(sum(Coordinates(Elements3(j,:),:))/3, (n-1)*dt)/6;
0137     end
0138 
0139     for j=1:size(Elements4,1)
0140         b(Elements4(j,:)) = b(Elements4(j,:))+ ...
0141             dt*det([1,1,1; Coordinates(Elements4(j,1:3),:)'])* ...
0142             f_t(sum(Coordinates(Elements4(j,:),:))/4, (n-1)*dt)/4;
0143     end
0144 
0145 
0146     % Neumann conditions
0147     for j = 1 : size(Neumann,1)
0148         b(Neumann(j,:))=b(Neumann(j,:)) + dt*norm(Coordinates(Neumann(j,1),:)- ...
0149             Coordinates(Neumann(j,2),:)) * g_t(sum(Coordinates(Neumann(j,:),:))/2, (n-1)*dt)/2;
0150     end
0151 
0152     % previous timestep
0153     b = b + B * U(:,n-1);
0154 
0155     % Dirichlet conditions
0156     u = sparse(size(Coordinates,1),1);
0157     u(unique(Dirichlet)) = u_d_t(Coordinates(unique(Dirichlet),:),(n-1)*dt);
0158     b = b - (dt * A + B) * u;
0159 
0160     % Computation of the solution
0161     u(FreeNodes) = (dt*A(FreeNodes,FreeNodes)+ ...
0162         B(FreeNodes,FreeNodes))\b(FreeNodes);
0163     U(:,n) = u;
0164 end
0165 
0166 show_parabolic(Elements3,Elements4,Coordinates,U,N);
0167 
0168 
0169 
0170 
0171 
0172 enorm=sqrt(U(:,N+1)'*A*U(:,N+1));
0173 disp(sprintf('energy norm is |u(1)|=%g',enorm));
0174 disp(sprintf('max         is |u(1)|=%g',max(U(:,N+1))));
0175 
0176 
0177 
0178 
0179 
0180 %
0181 % Appendix
0182 %
0183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0184 
0185 
0186 % % Assembly - faster version
0187 % % For meshes with many DOF an optimized way of assembling the stiffness
0188 % % matrix should be used, such as:
0189 %
0190 % indexi=zeros(3*3*size(Elements3,1),1);
0191 % indexj=zeros(3*3*size(Elements3,1),1);
0192 % value= zeros(3*3*size(Elements3,1),1);
0193 %
0194 % for j=1:size(Elements3,1)
0195 %     disp(j);
0196 %     M=hmm_stima_tri(NMicro,epsilon,MicroElements,...
0197 %         Coordinates(Elements3(j,:),:),Constraints, delta, bctype);
0198 %
0199 %     tmpi=repmat(Elements3(j,:),1,3);
0200 %     tmpj=repmat(Elements3(j,:),3,1);
0201 %
0202 %     indexi(1+(j-1)*9:j*9)=tmpi(:);
0203 %     indexj(1+(j-1)*9:j*9)=tmpj(:);
0204 %     value(1+(j-1)*9:j*9)=M(:);
0205 %
0206 %     %A(Elements3(j,:), Elements3(j,:)) = A(Elements3(j,:), Elements3(j,:))+M;
0207 % end
0208 % A=sparse(indexi,indexj,value,...
0209 %             size(Coordinates,1),size(Coordinates,1));
0210 %
0211 %
0212 %
0213 % indexi=zeros(4*4*size(Elements4,1),1);
0214 % indexj=zeros(4*4*size(Elements4,1),1);
0215 % value= zeros(4*4*size(Elements4,1),1);
0216 %
0217 % for j=1:size(Elements4,1)
0218 %     disp(j);
0219 %     M=hmm_stima_quad(NMicro,epsilon,MicroElements,...
0220 %         Coordinates(Elements4(j,:),:),Constraints, delta, bctype);
0221 %
0222 %     tmpi=repmat(Elements4(j,:),1,4);
0223 %     tmpj=repmat(Elements4(j,:),4,1);
0224 %
0225 %     indexi(1+(j-1)*16:j*16)=tmpi(:);
0226 %     indexj(1+(j-1)*16:j*16)=tmpj(:);
0227 %     value(1+(j-1)*16:j*16)=M(:);
0228 % end
0229 % A=A+sparse(indexi,indexj,value,...
0230 %             size(Coordinates,1),size(Coordinates,1));

Generated on Tue 21-Jul-2009 10:55:32 by m2html © 2003