1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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