Generating a 2D Composite Geometry with Subdomains for PDE Analysis in MATLAB

I am trying to create a composite geometry for a PDE problem in MATLAB. My geometry consists of two regions:

1. A general rectangular domain.
2. A subdomain with an arbitrary shape (even more complex than sawtooth function), which is entirely contained within the rectangular domain.

I need to combine these two geometries into a PDEModel while preserving both domains as distinct subregions so that I can apply separate boundary conditions or specifyCoefficients to each "face". I also need to generate a high-quality finite element mesh that respects both regions. I attempted to use polyshape because I could generate complex geometries, not just squares, circles like decsg. I put an example generating both geometries and trying to combine using subtraction (same with union), but just one region is generate.

Here the example:

close all; clear; clc;

%% Step 1: Define Custom Geometry (Subdomain)
% here I am using sawtooth, but in my real problem is more complex geometry 
T = 10 * (1 / 50); % Period of the sawtooth function
fs = 1000; % Sampling frequency
t = 0:1/fs:T-1/fs; % Time vector
fun = sawtooth(2 * pi * 50 * t); % Sawtooth function
ini_tickness = 1 * ones(size(fun)); % Base thickness
in_Domain = ini_tickness + normalize(fun, "range"); % Thickness profile normalized to range [1, 2]

%% Step 2: Define General Geometry (Domain)
Lx = 40; % Length of the rectangular domain [mm]
Ly = 20; % Height of the rectangular domain [mm]
pgon = polyshape([0, Lx, Lx, 0], [Ly, Ly, 0, 0]); % Rectangular domain

figure(1);
plot(pgon);
axis equal;
title('Rectangular Domain');

%% Step 3: Define the Subdomain Geometry 
x_holo = linspace(0, Ly - 1, length(in_Domain)); % Adjust x-coordinates to fit domain
y_holo = in_Domain; % Variable thickness as a function of x

% Create the internal shape shape as a polyshape
holo_polygon = polyshape([x_holo, fliplr(x_holo)], [zeros(size(x_holo)), fliplr(y_holo)]);

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

holo_polygon = rotate(holo_polygon, -90, [20, 0]); % Rotate the polyshape
holo_polygon = translate(holo_polygon, -19, -0.5); % Translate to fit inside the rectangular domain

figure(2);
plot(holo_polygon);
axis equal;
title('internal Geometry');

%% Step 4: Combine the Geometries
% Subtract the hologram from the rectangular domain
% also with union, etc I'll obtain just one "face" and I need two.
polyout = subtract(pgon, holo_polygon);

figure(3);
plot(polyout);
daspect([1, 1, 1]);
title('Composite Geometry');

%% Step 5: Generate Mesh and Define PDEModel
% Attempt to create a triangulation from the composite geometry
tr = triangulation(polyout);
polyout = boundaryshape(tr); % Convert to a boundaryshape object

% Define PDEModel and import the geometry
pdem = createpde;
pdem.geometryFromMesh(tr.Points', tr.ConnectivityList');

figure(4);
subplot(2, 1, 1);
pdegplot(pdem, 'EdgeLabels', "on");
title('Edges of the Geometry');

subplot(2, 1, 2);
pdegplot(pdem, 'FaceLabels', "on");
title('Faces of the Geometry');

%% Step 6: Generate Mesh % Please notice that there is only one face
generateMesh(pdem, "GeometricOrder", "linear");
figure(5);
pdeplot(pdem);
title('Generated Mesh');

Answer 1

[Divyam]

Hi @Darby Paez,

The "polyshape" method doesn't directly support multiple regions in the form of subdomains; it typically combines them into a single region, and hence it isn't possible to define your polygons as distinct regions.

Since you already tried the "decsg" method and found it incompatible with your use case, you could try defining a function handle for the location of the sawtooth model that will add the necessary properties you need to the polygon based on its location while treating the rectangle and the sawtooth as a single face.

% Pseudocode
applyBoundaryCondition(model, "dirichlet", Edge = 1, r = @myboundaryfun);
function bcMatrix = myboundaryfun(location, state)
    % add code here, can use the inpolygon method to check whether a point
    % lies in the sawtooth polygon or the rectangle
end

For more information regarding the use of function handles to define properties, refer to this documentation: https://www.mathworks.com/help/pde/ug/steps-to-specify-a-boundary-conditions-object.html?#buftxz0-1:~:text=Q%2Cg%3DG)%3B-,Nonconstant%20Boundary%20Conditions,-Use%20functions%20to

As an alternative, you can try out the PDE Modeler app to interact with the polygons and then define boundary conditions specifically for the selected polygon: https://www.mathworks.com/help/pde/ug/pdemodeler-app.html