fe_hmm2d

PURPOSE ^

FE_HMM2D two-dimensional heterogeneous multiscale finite element method method for elliptic problems.

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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