0001 function PLOT_SURFACE_TRIANGULATION(cell_shape,fname_cells,params_domain_geom)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 include_ECS = params_domain_geom.include_ECS;
0013 ECS_gap = params_domain_geom.ECS_gap;
0014 Rratio_IN = params_domain_geom.Rratio_IN;
0015 ndim = 3;
0016
0017 if (cell_shape == 2)
0018 [ncell,facets_cell,facets_labels_cell,nodes_cell,pt_in_cell,center,normal,Rcell,nslice_vec,nodes_ind_bottomring,nodes_ind_topring] ...
0019 = create_cylinders_geometry(fname_cells,Rratio_IN);
0020 elseif (cell_shape == 1)
0021 [ncell,facets_cell,facets_labels_cell,nodes_cell,pt_in_cell,center,normal,Rcell] = ...
0022 create_ellipses_geometry(fname_cells,Rratio_IN);
0023 end
0024
0025 facets = [];
0026 nodes = [];
0027 offset = 0;
0028
0029 boundary_attribute = [];
0030 nboundary = 0;
0031 for icell = 1:ncell
0032 facets = [facets,facets_cell{icell}+offset];
0033 nfacets = size(facets_cell{icell},2);
0034 nnodes = size(nodes_cell{icell},1);
0035 boundary_attribute = [boundary_attribute;facets_labels_cell{icell}'+nboundary];
0036 nboundary = nboundary + max(facets_labels_cell{icell}');
0037 nodes = [nodes;nodes_cell{icell}];
0038 offset = offset+nnodes;
0039 end
0040
0041 nodes_allcell = nodes;
0042 nnodes_allcell = size(nodes_allcell,1);
0043
0044 facets_allcell = facets;
0045 nfacets_allcell = size(facets_allcell,2);
0046
0047 pt_in_allcell = pt_in_cell;
0048 bd_attrib_allcell = boundary_attribute;
0049
0050 if (include_ECS == 2 & cell_shape == 2)
0051
0052 theta = linspace(0,2*pi,21);
0053 thetavec = theta(1:end-1);
0054 xvec=cos(thetavec);
0055 yvec=sin(thetavec);
0056
0057 nodes_all2d = [];
0058 for icell = 1:ncell
0059 center_one = center{icell};
0060 Rcell_one = Rcell{icell}*(1+ECS_gap);
0061 nodes_2d{icell}(:,1) = [xvec*Rcell_one(1,1)+center_one(1,1),center_one(1,1)]';
0062 nodes_2d{icell}(:,2) = [yvec*Rcell_one(1,1)+center_one(1,2),center_one(1,2)]';
0063
0064 nodes_all2d = [nodes_all2d;nodes_2d{icell}];
0065 Rvec(icell) = Rcell_one(1,1);
0066 end
0067
0068 Rmean = mean(Rvec(1:ncell));
0069 Rmin = min(Rvec(1:ncell));
0070 Rmax = max(Rvec(1:ncell));
0071
0072 shp = alphaShape(nodes_all2d(:,1),nodes_all2d(:,2),Rmax*2);
0073 vv = area(shp);
0074 shp = alphaShape(nodes_all2d(:,1),nodes_all2d(:,2),Rmax,...
0075 'holeThreshold',vv/2,'regionThreshold',vv/2);
0076
0077 if (numRegions(shp) > 1)
0078 disp('Can only have 1 region in ECS');
0079 stop
0080 end
0081 [tri, xy] = shp.boundaryFacets();
0082 shp_bd = polyshape([xy(:,1)],[xy(:,2)]);
0083 clear shp_cell;
0084 for icell = 1:ncell
0085 ii = nodes_ind_topring{icell};
0086 nodes_topring{icell} = [nodes_cell{icell}(ii,1:3)];
0087
0088 shp_cell(icell) = polyshape(nodes_topring{icell}(:,1),nodes_topring{icell}(:,2));
0089
0090 ii = nodes_ind_bottomring{icell};
0091 nodes_bottomring{icell} = [nodes_cell{icell}(ii,1:3)];
0092
0093 end
0094 shp_ecs = subtract(shp_bd,union(shp_cell(1:ncell)));
0095
0096 tri = shp_ecs.triangulation;
0097
0098
0099
0100
0101
0102 Pts_ind = zeros(size(tri.Points,1),2);
0103 for icell = 1:ncell
0104 ii = nodes_ind_topring{icell};
0105 for jj = 1:length(ii)
0106 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);
0107 Pts_ind(kk,1:2) = [icell,ii(jj)];
0108 end
0109 end
0110 ii = (find(Pts_ind(:,1)==0));
0111 Pts_ind(ii,1) = ncell+1;
0112 nodes_ecs = tri.Points(ii,1:2);
0113 Pts_ind(ii,2) = 1:length(ii);
0114
0115 Pts_ecs_top = zeros(size(nodes_ecs,1),3);
0116 Pts_ecs_top(:,1:2) = [nodes_ecs];
0117 Pts_ecs_top(:,3) = nodes_cell{1}(nodes_ind_topring{1}(1),3);
0118
0119 Pts_ecs_bottom = zeros(size(nodes_ecs,1),3);
0120 Pts_ecs_bottom(:,1:2) = [nodes_ecs];
0121 Pts_ecs_bottom(:,3) = nodes_cell{1}(nodes_ind_bottomring{1}(1),3);
0122
0123
0124 nodes = [nodes_allcell;Pts_ecs_bottom;Pts_ecs_top];
0125
0126 Cmat = tri.ConnectivityList;
0127 Pts_top = zeros(size(tri.Points,1),3);
0128 Pts_top(:,1:2) = [tri.Points];
0129 Pts_top(:,3) = nodes_cell{1}(nodes_ind_topring{1}(1),3);
0130
0131 Pts_bottom = zeros(size(tri.Points,1),3);
0132 Pts_bottom(:,1:2) = [tri.Points];
0133 Pts_bottom(:,3) = nodes_cell{1}(nodes_ind_bottomring{1}(1),3);
0134
0135 Pts_both = [Pts_bottom;Pts_top];
0136 nfacets = size(Pts_bottom,1);
0137 Cmat_both = [Cmat;Cmat+nfacets];
0138
0139
0140 tol = 1e-6;
0141 npt = size(Pts_both,1);
0142 rep_ind = zeros(npt,1);
0143 for ipt = 1:npt
0144 kk = find(abs(Pts_both(ipt,1)-nodes(:,1))<=tol & abs(Pts_both(ipt,2)...
0145 -nodes(:,2))<=tol & abs(Pts_both(ipt,3)-nodes(:,3))<=tol);
0146 rep_ind(ipt) = kk;
0147 end
0148 Cmat_new = zeros(size(Cmat_both));
0149 for ipt = 1:npt
0150 ii = find(Cmat_both==ipt);
0151 Cmat_new(ii) = rep_ind(ipt);
0152 end
0153 necs = size(nodes_ecs,1);
0154 [Cmat,Pts] = cylinder_connectivity(necs,necs);
0155 Cmat_out = Cmat+nnodes_allcell;
0156
0157
0158 nodes_box = [Pts_ecs_bottom; Pts_ecs_top];
0159 facets_box = [Cmat_new;Cmat_out']'-nnodes_allcell;
0160
0161 pt_in_box = Pts_ecs_bottom(1,1:3);
0162 pt_in_box(1,3) = 0;
0163
0164 bd_attrib_box = (nboundary+1)*ones(size(facets_box,2),1);
0165
0166 nboundary = nboundary+1;
0167
0168 nodes = [nodes_allcell;nodes_box];
0169
0170 facets = [facets_allcell,facets_box+nnodes_allcell];
0171
0172 bd_attrib = [bd_attrib_allcell; bd_attrib_box];
0173
0174 elseif (include_ECS == 2 & cell_shape == 1)
0175
0176 for icell = 1:ncell
0177 Rvec(icell) = mean(Rcell(icell,1:3));
0178 end
0179
0180 Rmean = mean(Rvec(1:ncell));
0181 Rmin = min(Rvec(1:ncell));
0182 Rmax = max(Rvec(1:ncell));
0183
0184 x1min = min(nodes(:,1)); x1max = max(nodes(:,1));
0185 x2min = min(nodes(:,2)); x2max = max(nodes(:,2));
0186 x3min = min(nodes(:,3)); x3max = max(nodes(:,3));
0187
0188 x1gap = Rmax*3;
0189 x2gap = Rmax*3;
0190 x3gap = Rmax*3;
0191
0192 x1min = x1min-x1gap;
0193 x1max = x1max+x1gap;
0194 x2min = x2min-x2gap;
0195 x2max = x2max+x2gap;
0196 x3min = x3min-x3gap;
0197 x3max = x3max+x3gap;
0198
0199 x1vec = linspace(x1min,x1max,80);
0200 x2vec = linspace(x2min,x2max,80);
0201 x3vec = linspace(x3min,x3max,80);
0202
0203 [X1,X2,X3] = ndgrid(x1vec,x2vec,x3vec);
0204
0205 xgap = ECS_gap*Rmax;
0206
0207 distmat = Inf*ones(size(X1));
0208 for icell = 1:ncell
0209 dd = max(0,sqrt((X1-center{icell}(1,1)).^2+(X2-center{icell}(1,2)).^2+...
0210 +(X3-center{icell}(1,3)).^2)-Rvec(icell));
0211 distmat = min(distmat,dd);
0212 end
0213
0214
0215
0216 [facets_box,nodes_box] = isosurface(X1,X2,X3,distmat,xgap);
0217
0218
0219 ii = find(distmat<=xgap*0.9 & distmat>=xgap*0.1);
0220 if (~isempty(ii))
0221 pt_in_box = [X1(ii(1)),X2(ii(1)),X3(ii(1))];
0222 else
0223 disp(['did not find pt in']);
0224 stop
0225 end
0226
0227 facets_box = facets_box';
0228
0229
0230
0231 bd_attrib_box = (nboundary+1)*ones(size(facets_box,2),1);
0232 nboundary = nboundary+1;
0233
0234 nodes = [nodes_allcell;nodes_box];
0235 facets = [facets_allcell,facets_box+nnodes_allcell];
0236 bd_attrib = [bd_attrib_allcell; bd_attrib_box];
0237
0238 elseif (include_ECS == 1)
0239 x1min = min(nodes(:,1)); x1max = max(nodes(:,1));
0240 x2min = min(nodes(:,2)); x2max = max(nodes(:,2));
0241 x3min = min(nodes(:,3)); x3max = max(nodes(:,3));
0242
0243 x1len = x1max-x1min;
0244 x2len = x2max-x2min;
0245 x3len = x3max-x3min;
0246
0247 x1gap = x1len*ECS_gap;
0248 x2gap = x2len*ECS_gap;
0249 x3gap = x3len*ECS_gap;
0250
0251 x1min = x1min-x1gap;
0252 x1max = x1max+x1gap;
0253 x2min = x2min-x2gap;
0254 x2max = x2max+x2gap;
0255 x3min = x3min-x3gap;
0256 x3max = x3max+x3gap;
0257
0258 [nodes_box,facets_box,pt_in_box] = create_box_geometry(x1min,x1max,x2min,x2max,x3min,x3max);
0259
0260 pt_in_box = [x1min+x1gap/2,x2min+x2gap/2,x3min+x3gap/2];
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 else
0269 nodes = [nodes_allcell];
0270 facets = [facets_allcell];
0271 bd_attrib = [bd_attrib_allcell];
0272 end
0273
0274 nfacets = size(facets,2);
0275 nnodes = size(nodes,1);
0276
0277 holes = [];
0278 nholes = size(holes,1);
0279
0280 regions = [pt_in_allcell];
0281 if (include_ECS ~= 0)
0282 regions = [regions;pt_in_box];
0283 end
0284
0285 nregions = size(regions,1);
0286
0287 nodesindex = [1:nnodes]';
0288
0289 figure; h=trisurf(facets',nodes(:,1),nodes(:,2),nodes(:,3));
0290 set(h,'facealpha',0.8); view(3); axis equal;
0291 title('Surface triangulation of canonical configuration');