MATLAB: How to generate a mesh with quadratic tetrahedral elements using function ‘generateMesh’ whose boundary nodes follow the underlying curved geometry

Partial Differential Equation Toolbox

It is sought for a workflow that allows to generate a mesh with quadratic tetrahedral elements from a curved three-dimensional shape, e.g. a sphere, where all the corresponding boundary nodes lie on the sphere. However, the additional nodes that are generated due to the quadratic nature of the elements are located in the median between the corner vertices of the tetrahedral elements. The latter can be observed using the following code snippet,
 
Radius = 5;
sp = multisphere(Radius);
m1 = createpde;
m1.Geometry = sp;
Hmax = 2;
mesh = generateMesh(m1,'Hmax',Hmax);
pdeplot3D(m1);
hold on
for i = 1:length(m1.Mesh.Elements(:, 1))
    nodeID = m1.Mesh.Elements(i, 1);
    plot3(m1.Mesh.Nodes(1, nodeID), m1.Mesh.Nodes(2, nodeID), m1.Mesh.Nodes(3,nodeID), 'ro');
end
which results in the following mesh,

Note that the additional nodes due to the quadratic tetrahedral elements are located in the median between the corner vertices of the tetrahedral elements and not on the exact geometry, as it would be desirable for some applications.
How can I generate a mesh with quadratic tetrahedral elements using function 'generateMesh' whose boundary nodes follow the underlying curved geometry?

Best Answer

This feature is unavailable as per MATLAB R2020b, as the high order nature of the elements is not based on the fact that the boundary faces of these elements match the geometry itself, but that they rather offer more nodes for a quadratic order approximation of the unkown field, which is standard in the industry.
Herein it is presented a simplified demo based on a tetrahedral discretization of a sphere. The workflow provided below allows for modifying the coordinates of these additional nodes that are located on the corresponding boundary face of the tetrahedral elements, such that all the underlying boundary nodes lie exactly on the spherical boundary shape. Moreover, this demo only shows how the corresponding nodal coordinates of the of the first element in the mesh that happens to be on the boundary of the sphere can be accordingly modified, such that its underlying nodes lie on the sperical boundary. This demo can be extended with a loop over all the elements, whereby the corresponding boundary elements also need to be identified, before their boundary nodes can be adjusted and it is only meant to serve as a demonstrator.
 
%% Create a sphere and a high order mesh
Radius = 5;
sp = multisphere(Radius);
m1 = createpde;
m1.Geometry = sp;
Hmax = 2;
mesh = generateMesh(m1,'Hmax',Hmax);
pdeplot3D(m1);
% Plot the original mesh that results when using function 'generateMesh'
% hold on
% for i = 1:length(m1.Mesh.Elements(:, 1))
%     nodeID = m1.Mesh.Elements(i, 1);
%     plot3(m1.Mesh.Nodes(1, nodeID), m1.Mesh.Nodes(2, nodeID), m1.Mesh.Nodes(3,nodeID), 'ro');
% end
%% Modify the nodal locations according to the spherical shape
% define the center of the sphere
center = [0; 0; 0];
% get the elements and the nodes of the mesh
els = m1.Mesh.Elements;
nds = m1.Mesh.Nodes;
% initialize auxiliary variables
displVct = zeros(3, length(els(:,1)));
displ = zeros(length(els(:,1)), 1);
% loop over all the nodes of the first element
for i = 1:length(els(:,1))
  % get the nodal ID
  nodeID = els(i,1);
  % compute the displacement vector and its magnitude from the center of
  % the sphere
  distCenterVct = nds(:,nodeID);
  distCenter = norm(distCenterVct);
  % compute the spherical coordinates of the node
  [azimuth, elevation, ~] = ...
    cart2sph(nds(1,nodeID), nds(2,nodeID), nds(3,nodeID));
  % compute the position vector of the point on the exact sphere
  vec = [Radius*cos(elevation)*cos(azimuth)
    Radius*cos(elevation)*sin(azimuth)
    Radius*sin(elevation)];
  % compute the necessary displacement to bring that node to the exact
  % spherical shape
  displVct(:, i) = vec - distCenterVct;
  displ(i, 1) = norm(displVct(:, i));
end
% sort the nodes by their displacement to the exact sphere
[displSorted, idx] = sort(displ);
% The three first idx(1:3) closest nodes are the vertex nodes of the 
% tetrahedral that lie on the sphere (due to approximation errors the 
% corresponding displacements won't be exactly zero)
%

% The three ones afterwards idx(4:6) are the high order nodes that lie
% within the edges of the tetrahedral at the boundary face of the 
% tetrahedral
%
% The rest of the nodes idx(7:10) are interior to the tetrahedral and have
% no influence on the outer shape
nodesToBeMoved = idx(4:6);
% loop over all nodes whose nodal locations are to be updated
for i = 1:length(nodesToBeMoved)
  % get the nodal IDs
  nodeID = els(nodesToBeMoved(i, 1),1);
  % update the nodal locations by the necessary displacements to achieve
  % that they lie on the sphere
  nds(:, nodeID) = nds(:, nodeID) + displVct(:, nodesToBeMoved(i, 1));
end
%% create a new pde problem and assign the new elements and nodes with modified locations to it
m2 = createpde;
geometryFromMesh(m2,nds,els);
figure
pdeplot3D(m2);
hold on
for i = 1:length(m2.Mesh.Elements(:, 1))
    nodeID = m2.Mesh.Elements(i, 1);
    plot3(m2.Mesh.Nodes(1, nodeID), m2.Mesh.Nodes(2, nodeID), m2.Mesh.Nodes(3,nodeID), 'ro');
end
Running the aforementioned code snippet results in having the boundary nodes of the first element, that also happens to lie on the boundary, lying on the exact spherical shape, please see the following screenshot,
Please note that the tetrahedrals are herein only visualized by their 4 corner nodes and not by the additional ones that are added by the quadratic approximation.
Related Question