importGeometry pde not closed volume H I Beam Partial Differential Equation Toolbox

Hello everyone. I'd appreciate if you can give me a hand.

I'm creating a H beam STL file in matlab and trying to load it in pde tool but it says:

Error using pde.DiscreteGeometry

Failed to create geometry, the input does not form a closed volume.

This may be due to missing faces or gaps in the input.

The thing is that the stl file is for sure okay. I check the geometry with two different software (3d builder from windows) and AutodeskSimstudio Tools 2016 R3. Both shows a closed volume mesh (empty of course). Even Simstudio can create a full volume (not empty) thanks to the fact that the mesh is closed.

I attach some pictures of the STL and the STL itself (.zip).

Inside view:

Outside view:

View from Autodesk (to take a look of the triangles) :

The code I used here (I attach the other two functions to run it):

clc

clear all

close all

a1=2;

b1=1;

c1=2;

d1=0.5;

k=20;

a2=2;

b2=1;

c2=2;

d2=0.5;

nodes = [...

-a1/2 0 0; ...

-a1/2 0 b1; ...

-d1/2 0 b1; ...

-d1/2 0 b1+c1; ...

-a1/2 0 b1+c1; ...

-a1/2 0 b1+b1+c1; ...

a1/2 0 b1+b1+c1; ...

a1/2 0 b1+c1; ...

d1/2 0 b1+c1; ...

d1/2 0 b1; ...

a1/2 0 b1; ...

a1/2 0 0; ...

a2/2 k 0+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k 0+(b1+c1/2-c2/2-b2) ];

%edges = [1 2;2 3;3 4;4 5;5 6;6 7;7 8;8 9;9 10;10 11;11 12];

% faces are oriented such that normals point outwards

f1 = [1 12 11 2; 3 10 9 4; 5 8 7 6]; %root face

f2 = f1+12; %tip face

f3 = [12 13 14 11; 10 11 14 15; 9 10 15 16; 8 9 16 17; 7 8 17 18; 6 7 18 19; 5 6 19 20; 4 5 20 21; 3 4 21 22; 2 3 22 23 ; 1 2 23 24 ; 12 1 24 13];

faces = [f1;f2;f3];

tri = triangulateFaces(faces);

fv.faces = tri;

fv.vertices = nodes;

stlwrite('Hbeam.stl',fv);

smodel = createpde('structural','static-solid');

importGeometry(smodel, 'Hbeam.stl');

Thank you!

[Cris LaPierre]

What is the value of n?

fv.vertices = n;

Does n = nodes?

[jose daniel hoyos giraldo]

Yes, a mistake sorry

Answer 1

[Cris LaPierre]

A quick check with an STL checker (Materialize MiniMagics) found 16 bad edges, which prevents the STL from being 3D printed (what the software is for). It is likely the same issue that is resulting in the warning you are seeing.

When I repair the stl (using an online tool, fixed stl is attached), it imports without error.

unzip('Hbeam_fixed.zip')
smodel = createpde('structural','static-solid');
gm = importGeometry(smodel, 'Hbeam_fixed.stl');

pdegplot(gm)

[Cris LaPierre]

First, MATLAB has triangulation and stlwrite functions already. You don't need to write your own. However, using these I still get the same error.

It does seem that creating an intermediate stl is unnecessary if your end goal is to import it back into MATLAB to perform a structural analysis. You could just go straight to creating the beam using the tools already in MATLAB.

I will caution that I am no expert in this space, but after a little playing around, here is what I came up with.

% Define the H shape
x = [-10 -10 -2.5 -2.5 -10 -10 10 10 2.5 2.5 10 10];
y = [0 10 10 30 30 40 40 30 30 10 10 0];
% 3 in first row specifies rectangle. 
% 12 in second row says how many points (12 x, 12 y)
gm = [3;12;x';y'];

% Create a sketch of the H
model = createpde;
g = decsg(gm);
pg = geometryFromEdges(model,g);

pdegplot(model,'EdgeLabels','on','FaceLabels','on')
axis equal

% Now extrude the H the desired amount
Hbm = extrude(pg,200);
pdegplot(Hbm,'FaceLabels','on')

% Create a structural, static solid model from this geometry
Hbeam_model = createpde("structural","static-solid");
Hbeam_model.Geometry = Hbm;

You can now add your forces, conditions, material properties, etc and run your analysis using Hbeam_model.

[jose daniel hoyos giraldo]

It really works for me. Any suggestion to create a non constant shape? for instance:

I'd appreciate your help because I dont understand your code very well, but I found that "decsg" is for 2D geometry.

Thank you.

[Cris LaPierre]

I'm not sure of a simple way to do something like this in MATLAB. It would involve adding and subtracting multiple objects. However, if you have access to Autodesk, I would create the object there (I use Inventor) and then export as an stl and import that stl into MATLAB using importGeometry.

[jose daniel hoyos giraldo]

I'd like do it in matlab bc I'll create an automatic routine.I'm trying something like:

clc

clear all

close all

a1=2;

b1=1;

c1=2;

d1=0.5;

k=20;

a2=1;

b2=0.5;

c2=1;

d2=0.3;

nodes = [...

-a1/2 0 0; ...

-a1/2 0 b1; ...

-d1/2 0 b1; ...

-d1/2 0 b1+c1; ...

-a1/2 0 b1+c1; ...

-a1/2 0 b1+b1+c1; ...

a1/2 0 b1+b1+c1; ...

a1/2 0 b1+c1; ...

d1/2 0 b1+c1; ...

d1/2 0 b1; ...

a1/2 0 b1; ...

a1/2 0 0; ...

a2/2 k 0+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k 0+(b1+c1/2-c2/2-b2) ];

shp = alphaShape(nodes(:,1),nodes(:,3),nodes(:,2));

[elements,nodes] = boundaryFacets(shp);

nodes = nodes';

elements = elements';

model = createpde();

geometryFromMesh(model,nodes,elements);

figure(1)

pdegplot(model,'VertexLabels','on','EdgeLabels','on','CellLabels','on','FaceLabels','on','FaceAlpha',0.5);

But there is something wrong with it:

[Cris LaPierre]

What shape is this supposed to be?

You way something is wrong with it. Could you please elaborate on what that is?

[jose daniel hoyos giraldo]

Those nodes create a stl with the following shape:

Then the shape created is not what I'm trying.

Thanks

[Cris LaPierre]

What are the dimensions supposed to be?

[jose daniel hoyos giraldo]

20 units onf length (i dont rlly care about if takes as mm or m but the shape).

% root. Height is b1+c1+b1

% a1=2; Flange width at root

% b1=1; Flange thickness at root

% c1=2; web height at root

% d1=0.2; web thickness at root

% tip. Height is b2+c2+b2

% a2=1; Flange width at tip

% b2=0.5; Flange thickness at tip

% c2=1; web height at tip

% d2=0.1; web thickness at tip

% k=20; beam length

Thanks

[Cris LaPierre]

Ok, so compare that to the nodes you have created.

a1=2;
b1=1;
c1=2;
d1=0.5;
k=20;
a2=1;
b2=0.5;
c2=1;
d2=0.3;
nodes = [...
    -a1/2 0 0; ...
    -a1/2 0 b1; ...
    -d1/2 0 b1; ...
    -d1/2 0 b1+c1; ...
    -a1/2 0 b1+c1; ...
    -a1/2 0 b1+b1+c1; ...
    a1/2 0 b1+b1+c1; ...
    a1/2 0 b1+c1; ...
    d1/2 0 b1+c1; ...
    d1/2 0 b1; ...
    a1/2 0 b1; ...
    a1/2 0 0; ...    
    a2/2 k 0+(b1+c1/2-c2/2-b2); ...  
    a2/2 k b2+(b1+c1/2-c2/2-b2); ...   
    d2/2 k b2+(b1+c1/2-c2/2-b2); ...  
    d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...  
    a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...  
    a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...
    -a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...
    -a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...
    -d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...
    -d2/2 k b2+(b1+c1/2-c2/2-b2); ... 
    -a2/2 k b2+(b1+c1/2-c2/2-b2); ...
    -a2/2 k 0+(b1+c1/2-c2/2-b2)  ];

figure
scatter3(nodes(:,1),nodes(:,2),nodes(:,3),'filled')

[jose daniel hoyos giraldo]

Yea if you take a look those points are the nodes or corners of the H beam at tip and root, but it doesnt works as I expected with these functions in order to use geometryFromMesh of pde

shp = alphaShape(nodes(:,1),nodes(:,3),nodes(:,2));

[elements,nodes] = boundaryFacets(shp);

Ty!

[jose daniel hoyos giraldo]

I added more points but i'm dealing with similar problem:

clc

clear all

close all

a1=2;

b1=1;

c1=2;

d1=0.5;

k=10;

a2=1;

b2=0.5;

c2=1;

d2=0.3;

