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

Last change on this file since 11027 was 11027, checked in by Eric.Larour, 13 years ago

merged trunk-jpl and trunk for revision 11026

File size: 2.3 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.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers);
28if ~isstruct(md.rifts.riftstruct),
29 error('TriMeshProcessRifts did not find any rift');
30end
31
32%Fill in rest of fields:
33numrifts=length(md.rifts.riftstruct);
34md.mesh.numberofelements=length(md.mesh.elements);
35md.mesh.numberofvertices=length(md.mesh.x);
36md.mesh.z=zeros(md.mesh.numberofvertices,1);
37md.mesh.vertexonboundary=zeros(length(md.mesh.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
38md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
39md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
40md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
41md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
42
43%get coordinates of rift tips
44for i=1:numrifts,
45 md.rifts.riftstruct(i).tip1coordinates=[md.mesh.x(md.rifts.riftstruct(i).tips(1)) md.mesh.y(md.rifts.riftstruct(i).tips(1))];
46 md.rifts.riftstruct(i).tip2coordinates=[md.mesh.x(md.rifts.riftstruct(i).tips(2)) md.mesh.y(md.rifts.riftstruct(i).tips(2))];
47end
48
49%In case we have rifts that open up the domain outline, we need to open them:
50flags=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,domainoutline,'node',0);
51found=0;
52for i=1:numrifts,
53 if flags(md.rifts.riftstruct(i).tips(1))==0,
54 found=1;
55 break;
56 end
57 if flags(md.rifts.riftstruct(i).tips(2))==0,
58 found=1;
59 break;
60 end
61end
62if found,
63 md=meshprocessoutsiderifts(md,domainoutline);
64end
65
66%get elements that are not correctly oriented in the correct direction:
67aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
68pos=find(aires<0);
69md.mesh.elements(pos,:)=[md.mesh.elements(pos,2) md.mesh.elements(pos,1) md.mesh.elements(pos,3)];
Note: See TracBrowser for help on using the repository browser.