Changeset 11004

Timestamp:

12/01/11 14:24:41 (13 years ago)

Author:

seroussi

Message:

extended better setflowequation to macayealstokes and pattynstokes coupling

File:

• issm/trunk-jpl/src/m/model/setflowequation.m (modified) (5 diffs)

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

@@ -67 +67 @@
 end
 
-%1: Hutter elements
+%Initialize node fields
 nodeonhutter=zeros(md.mesh.numberofvertices,1);
 nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
-
-%2: MacAyeal elements
 nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-
-%3: Pattyn elements
 nodeonpattyn=zeros(md.mesh.numberofvertices,1);
 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-
-%4: Stokes elements
+nodeonstokes=zeros(md.mesh.numberofvertices,1);
+nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+noneflag=zeros(md.mesh.numberofelements,1);
+
 %First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
 if any(stokesflag),
-        nodeonstokes=zeros(md.mesh.numberofvertices,1);
-        nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
         fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes));         %find all the nodes on the boundary of the domain without icefront
         fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
         stokesflag(find(fullspcelems))=0;
-end
-nodeonstokes=zeros(md.mesh.numberofvertices,1);
-nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-
-%5: None elements
-noneflag=zeros(md.mesh.numberofelements,1);
+        nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+end
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -139 +131 @@
                 macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
                 macayealpattynflag(find(commonelements))=1;
-                nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
+                nodeonmacayeal(:)=0;
                 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
 
@@ -148 +140 @@
                 list = list - any(sum(nodeonpattyn(md.mesh.elements(pos,:))  ,2),2);
                 pos1=find(list==1);
-                pattynflag(pos(pos1))=0;
                 macayealflag(pos(pos1))=1;
                 macayealpattynflag(pos(pos1))=0;
                 pos2=find(list==-1);
                 pattynflag(pos(pos2))=1;
-                macayealflag(pos(pos2))=0;
                 macayealpattynflag(pos(pos2))=0;
 
+                %Recompute nodes associated to these elements
                 nodeonmacayeal(:)=0;
                 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
@@ -175 +166 @@
                 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
 
-                %Now recreate nodeonpattynstokes
-                nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
+                %rule out elements that don't touch the 2 boundaries
+                pos=find(pattynstokesflag);
+                list=zeros(length(pos),1);
+                list = list + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2);
+                list = list - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2);
+                pos1=find(list==1);
+                stokesflag(pos(pos1))=1;
+                pattynstokesflag(pos(pos1))=0;
+                pos2=find(list==-1);
+                pattynflag(pos(pos2))=1;
+                pattynstokesflag(pos(pos2))=0;
+
+                %Recompute nodes associated to these elements
+                nodeonstokes(:)=0;
+                nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+                nodeonpattyn(:)=0;
+                nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
+                nodeonpattynstokes(:)=0;
                 nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1;
+
         elseif any(stokesflag) & any(macayealflag),
                 %Find node at the border
@@ -190 +198 @@
                 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
 
-                %Now recreate nodeonmacayealstokes
-                nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
+                %rule out elements that don't touch the 2 boundaries
+                pos=find(macayealstokesflag);
+                list=zeros(length(pos),1);
+                list = list + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
+                list = list - any(sum(nodeonstokes(md.mesh.elements(pos,:))  ,2),2);
+                pos1=find(list==1);
+                macayealflag(pos(pos1))=1;
+                macayealstokesflag(pos(pos1))=0;
+                pos2=find(list==-1);
+                stokesflag(pos(pos2))=1;
+                macayealstokesflag(pos(pos2))=0;
+
+                %Recompute nodes associated to these elements
+                nodeonmacayeal(:)=0;
+                nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+                nodeonstokes(:)=0;
+                nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+                nodeonmacayealstokes(:)=0;
                 nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1;
+
         elseif any(stokesflag) & any(hutterflag),
                 error('type of coupling not supported yet');