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.
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));