read_tetgen

 create FE mesh on canonical configuration; bend and twist the FE mesh nodes by analytical transformation
 
 Input:
     1. fname
     2. para_deform
     3. Ncmpt
     4. Nboundary
     
 Output:
     1. mymesh is a structure with 10 elements:
         Nnode
         Nele
         Nface
         Pts_cmpt_reorder
         Ele_cmpt_reorder
         Pts_ind
         Pts_boundary_reorder
         Fac_boundary_reorder
         Nboundary
         Ncmpt   
     2. cmpts_bdys_mat

0001 function [mymesh,cmpts_bdys_mat] = read_tetgen(fname,para_deform,Ncmpt,Nboundary)
0002 
0003 % create FE mesh on canonical configuration; bend and twist the FE mesh nodes by analytical transformation
0004 %
0005 % Input:
0006 %     1. fname
0007 %     2. para_deform
0008 %     3. Ncmpt
0009 %     4. Nboundary
0010 %
0011 % Output:
0012 %     1. mymesh is a structure with 10 elements:
0013 %         Nnode
0014 %         Nele
0015 %         Nface
0016 %         Pts_cmpt_reorder
0017 %         Ele_cmpt_reorder
0018 %         Pts_ind
0019 %         Pts_boundary_reorder
0020 %         Fac_boundary_reorder
0021 %         Nboundary
0022 %         Ncmpt
0023 %     2. cmpts_bdys_mat
0024 
0025 
0026 % [Pts_cmpt_reorder,Ele_cmpt_reorder,Pts_ind,Pts_boundary_reorder,Fac_boundary_reorder, Nboundary,Ncmpt]
0027 
0028 disp(['Reading from Tetgen FE mesh from ',fname]);
0029 
0030 filename = [fname,'.node'];
0031 fid = fopen(filename,'r');
0032 if (fid ~= -1)
0033     Nnode = fscanf(fid,'%d',[4,1]);
0034     Rnode = fscanf(fid, '%f', [4,inf]);
0035 end
0036 fclose(fid);
0037 mymesh.Nnode = Nnode;
0038 Pts_all = Rnode(2:end,:);
0039 
0040 filename = [fname,'.ele'];
0041 fid = fopen(filename,'r');
0042 if (fid ~= -1)
0043     Nele = fscanf(fid,'%d',[3,1]);
0044     Rele = fscanf(fid, '%f', [6,inf]);
0045 end
0046 fclose(fid);
0047 mymesh.Nele = Nele;
0048 
0049 Ele_all = Rele(2:5,:);
0050 Ele_attrib = Rele(6,:);
0051 
0052 filename = [fname,'.face'];
0053 fid = fopen(filename,'r');
0054 if (fid ~= -1)
0055     Nface = fscanf(fid,'%d',[2,1]);
0056     Rface = fscanf(fid, '%f', [5,inf]);
0057 end
0058 fclose(fid);
0059 mymesh.Nface = Nface;
0060 
0061 Fac_all = Rface(2:4,:);
0062 Fac_attrib = Rface(5,:);
0063 
0064 temp = [];
0065 Cmpt_attrib = unique(Ele_attrib);
0066 
0067 Ncmpt_tmp = length(Cmpt_attrib);
0068 
0069 disp(['Separating FE mesh into ',num2str(Ncmpt_tmp), ' compartments ']);
0070 
0071 if (Ncmpt_tmp ~= Ncmpt)
0072     disp(['FE mesh not good, use smaller hmax or change surf triangulation']);
0073     mymesh = [];
0074     cmpts_bdys_mat = [];
0075 else
0076 
0077 for icmpt = 1:Ncmpt
0078     jj = find(Ele_attrib == Cmpt_attrib(icmpt));
0079     Ele_cmpt{icmpt} = Ele_all(:,jj); 
0080     Pts_ind{icmpt} = unique(Ele_cmpt{icmpt}(:));
0081     Pts_cmpt_reorder{icmpt} = Pts_all(:,Pts_ind{icmpt});
0082 end
0083 
0084 if (~isempty(find(para_deform ~= 0)))
0085     for icmpt = 1:Ncmpt
0086         Pts_cmpt_reorder{icmpt} = deform_domain(Pts_cmpt_reorder{icmpt},para_deform);
0087     end
0088 end
0089 
0090 Boundary_attrib = unique(Fac_attrib);
0091 Nboundary_tmp = length(Boundary_attrib);
0092 
0093 disp(['Separating FE mesh with ',num2str(Nboundary_tmp), ' boundaries ']);
0094 
0095 if (Nboundary_tmp ~= Nboundary)
0096     disp(['FE mesh not good, use smaller hmax or change surf triangulation']);
0097     mymesh = [];
0098     cmpts_bdys_mat = [];
0099 else
0100 
0101 for iboundary = 1:Nboundary
0102     jj = find(Fac_attrib == Boundary_attrib(iboundary));
0103     Fac_boundary{iboundary} = Fac_all(:,jj);
0104     Pts_boundary{iboundary} = unique(Fac_all(:,jj));
0105 end
0106 
0107 for icmpt = 1:Ncmpt
0108     Ele_cmpt_reorder{icmpt} = Ele_cmpt{icmpt};
0109     for ii = 1:length(Pts_ind{icmpt})
0110         jj = find(Ele_cmpt{icmpt}==Pts_ind{icmpt}(ii));
0111         Ele_cmpt_reorder{icmpt}(jj) = ii;
0112     end
0113     for iboundary = 1:Nboundary
0114         onbd = 1;
0115         for ib = 1:size(Fac_boundary{iboundary},2)
0116             ind = find(Fac_boundary{iboundary}(ib)==Pts_ind{icmpt});
0117             if (isempty(ind))
0118                 onbd = 0;
0119             end
0120         end
0121         if (onbd == 1)
0122             Fac_boundary_reorder{icmpt}{iboundary} = Fac_boundary{iboundary};
0123             for ii = 1:length(Pts_ind{icmpt})
0124                 jj = find(Fac_boundary{iboundary}==Pts_ind{icmpt}(ii));
0125                 Fac_boundary_reorder{icmpt}{iboundary}(jj) = ii;
0126             end
0127             Pts_boundary_reorder{icmpt}{iboundary} = Pts_boundary{iboundary};
0128             for ii = 1:length(Pts_ind{icmpt})                
0129                 jj = find(Pts_boundary{iboundary}==Pts_ind{icmpt}(ii));
0130                 Pts_boundary_reorder{icmpt}{iboundary}(jj) = ii;                
0131             end
0132         else
0133             Fac_boundary_reorder{icmpt}{iboundary} = [];
0134             Pts_boundary_reorder{icmpt}{iboundary} = [];
0135         end
0136     end
0137 end
0138 
0139 mymesh.Pts_cmpt_reorder = Pts_cmpt_reorder;
0140 mymesh.Ele_cmpt_reorder = Ele_cmpt_reorder;
0141 mymesh.Pts_ind = Pts_ind;
0142 mymesh.Pts_boundary_reorder = Pts_boundary_reorder;
0143 mymesh.Fac_boundary_reorder = Fac_boundary_reorder;
0144 mymesh.Nboundary = Nboundary;
0145 mymesh.Ncmpt = Ncmpt;
0146 
0147 cmpts_bdys_mat = zeros(mymesh.Ncmpt,mymesh.Nboundary);
0148 for icmpt = 1:mymesh.Ncmpt
0149     for ibd = 1:mymesh.Nboundary
0150         cmpts_bdys_mat(icmpt,ibd)=~isempty(mymesh.Pts_boundary_reorder{icmpt}{ibd});
0151     end
0152 end
0153 % %%%%%%%%%%%%%%%%%%% Start checking the mesh quality %%%%%%%%%%%%%%%%%%%%%%%
0154 % aspect_ratio_lim = 0.05; % the worst -> [0, 1] <- the best
0155 % Qmesh=cell(1,mymesh.Ncmpt);
0156 % hmax=0;
0157 % for icmpt = 1:mymesh.Ncmpt
0158 %     Qmesh{icmpt} = mesh_quality(mymesh.Pts_cmpt_reorder{icmpt},mymesh.Ele_cmpt_reorder{icmpt});
0159 %     hmax = max(max(Qmesh{icmpt}.hout),hmax);
0160 % %    if Qmesh{icmpt}.quality1(1)<aspect_ratio_lim
0161 %         disp(['  Compartment ',num2str(icmpt),' - FE mesh with minimum aspect ratio of ',num2str(Qmesh{icmpt}.quality1(1),'%.1e')]);
0162 % %    end;
0163 % end;
0164 % %%%%%%%%%%%%%%%%%%% End of checking the mesh quality %%%%%%%%%%%%%%%%%%%%%%
0165 
0166 %%%%%%%%%%%%%%%%%%% Start checking the mesh quality %%%%%%%%%%%%%%%%%%%%%%%
0167 aspect_ratio_lim = 0.05; % the worst -> [0, 1] <- the best
0168 Qmesh = cell(1,mymesh.Ncmpt);
0169 hmax = 0;
0170 for icmpt = 1:mymesh.Ncmpt
0171     Qmesh{icmpt}=mesh_quality(mymesh.Pts_cmpt_reorder{icmpt},mymesh.Ele_cmpt_reorder{icmpt}) ;
0172     hmax = max(max(Qmesh{icmpt}.hout),hmax);
0173     disp(['  Compartment ',num2str(icmpt),' - FE mesh with minimum aspect ratio of ',num2str(Qmesh{icmpt}.quality(1),'%.1e')]);
0174 end
0175 %%%%%%%%%%%%%%%%%%% End of checking the mesh quality %%%%%%%%%%%%%%%%%%%%%%
0176 
0177 end
0178 end