highordertet

 [newnode,newelem]=highordertet(node,elem)

 generate high-order straight-edge tetrahedral mesh from
 the 1st order tetrahedral mesh

 author: Qianqian Fang, <q.fang at neu.edu>

 input:
    node: list of nodes
    elem: list of elements (each row are indices of nodes of each element)
    order: optional, the order of the generated mesh; if missing, order=2

 output:
    newnode: all new edge-nodes on the output mesh
    newelem: the indices of the edge nodes for each original tet element

    currently, this function only supports order=2
    to combine the newnode/newelem with the old mesh, one should use

    elemfull=[elem(:,1:4) newelem+size(node,1)]; % 10-node element
    nodefull=[node;newnode];

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [newnode,newelem]=highordertet(node,elem,order)
0002 %
0003 % [newnode,newelem]=highordertet(node,elem)
0004 %
0005 % generate high-order straight-edge tetrahedral mesh from
0006 % the 1st order tetrahedral mesh
0007 %
0008 % author: Qianqian Fang, <q.fang at neu.edu>
0009 %
0010 % input:
0011 %    node: list of nodes
0012 %    elem: list of elements (each row are indices of nodes of each element)
0013 %    order: optional, the order of the generated mesh; if missing, order=2
0014 %
0015 % output:
0016 %    newnode: all new edge-nodes on the output mesh
0017 %    newelem: the indices of the edge nodes for each original tet element
0018 %
0019 %    currently, this function only supports order=2
0020 %    to combine the newnode/newelem with the old mesh, one should use
0021 %
0022 %    elemfull=[elem(:,1:4) newelem+size(node,1)]; % 10-node element
0023 %    nodefull=[node;newnode];
0024 %
0025 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0026 %
0027 
0028 if(nargin<3)
0029     order=2;
0030 end
0031 if(order>=3 || order<=1)
0032     error('currently this function only supports order=2');
0033 end
0034 [edges,idx,newelem]=uniqedges(elem(:,1:4));
0035 newnode=node(edges',1:3);
0036 newnode=reshape(newnode',[3,2,size(edges,1)]);
0037 newnode=squeeze(mean(permute(newnode,[3 2 1]),2));