Changeset 2737

Timestamp:

01/04/10 14:26:01 (15 years ago)

Author:

Eric.Larour

Message:

meshnodensity: used to include rifts inside existing mesh: creates mesh without density specification. Densities of the PSLG domain
outlines are used.
meshplug: plug mesh into another one.
meshyamsrecreateriftsegments: plugin for yams meshing routine to refine around rifts.
rifttipsrefine,meshprocessrifts moved from upper directory mesh/
meshaddrifts: add rifts (specified in a rift outline) to an existing mesh! Kicks some serious ass!
updateriftindexing: used in meshaddrifts. reindexes rifts according to new mesh.

Location:

issm/trunk/src/m/classes/public/mesh

Files:

• meshnodensity.m (added)
• meshprocessrifts.m (deleted)
• meshyams.m (modified) (3 diffs)
• rifts/meshaddrifts.m (added)
• rifts/meshplug.m (added)
• rifts/meshprocessrifts.m (added)
• rifts/meshyamsrecreateriftsegments.m (added)
• rifts/rifttipsrefine.m (added)
• rifts/updateriftindexing.m (added)
• rifttipsrefine.m (deleted)

issm/trunk/src/m/classes/public/mesh/meshyams.m

@@ -54 +54 @@
 else
         md=mesh(md,domainoutline,riftoutline,resolution);
+        md=meshprocessrifts(md);
 end
 disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
@@ -90 +91 @@
         md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
 
+        %if we have rifts, we just messed them up, we need to recreate the segments that constitute those 
+        %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
+        if md.numrifts, 
+                md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+                md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+                md.segments=findsegments(md);
+                md=meshyamsrecreateriftsegments(md);
+        end
+
 end
         
@@ -123 +133 @@
 end
 md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
+
+%deal with rifts 
+if md.numrifts,
+        %first, recreate rift segments
+        md=meshyamsrecreateriftsegments(md);
+
+        %using the segments, recreate the penaltypairs
+        for j=1:md.numrifts,
+                rift=md.rifts(j);
+
+                %build normals and lengths of segments:
+                lengths=sqrt((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 );
+                normalsx=cos(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
+                normalsy=sin(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
+
+                %ok, build penaltypairs: 
+                numpenaltypairs=length(rift.segments)/2-1;
+                rift.penaltypairs=zeros(numpenaltypairs,7);
+
+                for i=1:numpenaltypairs,
+                        rift.penaltypairs(i,1)=rift.segments(i,2);
+                        rift.penaltypairs(i,2)=rift.segments(end-i,2);
+                        rift.penaltypairs(i,3)=rift.segments(i,3);
+                        rift.penaltypairs(i,4)=rift.segments(end-i,3);
+                        rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
+                        rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
+                        rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
+                end
+                %renormalize norms: 
+                norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
+                rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
+                rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
+
+                md.rifts(j)=rift;
+        end
+
+end
+
+
+
+
+