1## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo 2## 3## This file is part of: 4## MSH - Meshing Software Package for Octave 5## 6## MSH is free software; you can redistribute it and/or modify 7## it under the terms of the GNU General Public License as published by 8## the Free Software Foundation; either version 2 of the License, or 9## (at your option) any later version. 10## 11## MSH is distributed in the hope that it will be useful, 12## but WITHOUT ANY WARRANTY; without even the implied warranty of 13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14## GNU General Public License for more details. 15## 16## You should have received a copy of the GNU General Public License 17## along with MSH; If not, see <http://www.gnu.org/licenses/>. 18## 19## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> 20## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> 21 22## -*- texinfo -*- 23## @deftypefn {Function File} {[@var{varargout}]} = @ 24## msh3m_geometrical_properties(@var{mesh},[@var{string1},@var{string2},...]) 25## 26## 27## Compute @var{mesh} geometrical properties identified by input strings. 28## 29## Valid properties are: 30## @itemize @bullet 31## @item @code{"bar"}: return a matrix with size 3 times the number of mesh 32## elements containing the center of mass coordinates. 33## @item @code{"wjacdet"}: return the weigthed Jacobian determinant used 34## for the numerical integration with trapezoidal rule over an element. 35## @item @code{"shg"}: return a matrix of size 3 times the number of 36## elements matrix containing the gradient of P1 shape functions. 37## @item @code{"shp"}: return a matrix containing the the value of P1 shape 38## functions. 39## @item @code{"area"}: return a row vector containing the volume of each 40## element. 41## @end itemize 42## 43## The output will contain the geometrical properties requested in the 44## input in the same order specified in the function call. 45## 46## If an unexpected string is given as input, an empty vector is 47## returned in output. 48## 49## @seealso{msh2m_topological_properties, msh2m_geometrical_properties} 50## @end deftypefn 51 52 53function [varargout] = msh3m_geometrical_properties (imesh,varargin) 54 55 ## Check input 56 if (nargin < 2) # Number of input parameters 57 error ("msh3m_geometrical_properties: wrong number of input parameters."); 58 elseif (! (isstruct (imesh) && isfield (imesh,"p") && 59 isfield (imesh,"t") && isfield (imesh,"e"))) 60 error ("msh3m_geometrical_properties: first input is not a valid mesh structure."); 61 elseif (! iscellstr (varargin)) 62 error ("msh3m_geometrical_properties: only string value admitted for properties."); 63 endif 64 65 ## Compute properties 66 67 ## Extract tetrahedra node coordinates 68 x1 = imesh.p(1,imesh.t(1,:)); 69 y1 = imesh.p(2,imesh.t(1,:)); 70 z1 = imesh.p(3,imesh.t(1,:)); 71 x2 = imesh.p(1,imesh.t(2,:)); 72 y2 = imesh.p(2,imesh.t(2,:)); 73 z2 = imesh.p(3,imesh.t(2,:)); 74 x3 = imesh.p(1,imesh.t(3,:)); 75 y3 = imesh.p(2,imesh.t(3,:)); 76 z3 = imesh.p(3,imesh.t(3,:)); 77 x4 = imesh.p(1,imesh.t(4,:)); 78 y4 = imesh.p(2,imesh.t(4,:)); 79 z4 = imesh.p(3,imesh.t(4,:)); 80 81 nelem = columns(imesh.t); # Number of elements in the mesh 82 83 for nn = 1:length (varargin) 84 85 request = varargin{nn}; 86 87 switch request 88 89 case "bar" # Center of mass coordinates 90 if isfield (imesh,"bar") 91 varargout{nn} = imesh.bar; 92 else 93 b = zeros (3, nelem); 94 b(1,:) = ( x1 + x2 + x3 + x4 )/4; 95 b(2,:) = ( y1 + y2 + y3 + y4 )/4; 96 b(3,:) = ( z1 + z2 + z3 + z4 )/4; 97 varargout{nn} = b; 98 clear b; 99 endif 100 101 case "wjacdet" # Weighted Jacobian determinant 102 if isfield (imesh,"wjacdet") 103 varargout{nn} = imesh.wjacdet; 104 else 105 b = wjacdet (x1,y1,z1,... 106 x2,y2,z2,... 107 x3,y3,z3,... 108 x4,y4,z4); 109 varargout{nn} = b; 110 clear b 111 endif 112 113 case "area" # Element area 114 if isfield (imesh,"area") 115 varargout{nn} = imesh.area; 116 else 117 tmp = wjacdet (x1,y1,z1,... 118 x2,y2,z2,... 119 x3,y3,z3,... 120 x4,y4,z4); 121 b = sum (tmp,1); 122 varargout{nn} = b; 123 clear b; 124 endif 125 126 case "shg" # Gradient of shape functions 127 if isfield (imesh,"shg") 128 varargout{nn} = imesh.shg; 129 else 130 b = shg (x1,y1,z1,... 131 x2,y2,z2,... 132 x3,y3,z3,... 133 x4,y4,z4); 134 varargout{nn} = b; 135 clear b 136 endif 137 138 case "shp" # Value of shape functions 139 if isfield (imesh,"shp") 140 varargout{nn} = imesh.shp; 141 else 142 varargout{nn} = eye (4); 143 endif 144 145 otherwise 146 warning ("msh3m_geometrical_properties: unexpected value in property string. Empty vector passed as output.") 147 varargout{nn} = []; 148 149 endswitch 150 151 endfor 152 153endfunction 154 155function [b] = wjacdet(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) 156 157 ## Compute weighted yacobian determinant 158 159 weight = [1/4 1/4 1/4 1/4]'; 160 161 Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); 162 Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); 163 Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); 164 165 ## Determinant of the Jacobian of the 166 ## transformation from the base tetrahedron 167 ## to the tetrahedron K 168 detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; 169 ## Volume of the reference tetrahedron 170 Kkvolume = 1/6; 171 172 b(:,:) = Kkvolume * weight * detJ; 173 174endfunction 175 176function [b] = shg(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) 177 178 ## Compute gradient of shape functions 179 180 Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); 181 Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); 182 Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); 183 184 ## Determinant of the Jacobian of the 185 ## transformation from the base tetrahedron 186 ## to the tetrahedron K 187 detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; 188 189 ## Shape function gradients follow 190 ## First index represents space direction 191 ## Second index represents the shape function 192 ## Third index represents the tetrahedron number 193 b(1,1,:) = (y2.*(z4-z3) + y3.*(z2-z4) + y4.*(z3-z2))./ detJ; 194 b(2,1,:) = (x2.*(z3-z4) + x3.*(z4-z2) + x4.*(z2-z3))./ detJ; 195 b(3,1,:) = (x2.*(y4-y3) + x3.*(y2-y4) + x4.*(y3-y2))./ detJ; 196 197 b(1,2,:) = ( Nb2 ) ./ detJ; 198 b(2,2,:) = (x1.*(z4-z3) + x3.*(z1-z4) + x4.*(z3-z1)) ./ detJ; 199 b(3,2,:) = (x1.*(y3-y4) + x3.*(y4-y1) + x4.*(y1-y3)) ./ detJ; 200 201 b(1,3,:) = ( Nb3 ) ./ detJ; 202 b(2,3,:) = (x1.*(z2-z4) + x2.*(z4-z1) + x4.*(z1-z2)) ./ detJ; 203 b(3,3,:) = (x1.*(y4-y2) + x2.*(y1-y4) + x4.*(y2-y1)) ./ detJ; 204 205 b(1,4,:) = ( Nb4) ./ detJ; 206 b(2,4,:) = (x1.*(z3-z2) + x2.*(z1-z3) + x3.*(z2-z1)) ./ detJ; 207 b(3,4,:) = (x1.*(y2-y3) + x2.*(y3-y1) + x3.*(y1-y2)) ./ detJ; 208endfunction 209 210%!shared mesh,wjacdet,shg,shp 211% x = y = z = linspace(0,1,2); 212% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6) 213% [wjacdet] = msh3m_geometrical_properties(mesh,"wjacdet") 214% [shg] = msh3m_geometrical_properties(mesh,"shg") 215% [shp] = msh3m_geometrical_properties(mesh,"shp") 216%!test 217% assert(columns(mesh.t),columns(wjacdet)) 218%!test 219% assert(size(shg),[3 4 6]) 220%!test 221% assert(shp,eye(4)) 222%!test 223% fail(msh3m_geometrical_properties(mesh,"samanafattababbudoiu"),"warning","Unexpected value in passed string. Empty vector passed as output.") 224