How to find the (p,e,t) matrices for a regular meshgrid mesh ?

I have a regular mesh of simple geometry. I want to calculate the matrices p (nodes), e (edges) and t (elements) such that :

• p (points, the mesh nodes) is a 2-by-Np matrix of nodes, where Np is the number of nodes in the mesh. Each column p(:,k) consists of the x-coordinate of point k in p(1,k) and the y-coordinate of point k in p(2,k).
• e (edges) is a 7-by-Ne matrix of edges, where Ne is the number of edges in the mesh. The mesh edges in e and the edges of the geometry have a one-to-one correspondence. The e matrix represents the discrete edges of the geometry in the same manner as the t matrix represents the discrete faces. Each column in the e matrix represents one edge.
• e(1,k) is the index of the first point in mesh edge k.
• e(2,k) is the index of the second point in mesh edge k.
• e(3,k) is the parameter value at the first point of edge k. The parameter value is related to the arc length along the geometric edge.
• e(4,k) is the parameter value at the second point of edge k.
• e(5,k) is the ID of the geometric edge containing the mesh edge. You can see edge IDs by using the command pdegplot(geom,'EdgeLabels','on').
• e(6,k) is the subdomain number on the left side of the edge. The direction along the edge is given by increasing parameter values. The subdomain 0 is the exterior of the geometry.
• e(7,k) is the subdomain number on the right side of the edge.
• t (triangles) is a 4-by-Nt matrix of triangles or a 7-by-Nt matrix of triangles, depending on whether you call generateMesh with the GeometricOrder name-value pair set to 'quadratic' or 'linear', respectively. initmesh creates only 'linear' elements, which have size 4-by-Nt. Nt is the number of triangles in the mesh. Each column of t contains the indices of the points in p that form the triangle. The exception is the last entry in the column, which is the subdomain number.

The geometry is defined as follows:

% Plate
Ep_plaque = 5e-3;
R_plaque = 20e-2;

% Coil
R_int_ind   = 5e-2;
R_ext_ind   = 15e-2;
Ep_ind      = 2e-3;
Airgap      = 2e-3; 

% Airbox 
R_air = 5*R_plaque;
H_air = 100*Ep_plaque;

% zone 1  : air
% zone 10 : plate 
% zone 100 : coil

% Plate
x1 = 0;             y1 = 0;
x2 =  R_plaque;     y2 = 0;
x3 =  R_plaque;     y3 =  Ep_plaque;
x4 = 0;             y4 =  Ep_plaque;

% Coil
x100 =  R_int_ind;         y100 = - Airgap- Ep_ind;
x101 =  R_ext_ind;         y101 = - Airgap- Ep_ind;
x102 =  R_ext_ind;         y102 = - Airgap;
x103 =  R_int_ind;         y103 = - Airgap;

% Air
x20 = 0;            y20 = - H_air/2;
x21 =  R_air;        y21 = - H_air/2;
x22 =  R_air;        y22 =   H_air/2;
x23 = 0;            y23 =   H_air/2;

g = [ 2 x1   x2   y1   y2   10  1   0   0   0;...
    2 x2   x3   y2   y3   10  1   0   0   0;...
    2 x3   x4   y3   y4   10  1   0   0   0;...  
    2 x4   x1   y4   y1   10  0   0   0   0;...  
    2 x100 x101 y100 y101 100 1   0   0   0;...
    2 x101 x102 y101 y102 100 1   0   0   0;...
    2 x102 x103 y102 y103 100 1   0   0   0;...   
    2 x103 x100 y103 y100 100 1   0   0   0;... 
    2 x20  x21  y20  y21  1   0   0   0   0;...
    2 x21  x22  y21  y22  1   0   0   0   0;...
    2 x22  x23  y22  y23  1   0   0   0   0;...
    2 x23  x4   y23  y4   1   0   0   0   0;...
    2 x1   x20  y1   y20  1   0   0   0   0;...
    ].';

The mesh is defined as follows:

Nx = 100; % Nombre de points en x
Ny = 100; % Nombre de points en y

x_grid = linspace(0, R_air, Nx);
y_grid = linspace(-H_air/2, H_air/2, Ny);

[X, Y] = meshgrid(x_grid, y_grid);

x_nodes = X(:);
y_nodes = Y(:);

figure;
pdegplot(g); 
hold on;
plot(X,Y,'r',X',Y','r')
xlabel('X-axis');
ylabel('Y-axis');
legend('Geometry', 'Mesh');

For the node matrix, it's quite simple:

p = [x_nodes.';y_nodes.']; 

How do you calculate the other two matrices e and t, given that this is a rectangular mesh and the t element matrix contains 5 rows rather than 4? (4 nodes + zone number)

Answer 1

[Suraj Kumar]

Hi hakim,

As per my understanding, you are facing issues while creating the matrices for elements and edges. To generate the p (nodes), t (elements), and e (edges) matrices for a rectangular mesh in MATLAB you can go through the following steps:

1. Define the node matrix (p):

% Define the node matrix
p = [x_nodes.'; y_nodes.'];

2. Generate the elements matrix(t):

Each grid cell is divided into two triangles using its four corner nodes, forming triangles (n1, n2, n3) and (n2, n3, n4). The matrix t has 5 rows: three for node indices, one placeholder, and one for the zone number.

% Fill the element matrix t
index = 1;
for i = 1:Nx-1
    for j = 1:Ny-1
        n1 = (j-1)*Nx + i;
        n2 = n1 + 1;
        n3 = n1 + Nx;
        n4 = n3 + 1;

        % Define two triangles for each rectangle
        t(:, index) = [n1; n2; n3; 0; 1];
        index = index + 1;

        t(:, index) = [n2; n3; n4; 0; 1];
        index = index + 1;
    end
end

3. Generate the edge matrix(e) :

Each triangle's three edges in t are identified, sorted, and checked against a “map” to avoid duplicates. Unique edges are added to e with node indices, zone numbers, and subdomain assumptions. e is trimmed to represent unique edges accurately.

edge_index = 1;
edge_map = containers.Map('KeyType', 'char’,    'ValueType', 'any');

for i = 1:size(t, 2)
    tri = t(1:3, i);
    zone = t(5, i);
    edges = [tri(1), tri(2); tri(2), tri(3); tri(3),    tri(1)];
    for j = 1:3
        edge = sort(edges(j, :));
        edge_key = sprintf('%d-%d', edge(1), edge(2));
        if ~isKey(edge_map, edge_key)
            e(1:2, edge_index) = edge;
            e(5, edge_index) = zone
            e(6, edge_index) = 0; % For exterior
            e(7, edge_index) = zone; % For interior
            edge_map(edge_key) = edge_index;
            edge_index = edge_index + 1;
        end
    end
end

For more information on the “matrices” and “map” related functions in MATLAB, kindly go through the following documentations:

• Matrices: https://www.mathworks.com/help/matlab/math/creating-and-concatenating-matrices.html
• Map: https://www.mathworks.com/help/matlab/ref/containers.map.html

Hope it helps!