Changeset 11003

Timestamp:

12/01/11 14:03:21 (13 years ago)

Author:

Mathieu Morlighem

Message:

cleaning up setflowequations

File:

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

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

@@ -67 +67 @@
 end
 
-%add in model who is who
-md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
-
 %1: Hutter elements
 nodeonhutter=zeros(md.mesh.numberofvertices,1);
 nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
-md.flowequation.element_equation(find(hutterflag))=1;
 
 %2: MacAyeal elements
 nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-md.flowequation.bordermacayeal=nodeonmacayeal;
-md.flowequation.element_equation(find(macayealflag))=2;
 
 %3: Pattyn elements
 nodeonpattyn=zeros(md.mesh.numberofvertices,1);
 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-md.flowequation.borderpattyn=nodeonpattyn;
-md.flowequation.element_equation(find(pattynflag))=3;
 
 %4: Stokes elements
@@ -98 +90 @@
 nodeonstokes=zeros(md.mesh.numberofvertices,1);
 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-md.flowequation.borderstokes=nodeonstokes;
-md.flowequation.element_equation(find(stokesflag))=4;
+
+%5: None elements
+noneflag=zeros(md.mesh.numberofelements,1);
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -106 +99 @@
                 pattynflag(~stokesflag)=1;
                 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-                md.flowequation.borderpattyn=nodeonpattyn;
-                md.flowequation.element_equation(find(~stokesflag))=3;
         elseif any(macayealflag), %fill with macayeal
                 macayealflag(~stokesflag)=1;
                 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-                md.flowequation.bordermacayeal=nodeonmacayeal;
-                md.flowequation.element_equation(find(~stokesflag))=2;
         else %fill with none 
-        %5: None elements (non Stokes)
-                md.flowequation.element_equation(find(~stokesflag))=0;
+                noneflag(find(~stokesflag))=1;
         end
 end
@@ -124 +112 @@
 nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
 nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
+macayealpattynflag=zeros(md.mesh.numberofelements,1);
+macayealstokesflag=zeros(md.mesh.numberofelements,1);
+pattynstokesflag=zeros(md.mesh.numberofelements,1);
 if strcmpi(coupling_method,'penalties'),
         %Create the border nodes between Pattyn and MacAyeal and extrude them
@@ -147 +138 @@
                 commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
                 macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
-                macayealpattynflag=zeros(md.mesh.numberofelements,1);
                 macayealpattynflag(find(commonelements))=1;
                 nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
                 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-                md.flowequation.bordermacayeal=nodeonmacayeal;
-
-                %Create MacaAyealPattynApproximation where needed
-                md.flowequation.element_equation(find(macayealpattynflag))=5;
-
-                %Now recreate nodeonmacayeal and nodeonmacayealpattyn
+
+                %rule out elements that don't touch the 2 boundaries
+                pos=find(macayealpattynflag);
+                list=zeros(length(pos),1);
+                list = list + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
+                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;
+
+                nodeonmacayeal(:)=0;
+                nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+                nodeonpattyn(:)=0;
+                nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
+                nodeonmacayealpattyn(:)=0;
                 nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1;
+
         elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
                 %Find node at the border
@@ -166 +171 @@
                 commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
                 stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
-                pattynstokesflag=zeros(md.mesh.numberofelements,1);
                 pattynstokesflag(find(commonelements))=1;
                 nodeonstokes=zeros(md.mesh.numberofvertices,1);
                 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-                md.flowequation.borderstokes=nodeonstokes;
-
-                %Create MacaAyealPattynApproximation where needed
-                md.flowequation.element_equation(find(pattynstokesflag))=7;
 
                 %Now recreate nodeonpattynstokes
@@ -186 +186 @@
                 commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal
                 stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements
-                macayealstokesflag=zeros(md.mesh.numberofelements,1);
                 macayealstokesflag(find(commonelements))=1;
                 nodeonstokes=zeros(md.mesh.numberofvertices,1);
                 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-                md.flowequation.borderstokes=nodeonstokes;
-
-                %Create MacaAyeal Approximation where needed
-                md.flowequation.element_equation(find(macayealstokesflag))=6;
 
                 %Now recreate nodeonmacayealstokes
@@ -202 +197 @@
         end
 end
+
+%Create MacaAyealPattynApproximation where needed
+md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
+md.flowequation.element_equation(find(noneflag))=0;
+md.flowequation.element_equation(find(hutterflag))=1;
+md.flowequation.element_equation(find(macayealflag))=2;
+md.flowequation.element_equation(find(pattynflag))=3;
+md.flowequation.element_equation(find(stokesflag))=4;
+md.flowequation.element_equation(find(macayealpattynflag))=5;
+md.flowequation.element_equation(find(macayealstokesflag))=6;
+md.flowequation.element_equation(find(pattynstokesflag))=7;
+
+%border
+md.flowequation.borderpattyn=nodeonpattyn;
+md.flowequation.bordermacayeal=nodeonmacayeal;
+md.flowequation.borderstokes=nodeonstokes;
 
 %Create vertices_type
@@ -213 +224 @@
 pos=find(nodeonhutter);
 md.flowequation.vertex_equation(pos)=1;
-pos=find(nodeonpattyn & nodeonmacayeal);
-md.flowequation.vertex_equation(pos)=3;
 pos=find(nodeonmacayealpattyn);
 md.flowequation.vertex_equation(pos)=5;
@@ -235 +244 @@
 md.flowequation.isstokes=double(any(md.flowequation.element_equation==4));
 
-end
+return
+
+%Check that tiling can work:
+if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1),
+        error('error coupling domain too irregular');
+end
+if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1),
+        error('error coupling domain too irregular');
+end
+if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1),
+        error('error coupling domain too irregular');
+end