Changeset 11027 for issm/trunk/src/m/model/setflowequation.m

Timestamp:

12/09/11 20:36:51 (13 years ago)

Author:

Eric.Larour

Message:

merged trunk-jpl and trunk for revision 11026

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src/m/model/setflowequation.m (modified) (8 diffs)

issm/trunk/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
+%Initialize node fields
 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
+nodeonstokes=zeros(md.mesh.numberofvertices,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;
-md.flowequation.borderstokes=nodeonstokes;
-md.flowequation.element_equation(find(stokesflag))=4;
+        nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+end
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -106 +90 @@
                 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 +103 @@
 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 +129 @@
                 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
+                nodeonmacayeal(:)=0;
+                nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+
+                %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);
+                macayealflag(pos(pos1))=1;
+                macayealpattynflag(pos(pos1))=0;
+                pos2=find(list==-1);
+                pattynflag(pos(pos2))=1;
+                macayealpattynflag(pos(pos2))=0;
+
+                %Recompute nodes associated to these elements
+                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 +161 @@
                 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
-                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
@@ -186 +193 @@
                 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
-                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');
         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 +248 @@
 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 +268 @@
 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