Why doesn't the subtract function work with 3D DiscreteGeometry in R2024b?

%  Why doesn't the subtract function work with 3D DiscreteGeometry in R2024b?

model = femodel();  % Create a Finite element analysis model object

% Define dimensions of the rectangular plate (width, height, thickness)
plateWidth = 10;
plateHeight = 6;
plateThickness = 1;

% Create the rectangular plate
rectPlate = multicuboid(plateWidth, plateHeight, plateThickness);

% Define dimensions of the circular hole (radius, height)
holeRadius = 1;
holeHeight = plateThickness;

% Create the cylinder for the hole
cylinderHole = multicylinder(holeRadius, holeHeight);

% Subtract the cylinder from the plate
plateWithHole = subtract(rectPlate, cylinderHole);

% The error message follows...

'subtract' is used in the following examples:

Outputs Defined for a Function

Outputs for Current Function

-----------------------------------------------------------

I have also tried to use:

plateWithHole = addVoid(rectPlate, cylinderHole);

and I get the following error message.

Error using pde.DiscreteGeometry/addVoid

Added cells must be located strictly inside a single cell of the original geometry.

In the error-generating case, the cylindrical void is the same thickness as the rectangular plate. If I make the cylindrical void ever so slightly shorter than the thickness of the plate, it works without error, but then I don't get a hole that goes completely through the plate.

-----------------------------------------------------------

What other Boolean functions will work with these DiscreteGeometry types that result in a hole going completely through the flat plate?

-----------------------------------------------------------

I’ve heard it said that subtract works in R2023, but I don't want to revert to that version.

Thank you,

Tim

[DGM]

The only subtract() function/method I can find is the one for 2D polyshape objects. I wouldn't expect it to work with other objects.

which subtract -all

/MATLAB/toolbox/matlab/polyfun/@polyshape/subtract.m  % polyshape method

I don't see that DiscreteGeometry objects have any subtract method, but I don't have PDE toolbox, so I'm not going to go digging under the hood.

That said, I have no idea what the error message is supposed to mean. When I run the same code, it doesn't give links to anything. None of the lines in the error message yield any useful search results (but most things don't these days). I'm probably missing something, but this isn't a very helpful error message.

% Subtract the cylinder from the plate
plateWithHole = subtract(rectPlate, cylinderHole);

'subtract' is used in the following examples:
Outputs Defined for a Function
Outputs for Current Function

[Paul]

If there was a "subtract" method for DiscreteGeometry objects in R2023 and there isn't now, perhaps that decision is discussed in the release notes for whatever version removed that method, perhaps pointing to an alternative method.

