hmm_stima_quad

PURPOSE ^

HMM_STIMA_QUAD Computes element stiffness matrix for macro parallelograms.

SYNOPSIS ^

function [A_macro]=hmm_stima_quad(NMicro, epsilon, MicroElements,MacroVertices, Constraints, delta, bctype)

DESCRIPTION ^

HMM_STIMA_QUAD   Computes element stiffness matrix for macro parallelograms.
   A_MACRO = HMM_STIMA_QUAD(VERTICES, MACROQUADNODE, EPSILON) 
   computes the entry to the macro element stiffness matrix for macro
   parallelograms using parallelograms on the micro sampling domain and
   specified boundary conditions on the micro problem.

   NMICRO  is the number of degrees of freedom per space dimension in the
   micro domain

   EPSILON is the periodicity of the tensor a.

   MICROELEMENTS has dimension (NMicro-1)^2 x 4 and contain the node
   numbers of the triangulation of the micro sampling domain.

   MACROVERTICES contains the node coordinates of the macro parallelograms.

   CONSTRAINTS contains information about the micro boundary problem. The
   format depends on whether dirichlet or periodic boundary conditions are
   used. See make_contraints.m for details.

   DELTA is the width or height of the square sampling domain.

   BCTYPE is either "dirichlet" or "periodic" and determines the boundary
   condition of the micro problem.
   
   A_MACRO has dimension 4 x 4. The vertices are numbered 
   anti-clockwise. 


   This function should not be modified.


   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.


   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 function [A_macro]=hmm_stima_quad(NMicro, epsilon, MicroElements,...
0002     MacroVertices, Constraints, delta, bctype)
0003 %HMM_STIMA_QUAD   Computes element stiffness matrix for macro parallelograms.
0004 %   A_MACRO = HMM_STIMA_QUAD(VERTICES, MACROQUADNODE, EPSILON)
0005 %   computes the entry to the macro element stiffness matrix for macro
0006 %   parallelograms using parallelograms on the micro sampling domain and
0007 %   specified boundary conditions on the micro problem.
0008 %
0009 %   NMICRO  is the number of degrees of freedom per space dimension in the
0010 %   micro domain
0011 %
0012 %   EPSILON is the periodicity of the tensor a.
0013 %
0014 %   MICROELEMENTS has dimension (NMicro-1)^2 x 4 and contain the node
0015 %   numbers of the triangulation of the micro sampling domain.
0016 %
0017 %   MACROVERTICES contains the node coordinates of the macro parallelograms.
0018 %
0019 %   CONSTRAINTS contains information about the micro boundary problem. The
0020 %   format depends on whether dirichlet or periodic boundary conditions are
0021 %   used. See make_contraints.m for details.
0022 %
0023 %   DELTA is the width or height of the square sampling domain.
0024 %
0025 %   BCTYPE is either "dirichlet" or "periodic" and determines the boundary
0026 %   condition of the micro problem.
0027 %
0028 %   A_MACRO has dimension 4 x 4. The vertices are numbered
0029 %   anti-clockwise.
0030 %
0031 %
0032 %   This function should not be modified.
0033 %
0034 %
0035 %   The code is available at http://anmc.epfl.ch/ and described in
0036 %   further detail in
0037 %
0038 %   A. Abdulle and A. Nonnenmacher
0039 %   "A short and versatile finite element multiscale code for
0040 %   homogenization problems"
0041 %   Computer Methods in Applied Mechanics and Engineering,
0042 %   http://dx.doi.org/10.1016/j.cma.2009.03.019
0043 %
0044 %   Please cite this article in any publication describing research
0045 %   performed using the software.
0046 %
0047 %
0048 %   Email           : assyr.abdulle@epfl.ch and achim.nonnenmacher@epfl.ch
0049 %   Last updated    : 04/29/2009 with MATLAB 7.4
0050 %
0051 %   FE_HMM2D is Copyright (C) 2009 A. Abdulle and A. Nonnenmacher.
0052 %   The software is provided free for non-commercial use unter the terms of
0053 %   the GNU General Public License. See "copyright.m" for full details.
0054 
0055 %
0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 
0058 % initialization of variables
0059 A_macro=zeros(4,4);
0060 
0061 % transformation matrix for affine map to reference element
0062 Map=[MacroVertices(2,:)-MacroVertices(1,:);MacroVertices(4,:)-MacroVertices(1,:)]';
0063 % reference and physical quadrilateral quadrature nodes
0064 quadnodes_ref=[.5-sqrt(3)/6, .5+sqrt(3)/6, .5-sqrt(3)/6, .5+sqrt(3)/6;...
0065                .5-sqrt(3)/6, .5-sqrt(3)/6, .5+sqrt(3)/6, .5+sqrt(3)/6];
0066 quadnodes= Map*quadnodes_ref+repmat(MacroVertices(1,:)',1,4);
0067 
0068 for node_no=1:4
0069     % select current quadrature nodes
0070     quadnode=quadnodes(:,node_no);
0071     refquadnode=quadnodes_ref(:,node_no);
0072 
0073     % Generate micro-coordinates
0074     [MicroCoordinates]=micromesh_coords(NMicro, quadnode, delta);
0075 
0076     % Initialization
0077     NoOfNodes = size(MicroCoordinates,1);
0078     A = sparse(NoOfNodes,NoOfNodes);
0079 
0080     %Stiffness matrix Assembly
0081     % For meshes with many DOF this can be very slow due to Matlab's way of
0082     % handling sparse matrices. An optimized code can be found at the very
0083     % bottom of this file.
0084     for j=1:size(MicroElements,1)
0085         A(MicroElements(j,:), MicroElements(j,:)) = ...
0086             A(MicroElements(j,:), MicroElements(j,:)) + ...
0087             hmm_microstima_quad(MicroCoordinates(MicroElements(j,:),:),...
0088             quadnode, epsilon);
0089     end
0090         
0091 
0092     % Get constraint matrix and corresponding rhs for the lagrange multiplier
0093     % micro boundary conditions.
0094     switch(lower(bctype))
0095         case{'dirichlet'}
0096             [ConstraintMat, ConstraintRhs]=...
0097                 micro_constraints_dirichlet(Constraints, Map, MacroVertices ,...
0098                 MicroElements, MicroCoordinates, refquadnode);
0099         case{'periodic'}
0100             [ConstraintMat, ConstraintRhs]=...
0101                 micro_constraints_periodic(Constraints, Map, MacroVertices ,...
0102                 MicroElements, MicroCoordinates, refquadnode);
0103     end
0104 
0105     % Assemble matrix AConstr for constrained system and corresponding rhs
0106     AConstr=[A ConstraintMat';...
0107         ConstraintMat sparse(size(ConstraintMat,1),size(ConstraintMat,1))];
0108     RhsConstr=[zeros(NoOfNodes,4);ConstraintRhs];
0109 
0110     % solve constrained microproblem
0111     x=AConstr\RhsConstr;
0112     alpha=x(1:size(MicroCoordinates,1),:);
0113 
0114     % area of macro-element and micro-domain
0115     K_macro=det([1,1,1; MacroVertices(1:3,:)']);
0116     K_micro=delta^2;
0117 
0118     %The contribution to the macro stiffness matrix
0119     A_macro=A_macro+ .25 * K_macro/K_micro *alpha'*A*alpha;
0120 end
0121 
0122 end
0123 
0124 
0125 
0126 
0127 
0128 
0129 
0130 %
0131 % Appendix
0132 %
0133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0134 
0135 %     % stiffness matrix assembly - fast version
0136 %
0137 %     indexx=zeros(4*4*size(MicroElements,1),1);
0138 %     indexy=zeros(4*4*size(MicroElements,1),1);
0139 %     value= zeros(4*4*size(MicroElements,1),1);
0140 %
0141 %     for j=1:size(MicroElements,1)
0142 %         M=hmm_microstima_quad(MicroCoordinates(MicroElements(j,:),:),...
0143 %              quadnode, epsilon);
0144 %
0145 %         tmpa=repmat(MicroElements(j,:),1,4);
0146 %         tmpb=repmat(MicroElements(j,:),4,1);
0147 %         tmpM(:,:)=M;
0148 %         indexx(1+(j-1)*16:j*16)=tmpa(:);
0149 %         indexy(1+(j-1)*16:j*16)=tmpb(:);
0150 %         value(1+(j-1)*16:j*16)=tmpM(:);
0151 %     end
0152 %     A=sparse(indexx,indexy,value,...
0153 %         size(MicroCoordinates,1),size(MicroCoordinates,1));

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