Boundary Nodes, Elements and Edges: complex mesh boundaries

When the geometry of a mesh is not regular or its boundary is decomposed of many different parts, we can use the edges of the mesh to find the nodes that belongs to the boundary.

The edges corresponding to the boundary are the only ones that belongs only to a unique element. The internal ones are always sheared by two elements. Using this property, we have build the boundaryNodes function (at the end of this document) to find all the nodes and all the edges belonging to the boundary.

Example:

eval('meshDataHoleS5');

[numElem, ndim] = size(elem);

figure()

plotElements(nodes, elem, 0);

[indNodBd, indElemBd, indLocalEdgBd, edges] = boundaryNodes(nodes, elem);

hold on;

for i=1:length(indElemBd)

fill(nodes(elem(indElemBd(i),:),1),nodes(elem(indElemBd(i),:),2),'y')

end

plot(nodes(indNodBd,1), nodes(indNodBd,2),'rs');

hold off;

Plot only the corresponding edges of the boundary elements

figure()

for i=1:length(indElemBd)

k1=indLocalEdgBd(i);

k2=k1+1;

if (ndim==3)

if (k1==3)

k2=1;

end

elseif(ndim==4)

if (k1==4)

k2=1;

end

v1=nodes(elem(indElemBd(i),k1),:);

v2=nodes(elem(indElemBd(i),k2),:);

X=[v1(1); v2(1)];

Y=[v1(2); v2(2)];

plot(X,Y,'m')

hold on;

end

axis equal

Only one part of the boundary

We can mix the find function together with boundaryNodes function to obtain complementary boundaries using setdiff.

Example: Let's assume that at the boundary we have the same BC and it is only different on the left side of the boundary and in the circle.

We must have three different node list: indLeft, indCircle and indEqual (all the rest).

indLeft = find(nodes(:,1) < 0.001);

indCircle = find(((nodes(:,1)-4).^2+(nodes(:,2)-1).^2 - 0.5^2) < 0.01);

indEqual = setdiff(indNodBd,indLeft); %substract left nodes from all boundary nodes

indEqual = setdiff(indEqual,indCircle); %now substract circle nodes

figure()

plot(nodes(indLeft,1),nodes(indLeft,2),'ro'); %plot on the previous figure

hold on;

plot(nodes(indCircle,1),nodes(indCircle,2),'bo');

plot(nodes(indEqual,1),nodes(indEqual,2),'gs');

axis equal

hold off;

Exercise 1: (easy)

From the mesh file meshUTriang.m and using find and boundaryNodes, select all the boundary nodes according to the figure:

Hint: Total number of nodes in the complete boundary 76. Final number of selected nodes 60.

Exercise 2: (middle)

From the mesh file meshUTriang.m and using find and boundaryNodes, select all the boundary nodes according to the figure:

Hint: Final number of selected nodes 26.

Exercise 3: (difficult)

From the mesh file meshClipQuad.m and using find and boundaryNodes, select all the boundary nodes according to the figure:

function boundaryNodes: All nodes belonging to the boundary

function [indNodBd, indElemBd, indLocalEdgBd, edges] = boundaryNodes(nodes, elem)
% output:
% indNodBd  -->  list of all nodes on the boundary
% indElemBd -->  list of indices of the Elements on the boundary
% indLocalEdgBd  -->  list of indices of the local element edges on the boundary
% edges     -->  list of all edges defined in the total mesh
% 
    numNod  = size(nodes,1);
    [numElem, ndim] = size(elem);
    if (ndim == 3) %triangle edges
        edges = unique(sort([elem(:,[1,2]);elem(:,[2,3]);elem(:,[3,1])],2),'rows');
    elseif (ndim == 4) %quadrilateral edges
        edges = unique(sort([elem(:,[1,2]);elem(:,[2,3]);elem(:,[3,4]);elem(:,[4,1])],2),'rows');
    end    
    indNodBd=[];
    indLocalEdgBd=[];
    indElemBd=[];
    % look for the edges belonging only to one element
    for i=1:size(edges,1)
        n1=edges(i,1);
        n2=edges(i,2);
        [indRow,indCol]=find(elem == n1); %find elements owning the first node
        [indElem,col]=find(elem(indRow,:) == n2); % owing also the second one
        if (length(indElem) == 1) %boundary edges
            indNodBd=[indNodBd; n1; n2];
            indElemBd=[indElemBd; indRow(indElem)];
            lloc1=find(elem(indRow(indElem),:)==n1);
            lloc2=find(elem(indRow(indElem),:)==n2);
            if (ndim == 3) %triangle edges
                aux=[0,0,0];
                aux(lloc1)=1;
                aux(lloc2)=1;
                number = aux(1)+2*aux(2)+4*aux(3); 
                switch (number) %identify the appropriate element edge
                    case 3
                        edgeBd=1; 
                    case 5
                        edgeBd=3; 
                    case 6
                        edgeBd=2; 
                    otherwise, error('edge not allowed');
                end
            elseif (ndim == 4) %quadrilateral edges
                aux=[0,0,0,0];
                aux(lloc1)=1;
                aux(lloc2)=1;
                number = aux(1)+2*aux(2)+4*aux(3)+8*aux(4); 
                switch (number) %identify the appropriate element edge
                    case 3
                        edgeBd=1; 
                    case 6
                        edgeBd=2; 
                    case 9
                        edgeBd=4; 
                    case 12
                        edgeBd=3; 
                    otherwise, error('edge not allowed');
                end
            end
            indLocalEdgBd=[indLocalEdgBd; edgeBd];
        end
    end
    indNodBd=unique(indNodBd);
end

(c)Numerical Factory 2020