From: <jpi...@us...> - 2012-08-24 09:33:12
|
Revision: 10906 http://octave.svn.sourceforge.net/octave/?rev=10906&view=rev Author: jpicarbajal Date: 2012-08-24 09:33:05 +0000 (Fri, 24 Aug 2012) Log Message: ----------- mechanics: Doc improvement. Bug fixes. Adding new functions Modified Paths: -------------- trunk/octave-forge/main/mechanics/DESCRIPTION trunk/octave-forge/main/mechanics/INDEX trunk/octave-forge/main/mechanics/NEWS trunk/octave-forge/main/mechanics/inst/core/forcematrix.m trunk/octave-forge/main/mechanics/inst/core/inertiamoment.m trunk/octave-forge/main/mechanics/inst/core/masscenter.m trunk/octave-forge/main/mechanics/inst/core/principalaxes.m Added Paths: ----------- trunk/octave-forge/main/mechanics/devel/ trunk/octave-forge/main/mechanics/inst/core/private/drawAxis3D.m trunk/octave-forge/main/mechanics/inst/core/rodmassmatrix.m Removed Paths: ------------- trunk/octave-forge/main/mechanics/inst/core/drawAxis3D.m Modified: trunk/octave-forge/main/mechanics/DESCRIPTION =================================================================== --- trunk/octave-forge/main/mechanics/DESCRIPTION 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/DESCRIPTION 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,6 +1,6 @@ Name: Mechanics -Version: 1.2.0 -Date: 2011-12-09 +Version: 1.3.0 +Date: 2012-24-08 Author: Juan Pablo Carbajal <car...@if...>, Johan Beke <joh...@ho...> Maintainer: Juan Pablo Carbajal <car...@if...>, Johan Beke <joh...@ho...> Title: Classical Mechanics & Structural Analysis Modified: trunk/octave-forge/main/mechanics/INDEX =================================================================== --- trunk/octave-forge/main/mechanics/INDEX 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/INDEX 2012-08-24 09:33:05 UTC (rev 10906) @@ -10,19 +10,19 @@ principalaxes.m 3D Mechanics EAmatrix.m - RBequations_rot.m - RBexample.m - drawAxis3D.m mat2quat.m quat2mat.m quatconj.m quatprod.m quatvrot.m + rodmassmatrix.m + RBequations_rot.m + RBexample.m Molecular Dynamics + lennardjones + mdsim verletstep verletstep_boxed - mdsim - lennardjones 2D Structural analysis (statics) MSNForces PlotDiagrams Modified: trunk/octave-forge/main/mechanics/NEWS =================================================================== --- trunk/octave-forge/main/mechanics/NEWS 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/NEWS 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,13 +1,18 @@ Summary of important user-visible changes for releases of the mechanics package =============================================================================== -mechanics-1.2.X Release Date: XX Release Manager: Juan Pablo Carbajal +mechanics-1.3.0 Release Date: 2012-24-08 Release Manager: Juan Pablo Carbajal =============================================================================== -* Bug fixes - - masscenter.m +* Functions with bug fixes + masscenter +* Added functions + rodmassmatrix +* Removed functions + drawAxes3D is now a private function. + =============================================================================== mechanics-1.2.0 Release Date: 2011-10-16 Release Manager: Juan Pablo Carbajal =============================================================================== Deleted: trunk/octave-forge/main/mechanics/inst/core/drawAxis3D.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/drawAxis3D.m 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/inst/core/drawAxis3D.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,21 +0,0 @@ -%%copyright (c) 2011 Juan Pablocarbajal <car...@if...> -%% -%% This program is free software: youcan redistribute itand/or modify -%% it under the terms of the GNU General Public Licenseas publishedby -%% the Free Software Foundation, either version 3 of the License, or -%% any later version. -%% -%% This program is distributed in the hope that it willbe useful, -%% but WITHOUTaNY WARRANTY; without even the implied warranty of -%% MERCHANTABILITY or FITNESS FORa PARTICULAR PURPOSE. See the -%% GNU General Public License for more details. -%% -%% You should have receivedacopy of the GNU General Public License -%% along with this program. If not, see <http://www.gnu.org/licenses/>. - -function drawAxis3D (p,E,c) - for i=1:3 - l =[p; p+E(i,:)]; - line(l(:,1),l(:,2),l(:,3),'color',c(i,:),'linewidth',2); - end -end Modified: trunk/octave-forge/main/mechanics/inst/core/forcematrix.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/forcematrix.m 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/inst/core/forcematrix.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,5 +1,5 @@ %% Copyright (c) 2010 Juan Pablo Carbajal <car...@if...> -%% +%% %% This program is free software: you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published by %% the Free Software Foundation, either version 3 of the License, or @@ -26,7 +26,7 @@ %% %% @item %% @var{A} is the a connectivity (adjacency) matrix expressed as a vector. That is, if -%% @var{M} is the complete NxN symmetric connectivity matrix, then +%% @var{M} is the complete NxN symmetric connectivity matrix, then %% @code{@var{A}= vech(@var{M})}. The diagonal of @var{M} is not used. The elements %% of A are indexes to the corresponding interaction force in @var{funcs}. For %% example, there are 3 points and points 2 and 3 interact with a force @@ -45,7 +45,7 @@ %% @itemize %% @item %% @var{F} is a Nxdim array of forces acting on each point. -%% +%% %% @item %% @var{D} is an array in the same format as @var{A} containing the euclidean %% distances between the points. @@ -57,17 +57,18 @@ %% @seealso{pointmassmesh, polyval, vech, sub2ind} %% @end deftypefn -function [F D] = forcematrix (pos, A, funcs) +function [F D dr] = forcematrix (pos, A, funcs) %% Argument parsing [N dim] = size(pos); - nA = numel(A); - + nA = numel(A); + % Relative position matrix dr = zeros (nA,dim); for ic = 1:dim dr(:,ic) = vech (bsxfun (@minus, pos(:,ic)', pos(:,ic))); end + % Distance vector D = sqrt (sum (dr.^2, 2)); % Versor @@ -77,15 +78,14 @@ funcs{end+1} = 0; nf = numel(funcs); A(A==0) = nf; - + %% pair to pair forces f = cellfun (@handleOpoly, {funcs{A}}.', mat2cell (D, ones(nA,1), 1)); f = repmat(f,1,dim) .* drhat; - + %% net force on p1 F = zeros (N,dim); - - %% FIXME avoid looping + aux = ones(1,N-1); p2 = sub2ind_tril (N, 2:N, aux); F(1,:) = sum (f(p2,:)); @@ -99,7 +99,7 @@ ndhalf = idx((p1+1):N); F(p1,:) = sum (f(ndhalf,:)) - sum (f(sthalf,:)); end - + endfunction function f = handleOpoly (y,x) @@ -154,7 +154,7 @@ % 3D %!test % Cube %! pos = [0 0 0; 1 0 0; 1 1 0; 0 1 0; ... -%! 0 0 1; 1 0 1; 1 1 1; 0 1 1;]; +%! 0 0 1; 1 0 1; 1 1 1; 0 1 1;]; %! A= [0;1;0;1;1;0;0;0; ... %! 0;1;0;0;1;0;0; ... %! 0;1;0;0;1;0; ... @@ -162,7 +162,7 @@ %! 0;1;0;1; ... %! 0;1;0; ... %! 0;1; ... -%! 0; ]; +%! 0; ]; %! funcs={[-1 1]}; %! [F D] = forcematrix (pos, A, funcs); %! assert (F, zeros(8,3), 1e-8); @@ -178,7 +178,7 @@ %!test % Cube 2 funcs %! pos = [0 0 0; 1 0 0; 1 1 0; 0 1 0; ... -%! 0 0 1; 1 0 1; 1 1 1; 0 1 1;]; +%! 0 0 1; 1 0 1; 1 1 1; 0 1 1;]; %! A= [0;1;0;1;2;0;0;0; ... %! 0;1;0;0;2;0;0; ... %! 0;1;0;0;2;0; ... @@ -186,7 +186,7 @@ %! 0;1;0;1; ... %! 0;1;0; ... %! 0;1; ... -%! 0; ]; +%! 0; ]; %! funcs={[-1 1], @(x)-x}; %! [F D] = forcematrix (pos, A, funcs); %! F_sb = [zeros(4,2) ones(4,1); ... @@ -201,4 +201,3 @@ %! 0;sqrt(1); ... %! 0]; %! assert (D, D_sb, 1e-8); - Modified: trunk/octave-forge/main/mechanics/inst/core/inertiamoment.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/inertiamoment.m 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/inst/core/inertiamoment.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,4 +1,4 @@ -%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...> +%% Copyright (c) 2011-2012 Juan Pablo Carbajal <car...@if...> %% %% This program is free software; you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published by Modified: trunk/octave-forge/main/mechanics/inst/core/masscenter.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/masscenter.m 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/inst/core/masscenter.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -1,4 +1,4 @@ -%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...> +%% Copyright (c) 2011-2012 Juan Pablo Carbajal <car...@if...> %% %% This program is free software; you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published by Modified: trunk/octave-forge/main/mechanics/inst/core/principalaxes.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/principalaxes.m 2012-08-24 02:56:21 UTC (rev 10905) +++ trunk/octave-forge/main/mechanics/inst/core/principalaxes.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -21,6 +21,8 @@ ## axes of the shape. @var{l} is the second moment of area around the correspoding ## principal axis. @var{axes} is order from lower to higher @var{l}. ## +## @var{shape} can be defined by a polygon or by a piece-wise smooth shape. +## ## @seealso{inertiamoment, masscenter} ## @end deftypefn @@ -68,7 +70,7 @@ %! h = 1; b = 2; %! rectangle = [-b/2 -h/2; b/2 -h/2; b/2 h/2; -b/2 h/2]; %! [PA l] = principalaxes(rectangle); -%! assert ( [0 1; 1 0], PA, 1e-6); +%! assert ( [1 0; 0 1], PA, 1e-6); %! assert ([b*h^3; h*b^3]/12, l); %!demo @@ -78,7 +80,8 @@ %! [PAr l] = principalaxes(shapeR); %! [PA l] = principalaxes(shape); %! -%! cla +%! figure (1) +%! clf %! plot(shape(:,1),shape(:,2),'-k'); %! line([0 PA(1,1)],[0 PA(1,2)],'color','r'); %! line([0 PA(2,1)],[0 PA(2,2)],'color','b'); Copied: trunk/octave-forge/main/mechanics/inst/core/private/drawAxis3D.m (from rev 10905, trunk/octave-forge/main/mechanics/inst/core/drawAxis3D.m) =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/private/drawAxis3D.m (rev 0) +++ trunk/octave-forge/main/mechanics/inst/core/private/drawAxis3D.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -0,0 +1,21 @@ +%%copyright (c) 2011 Juan Pablocarbajal <car...@if...> +%% +%% This program is free software: youcan redistribute itand/or modify +%% it under the terms of the GNU General Public Licenseas publishedby +%% the Free Software Foundation, either version 3 of the License, or +%% any later version. +%% +%% This program is distributed in the hope that it willbe useful, +%% but WITHOUTaNY WARRANTY; without even the implied warranty of +%% MERCHANTABILITY or FITNESS FORa PARTICULAR PURPOSE. See the +%% GNU General Public License for more details. +%% +%% You should have receivedacopy of the GNU General Public License +%% along with this program. If not, see <http://www.gnu.org/licenses/>. + +function drawAxis3D (p,E,c) + for i=1:3 + l =[p; p+E(i,:)]; + line(l(:,1),l(:,2),l(:,3),'color',c(i,:),'linewidth',2); + end +end Added: trunk/octave-forge/main/mechanics/inst/core/rodmassmatrix.m =================================================================== --- trunk/octave-forge/main/mechanics/inst/core/rodmassmatrix.m (rev 0) +++ trunk/octave-forge/main/mechanics/inst/core/rodmassmatrix.m 2012-08-24 09:33:05 UTC (rev 10906) @@ -0,0 +1,113 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% This program is free software: you can redistribute it and/or modify +%% it under the terms of the GNU General Public License as published by +%% the Free Software Foundation, either version 3 of the License, or +%% any later version. +%% +%% This program is distributed in the hope that it will be useful, +%% but WITHOUT ANY WARRANTY; without even the implied warranty of +%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%% GNU General Public License for more details. +%% +%% You should have received a copy of the GNU General Public License +%% along with this program. If not, see <http://www.gnu.org/licenses/>. + +%% -*- texinfo -*- +%% @deftypefn {Function File} {[@var{J} @var{sigma}]= } rodmassmatrix (@var{sigma},@var{l}, @var{rho}) +%% Mass matrix of one dimensional rod in 3D. +%% +%% Let q be the configuration vector of the rod, with the first three elements of +%% q being the spatial coordinates (e.g. x,y,z) and the second three elements of +%% q the rotiational coordinates (e.g. Euler angles), then the kinetical energy +%% of the rod is given by +%% T = 1/2 (dqdt)^T kron(J,eye(3)) dqdt +%% +%% @var{sigma} is between 0 and 1. Corresponds to the point in the rod that is +%% being used to indicate the position of the rod in space. +%% If @var{sigma} is a string then the value corresponding to the center of mass +%% of the rod. This makes @var{J} a diagonal matrix. If @var{sigma} is a string +%% the return value of @var{sigma} corresponds to the value pointing to the +%% center of mass. +%% +%% @var{l} is the length of the rod. If omitted the rod has unit length. +%% +%% @var{rho} is a function handle to the density of the rod defined in the +%% interval 0,1. The integral of this density equals the mass and is stored in +%% @code{@var{J}(1,1)}. If omitted, the default is a uniform rod with unit mass. +%% +%% Run @code{demo rodmassmatrix} to see some examples. +%% +%% @end deftypefn + +function [J varargout] = rodmassmatrix(sigma, l = 1, dens = @(x)1) + + if ischar (sigma) + sigma = quadgk (@(x)x.*dens(x), 0,1); + end + + u = [-sigma*l (1-sigma)*l]; + + m = quadgk (@(x)dens(sigma+x/l), u(1),u(2))/l; + f = quadgk (@(x)dens(sigma+x/l).*x, u(1),u(2))/l; + iner_m = quadgk (@(x)dens(sigma+x/l).*x.^2, u(1),u(2))/l; + + J = [m f; f iner_m]; + + if nargout > 0 + varargout{1} = quadgk (@(x)x.*dens(x), 0,1); + end + +endfunction + +%!demo +%! barlen = 2; +%! [Jc, s] = rodmassmatrix (0, barlen); +%! +%! printf ("Inertia matrix from the extrema : \n") +%! disp (Jc) +%! printf ("Sigma value to calculate from center of mass : %g \n",s) +%! +%! J = rodmassmatrix (s); +%! printf ("Inertia matrix from the CoM : \n") +%! disp (J) +%! +%! J2 = rodmassmatrix ("com"); +%! tf = all((J2 == J)(:)); +%! disp (["Are J and J2 equal? " "no"*not(tf) "yes"*tf]) +%! +%! % ---------------------------------------------------------------------------- +%! % This example shows the calculations for rod of length 2. First we place one +%! % of its extrema in the origin. Then we use the value of sigma provided by +%! % the function to do the same calculation form the center of mass. + +%!demo +%! % A normalized density function +%! density = @(x) (0.5*ones(size(x)) + 10*(x<0.1)).*(x>=0 & x<=1)/1.5; +%! [Jc, s] = rodmassmatrix (0,1,density); +%! +%! printf ("Inertia matrix from the extrema : \n") +%! disp (Jc) +%! printf ("Sigma value to calculate from center of mass : %g \n",s) +%! J = rodmassmatrix (s,1,density); +%! +%! printf ("Inertia matrix from the CoM : \n") +%! disp (J) +%! +%! figure (1) +%! clf +%! x = linspace (0,1,100)'; +%! h = plot (x,density(x),'b-;density;'); +%! set (h,'linewidth',2) +%! axis tight +%! v = axis(); +%! hold on +%! h = plot ([s s],v([3 4]),'k-;CoM;'); +%! set (h, 'linewidth', 2); +%! hold off +%! axis auto +%! % ---------------------------------------------------------------------------- +%! % This example defines a density function with an accumulation of mass near +%! % one end of the rod. First we place one of its extrema in the origin. Then +%! % we use the value of sigma provided by the function to do the same +%! % calculation form the center of mass. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |