How to create a geometry consisting of a prism, with an inner inclined plane, for a diffusion pde problem?

For example, the nodes could be placed at (0, 0 ,0), (1, 0 ,0), (1, 1 ,0), (0, 1 ,0), (0, 0 ,1), (1, 0 ,1), (1, 1 ,1), (0, 1 ,1). The plane could go from the upper left edge to the lower right edge, so its nodes would be (0, 0 ,1), (0, 1 ,0), (1, 0, 0), (1, 1, 0).

Everytime I tried

model = createpde(); geometryFromMesh(model,nodes,elements); or importGeometry(model,'PrismAndPlane.stl');

I have gotten all kinds of error messages (no matter how I generated the nodes and elements), including

'The input mesh does not form a closed volume. This may be due to missing faces or gaps in the input.' and 'Failed to create geometry. The stl file is invalid, more than two facets share an edge.'

I cannot include specific code, because I have literally tried a dozen different solutions (alpha shapes, geometry from meshes, create and import stl files, etc). The only one that worked provided a double inner plane, which is not ok to solve a diffusion pde.

So take it as a quiz. Glory for the one who solves it!

Thanks in advance,

Gabriel

[Gabriel Madirolas]

Ah, I also have gotten many times the error message

'The input mesh is invalid. Neighboring elements in the mesh are not properly connected.'

[Alan Weiss]

Can you please sketch the geometry you want, specifying holes and solid portions? I suspect that, with the geometry you seem to be specifying, the resulting solid is not one contiguous piece, which is the reason for the errors you observe. But a sketch would help confirm or disprove that suspicion.

Alan Weiss

MATLAB mathematical toolbox documentation

[Gabriel Madirolas]

Thanks Alan,

The code below, for example, generates two contiguous regions:

% Define the two volumes (wedges) that will form the geometry

xUpperWedge = [0 0 1 1 1 1]';
yUpperWedge = [0 1 0 1 0 1]';
zUpperWedge = [1 1 1 1 0 0]';

xLowerWedge = [0 0 0 0 1 1]';
yLowerWedge = [0 1 0 1 0 1]';
zLowerWedge = [1 1 0 0 0 0]';

% Create an alpha shape for each wedge,
% and obtain a surface mesh from each

upperWedge = alphaShape(xUpperWedge,yUpperWedge,zUpperWedge);
[upperElements,upperNodes] = boundaryFacets(upperWedge);

lowerWedge = alphaShape(xLowerWedge,yLowerWedge,zLowerWedge);
[lowerElements,lowerNodes] = boundaryFacets(lowerWedge);
lowerElements = lowerElements +6; % this is to concatenate with upper 

% Put the data in the correct shape for geometryFromMesh.

upperElements = upperElements';
upperNodes = upperNodes';
lowerElements = lowerElements';
lowerNodes = lowerNodes';

elements = [upperElements lowerElements];
nodes = [upperNodes lowerNodes];

% Create a PDE model and import the surface mesh.

plate = createpde();
geometryFromMesh(plate,nodes,elements);

% Plot for visualzation purposes

figure
plot(upperWedge,'facecolor','g','facealpha',0.3)
hold on
plot(lowerWedge,'facecolor','y','facealpha',0.3)

% View the geometry and face numbers.

figure
pdegplot(plate,'FaceLabels','on','FaceAlpha',0.5)

I want to study diffusion from one to the other. But the problem in this case is that I get a double inner face (labeled F3 and F9, at least when I run it in my computer). I would like the boundary between the two regions to be a single face, because otherwise the pde solver won't consider both regions to be connected (as far as I know).

I found an example of a multiple region in the documentation of geometryFromMesh:

load MultidomainMesh3D
model = createpde;
geometryFromMesh(model,nodes,elements,ElementIdToRegionId);
pdegplot(model,'CellLabels','on','facelabels','on','facealpha',0.5)

This looks similar, but I would want the inner plane to be inclined. Unfortunately, it is not specified how the MultidomainMesh3D geometry was created, so I cannot custom my own geometry inspired by this one.

Again, I have tried many things, including create and import an stl file, but nothing worked.

Thanks in advance,

Gabriel Madirolas

[Simone Galante]

Hello Gabriel,

Have you been able to solve your problem? Have you found the method to merge your inner face F3 and F9? I have found a similar problem with the edges about two cyilindical solids that I show you in the attached immage file. Thanks in advance, Simone

[Gabriel Madirolas]

Hi Simone,

Unfortunately, I couldn't manage a solution to my problem. I think Matlab's pde is still not working smoothly enough with multi-region 3D geometries.

What I finally did is to write my own pde solver code, using finite differences method.

Gabriel Madirolas

Answer 1

[Konstantin Kovalev]

Hello Gabriel,

Please find below a script that does the job. Feel free to comment back should you have any questions.

Best Regards,

Konstantin

close all
clear variables

% columns 5 and 6 belong to upper wedge only, 7 and 8 - lower
x = [0 0 1 1 1 1 0 0];
y = [0 1 0 1 0 1 0 1];
z = [1 1 0 0 1 1 0 0];

% create delaunay triangulation for upper wedge
w1 = delaunayTriangulation(x(1:6)',y(1:6)',z(1:6)');

% copy elements (the copy will be lower wedge)
e1 = w1.ConnectivityList;
e2 = e1;

% update nodes 5-6 with nodes 7-8 (lower elements are inverted, V<0)
e2(e2==5) = 7;
e2(e2==6) = 8;

% combine both wedges' elements
% (1st two nodes of lower elements are swapped to undo inversion)
e = [e1; e2(:,[2 1 3:end])]';

% domain IDs - elements 1-3 are upper wedge, 4-6 lower
d = [1 1 1 2 2 2]';

pdem = createpde;
g = geometryFromMesh(pdem, [x; y; z], e, d);
pdegplot(pdem,'FaceLabels','on','CellLabels','on','FaceAlpha',0.75);

return

Answer 2

[Gabriel Madirolas]

Hi Konstantin,

Sorry for the late reply. Thanks a lot for the code. It works nicely!

Best regards,

Gabriel