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