Don't know why the error message doesn't come with a link. But it looks like if "subtract" isn't found then Matlab does search through the documentation to find potential use-cases. I've never seen this type of message before (I've always seen "function not found" or something like that), perhaps it's a new feature.

Anyway, I did a search for "Outputs for Current Function" and found nargout, which does seem to be relevant to the error message but nothing the OP would care about (Outputs Defined for a Function is the other example at nargout).

[DGM]

I'm not sure, but the lack of links in the error message might just be a forum-only thing. I don't have the PDE toolbox, so I don't know.

Certainly, I found a lot of function() related docs links, but as much as those are relevant to function(), none are relevant to subtract().

Answer 1

[Timothy]

Walter et al,

I would like to thank all of you who helped me get to the answer for the problem I ran into. Your insights and suggestions really helped me to learn more about MATLAB and discover more of its nooks and crannies.

I will publish the complete Live Script for how I was able to solve the thermal FEA problem. It’s probably longer than the experts need, but it may be helpful to another non-expert like me.

Unfortunately for people like “DGM” I am using the PDE toolbox. There I generate a 2D shape and then within the Live Script I extrude it to the 3D block with holes that I need.

Again, thank you for your help.

-Tim

-----------------------------------------------------------------------------------------------------------------------------------------

Thermal FEA of a Three Dimensional (3-D) Object

Create the FEA Model Object

model = femodel();

Define Analysis Type

%model.AnalysisType = "thermalSteady";
model.AnalysisType = "thermalTransient";

Geometry 2-D: Get 2-D Geometry from PDE the Modeler App

%Note:  "g" can be found in the attached file.
load g_b
model.Geometry = g;  % This variable "g" comes from the PDE Modeler App using the "Boundary -> Export Decomposed Geometry..." command
pdegplot(g,"FaceLabels","on")

Geometry 3-D: Extrude 2-D Geometry to Make 3-D Object

%model = generateMesh(model,'Hmax',0.025,'GeometricOrder','quadratic','Hgrad',1.1);  % 'Hmax' defines the maximum element size
extrudeDistance = 0.1;
model.Geometry = extrude(model.Geometry,extrudeDistance);
pdegplot(model.Geometry,"FaceLabels","on"); % Now I have 3-D geometry that works in an FEA model

xlim([-0.138 0.941])
ylim([-0.219 0.281])
zlim([-0.231 0.100])
view([30.300 49.623])

Face_Hot  = 4;
Face_Cold = 6;

Generate Mesh:

model = generateMesh(model,'Hmax',0.05,'GeometricOrder','quadratic','Hgrad',1.1);  % 'Hmax' defines the maximum element size
pdemesh(model);  
xlim([-0.138 0.941])
ylim([-0.219 0.281])
zlim([-0.231 0.100])
view([30.300 49.623])

Define Physics: Material Properties:

E   = 210e9;    % Young's modulus: N/m^2
nu  = 0.3;      % Poisson's ratio: unitless
rho = 7850.0;   % Density:kg/m³
k   = 60.5;     % Thermal Conductivity: W/m/K
Cp  =  0.49;    % Specific Heat Capacity:  kJ/kg.K
cte = 12*10^-6; % Coefficient of Thermal Expansion : m/m/K

model.MaterialProperties = materialProperties(YoungsModulus=E, ...
    PoissonsRatio=nu, ...
    MassDensity=rho, ...
    CTE=cte, ...
    SpecificHeat=Cp, ...
    ThermalConductivity=k);

Define Physics: Boundary Conditions & Loads:

model.FaceLoad(1,1:max(model.Geometry.NumFaces))   = faceLoad(ConvectionCoefficient=5.0, AmbientTemperature=350); % All the other faces

model.FaceLoad(1, Face_Hot).ConvectionCoefficient   = [];
model.FaceLoad(1, Face_Hot).AmbientTemperature      = [];
model.FaceLoad(1, Face_Cold).ConvectionCoefficient  = [];
model.FaceLoad(1, Face_Cold).AmbientTemperature     = [];

model.FaceBC(Face_Hot)                              = faceBC(Temperature=500);  
model.FaceBC(Face_Cold)                             = faceBC(Temperature=200);  

model.CellIC = cellIC(Temperature=300);

Solve FEA Problem:

time_start  = -2;    % starting exponent (e.g., 10^-1)
time_end    = 1;     % ending exponent (e.g., 10^2) 
num_points  = 10;   % number of points in the vector 
tlist       = logspace(time_start, time_end, num_points);


%tlist = 0:0.1:10;  %Time steps
thermalresults = solve(model,tlist);  % Solve the problem for times in tlist.

%thermalresults = solve(model);  % Solve the problem.

Postprocess Results: Plotting the Temperature Results:

% Plot the solution at different time steps 
for i = 1:length(tlist) 
    pdeplot3D(thermalresults.Mesh.Nodes,thermalresults.Mesh.Elements,'ColorMapData',thermalresults.Temperature(:,i));
    xlim([-0.138 0.941])
    ylim([-0.219 0.281])
    zlim([-0.231 0.100])
    view([30.300 49.623])
    title(['Temperature at t = ', num2str(tlist(i)), ' seconds.  Steps to go: ',num2str(length(tlist) -i)]); 
    pause(0.001); % Pause to visualize each time step 
end

Results: Plotting the heat flow results:

QnHot       = evaluateHeatRate(thermalresults,Face=Face_Hot);  % if the problem is specified in SI, the units of Qn are Watts/
QnCold    = evaluateHeatRate(thermalresults,Face=Face_Cold);

semilogx(tlist,QnHot);
hold on;
semilogx(tlist,QnCold);
hold off;
% Set y-axis limits 
%ylim([-50 50]);
xlabel("Time")
ylabel("Heat Flow Rate")
grid on
legend(["QnHotFace", "QnColdFace"], "Position", [0.7077 0.6853 0.1893, 0.0774])

[DGM]

I appreciate a good post either way.

Answer 2

[Walter Roberson]

In R2023a, the methods available for DiscreteGeometry are

addCell addVertex cellEdges extrude facesAttachedToEdges nearestFace scale

addFace addVoid cellFaces faceEdges nearestEdge rotate translate

You can subtract() certain kinds of Antenna Toolbox shapes; https://www.mathworks.com/matlabcentral/answers/2092471-creating-solid-geometry-with-pde or some RF Toolbox shapes

For subtracting geometries with the PDE toolbox, you use decsg

[Paul]

But why isn't the shape.Box method found by which -all?

which subtract -all

/MATLAB/toolbox/matlab/polyfun/@polyshape/subtract.m  % polyshape method

box = shape.Box("Height",10 ,"Length",30,"Color","r","Width",10);

Now which finds the Shape.Box method because the namespace has now been imported?

which subtract -all

/MATLAB/toolbox/matlab/polyfun/@polyshape/subtract.m    % polyshape method
/MATLAB/toolbox/shared/em_cad/authoring/+shape/Shape.m  % shape.Box method

Why did we have to instantiate a Shape.Box for its subtract method to be found? Is that because the method is implemented in Shape.m and not subtract.m (as for polyshape)? Or maybe it has something to do with Shape.m being in a namespace? If so, does the same apply for ordinary m-functions in a namespace, i.e., an m-function in a namespace won't be fround by which until after its been used, or maybe not until after something in the namespace has been used?

[Timothy]

Is the addVoid method specifically designed so that it will not work if the thicknesses of the two items (e.g., rectPlate, cylinderHole in this case) are identical?

plateWithHole = addVoid(rectPlate, cylinderHole);

[Timothy]

This MATLAB Answer is where I got the idea that the subtract method worked in R2003b. It's the same answer referenced above.

https://www.mathworks.com/matlabcentral/answers/2092471-creating-solid-geometry-with-pde

The given solution is to a problem that is esentially the same as mine. Unfortunately, it doesn't work in R2024b Update 3.

[Stephen23]

"Unfortunately, it doesn't work in R2024b Update 3."

It seems to work on this forum:

% Define the start and end points
xStart = -13;
xEnd = 13;

% Number of cylinders to create
numCylinders = 5;

% Calculate the spacing between each cylinder
spacing = (xEnd - xStart) / (numCylinders - 1);

% Initialize an array to hold the Cylinder objects
cylinders = cell(1, numCylinders);

% Loop to create each Cylinder and set its center
for i = 1:numCylinders
    xCenter = xStart + (i-1) * spacing;
    cylinders{i} = shape.Cylinder("Radius", 1, "Height", 10, ...
        "Color", [0.9290, 0.6940, 0.1250], ...
        "Center", [xCenter, 0, 0]);
end

box = shape.Box("Height",10 ,"Length",30,"Color","r","Width",10);
show(box)

for i = 1:numCylinders
    box = subtract(box, cylinders{i});
end
show(box)

which subtract(box, cylinders{i});

/MATLAB/toolbox/shared/em_cad/authoring/+shape/Shape.m  % shape.Custom3D method

[Walter Roberson]

Note that this subtract() is of Antenna Toolbox shapes, not of DiscreteGeometry

[Timothy]

Walter,

Thank you very much. Your example surely does work.

Given that "subtract" does work in R2024b Update 3, what is wrong with the following code where "subtract" does not work? Is the problem related to using "multicuboid" and/or "multicylinder"?

%  Why doesn't the subtract function work with 3D DiscreteGeometry in R2024b?
model = femodel();  % Create a Finite element analysis model object
% Define dimensions of the rectangular plate (width, height, thickness)
plateWidth = 10;
plateHeight = 6;
plateThickness = 1;
% Create the rectangular plate
rectPlate = multicuboid(plateWidth, plateHeight, plateThickness);
% Define dimensions of the circular hole (radius, height)
holeRadius = 1;
holeHeight = plateThickness;
% Create the cylinder for the hole
cylinderHole = multicylinder(holeRadius, holeHeight);
% Subtract the cylinder from the plate
plateWithHole = subtract(rectPlate, cylinderHole);

I should make it clear that I am not dedicated to discreteGeometry. I just need a way to make 3D geometry for use in finite element projects.

-Tim

[Walter Roberson]

The example is by @Stephen23

[Walter Roberson]

%  Why doesn't the subtract function work with 3D DiscreteGeometry in R2024b?
model = femodel();  % Create a Finite element analysis model object
% Define dimensions of the rectangular plate (width, height, thickness)
plateWidth = 10;
plateHeight = 6;
plateThickness = 1;
% Create the rectangular plate
rectPlate = multicuboid(plateWidth, plateHeight, plateThickness);
% Define dimensions of the circular hole (radius, height)
holeRadius = 1;
holeHeight = plateThickness;
% Create the cylinder for the hole
cylinderHole = multicylinder(holeRadius, holeHeight);
% Subtract the cylinder from the plate
whos rectPlate cylinderHole

  Name              Size            Bytes  Class                   Attributes

  cylinderHole      1x1                 8  pde.DiscreteGeometry              
  rectPlate         1x1                 8  pde.DiscreteGeometry              

methods(rectPlate)

Methods for class pde.DiscreteGeometry:

addCell               addVertex             cellEdges             extrude               facesAttachedToEdges  nearestEdge           rotate                translate             
addFace               addVoid               cellFaces             faceEdges             mergeCells            nearestFace           scale                 

Call "methods('handle')" for methods of pde.DiscreteGeometry inherited from handle.

rectPlate is a DiscreteGeometry. subtract() is not defined for DiscreteGeometry. subtract() is defined for shape.* objects -- which are objects from the Antenna Toolbox, not from the PDE Toolbox.

[DGM]

I'm the tangential simpleton here, but decsg() sounds terribly obtuse for object-level CSG, especially considering that I wouldn't have found out about it even after an hour's worth of dredging the docs for a toolbox I don't have. I'm going to re-suggest a comment I deleted earlier.

Is it possible or practical with your expected workflow to develop the composite geometry as a mesh and then convert it to a DiscreteGeometry object afterwards? I have seen that GPToolbox (FEX) has some CSG tools which I haven't tested, and there are any number of external tools that could do this. Even external tools can be exploited programmatically via system().

https://www.mathworks.com/matlabcentral/answers/2169008-are-there-matlab-functions-or-scripts-to-perform-boolean-operations-on-triangulated-solids

I don't have PDE toolbox and don't use DiscreteGeometry for anything, so I don't know if this is helpful. I'm just throwing this out there, since none of this seems like it's a satisfying resolution to the problem.

[Timothy]

I will give your suggestions a try.