Initial commit R2021a

Models/Bouncing_Ball/sm_contact_ball_plot1time.m

@@ -3,7 +3,7 @@
 %
 % <enter plot description here if desired>
 %
-% Copyright 2016-2020 The MathWorks, Inc.
+% Copyright 2016-2021 The MathWorks, Inc.
 
 % Reuse figure if it exists, else create new figure
 if ~exist('h1_sm_contact_ball', 'var') || ...

Models/Point_Cloud_StrWhl/sm_point_cloud_steering_wheel_adjust.m

@@ -0,0 +1,14 @@
+%% Read in Steering Wheel STL file
+% Copyright 2021 The MathWorks, Inc.
+steering_wheel_ptCloud = stlread('steering_wheel_orig.STL');
+
+%% Move Steering wheel to center of solid
+steering_wheel_ctr = ...
+    triangulation(...
+    steering_wheel_ptCloud.ConnectivityList,...
+    steering_wheel_ptCloud.Points-((max(steering_wheel_ptCloud.Points)-min(steering_wheel_ptCloud.Points))/2)-min(steering_wheel_ptCloud.Points));
+
+%% Write new steering wheel file
+stlwrite(steering_wheel_ctr,'steering_wheel_ctr.stl')
+
+

README.md

@@ -1,11 +1,12 @@
 # **Spatial Contact Force Block Examples in Simscape Multibody™**
