[newnode,newelem,newface]=meshrefine(node,elem,face,opt)
refine a tetrahedral mesh by adding new nodes or constraints
author: Qianqian Fang, <q.fang at neu.edu>
input parameters:
node: existing tetrahedral mesh node list
elem: existing tetrahedral element list
face: (optional) existing tetrahedral mesh surface triangle list
opt: options for mesh refinement:
if opt is a Nx3 array, opt is treated as a list of new nodes to
be inserted into the mesh. the new nodes must be located on the
surface or inside the original mesh. external nodes are
discarded, unless the opt.extcmdopt is specified.
if opt is a vector with a length that equals to that of node,
it will be used to specify the desired edge-length at each node;
setting a node value to 0 will by-pass the refinement at this node
if opt is a vector with a length that equals to that of elem,
it will be used as the desired maximum element volume of each
tetrahedron; setting to 0 will by-pass the refinement of that element.
if opt is a struct, it can have the following fields:
opt.newnode: same as setting opt to an Nx3 array
opt.reratio: radius-edge ratio, by default, iso2mesh uses 1.414
opt.maxvol: maximum element volume
opt.sizefield: a vector specifying either the desired edge-length
at each node, or the maximum volume constraint within each
tetrahedron, see above for details.
opt.extcmdopt: by default, meshrefine can only insert nodes
that are inside the original mesh. if one prefers to insert
nodes that are outside of the original mesh, one can define
this parameter to specify the meshing option (for tetgen)
for the extended domain, i.e. the convex hull including
both the original and the external nodes. If not defined,
'-Y' option is used by default (prevent tetgen from
inserting new nodes on the surface).
opt.extlabel: when external nodes are inserted, the new elements
will be assigned with an element label to group them
together, by default, this label is 0, unless opt.extlabel
is given
opt.extcorelabel: when external nodes are inserted, par of the
new elements share the polyhedra between the inserted nodes,
these special elements will be marked by opt.extcorelabel,
otherwise the label will be set to -1
outputs:
newnode: node coordinates of the tetrahedral mesh
newelem: element list of the tetrahedral mesh
newface: mesh surface element list of the tetrahedral mesh
the last column denotes the boundary ID
examples:
[node,face,elem]=meshasphere([0 0 0],24,1,2);
elem(:,5)=1;
% inserting nodes that are inside the original mesh
innernodes=double([1 1 1; 2 2 2; 3 3 3]);
[newno,newel]=meshrefine(node,elem,innernodes);
all(ismember(round(innernodes*1e10)*1e-10,round(newno*1e10)*1e-10,'rows'))
plotmesh(newno,[],newel,'x>-3')
% inserting nodes that are external to the original mesh
extnodes=double([-5 -5 25;-5 5 25;5 5 25;5 -5 25]);
[newno,newel]=meshrefine(node,elem,struct('newnode',extnodes,'extcmdopt','-Y'));
figure;
plotmesh(newno,[],newel,'x>-3')
-- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)0001 function [newnode,newelem,newface]=meshrefine(node,elem,varargin) 0002 % 0003 % [newnode,newelem,newface]=meshrefine(node,elem,face,opt) 0004 % 0005 % refine a tetrahedral mesh by adding new nodes or constraints 0006 % 0007 % author: Qianqian Fang, <q.fang at neu.edu> 0008 % 0009 % input parameters: 0010 % node: existing tetrahedral mesh node list 0011 % elem: existing tetrahedral element list 0012 % face: (optional) existing tetrahedral mesh surface triangle list 0013 % opt: options for mesh refinement: 0014 % if opt is a Nx3 array, opt is treated as a list of new nodes to 0015 % be inserted into the mesh. the new nodes must be located on the 0016 % surface or inside the original mesh. external nodes are 0017 % discarded, unless the opt.extcmdopt is specified. 0018 % if opt is a vector with a length that equals to that of node, 0019 % it will be used to specify the desired edge-length at each node; 0020 % setting a node value to 0 will by-pass the refinement at this node 0021 % if opt is a vector with a length that equals to that of elem, 0022 % it will be used as the desired maximum element volume of each 0023 % tetrahedron; setting to 0 will by-pass the refinement of that element. 0024 % if opt is a struct, it can have the following fields: 0025 % opt.newnode: same as setting opt to an Nx3 array 0026 % opt.reratio: radius-edge ratio, by default, iso2mesh uses 1.414 0027 % opt.maxvol: maximum element volume 0028 % opt.sizefield: a vector specifying either the desired edge-length 0029 % at each node, or the maximum volume constraint within each 0030 % tetrahedron, see above for details. 0031 % opt.extcmdopt: by default, meshrefine can only insert nodes 0032 % that are inside the original mesh. if one prefers to insert 0033 % nodes that are outside of the original mesh, one can define 0034 % this parameter to specify the meshing option (for tetgen) 0035 % for the extended domain, i.e. the convex hull including 0036 % both the original and the external nodes. If not defined, 0037 % '-Y' option is used by default (prevent tetgen from 0038 % inserting new nodes on the surface). 0039 % opt.extlabel: when external nodes are inserted, the new elements 0040 % will be assigned with an element label to group them 0041 % together, by default, this label is 0, unless opt.extlabel 0042 % is given 0043 % opt.extcorelabel: when external nodes are inserted, par of the 0044 % new elements share the polyhedra between the inserted nodes, 0045 % these special elements will be marked by opt.extcorelabel, 0046 % otherwise the label will be set to -1 0047 % 0048 % outputs: 0049 % newnode: node coordinates of the tetrahedral mesh 0050 % newelem: element list of the tetrahedral mesh 0051 % newface: mesh surface element list of the tetrahedral mesh 0052 % the last column denotes the boundary ID 0053 % 0054 % examples: 0055 % 0056 % [node,face,elem]=meshasphere([0 0 0],24,1,2); 0057 % elem(:,5)=1; 0058 % 0059 % % inserting nodes that are inside the original mesh 0060 % innernodes=double([1 1 1; 2 2 2; 3 3 3]); 0061 % [newno,newel]=meshrefine(node,elem,innernodes); 0062 % all(ismember(round(innernodes*1e10)*1e-10,round(newno*1e10)*1e-10,'rows')) 0063 % plotmesh(newno,[],newel,'x>-3') 0064 % 0065 % % inserting nodes that are external to the original mesh 0066 % extnodes=double([-5 -5 25;-5 5 25;5 5 25;5 -5 25]); 0067 % [newno,newel]=meshrefine(node,elem,struct('newnode',extnodes,'extcmdopt','-Y')); 0068 % figure; 0069 % plotmesh(newno,[],newel,'x>-3') 0070 % 0071 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net) 0072 % 0073 0074 exesuff=getexeext; 0075 exesuff=fallbackexeext(exesuff,'tetgen'); 0076 0077 newpt=[]; 0078 sizefield=[]; 0079 if(size(node,2)==4) 0080 sizefield=node(:,4); 0081 node=node(:,1:3); 0082 end 0083 opt=struct; 0084 if(length(varargin)==1) 0085 face=[]; 0086 if(isstruct(varargin{1})) 0087 opt=varargin{1}; 0088 else 0089 if(length(varargin{1})==size(node,1) || length(varargin{1})==size(elem,1)) 0090 sizefield=varargin{1}; 0091 else 0092 newpt=varargin{1}; 0093 end 0094 end 0095 elseif(length(varargin)>=2) 0096 face=varargin{1}; 0097 if(isstruct(varargin{2})) 0098 opt=varargin{2}; 0099 else 0100 if(length(varargin{2})==size(node,1) || length(varargin{1})==size(elem,1)) 0101 sizefield=varargin{2}; 0102 else 0103 newpt=varargin{2}; 0104 end 0105 end 0106 else 0107 error('meshrefine requires at least 3 inputs'); 0108 end 0109 if(isstruct(opt) && isfield(opt,'newnode')) 0110 newpt=opt.newnode; 0111 end 0112 if(isstruct(opt) && isfield(opt,'sizefield')) 0113 sizefield=opt.sizefield; 0114 end 0115 0116 % call tetgen to create volumetric mesh 0117 deletemeshfile(mwpath('pre_refine.*')); 0118 deletemeshfile(mwpath('post_refine.*')); 0119 0120 moreopt=''; 0121 setquality=0; 0122 if(isstruct(opt) && isfield(opt,'reratio')) 0123 moreopt=[moreopt sprintf(' -q %.10f ',opt.reratio)]; 0124 setquality=1; 0125 end 0126 if(isstruct(opt) && isfield(opt,'maxvol')) 0127 moreopt=[moreopt sprintf(' -a %.10f ',opt.maxvol)]; 0128 end 0129 0130 externalpt=[]; 0131 if(isstruct(opt) && isfield(opt,'extcmdopt')) 0132 isinside=tsearchn(node(:,1:3),elem(:,1:4),newpt(:,1:3)); 0133 externalpt=newpt(isnan(isinside),:); 0134 newpt=newpt(find(isnan(isinside)==0),:); 0135 end 0136 0137 if(~isempty(newpt)) 0138 if(size(newpt,1)<4) 0139 newpt=[newpt; repmat(newpt(1,:),4-size(newpt,1),1)]; 0140 end 0141 savetetgennode(newpt,mwpath('pre_refine.1.a.node')); 0142 moreopt=' -i '; 0143 end 0144 0145 if(~isempty(sizefield)) 0146 if(length(sizefield)==size(node,1)) 0147 fid=fopen(mwpath('pre_refine.1.mtr'),'wt'); 0148 fprintf(fid,'%d 1\n',size(sizefield,1)); 0149 fprintf(fid,'%.16f\n',sizefield); 0150 fclose(fid); 0151 moreopt=[moreopt ' -qa ']; 0152 else 0153 fid=fopen(mwpath('pre_refine.1.vol'),'wt'); 0154 fprintf(fid,'%d\n',size(sizefield,1)); 0155 fprintf(fid,'%d\t%.16f\n',[(1:size(sizefield,1))' sizefield]'); 0156 fclose(fid); 0157 moreopt=[moreopt ' -qa ']; 0158 end 0159 end 0160 0161 if(size(elem,2)==3 && setquality==0) 0162 if(~isempty(newpt)) 0163 error('inserting new point can not be used for surfaces'); 0164 end 0165 nedge=savegts(node, elem,mwpath('pre_refine.gts')); 0166 exesuff=fallbackexeext(getexeext,'gtsrefine'); 0167 elseif(size(elem,2)==3) 0168 savesurfpoly(node,elem,[],[],[],[],mwpath('pre_refine.poly')); 0169 else 0170 savetetgennode(node, mwpath('pre_refine.1.node')); 0171 savetetgenele (elem, mwpath('pre_refine.1.ele')); 0172 end 0173 0174 fprintf(1,'refining the input mesh ...\n'); 0175 0176 if(size(elem,2)==3 && setquality==0) 0177 if(isstruct(opt) && isfield(opt,'scale')) 0178 moreopt=sprintf('%s -n %d ',moreopt,round(nedge*opt.scale)); 0179 else 0180 error('you must give opt.scale value for refining a surface'); 0181 end 0182 end 0183 if(isstruct(opt) && isfield(opt,'moreopt')) 0184 moreopt=[moreopt opt.moreopt]; 0185 end 0186 0187 if(size(elem,2)==3 && setquality==0) 0188 system([' "' mcpath('gtsrefine') exesuff '" ' moreopt ' < "' ... 0189 mwpath('pre_refine.gts') '" > "' mwpath('post_refine.gts') '"']); 0190 [newnode,newelem]=readgts(mwpath('post_refine.gts')); 0191 newface=newelem; 0192 elseif(size(elem,2)==3) 0193 system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -p -A "' mwpath('pre_refine.poly') '"']); 0194 [newnode,newelem,newface]=readtetgen(mwpath('pre_refine.1')); 0195 elseif(~isempty(moreopt)) 0196 system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -r "' mwpath('pre_refine.1') '"']); 0197 [newnode,newelem,newface]=readtetgen(mwpath('pre_refine.2')); 0198 else 0199 newnode=node; 0200 newelem=elem; 0201 newface=face; 0202 end 0203 0204 if(~isempty(externalpt)) % user request to insert nodes that are outside of the original mesh 0205 % create a mesh including the external points 0206 externalpt=unique(externalpt,'rows'); 0207 allnode=[newnode;externalpt]; 0208 0209 % define the convexhull as the external surface 0210 try 0211 outface=convhull(allnode,'simplify',false); 0212 catch 0213 outface=convhulln(allnode); 0214 end 0215 outface=sort(outface,2); 0216 face=volface(newelem(:,1:4)); 0217 inface=sort(face(:,1:3),2); 0218 0219 % define the surface that bounds the newly extended convex hull space 0220 bothsides=removedupelem([outface;inface]); 0221 0222 % define a seed point to avoid meshing the interior space 0223 holelist=surfseeds(node,face(:,1:3)); 0224 0225 % mesh the extended space 0226 ISO2MESH_TETGENOPT=jsonopt('extcmdopt','-Y',opt); 0227 if(size(bothsides,1)>=size(inface,1)) 0228 [no,el]=surf2mesh(allnode,bothsides,[],[],1,10,[],holelist); 0229 else 0230 [no,el]=surf2mesh(allnode,bothsides,[],[],1,10); 0231 end 0232 0233 [isinside,map]=ismember(round(no*1e10)*1e-10,round(allnode*1e10)*1e-10,'rows'); 0234 snid=[length(newnode)+1:length(allnode)]; 0235 0236 if(~isempty(map==0)) 0237 oldsize=size(allnode,1); 0238 allnode=[allnode;no(map==0,:)]; 0239 map(map==0)=oldsize+1:size(allnode,1); 0240 end 0241 % merge the external space with the original mesh 0242 el2=map(el(:,1:4)); 0243 0244 % label all new elements with -1 0245 if(size(newelem,2)==5) 0246 el2(:,5)=jsonopt('extlabel',0,opt); 0247 % search elements that contain source(s) and save their id(s) 0248 iselm=ismember(el2(:,1:4),snid); 0249 el2(sum(iselm,2)>=3,5)=jsonopt('extcorelabel',-1,opt); 0250 end 0251 0252 % merge nodes/elements and replace the original ones 0253 newnode=allnode; 0254 newelem=[newelem;el2]; 0255 end 0256 0257 % read in the generated mesh 0258 0259 fprintf(1,'mesh refinement is complete\n'); 0260