Changeset 12612

Timestamp:

07/05/12 15:28:36 (13 years ago)

Author:

seroussi

Message:

changes modelextract and contourenvelope for 3d meshes (problem when redoing BC after extraction)

Location:

issm/trunk-jpl/src/m/model

Files:

• contourenvelope.m (modified) (5 diffs)
• modelextract.m (modified) (1 diff)

issm/trunk-jpl/src/m/model/contourenvelope.m

@@ -35 +35 @@
 %Now, build the connectivity tables for this mesh.
 %Computing connectivity
-if size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices,
+if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d),
         md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 end
-if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements,
+if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d),
         md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
 end
@@ -44 +44 @@
 %get nodes inside profile
 mesh.elementconnectivity=md.mesh.elementconnectivity;
+if md.mesh.dimension==2;
+        mesh.elements=md.mesh.elements;
+        mesh.x=md.mesh.x;
+        mesh.y=md.mesh.y;
+        mesh.numberofvertices=md.mesh.numberofvertices;
+        mesh.numberofelements=md.mesh.numberofelements;
+else
+        mesh.elements=md.mesh.elements2d;
+        mesh.x=md.mesh.x2d;
+        mesh.y=md.mesh.y2d;
+        mesh.numberofvertices=md.mesh.numberofvertices2d;
+        mesh.numberofelements=md.mesh.numberofelements2d;
+end
+
 if nargin==2,
+
         if isfile,
                 %get flag list of elements and nodes inside the contour
-                nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1);
-                elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
+                nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1);
+                elemin=(sum(nodein(mesh.elements),2)==size(mesh.elements,2));
                 %modify element connectivity
                 elemout=find(~elemin);
@@ -55 +70 @@
         else
                 %get flag list of elements and nodes inside the contour
-                nodein=zeros(md.mesh.numberofvertices,1); 
-                elemin=zeros(md.mesh.numberofelements,1); 
+                nodein=zeros(mesh.numberofvertices,1); 
+                elemin=zeros(mesh.numberofelements,1); 
                 
                 pos=find(flags); 
                 elemin(pos)=1;
-                nodein(md.mesh.elements(pos,:))=1;
+                nodein(mesh.elements(pos,:))=1;
 
                 %modify element connectivity
@@ -87 +102 @@
         els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
         if length(els2)>1,
-                flag=intersect(md.mesh.elements(els2(1),:),md.mesh.elements(els2(2),:));
-                nods1=md.mesh.elements(el1,:);
+                flag=intersect(mesh.elements(els2(1),:),mesh.elements(els2(2),:));
+                nods1=mesh.elements(el1,:);
                 nods1(find(nods1==flag))=[];
                 segments(count,:)=[nods1 el1];
 
-                ord1=find(nods1(1)==md.mesh.elements(el1,:));
-                ord2=find(nods1(2)==md.mesh.elements(el1,:));
+                ord1=find(nods1(1)==mesh.elements(el1,:));
+                ord2=find(nods1(2)==mesh.elements(el1,:));
 
                 %swap segment nodes if necessary
@@ -104 +119 @@
                 count=count+1;
         else
-                nods1=md.mesh.elements(el1,:);
-                flag=setdiff(nods1,md.mesh.elements(els2,:));
+                nods1=mesh.elements(el1,:);
+                flag=setdiff(nods1,mesh.elements(els2,:));
                 for j=1:3,
                         nods=nods1; nods(j)=[];
                         if any(ismember(flag,nods)),
                                 segments(count,:)=[nods el1];
-                                ord1=find(nods(1)==md.mesh.elements(el1,:));
-                                ord2=find(nods(2)==md.mesh.elements(el1,:));
+                                ord1=find(nods(1)==mesh.elements(el1,:));
+                                ord2=find(nods(2)==mesh.elements(el1,:));
                                 if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
                                         temp=segments(count,1);

issm/trunk-jpl/src/m/model/modelextract.m

@@ -210 +210 @@
                 md2.mesh.segments=contourenvelope(md2);
                 md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+        else
+                %First do the connectivity for the contourenvelope in 2d
+                md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
+                md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
+                md2.mesh.segments=contourenvelope(md2);
+                md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+                md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
+                %Then do it for 3d as usual
+                md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
+                md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
         end