HMM_STIMA_TRI Computes element stiffness matrix for macro triangles. A_MACRO = HMM_STIMA_TRI(VERTICES, MACROQUADNODE, EPSILON) computes the entry to the macro element stiffness matrix for macro triangles 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 triangles. 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_tri(NMicro, epsilon, MicroElements,... 0002 MacroVertices, Constraints, delta, bctype) 0003 %HMM_STIMA_TRI Computes element stiffness matrix for macro triangles. 0004 % A_MACRO = HMM_STIMA_TRI(VERTICES, MACROQUADNODE, EPSILON) 0005 % computes the entry to the macro element stiffness matrix for macro 0006 % triangles 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 triangles. 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 % This function should not be modified. 0032 % 0033 % 0034 % The code is available at http://anmc.epfl.ch/ and described in 0035 % further detail in 0036 % 0037 % A. Abdulle and A. Nonnenmacher 0038 % "A short and versatile finite element multiscale code for 0039 % homogenization problems" 0040 % Computer Methods in Applied Mechanics and Engineering, 0041 % http://dx.doi.org/10.1016/j.cma.2009.03.019 0042 % 0043 % Please cite this article in any publication describing research 0044 % performed using the software. 0045 % 0046 % 0047 % Email : assyr.abdulle@epfl.ch and achim.nonnenmacher@epfl.ch 0048 % Last updated : 04/29/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 % generate micro-coordinates 0059 bary=sum(MacroVertices(:,:))/3; 0060 [MicroCoordinates]=micromesh_coords(NMicro, bary, delta); 0061 0062 % initialization of variables 0063 NoOfNodes = size(MicroCoordinates,1); 0064 A = sparse(NoOfNodes,NoOfNodes); 0065 0066 % stiffness matrix assembly 0067 % For meshes with many DOF this can be very slow due to Matlab's way of 0068 % handling sparse matrices. An optimized code can be found at the very 0069 % bottom of this file. 0070 for j=1:size(MicroElements,1) 0071 A(MicroElements(j,:), MicroElements(j,:)) = ... 0072 A(MicroElements(j,:), MicroElements(j,:)) + ... 0073 hmm_microstima_quad(MicroCoordinates(MicroElements(j,:),:),... 0074 bary, epsilon); 0075 end 0076 0077 0078 0079 0080 % solution of cell problems dilatation for the macro element 0081 D=[MacroVertices(2,:)-MacroVertices(1,:);MacroVertices(3,:)-MacroVertices(1,:)]'; 0082 0083 % get constraint matrix and corresponding rhs for the lagrange multiplier 0084 0085 switch(lower(bctype)) 0086 case{'dirichlet'} 0087 [ConstraintMat, ConstraintRhs]=... 0088 micro_constraints_dirichlet(Constraints, D, MacroVertices ,... 0089 MicroElements, MicroCoordinates, 0); 0090 case{'periodic'} 0091 [ConstraintMat, ConstraintRhs]=... 0092 micro_constraints_periodic(Constraints, D, MacroVertices ,... 0093 MicroElements, MicroCoordinates, 0); 0094 end 0095 0096 0097 % assemble matrix for constraint system and corresponding rhs 0098 AConstr=[A ConstraintMat';... 0099 ConstraintMat sparse(size(ConstraintMat,1),size(ConstraintMat,1))]; 0100 RhsConstr=[zeros(NoOfNodes,3);ConstraintRhs]; 0101 0102 % finally solve 0103 x=AConstr\RhsConstr; 0104 alpha=x(1:size(MicroCoordinates,1),:); 0105 0106 % area of macro-element and micro-domain 0107 K_macro=.5*det([1,1,1; MacroVertices']); 0108 K_micro=delta^2; 0109 0110 % the entries for the macro stiffness matrix 0111 A_macro=K_macro/K_micro *alpha'*A*alpha; 0112 0113 end 0114 0115 0116 0117 0118 % 0119 % Appendix 0120 % 0121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0122 0123 % % stiffness matrix assembly - fast version 0124 % indexx=zeros(4*4*size(MicroElements,1),1); 0125 % indexy=zeros(4*4*size(MicroElements,1),1); 0126 % value= zeros(4*4*size(MicroElements,1),1); 0127 % 0128 % for j=1:size(MicroElements,1) 0129 % M=hmm_microstima_quad(MicroCoordinates(MicroElements(j,:),:),... 0130 % bary, epsilon); 0131 % 0132 % tmpa=repmat(MicroElements(j,:),1,4); 0133 % tmpb=repmat(MicroElements(j,:),4,1); 0134 % tmpM(:,:)=M; 0135 % indexx(1+(j-1)*16:j*16)=tmpa(:); 0136 % indexy(1+(j-1)*16:j*16)=tmpb(:); 0137 % value(1+(j-1)*16:j*16)=tmpM(:); 0138 % end 0139 % A=sparse(indexx,indexy,value,... 0140 % size(MicroCoordinates,1),size(MicroCoordinates,1));