-Copyright 2020 The MathWorks, Inc.
+Copyright 2021 The MathWorks, Inc.
 
 This repository contains a set of examples that demonstrate the 
 [Spatial Contact Force block](https://www.mathworks.com/help/physmod/sm/ref/spatialcontactforce.html).
 The examples cover:
 * **Solid-to-solid contact** between parameterized solids like bricks, spheres, and cylinders
 * **Contact between CAD geometry** including STEP files.
+* **[Point Cloud](https://www.mathworks.com/help/physmod/sm/ref/pointcloud.html) definition** for efficient contact modeling of complex shapes.
 
 Open the project file Spatial_Contact_Force_Examples.prj to get started.
 
@@ -19,3 +20,6 @@ To learn more about contact modeling with Simscape Multibody, please visit:
    * [Simscape Electrical™](https://www.mathworks.com/products/simscape-electrical.html)
    * [Simscape Fluids™](https://www.mathworks.com/products/simscape-fluids.html)
    * [Simscape Multibody™](https://www.mathworks.com/products/simscape-multibody.html)
+
+This submission uses the [S^2 Sampling Toolbox](https://www.mathworks.com/matlabcentral/fileexchange/37004-suite-of-functions-to-perform-uniform-sampling-of-a-sphere)
+Anton Semechko (2022). Suite of functions to perform uniform sampling of a sphere (https://github.com/AntonSemechko/S2-Sampling-Toolbox), GitHub. Retrieved March 24, 2022.

Scripts_Data/Point_Cloud/Point_Cloud_Data_Brick.m

@@ -0,0 +1,142 @@
+function ptcld = Point_Cloud_Data_Brick(x,y,z,pointset,varargin)
+%Point_Cloud_Data_Brick Produce point cloud for exterior surface of a brick
+%   [ptcld] = Point_Cloud_Data_Brick(x,y,z)
+%
+%   You can specify:
+%       x         length along (first column of ptcld)
+%       y         length along (second column of ptcld)
+%       z         length along (third column of ptcld)
+%       pointset  'full':    corners, edges, and faces
+%                 'corners': corners only
+%   Point cloud data will be centered at [0,0,0]
+%
+% Copyright 2021-2022 The MathWorks, Inc.
+
+% Default data to show diagram
+if (nargin == 0)
+    x = 3;
+    y = 2;
+    z = 1;
+    pointset = 'full';
+end
+
+% Check if plot should be produced
+if (isempty(varargin))
+    showplot = 'n';
+else
+    showplot = varargin;
+end
+
+ptcld_corners = [...
+    1  -1  -1;
+    1   1  -1;
+    -1  -1  -1;
+    -1   1  -1;
+    1  -1   1;
+    1   1   1;
+    -1  -1   1;
+    -1   1   1;
+    ];
+
+ptcld_face_ctr = [...
+    1   0   0;
+    -1   0   0;
+    0  -1   0;
+    0   1   0;
+    0   0  -1;
+    0   0   1;
+    ];
+
+ptcld_edges = [...
+    1  -1   0;
+    1   1   0;
+    -1  -1   0;
+    -1   1   0;
+    1   0   1;
+    -1   0   1;
+    1   0  -1;
+    -1   0  -1;
+    0  -1  -1;
+    0   1  -1;
+    0  -1   1;
+    0   1   1;
+    ];
+
+ptcld_faces = [...
+    1    0.5   0.5;
+    1    0.5  -0.5;
+    1   -0.5   0.5;
+    1   -0.5  -0.5;
+    -1    0.5   0.5;
+    -1    0.5  -0.5;
+    -1   -0.5   0.5;
+    -1   -0.5  -0.5;
+    0.5  1     0.5;
+    0.5  1    -0.5;
+    -0.5  1     0.5;
+    -0.5  1    -0.5;
+    0.5 -1     0.5;
+    0.5 -1    -0.5;
+    -0.5 -1     0.5;
+    -0.5 -1    -0.5;
+    0.5   0.5  1;
+    0.5  -0.5  1;
+    -0.5   0.5  1;
+    -0.5  -0.5  1;
+    0.5   0.5 -1;
+    0.5  -0.5 -1;
+    -0.5   0.5 -1;
+    -0.5  -0.5 -1;
+    ];
+
+switch pointset
+    case 'full'
+        ptcld = [...
+            ptcld_corners
+            ptcld_face_ctr
+            ptcld_edges
+            ptcld_faces].*[x y z]/2;
+    case 'corners'
+        ptcld = [...
+            ptcld_corners].*[x y z]/2;
+end
+
+% Plot diagram to show parameters and extrusion
+if (nargin == 0 || strcmpi(showplot,'plot'))
+
+    % Figure name
+    figString = ['h1_' mfilename];
+    % Only create a figure if no figure exists
+    figExist = 0;
+    fig_hExist = evalin('base',['exist(''' figString ''')']);
+    if (fig_hExist)
+        figExist = evalin('base',['ishandle(' figString ') && strcmp(get(' figString ', ''type''), ''figure'')']);
+    end
+    if ~figExist
+        fig_h = figure('Name',figString);
+        assignin('base',figString,fig_h);
+    else
+        fig_h = evalin('base',figString);
+    end
+    figure(fig_h)
+    clf(fig_h)
+
+    temp_colororder = get(gca,'defaultAxesColorOrder');
+
+    plot3(ptcld(:,1),ptcld(:,2),ptcld(:,3),'o','MarkerFaceColor',temp_colororder(2,:))
+    hold on
+
+    xlabel(['x = ' num2str(x)],'Color','r');
+    ylabel(['y = ' num2str(y)],'Color','g');
+    zlabel(['z = ' num2str(z)],'Color','b');
+
+    DT = delaunayTriangulation(ptcld);
+    [K,~] = convexHull(DT);
+    trisurf(K,DT.Points(:,1),DT.Points(:,2),DT.Points(:,3),'FaceColor',[0.6 0.6 0.9],'EdgeColor',[0.6 0.6 0.9],'FaceAlpha',0.3)
+
+    title(['[ptcld] = Point\_Cloud\_Data\_Brick(x, y, z, pointset);']);
+    hold off
+    box on
+    axis equal
+    grid off
+end

Scripts_Data/Point_Cloud/Point_Cloud_Data_Circle.m

@@ -0,0 +1,72 @@
+function ptcld = Point_Cloud_Data_Circle(radius,numpts,varargin)
+%Point_Cloud_Data_Circle Produce point cloud for area within a circle
+%   [ptcld] = Point_Cloud_Data_Circle(radius,numpts)
+%   
+%   You can specify:
+%       radius    circle radius
+%       numpts    number of points
+%
+%   Point cloud data will be centered at [0,0,0]
+%
+% Copyright 2021-2022 The MathWorks, Inc.
+
+% Default data to show diagram
+if (nargin == 0)
+    radius = 5;
+    numpts = 100;
+    showplot = 'plot';
+end
+
+% Check if plot should be produced
+if (isempty(varargin))
+    showplot = 'n';
+else
+    showplot = varargin;
+end
+
+ptcld=[0 0 0];
+golden_angle = pi * (3 - sqrt(5));
+
+for i = 1:numpts
+    theta = i*golden_angle;
+    r = sqrt(i)/sqrt(numpts);
+    ptcld(i,:) = [r*cos(theta) r*sin(theta) 0]*radius;
+end
+
+
+% Plot diagram to show parameters and extrusion
+if (nargin == 0 || strcmpi(showplot,'plot'))
+
+    % Figure name
+    figString = ['h1_' mfilename];
+    % Only create a figure if no figure exists
+    figExist = 0;
+    fig_hExist = evalin('base',['exist(''' figString ''')']);
+    if (fig_hExist)
+        figExist = evalin('base',['ishandle(' figString ') && strcmp(get(' figString ', ''type''), ''figure'')']);
+    end
+    if ~figExist
+        fig_h = figure('Name',figString);
+        assignin('base',figString,fig_h);
+    else
+        fig_h = evalin('base',figString);
+    end
+    figure(fig_h)
+    clf(fig_h)
+
+    temp_colororder = get(gca,'defaultAxesColorOrder');
+
+    plot3(ptcld(:,1),ptcld(:,2),ptcld(:,3),'o','MarkerFaceColor',temp_colororder(2,:))
+    hold on
+
+    angles = 0:0.1:2*pi;
+    plot3(sin(angles)*radius,cos(angles)*radius,zeros(size(angles)),'b')
+    plot3([0 radius],[0 0],[0 0],'b-d','MarkerFaceColor','b','LineWidth',1)
+    text(radius/2,0,0,'radius','Color','blue')
+
+    title(['[ptcld] = Point\_Cloud\_Data\_Circle(radius, height, numpts);']);
+    hold off
+    box on
+    axis equal
+    grid off
+end

Scripts_Data/Point_Cloud/Point_Cloud_Data_Cylinder.m

@@ -0,0 +1,121 @@
+function ptcld = Point_Cloud_Data_Cylinder(radius,height,numpts,caps,varargin)
+%Point_Cloud_Data_Cylinder Produce point cloud for exterior surface of a cylinder
+%   [ptcld] = Point_Cloud_Data_Cylinder(radius,height,numpts,caps)
+%
+%   You can specify:
+%       radius    Cylinder radius
+%       height    Cylinder height
+%       caps      Add points on ends (true/false)
+%       numpts    Estimated number of points about circumference
+%                 Points on cylinder ends will be in addition
+%                 to this estimate, placed at the same points/area
+%                 as on the cylinder circumference
+%
+%   Point cloud data will be centered at [0,0,0]
+%
+% Copyright 2022 The MathWorks, Inc.
+
+% Default data to show diagram
+if (nargin == 0)
+    radius = 5;
+    height = 5;
+    caps   = true;
+    numpts = 100;
+    showplot = 'plot';
+end
+
+% Check if plot should be produced
+if (isempty(varargin))
+    showplot = 'n';
+else
+    showplot = varargin;
+end
+
+% Calculate cylinder circumference
+cyl_cir = 2*pi*radius;
+
+% Starting point for search
+n_h = 2;
+n_c = 6;
+
+% Increase number of points along height and about circumference
+% until number of points meets or exceeds request
+while (n_h*n_c<=numpts)
+    n_h   = n_h+2;      % Always add 2 lines of points about circumference
+    box_h = height/n_h; % Calculate separate of points about circumference
+    n_c   =  floor(cyl_cir/(box_h*2));  % Calculate vertical point separation
+end
+
+% Calculate angle vectors for staggered rings about circumference
+ang_ca = linspace(0,2*pi-(2*pi/n_c),n_c);
+ang_cb = ang_ca+pi/n_c;
+
+% Assemble point cloud, ring by ring
+ptcld = [];
+for i = 0:n_h
+    if(rem(i,2)==0)
+        pts_x = cos(ang_ca)*radius;
+        pts_y = sin(ang_ca)*radius;
+    else
+        pts_x = cos(ang_cb)*radius;
+        pts_y = sin(ang_cb)*radius;
+    end
+    pts_z = i/n_h*height*ones(size(pts_x))-height/2;
+    ptcld = [ptcld;pts_x' pts_y' pts_z'];
+end
+
+if(caps)
+    pts_per_area = size(ptcld,1)/(cyl_cir*height);
+    area_end = pi*(radius-box_h/2)^2;
+    pts_end  = ceil(area_end*pts_per_area);
+
+    pts_end = Point_Cloud_Data_Circle(radius-box_h/2,pts_end);
+
+    pts_end(:,3) = height/2;
+    ptcld = [ptcld; pts_end];
+    pts_end(:,3) = -height/2;
+    ptcld = [ptcld; pts_end];
+end
+
+
+% Plot diagram to show parameters and extrusion
+if (nargin == 0 || strcmpi(showplot,'plot'))
+
+    % Figure name
+    figString = ['h1_' mfilename];
+    % Only create a figure if no figure exists
+    figExist = 0;
+    fig_hExist = evalin('base',['exist(''' figString ''')']);
+    if (fig_hExist)
+        figExist = evalin('base',['ishandle(' figString ') && strcmp(get(' figString ', ''type''), ''figure'')']);
+    end
+    if ~figExist
+        fig_h = figure('Name',figString);
+        assignin('base',figString,fig_h);
+    else
+        fig_h = evalin('base',figString);
+    end
+    figure(fig_h)
+    clf(fig_h)
+
+    temp_colororder = get(gca,'defaultAxesColorOrder');
+
+    plot3(ptcld(:,1),ptcld(:,2),ptcld(:,3),'o','MarkerFaceColor',temp_colororder(2,:))
+    hold on
+
+    plot3([0 0],[0 0],[-1 1]*height/2,'b-d','MarkerFaceColor','b','LineWidth',1)
+    text(radius*0.05,radius*0.05,0,'height','Color','blue')
+
+    plot3([0 radius],[0 0],[1 1]*height/2,'k-d','MarkerFaceColor','k','LineWidth',1)
+    text(radius/2,0,height*0.45,'radius','Color','black')
+
+    DT = delaunayTriangulation(ptcld);
+    [K,~] = convexHull(DT);
+    trisurf(K,DT.Points(:,1),DT.Points(:,2),DT.Points(:,3),'FaceColor',[0.6 0.6 0.9],'EdgeColor',[0.6 0.6 0.9],'FaceAlpha',0.3)
+
+    title(['[ptcld] = Point\_Cloud\_Data\_Cylinder(radius, height, caps, numpts);']);
+    hold off
+    box on
+    axis equal
+    grid off
+end