Hello,
I understand that you want to get the stress and strain for each mesh part of structuralmodel into the upper or the lower cell and classify the part into the cell it belongs to. I assume parts of mesh mean the Elements in a FEMesh object. Please find below the code to achieve the same.
in = 25.4e-3;
F0 = 100; % Assuming a value for F0
w0 = 1; % Assuming a value for w0
structuralModel = createpde('structural','static-solid');
cylinderGeom = multicylinder(25.4e-3,[0.25*in,0.25*in], "Zoffset",[0,0.25*25.4e-3]);
% Assign the geometry to the structural model
structuralModel.Geometry = cylinderGeom;
% Set the mesh size
generateMesh(structuralModel, 'Hmax', 2.0e-3);
E1 = 6.8e10;
sigma1 = 0.18;
E2 = 6.9e10;
sigma2 = 0.19;
structuralProperties(structuralModel, "Cell", 1, ...
'YoungsModulus', E1, ...
'PoissonsRatio', sigma1);
structuralProperties(structuralModel, "Cell", 2, ...
'YoungsModulus', E2, ...
'PoissonsRatio', sigma2);
structuralBC(structuralModel, 'Face', 1, 'Constraint', 'fixed'); % Fix back face of mirror
structFunForce = @(location,state) HGforce(location, state, 0, 0, F0, w0);
structuralBoundaryLoad(structuralModel, ...
'Face', 4, ...
'SurfaceTraction', structFunForce);
structuralResults = solve(structuralModel);
% Cell arrays to store szz component of stress and index of every mesh
% element based on whether they belong to upper cylindrical cell or lower
upperElements = {};
lowerElements = {};
% % Loop over each element and classify based on z-coordinate
numElements = size(structuralModel.Mesh.Elements, 2);
for i = 1:numElements
% For each element, find the corresponding node indices
elementNodes = structuralModel.Mesh.Elements(:, i);
% Find the coordinates of each node in the element
elementCoordinates = structuralModel.Mesh.Nodes(:, elementNodes);
% Find the szz component of stress for the element under consideration
stress_element_z = mean(structuralResults.Stress.szz(elementNodes));
% Find the centroid of the current element by taking mean
centroid = mean(elementCoordinates, 2);
% If the z-coordinate of the centroid is above 0.25*in, we classify the
% element to come from upper cell, and vice versa
zCoord = centroid(3);
if zCoord >= 0.25*in
upperElements{end+1} = {i, stress_element_z};
else
lowerElements{end+1} = {i, stress_element_z};
end
end
% Dummy HGForce
function force = HGforce(location, state, ~, ~, F0, w0)
force = [0; 0; F0 * cos(w0 * state.time)];
end
The above code finds the szz component of stress (FEStruct object) for each mesh element. We can similarly find the other components of stress and repeat the process for strain (FeStruct Object).
You can refer to the following documentations for more information on the code above:
- https://www.mathworks.com/help/pde/ug/pde.femesh-properties.html
- https://www.mathworks.com/help/pde/ug/pde.staticstructuralresults.html
Hope this helps.
