create_femesh_fromcells

 create FE mesh from cells
 
 Input:
     1. params_cells is a structure with
         a. 7 elements for spheres (cell_shape = 1):
             cell_shape
             ncell
             Rmin
             Rmax
             dmin
             dmax
             para_deform 
         b. 8 elements for cylinders (cell_shape = 2):
             cell_shape
             ncell
             Rmin
             Rmax
             dmin
             dmax
             para_deform
             Hcyl  

     2. fname_cells

     3. params_domain_geom is a structure with 3 elements:
         Rratio_IN
         include_ECS
         ECS_gap

     4. params_domain_femesh is a structure with 2 elements:
         Htetgen
         tetgen_cmd

     5. fname_tetgen
 
 Output: 
     fname_tetgen_femesh

0001 function [fname_tetgen_femesh] = ...
0002    create_femesh_fromcells(params_cells,fname_cells,params_domain_geom,params_domain_femesh,fname_tetgen)
0003 
0004 % create FE mesh from cells
0005 %
0006 % Input:
0007 %     1. params_cells is a structure with
0008 %         a. 7 elements for spheres (cell_shape = 1):
0009 %             cell_shape
0010 %             ncell
0011 %             Rmin
0012 %             Rmax
0013 %             dmin
0014 %             dmax
0015 %             para_deform
0016 %         b. 8 elements for cylinders (cell_shape = 2):
0017 %             cell_shape
0018 %             ncell
0019 %             Rmin
0020 %             Rmax
0021 %             dmin
0022 %             dmax
0023 %             para_deform
0024 %             Hcyl
0025 %
0026 %     2. fname_cells
0027 %
0028 %     3. params_domain_geom is a structure with 3 elements:
0029 %         Rratio_IN
0030 %         include_ECS
0031 %         ECS_gap
0032 %
0033 %     4. params_domain_femesh is a structure with 2 elements:
0034 %         Htetgen
0035 %         tetgen_cmd
0036 %
0037 %     5. fname_tetgen
0038 %
0039 % Output:
0040 %     fname_tetgen_femesh
0041 
0042 cell_shape = params_cells.cell_shape;
0043 hreq = params_domain_femesh.Htetgen;
0044 tetgen_cmd = params_domain_femesh.tetgen_cmd;
0045 include_ECS = params_domain_geom.include_ECS;
0046 ECS_gap = params_domain_geom.ECS_gap;
0047 Rratio_IN = params_domain_geom.Rratio_IN;
0048 
0049 ndim = 3;
0050 
0051 if (cell_shape == 2)
0052     [ncell,facets_cell,facets_labels_cell,nodes_cell,pt_in_cell,center,normal,Rcell,nslice_vec,nodes_ind_bottomring,nodes_ind_topring] ...
0053         = create_cylinders_geometry(fname_cells,Rratio_IN);    
0054 elseif (cell_shape == 1)
0055     [ncell,facets_cell,facets_labels_cell,nodes_cell,pt_in_cell,center,normal,Rcell] = ...
0056         create_ellipses_geometry(fname_cells,Rratio_IN);
0057 end
0058 
0059 facets = [];
0060 nodes = [];
0061 offset = 0;
0062 
0063 boundary_attribute = [];
0064 nboundary = 0;
0065 for icell = 1:ncell
0066     facets = [facets,facets_cell{icell}+offset];    
0067     nfacets = size(facets_cell{icell},2);
0068     nnodes = size(nodes_cell{icell},1);            
0069     boundary_attribute = [boundary_attribute;facets_labels_cell{icell}'+nboundary];    
0070     nboundary = nboundary + max(facets_labels_cell{icell}');
0071     nodes = [nodes;nodes_cell{icell}];
0072     offset = offset+nnodes;
0073 end
0074 
0075 
0076 nodes_allcell = nodes;
0077 nnodes_allcell = size(nodes_allcell,1);
0078 
0079 facets_allcell = facets;
0080 nfacets_allcell = size(facets_allcell,2);
0081 
0082 pt_in_allcell = pt_in_cell;
0083 bd_attrib_allcell = boundary_attribute;
0084 
0085 if (include_ECS == 2 & cell_shape == 2)
0086     theta = linspace(0,2*pi,21);
0087     thetavec = theta(1:end-1);
0088     xvec=cos(thetavec);
0089     yvec=sin(thetavec);
0090     
0091     nodes_all2d = [];
0092     for icell = 1:ncell                  
0093         center_one = center{icell};                   
0094         Rcell_one = Rcell{icell}*(1+ECS_gap);
0095         nodes_2d{icell}(:,1) = [xvec*Rcell_one(1,1)+center_one(1,1),center_one(1,1)]';
0096         nodes_2d{icell}(:,2) = [yvec*Rcell_one(1,1)+center_one(1,2),center_one(1,2)]';
0097         
0098         nodes_all2d = [nodes_all2d;nodes_2d{icell}];
0099         Rvec(icell) = Rcell_one(1,1);
0100     end
0101     
0102     Rmean = mean(Rvec(1:ncell));
0103     Rmin = min(Rvec(1:ncell));
0104     Rmax = max(Rvec(1:ncell));
0105     
0106     shp = alphaShape(nodes_all2d(:,1),nodes_all2d(:,2),Rmax*2);
0107     vv = area(shp);
0108     shp = alphaShape(nodes_all2d(:,1),nodes_all2d(:,2),Rmax,...
0109         'holeThreshold',vv/2,'regionThreshold',vv/2);
0110   
0111     if (numRegions(shp) > 1)
0112         disp('Can only have 1 region in ECS');
0113         stop
0114     end    
0115     [tri, xy] = shp.boundaryFacets();
0116     shp_bd = polyshape([xy(:,1)],[xy(:,2)]);
0117     clear shp_cell;
0118     for icell = 1:ncell
0119         ii = nodes_ind_topring{icell};
0120         nodes_topring{icell} = [nodes_cell{icell}(ii,1:3)];
0121          
0122         shp_cell(icell) = polyshape(nodes_topring{icell}(:,1),nodes_topring{icell}(:,2));
0123 
0124         ii = nodes_ind_bottomring{icell};
0125         nodes_bottomring{icell} = [nodes_cell{icell}(ii,1:3)];        
0126 
0127     end
0128     shp_ecs = subtract(shp_bd,union(shp_cell(1:ncell)));
0129 
0130     tri = shp_ecs.triangulation;
0131     % figure;
0132     % triplot(tri.ConnectivityList,tri.Points(:,1),tri.Points(:,2));
0133     % title('ECS');
0134     % axis equal;
0135     
0136     Pts_ind = zeros(size(tri.Points,1),2);
0137     for icell = 1:ncell
0138         ii = nodes_ind_topring{icell};              
0139         for jj = 1:length(ii)
0140             kk = find(abs(nodes_cell{icell}(ii(jj),1)-tri.Points(:,1))<1e-6 & abs(nodes_cell{icell}(ii(jj),2)-tri.Points(:,2))<1e-6);  
0141             Pts_ind(kk,1:2) = [icell,ii(jj)];
0142         end      
0143     end
0144     ii = (find(Pts_ind(:,1)==0));
0145     Pts_ind(ii,1) = ncell+1;
0146     nodes_ecs = tri.Points(ii,1:2);
0147     Pts_ind(ii,2) = 1:length(ii);
0148 
0149     Pts_ecs_top = zeros(size(nodes_ecs,1),3);
0150     Pts_ecs_top(:,1:2) = [nodes_ecs];
0151     Pts_ecs_top(:,3) = nodes_cell{1}(nodes_ind_topring{1}(1),3);
0152 
0153     Pts_ecs_bottom = zeros(size(nodes_ecs,1),3);
0154     Pts_ecs_bottom(:,1:2) = [nodes_ecs];
0155     Pts_ecs_bottom(:,3) = nodes_cell{1}(nodes_ind_bottomring{1}(1),3);
0156 
0157 
0158     nodes = [nodes_allcell;Pts_ecs_bottom;Pts_ecs_top];
0159 
0160     Cmat = tri.ConnectivityList;
0161     Pts_top = zeros(size(tri.Points,1),3);
0162     Pts_top(:,1:2) = [tri.Points];
0163     Pts_top(:,3) = nodes_cell{1}(nodes_ind_topring{1}(1),3);
0164 
0165     Pts_bottom = zeros(size(tri.Points,1),3);
0166     Pts_bottom(:,1:2) = [tri.Points];
0167     Pts_bottom(:,3) = nodes_cell{1}(nodes_ind_bottomring{1}(1),3);
0168 
0169     Pts_both = [Pts_bottom;Pts_top];
0170     nfacets = size(Pts_bottom,1);
0171     Cmat_both = [Cmat;Cmat+nfacets];    
0172 
0173     
0174     tol = 1e-6;
0175     npt = size(Pts_both,1);
0176     rep_ind = zeros(npt,1);
0177     for ipt = 1:npt
0178         kk = find(abs(Pts_both(ipt,1)-nodes(:,1))<=tol & abs(Pts_both(ipt,2)...
0179                 -nodes(:,2))<=tol & abs(Pts_both(ipt,3)-nodes(:,3))<=tol);
0180         rep_ind(ipt) = kk;
0181     end
0182     Cmat_new = zeros(size(Cmat_both));
0183     for ipt = 1:npt
0184         ii = find(Cmat_both==ipt);
0185         Cmat_new(ii) = rep_ind(ipt);        
0186     end    
0187     necs = size(nodes_ecs,1);
0188     [Cmat,Pts] = cylinder_connectivity(necs,necs);
0189     Cmat_out = Cmat+nnodes_allcell;
0190 
0191     
0192     nodes_box = [Pts_ecs_bottom; Pts_ecs_top];
0193     facets_box = [Cmat_new;Cmat_out']'-nnodes_allcell;
0194 
0195     pt_in_box = Pts_ecs_bottom(1,1:3);
0196     pt_in_box(1,3) = 0;
0197 
0198     bd_attrib_box = (nboundary+1)*ones(size(facets_box,2),1);
0199 
0200     nboundary = nboundary+1;
0201     
0202     nodes = [nodes_allcell;nodes_box];
0203     
0204     facets = [facets_allcell,facets_box+nnodes_allcell];
0205 
0206     bd_attrib = [bd_attrib_allcell; bd_attrib_box];    
0207   
0208 elseif (include_ECS == 2 & cell_shape == 1)
0209 
0210     for icell = 1:ncell
0211         Rvec(icell) = mean(Rcell(icell,1:3));
0212     end
0213     
0214     Rmean = mean(Rvec(1:ncell));
0215     Rmin = min(Rvec(1:ncell));
0216     Rmax = max(Rvec(1:ncell));
0217         
0218     x1min = min(nodes(:,1)); x1max = max(nodes(:,1));
0219     x2min = min(nodes(:,2)); x2max = max(nodes(:,2));
0220     x3min = min(nodes(:,3)); x3max = max(nodes(:,3));  
0221     
0222     x1gap = Rmax*3;
0223     x2gap = Rmax*3;
0224     x3gap = Rmax*3;
0225     
0226     x1min = x1min-x1gap;
0227     x1max = x1max+x1gap;
0228     x2min = x2min-x2gap;
0229     x2max = x2max+x2gap;
0230     x3min = x3min-x3gap;
0231     x3max = x3max+x3gap;
0232     
0233     x1vec = linspace(x1min,x1max,80);
0234     x2vec = linspace(x2min,x2max,80);
0235     x3vec = linspace(x3min,x3max,80);
0236     
0237     [X1,X2,X3] = ndgrid(x1vec,x2vec,x3vec);
0238     
0239     xgap = ECS_gap*Rmax;
0240     
0241     distmat = Inf*ones(size(X1));
0242     for icell = 1:ncell
0243         dd = max(0,sqrt((X1-center{icell}(1,1)).^2+(X2-center{icell}(1,2)).^2+...
0244             +(X3-center{icell}(1,3)).^2)-Rvec(icell));
0245         distmat = min(distmat,dd);
0246     end
0247  
0248     [facets_box,nodes_box] = isosurface(X1,X2,X3,distmat,xgap);
0249 
0250     ii = find(distmat<=xgap*0.9 & distmat>=xgap*0.1);
0251     if (~isempty(ii))
0252         pt_in_box = [X1(ii(1)),X2(ii(1)),X3(ii(1))];
0253     else
0254         disp(['did not find pt in']);
0255         stop
0256     end
0257     
0258     facets_box = facets_box';
0259 
0260     %figure; h=trisurf(facets_box',nodes_box(:,1),nodes_box(:,2),nodes_box(:,3)); set(h,'facealpha',0.1); view(3); axis equal;
0261  
0262     bd_attrib_box = (nboundary+1)*ones(size(facets_box,2),1);
0263     nboundary = nboundary+1;
0264     
0265     nodes = [nodes_allcell;nodes_box];
0266     facets = [facets_allcell,facets_box+nnodes_allcell];
0267     bd_attrib = [bd_attrib_allcell; bd_attrib_box];
0268 
0269 elseif (include_ECS == 1)
0270     x1min = min(nodes(:,1)); x1max = max(nodes(:,1));
0271     x2min = min(nodes(:,2)); x2max = max(nodes(:,2));
0272     x3min = min(nodes(:,3)); x3max = max(nodes(:,3));
0273     
0274     x1len = x1max-x1min;
0275     x2len = x2max-x2min;
0276     x3len = x3max-x3min;
0277     
0278     x1gap = x1len*ECS_gap;
0279     x2gap = x2len*ECS_gap;
0280     x3gap = x3len*ECS_gap;
0281     
0282     x1min = x1min-x1gap;
0283     x1max = x1max+x1gap;
0284     x2min = x2min-x2gap;
0285     x2max = x2max+x2gap;
0286     x3min = x3min-x3gap;
0287     x3max = x3max+x3gap;
0288     
0289     [nodes_box,facets_box,pt_in_box] = create_box_geometry(x1min,x1max,x2min,x2max,x3min,x3max);
0290     
0291     pt_in_box = [x1min+x1gap/2,x2min+x2gap/2,x3min+x3gap/2];
0292     
0293     bd_attrib_box = (nboundary+1)*ones(size(facets_box,2),1);
0294     nboundary = nboundary+1;
0295     
0296     nodes = [nodes_allcell;nodes_box];
0297     facets = [facets_allcell,facets_box+nnodes_allcell];
0298     bd_attrib = [bd_attrib_allcell; bd_attrib_box];
0299 else
0300     nodes = [nodes_allcell];
0301     facets = [facets_allcell];
0302     bd_attrib = [bd_attrib_allcell];
0303 end
0304 
0305 nfacets = size(facets,2);
0306 nnodes = size(nodes,1);
0307 
0308 holes = [];
0309 nholes = size(holes,1);
0310 
0311 regions = [pt_in_allcell];
0312 if (include_ECS ~= 0)
0313     regions = [regions;pt_in_box];
0314 end
0315     
0316 nregions = size(regions,1);
0317 
0318 nodesindex = [1:nnodes]';
0319 
0320 filename = [fname_tetgen,'.node'];
0321 fid = fopen(filename,'w');
0322 fprintf(fid, '%s\n', '# Part 1 - node list');
0323 fprintf(fid, '%s\n', '# node count, 3 dim, no attribute, no boundary marker');
0324 fprintf(fid, '%d %d %d %d\n', nnodes,ndim,0,0);
0325 fprintf(fid, '%s\n', '# Node index, node coordinates');
0326 if (fid ~= -1) 
0327   fprintf(fid, '%d %26.16f %26.16f %26.16f\n', [nodesindex,nodes]');
0328 end
0329 fclose(fid);
0330 
0331 filename = [fname_tetgen,'.poly'];
0332 fid = fopen(filename,'w');
0333 
0334 fprintf(fid, '%s\n', '# Part 1 - node list');
0335 fprintf(fid, '%s\n', '#  0 indicates the node list is stored in file .node');
0336 fprintf(fid, '%d\n', 0);
0337 fprintf(fid, '%s\n', '# Part 2 - facet list');
0338 fprintf(fid, '%s\n', '# facet count, yes boundary marker');
0339 fprintf(fid, '%d %d\n', nfacets,1);
0340 
0341 fprintf(fid, '%s\n', '# Node index, node coordinates');
0342 if (fid ~= -1) 
0343   for ii = 1:nfacets
0344     fprintf(fid, '%d %d %d\n', 1,0,bd_attrib(ii));  
0345     fprintf(fid, '%d ', 3);
0346     for jj = 1:3
0347       fprintf(fid, '%d ', facets(jj,ii));
0348     end
0349     fprintf(fid, '\n');
0350   end
0351 end
0352 
0353 
0354 
0355 %%%Part 3 - hole list
0356 %%%Holes in the volume are specified by identifying a point inside each hole.
0357 %%%  One line: <# of holes>
0358 %%%  Following lines list # of holes:
0359 %%%    <hole #> <x> <y> <z>
0360 %%%    ...
0361 
0362 fprintf(fid, '%s\n', ' # Part 3 - hole list');
0363 fprintf(fid, '%d\n', nholes);
0364 if (nholes>=1) 
0365     for ih = 1:nholes
0366         fprintf(fid, '%d %f %f %f\n', ih,holes(ih,1),holes(ih,2),holes(ih,3));
0367     end
0368 end
0369 
0370 
0371 fprintf(fid, '%s\n', ' # Part 4 - region list');
0372 fprintf(fid, '%d\n', nregions);
0373 if (nregions >= 1) 
0374     for ir = 1:nregions
0375         fprintf(fid, '%d %f %f %f %d %f\n', ir,regions(ir,1),regions(ir,2),regions(ir,3),ir,0.1);
0376     end
0377 end
0378 fclose(fid);
0379 
0380 if (hreq > 0) 
0381     tetgen_options = ['-pqA','a',num2str(hreq)];
0382 else
0383     tetgen_options = ['-pqA'];
0384 end
0385 disp(['Running command "',tetgen_cmd,' ',tetgen_options,' ',filename,'"']);
0386 disp(['*****Start Tetgen ']);
0387 system([tetgen_cmd,' ',tetgen_options,' ',filename]);
0388 disp(['*****End Tetgen']);
0389 
0390 fname_tetgen_femesh = [fname_tetgen,'.1'];