function [face_nodes,face_elts,crackn]=TreatmentofEnrichment3D(ic,nmod,crack,front,zone) iscale=1; pscale=1; load(fullfile('TMP',sprintf('sample0_%d',0)),'sizeim'); load(fullfile('TMP','params'),'param'); roi=param.roi; check=1; load(fullfile('TMP',sprintf('%d_3d_mesh_%d',nmod,iscale-1)),'xo','yo','zo','Nnodes','Nelems','conn','elt','Smesh'); xo=xo+roi(1)-zone(1); yo=yo+roi(3)-zone(3); zo=zo+roi(5)-zone(5); sizec=size(crack); inbox=~((xo<1)|(xo>sizec(1))|(yo<1)|(yo>sizec(2))|(zo<1)|(zo>sizec(3))); inside_nodes=find(inbox); ind_inside=zeros(Nnodes); ind_inside(inside_nodes)=1:length(inside_nodes); xn=xo(inside_nodes);yn=yo(inside_nodes);zn=zo(inside_nodes); crackn=mexInterpLinear3D(xn,yn,zn,crack); frontn=mexInterpLinear3D(xn,yn,zn,front); hnodes=2*double(crackn>=0)-1; face_elts=[]; selected_nodes=zeros(Nnodes); for ie=1:prod(Nelems) inods=conn(ie,1:elt(ie)); if ~any(~inbox(inods)) xn=xo(inods);yn=yo(inods);zn=zo(inods); iin=ind_inside(inods); above_below=abs(mean(hnodes(iin))); if (above_below<1) infront_behind=frontn(iin)<0; if any(infront_behind) face_elts=[face_elts,ie]; selected_nodes(inods(infront_behind))=1; end end end end face_nodes=find(selected_nodes); xn=xo(face_nodes);yn=yo(face_nodes);zn=zo(face_nodes); crackn=mexInterpLinear3D(xn,yn,zn,crack); frontn=mexInterpLinear3D(xn,yn,zn,front); mesh_file=fullfile('TMP',sprintf('%d_3d_mesh_%d',nmod,iscale-1)); enriched_face_nodes=AddOneNodeLayer(conn,face_nodes); phi=CreateFiniteElementBasis3D(mesh_file,sizeim,pscale,enriched_face_nodes,'sub_cells',false,face_elts); xp=phi*xo(enriched_face_nodes);yp=phi*yo(enriched_face_nodes);zp=phi*zo(enriched_face_nodes); crackp=mexInterpLinear3D(xp,yp,zp,crack); frontp=mexInterpLinear3D(xp,yp,zp,front); if check figure plot3(xo,yo,zo,'k+','MarkerSize',7) hold on % plot3(xo(inside_nodes),yo(inside_nodes),zo(inside_nodes),'rx','MarkerSize',7) plot3(xo(face_nodes),yo(face_nodes),zo(face_nodes),'bo','MarkerSize',7) xlabel('X');ylabel('Y');zlabel('Z'); axis equal print('-djpeg',[param.result_file,'-enrichment.jpg']); print('-depsc2',[param.result_file,'-enrichment.eps']); pause(0.1) end enriched_nodes=face_nodes; save(fullfile('TMP',sprintf('%d_enrichment_%d',nmod,ic)),'enriched_nodes','face_elts','face_nodes','frontn','crackn','zone','crackp','frontp'); end