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