source: issm/trunk/src/m/model/mesh/rifts/meshprocessrifts.m@ 9616

Last change on this file since 9616 was 9616, checked in by Mathieu Morlighem, 14 years ago

Added check

File size: 2.1 KB
Line 
1function md=meshprocessrifts(md,domainoutline)
2%MESHPROCESSRIFTS - process mesh when rifts are present
3%
4% split rifts inside mesh (rifts are defined by presence of
5% segments inside the domain outline)
6% if domain outline is provided, check for rifts that could touch it, and open them up.
7%
8% Usage:
9% md=meshprocessrifts(md,domainoutline)
10%
11% Ex:
12% md=meshprocessrifts(md,'DomainOutline.exp');
13%
14
15%some checks on arguments:
16if nargout~=1,
17 help meshprocessrifts
18 error('meshprocessrifts usage error:');
19end
20
21if nargin~=2,
22 help meshprocessrifts
23 error('meshprocessrifts usage error:');
24end
25
26%Call MEX file
27[md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);
28if ~isstruct(md.rifts),
29 error('TriMeshProcessRifts did not find any rift');
30end
31
32%Fill in rest of fields:
33md.numberofelements=length(md.elements);
34md.numberofnodes=length(md.x);
35md.z=zeros(md.numberofnodes,1);
36md.nodeonboundary=zeros(length(md.x),1); md.nodeonboundary(md.segments(:,1:2))=1;
37md.numrifts=length(md.rifts);
38md.elements_type=3*ones(md.numberofelements,1);
39md.nodeonbed=ones(md.numberofnodes,1);
40md.nodeonsurface=ones(md.numberofnodes,1);
41md.elementonbed=ones(md.numberofelements,1);
42md.elementonsurface=ones(md.numberofelements,1);
43
44%get coordinates of rift tips
45for i=1:md.numrifts,
46 md.rifts(i).tip1coordinates=[md.x(md.rifts(i).tips(1)) md.y(md.rifts(i).tips(1))];
47 md.rifts(i).tip2coordinates=[md.x(md.rifts(i).tips(2)) md.y(md.rifts(i).tips(2))];
48end
49
50%In case we have rifts that open up the domain outline, we need to open them:
51flags=ContourToMesh(md.elements,md.x,md.y,domainoutline,'node',0);
52found=0;
53for i=1:md.numrifts,
54 if flags(md.rifts(i).tips(1))==0,
55 found=1;
56 break;
57 end
58 if flags(md.rifts(i).tips(2))==0,
59 found=1;
60 break;
61 end
62end
63if found,
64 md=meshprocessoutsiderifts(md,domainoutline);
65end
66
67%get elements that are not correctly oriented in the correct direction:
68aires=GetAreas(md.elements,md.x,md.y);
69pos=find(aires<0);
70md.elements(pos,:)=[md.elements(pos,2) md.elements(pos,1) md.elements(pos,3)];
Note: See TracBrowser for help on using the repository browser.