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