Changeset 5613

Timestamp:

08/27/10 11:56:03 (15 years ago)

Author:

seroussi

Message:

added pattyn/stokes coupling in setelementstype

Location:

issm/trunk/src/m/classes/public

Files:

• ismodelselfconsistent.m (modified) (1 diff)
• plot/plot_elementstype.m (modified) (2 diffs)
• setelementstype.m (modified) (7 diffs)

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -58 +58 @@
 checksize(md,fields,[md.numberofelements 1]);
 %Check the values of elements_type
-checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
+checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() PattynStokesApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
 if (md.dim==2),
         checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]);

issm/trunk/src/m/classes/public/plot/plot_elementstype.m

@@ -69 +69 @@
         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        %Stokes elements
+        posS=find(data==StokesApproximationEnum);
+        A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
+        p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
         %MacAyealPattyn elements
         posP=find(data==MacAyealPattynApproximationEnum);
@@ -77 +85 @@
         patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
         patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
-        %Stokes elements
-        alpha=0.35;
-        posS=find(data==StokesApproximationEnum);
-        if ~isempty(posS)
-                A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
-                p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-                patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-                patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-                patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-                patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-                legend([p1 p2 p3 p5 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements','Stokes''s elements');
-        else
-                legend([p1 p2 p3 p5],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements');
-        end
+        %PattynStokes elements
+        PosPS=find(data==PattynStokesApproximationEnum);
+        A=elements(PosPS,1); B=elements(PosPS,2); C=elements(PosPS,3); D=elements(PosPS,4); E=elements(PosPS,5); F=elements(PosPS,6);
+        p6=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+        legend([p1 p2 p3 p4 p5 p6],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements','MacAyealPattyn''s elements','PattynStokes''s elements');
 end
 

issm/trunk/src/m/classes/public/setelementstype.m

@@ -43 +43 @@
         macayealflag(find(~hutterflag & ~pattynflag))=1;
 elseif strcmpi(filltype,'pattyn'),
-        pattynflag(find(~hutterflag & ~macayealflag))=1;
+        pattynflag(find(~hutterflag & ~macayealflag & ~stokesflag))=1;
 end
 
@@ -52 +52 @@
 
 %check that each element has only one flag
-if any(hutterflag+ macayealflag+pattynflag>1),
+if any(hutterflag+ macayealflag+pattynflag+stokesflag>1),
         disp('setelementstype warning message: some elements have several types, higher order type is used for them')
         hutterflag(find(hutterflag & macayealflag))=0;
+        hutterflag(find(hutterflag & pattynflag))=0;
         macayealflag(find(macayealflag & pattynflag))=0;
-        pattynflag(find(pattynflag & macayealflag))=0;
 end
 
@@ -67 +67 @@
 
 %Stokes can only be used alone for now:
-if any(stokesflag) &any(hutterflag+pattynflag+macayealflag),
+if any(stokesflag) &any(hutterflag+macayealflag),
         error('setelementstype error message: stokes cannot be used with any other model for now, put stokes everywhere')
 end
 
-if any(stokesflag),
-        %Modify stokesflag to get rid of elements contrained everywhere
-        fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3);         %find all the grids on the boundary of the domain without icefront
-        fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
-        stokesflag(find(fullspcelems))=0;
-end
-        
 %add in model who is who
 md.elements_type=zeros(md.numberofelements,1);
@@ -100 +93 @@
 
 %4: Stokes elements
+%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
+if any(stokesflag),
+        gridonstokes=zeros(md.numberofgrids,1);
+        gridonstokes(md.elements(find(stokesflag),:))=1;
+        fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3 | (gridonpattyn & gridonstokes));         %find all the grids on the boundary of the domain without icefront
+        fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
+        stokesflag(find(fullspcelems))=0;
+end
 gridonstokes=zeros(md.numberofgrids,1);
 gridonstokes(md.elements(find(stokesflag),:))=1;
@@ -105 +106 @@
 md.elements_type(find(stokesflag))=StokesApproximationEnum();
 
-%Complete with NoneApproximation if there is no stokes
+%Then complete with NoneApproximation or the other model used if there is no stokes
 if any(stokesflag), 
+        if any(pattynflag), %fill with pattyn
+                pattynflag(~stokesflag)=1;
+                gridonpattyn(md.elements(find(pattynflag),:))=1;
+                md.gridonpattyn=gridonpattyn;
+                md.elements_type(find(~stokesflag))=PattynApproximationEnum();
+        else %fill with none 
         %5: None elements (non Stokes)
-        md.elements_type(find(~stokesflag))=NoneApproximationEnum();
+                md.elements_type(find(~stokesflag))=NoneApproximationEnum();
+        end
 end
 
 %Now take care of the coupling between MacAyeal and Pattyn
 md.penalties=[];
+gridonmacayealpattyn=zeros(md.numberofgrids,1);
+gridonpattynstokes=zeros(md.numberofgrids,1);
 if strcmpi(coupling_method,'penalties'),
-        gridonmacayealpattyn=zeros(md.numberofgrids,1);
         %Create the border grids between Pattyn and MacAyeal and extrude them
         numgrids2d=md.numberofgrids2d;
@@ -129 +138 @@
         end
 elseif strcmpi(coupling_method,'tiling'),
-        %Find grid at the border
-        gridonmacayealpattyn=zeros(md.numberofgrids,1);
-        gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
-        %Macayeal elements in contact with this layer become MacAyealPattyn elements
-        matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
-        commonelements=sum(matrixelements,2)~=0;
-        commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
-        macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
-        macayealpattynflag=zeros(md.numberofelements,1);
-        macayealpattynflag(find(commonelements))=1;
-        gridonmacayeal=zeros(md.numberofgrids,1);
-        gridonmacayeal(md.elements(find(macayealflag),:))=1;
-        md.gridonmacayeal=gridonmacayeal;
-
-        %Create MacaAyealPattynApproximation where needed
-        md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
-
-        %Now recreate gridonmacayeal and gridonmacayealpattyn
-        gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
+        if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn
+                %Find grid at the border
+                gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
+                %Macayeal elements in contact with this layer become MacAyealPattyn elements
+                matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
+                commonelements=sum(matrixelements,2)~=0;
+                commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
+                macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
+                macayealpattynflag=zeros(md.numberofelements,1);
+                macayealpattynflag(find(commonelements))=1;
+                gridonmacayeal=zeros(md.numberofgrids,1);
+                gridonmacayeal(md.elements(find(macayealflag),:))=1;
+                md.gridonmacayeal=gridonmacayeal;
+
+                %Create MacaAyealPattynApproximation where needed
+                md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
+
+                %Now recreate gridonmacayeal and gridonmacayealpattyn
+                gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
+        elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
+                %Find grid at the border
+                gridonpattynstokes(find(gridonpattyn & gridonstokes))=1;
+                %Stokes elements in contact with this layer become PattynStokes elements
+                matrixelements=ismember(md.elements,find(gridonpattynstokes));
+                commonelements=sum(matrixelements,2)~=0;
+                commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
+                stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
+                pattynstokesflag=zeros(md.numberofelements,1);
+                pattynstokesflag(find(commonelements))=1;
+                gridonstokes=zeros(md.numberofgrids,1);
+                gridonstokes(md.elements(find(stokesflag),:))=1;
+                md.gridonstokes=gridonstokes;
+
+                %Create MacaAyealPattynApproximation where needed
+                md.elements_type(find(pattynstokesflag))=PattynStokesApproximationEnum();
+
+                %Now recreate gridonmacayeal and gridonmacayealpattyn
+                gridonpattynstokes(md.elements(find(pattynstokesflag),:))=1;
+        elseif any(stokesflag) & any(macayealflag+hutterflag),
+                error('type of coupling not supported yet');
+        end
 end
 
@@ -164 +195 @@
 pos=find(gridonmacayealpattyn);
 md.vertices_type(pos)=MacAyealPattynApproximationEnum();
+pos=find(gridonpattynstokes);
+md.vertices_type(pos)=PattynStokesApproximationEnum();
 pos=find(gridonstokes);
 md.vertices_type(pos)=StokesApproximationEnum();