nodes = [...

-a1/2 0 0; ...

-a1/2 0 b1; ...

%-a1/4 - d1/4 0 b1; ...

-d1/2 0 b1; ...

-d1/2 0 b1+c1/4; ...

-d1/2 0 b1+c1/2; ...

-d1/2 0 b1+3*c1/4; ...

-d1/2 0 b1+c1; ...

-a1/2 0 b1+c1; ...

-a1/2 0 b1+b1+c1; ...

-a1/4 0 b1+b1+c1; ...

-0 0 b1+b1+c1; ...

a1/4 0 b1+b1+c1; ...

a1/2 0 b1+b1+c1; ...

a1/2 0 b1+c1; ...

d1/2 0 b1+c1; ...

d1/2 0 b1+3*c1/4; ...

d1/2 0 b1+c1/2; ...

d1/2 0 b1+c1/4; ...

d1/2 0 b1; ...

a1/2 0 b1; ...

a1/2 0 0; ...

a1/4 0 0; ...

0 0 0; ...

-a1/4 0 0; ...

a2/2 k 0+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+c2/4+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+c2/2+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+3*c2/4+(b1+c1/2-c2/2-b2); ...

d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

a2/4 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

0 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/4 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+c2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+3*c2/4+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+c2/2+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+c2/4+(b1+c1/2-c2/2-b2); ...

-d2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k b2+(b1+c1/2-c2/2-b2); ...

-a2/2 k 0+(b1+c1/2-c2/2-b2); ...

-a2/4 k 0+(b1+c1/2-c2/2-b2); ...

0 k 0+(b1+c1/2-c2/2-b2); ...

a2/4 k 0+(b1+c1/2-c2/2-b2)];

figure(1)

scatter3(nodes(:,1),nodes(:,2),nodes(:,3),'filled')

axis equal

shp = alphaShape(nodes(:,3)*100,nodes(:,1)*100,nodes(:,2)*100);

[elements,nodes] = boundaryFacets(shp);

nodes = nodes';

elements = elements';

model = createpde();

geometryFromMesh(model,nodes,elements);

figure(2)

%pdegplot(model,'VertexLabels','on','EdgeLabels','on','CellLabels','on','FaceLabels','on');

pdegplot(model);

[Cris LaPierre]

I can't find a way to complete the workflow in MATLAB, but feel free to keep trying. Here is an approach that creates the beam as a triangulation that can then be saved as an stl and then imported using importGeometry.

% Define parameters
a1=2;
b1=1;
c1=2;
d1=0.5;
k=20;
a2=1;
b2=0.5;
c2=1;
d2=0.3;

% Define points/vertices
nodes = [...
 -a1/2 0 -b1-c1/2; ...
 -a1/2 0 -c1/2; ...
 -d1/2 0 -c1/2; ...
 -d1/2 0 c1/2; ...
 -a1/2 0 c1/2; ...
 -a1/2 0 b1+c1/2; ...
 a1/2 0 b1+c1/2; ...
 a1/2 0 c1/2; ...
 d1/2 0 c1/2; ...
 d1/2 0 -c1/2; ...
 a1/2 0 -c1/2; ...
 a1/2 0 -b1-c1/2; ...
 -a2/2 k -b2-c2/2; ...
 -a2/2 k -c2/2; ...
 -d2/2 k -c2/2; ...
 -d2/2 k c2/2; ...
 -a2/2 k c2/2; ...
 -a2/2 k b2+c2/2; ...
 a2/2 k b2+c2/2; ...
 a2/2 k c2/2; ...
 d2/2 k c2/2; ...
 d2/2 k -c2/2; ...
 a2/2 k -c2/2; ...
 a2/2 k -b2-c2/2];

% Define triangulation connectivity
T = [1 2 3;
 1 3 12;
 3 10 12;
 10 11 12;
 3 4 10;
 4 9 10;
 4 5 6;
 4 6 7;
 4 7 9;
 7 8 9;
 13 14 15;
 13 15 24;
 15 22 24;
 22 23 24;
 15 22 16;
 16 21 22;
 16 17 18;
 16 18 19;
 16 19 21;
 19 20 21;

 1 13 12;
 12 13 24;
 1 13 2;
 2 13 14;
 2 14 3;
 3 14 15;
 3 4 15;
 4 15 16;
 4 5 16;
 5 16 17;
 5 6 17;
 6 17 18;
 6 7 18;
 7 18 19;
 7 19 8;
 8 19 20;
 8 9 20;
 9 20 21;
 9 10 21;
 10 21 22;
 10 11 22;
 11 22 23;
 11 12 23;
 12 23 24];

x=nodes(:,1);
y=nodes(:,2);
z=nodes(:,3);

% For reference only, nodel numbering
plot3(x,y,z,'k')
text(x,y,z, cellstr(num2str((1:length(nodes))')));
axis off

% Create triangulation
TR = triangulation(T,x,y,z);
figure
trisurf(TR)

% Save triangulation to stl
stlwrite(TR,"Hbeam2.stl")

% Not ideal, but temporary solution using stl works
model = createpde("structural","static-solid");
gm = importGeometry(model,"Hbeam2.stl");
figure
pdegplot(model,'FaceLabels','on')

[jose daniel hoyos giraldo]

Thank you, definitely it is a solution. I'll take a look of the code to uderstand it.

[jose daniel hoyos giraldo]

Could you take a look of this one? https://la.mathworks.com/matlabcentral/answers/1625710-momentum-as-boundary-fea-structuralboundaryload-partial-differential-equation-toolbox

Thanks