Changeset 5613


Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r5575 r5613  
    5858checksize(md,fields,[md.numberofelements 1]);
    5959%Check the values of elements_type
    60 checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
     60checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() PattynStokesApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
    6161if (md.dim==2),
    6262        checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]);
  • issm/trunk/src/m/classes/public/plot/plot_elementstype.m

    r5433 r5613  
    6969        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
    7070        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     71        %Stokes elements
     72        posS=find(data==StokesApproximationEnum);
     73        A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
     74        p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     75        patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     76        patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     77        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     78        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
    7179        %MacAyealPattyn elements
    7280        posP=find(data==MacAyealPattynApproximationEnum);
     
    7785        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
    7886        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
    79         %Stokes elements
    80         alpha=0.35;
    81         posS=find(data==StokesApproximationEnum);
    82         if ~isempty(posS)
    83                 A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
    84                 p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
    85                 patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
    86                 patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
    87                 patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
    88                 patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
    89                 legend([p1 p2 p3 p5 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements','Stokes''s elements');
    90         else
    91                 legend([p1 p2 p3 p5],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements');
    92         end
     87        %PattynStokes elements
     88        PosPS=find(data==PattynStokesApproximationEnum);
     89        A=elements(PosPS,1); B=elements(PosPS,2); C=elements(PosPS,3); D=elements(PosPS,4); E=elements(PosPS,5); F=elements(PosPS,6);
     90        p6=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     91        patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     92        patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     93        patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     94        patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
     95        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');
    9396end
    9497
  • issm/trunk/src/m/classes/public/setelementstype.m

    r5534 r5613  
    4343        macayealflag(find(~hutterflag & ~pattynflag))=1;
    4444elseif strcmpi(filltype,'pattyn'),
    45         pattynflag(find(~hutterflag & ~macayealflag))=1;
     45        pattynflag(find(~hutterflag & ~macayealflag & ~stokesflag))=1;
    4646end
    4747
     
    5252
    5353%check that each element has only one flag
    54 if any(hutterflag+ macayealflag+pattynflag>1),
     54if any(hutterflag+ macayealflag+pattynflag+stokesflag>1),
    5555        disp('setelementstype warning message: some elements have several types, higher order type is used for them')
    5656        hutterflag(find(hutterflag & macayealflag))=0;
     57        hutterflag(find(hutterflag & pattynflag))=0;
    5758        macayealflag(find(macayealflag & pattynflag))=0;
    58         pattynflag(find(pattynflag & macayealflag))=0;
    5959end
    6060
     
    6767
    6868%Stokes can only be used alone for now:
    69 if any(stokesflag) &any(hutterflag+pattynflag+macayealflag),
     69if any(stokesflag) &any(hutterflag+macayealflag),
    7070        error('setelementstype error message: stokes cannot be used with any other model for now, put stokes everywhere')
    7171end
    7272
    73 if any(stokesflag),
    74         %Modify stokesflag to get rid of elements contrained everywhere
    75         fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3);         %find all the grids on the boundary of the domain without icefront
    76         fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
    77         stokesflag(find(fullspcelems))=0;
    78 end
    79        
    8073%add in model who is who
    8174md.elements_type=zeros(md.numberofelements,1);
     
    10093
    10194%4: Stokes elements
     95%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
     96if any(stokesflag),
     97        gridonstokes=zeros(md.numberofgrids,1);
     98        gridonstokes(md.elements(find(stokesflag),:))=1;
     99        fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3 | (gridonpattyn & gridonstokes));         %find all the grids on the boundary of the domain without icefront
     100        fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
     101        stokesflag(find(fullspcelems))=0;
     102end
    102103gridonstokes=zeros(md.numberofgrids,1);
    103104gridonstokes(md.elements(find(stokesflag),:))=1;
     
    105106md.elements_type(find(stokesflag))=StokesApproximationEnum();
    106107
    107 %Complete with NoneApproximation if there is no stokes
     108%Then complete with NoneApproximation or the other model used if there is no stokes
    108109if any(stokesflag),
     110        if any(pattynflag), %fill with pattyn
     111                pattynflag(~stokesflag)=1;
     112                gridonpattyn(md.elements(find(pattynflag),:))=1;
     113                md.gridonpattyn=gridonpattyn;
     114                md.elements_type(find(~stokesflag))=PattynApproximationEnum();
     115        else %fill with none
    109116        %5: None elements (non Stokes)
    110         md.elements_type(find(~stokesflag))=NoneApproximationEnum();
     117                md.elements_type(find(~stokesflag))=NoneApproximationEnum();
     118        end
    111119end
    112120
    113121%Now take care of the coupling between MacAyeal and Pattyn
    114122md.penalties=[];
     123gridonmacayealpattyn=zeros(md.numberofgrids,1);
     124gridonpattynstokes=zeros(md.numberofgrids,1);
    115125if strcmpi(coupling_method,'penalties'),
    116         gridonmacayealpattyn=zeros(md.numberofgrids,1);
    117126        %Create the border grids between Pattyn and MacAyeal and extrude them
    118127        numgrids2d=md.numberofgrids2d;
     
    129138        end
    130139elseif strcmpi(coupling_method,'tiling'),
    131         %Find grid at the border
    132         gridonmacayealpattyn=zeros(md.numberofgrids,1);
    133         gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
    134         %Macayeal elements in contact with this layer become MacAyealPattyn elements
    135         matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
    136         commonelements=sum(matrixelements,2)~=0;
    137         commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
    138         macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
    139         macayealpattynflag=zeros(md.numberofelements,1);
    140         macayealpattynflag(find(commonelements))=1;
    141         gridonmacayeal=zeros(md.numberofgrids,1);
    142         gridonmacayeal(md.elements(find(macayealflag),:))=1;
    143         md.gridonmacayeal=gridonmacayeal;
    144 
    145         %Create MacaAyealPattynApproximation where needed
    146         md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
    147 
    148         %Now recreate gridonmacayeal and gridonmacayealpattyn
    149         gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
     140        if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn
     141                %Find grid at the border
     142                gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
     143                %Macayeal elements in contact with this layer become MacAyealPattyn elements
     144                matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
     145                commonelements=sum(matrixelements,2)~=0;
     146                commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
     147                macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
     148                macayealpattynflag=zeros(md.numberofelements,1);
     149                macayealpattynflag(find(commonelements))=1;
     150                gridonmacayeal=zeros(md.numberofgrids,1);
     151                gridonmacayeal(md.elements(find(macayealflag),:))=1;
     152                md.gridonmacayeal=gridonmacayeal;
     153
     154                %Create MacaAyealPattynApproximation where needed
     155                md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
     156
     157                %Now recreate gridonmacayeal and gridonmacayealpattyn
     158                gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
     159        elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
     160                %Find grid at the border
     161                gridonpattynstokes(find(gridonpattyn & gridonstokes))=1;
     162                %Stokes elements in contact with this layer become PattynStokes elements
     163                matrixelements=ismember(md.elements,find(gridonpattynstokes));
     164                commonelements=sum(matrixelements,2)~=0;
     165                commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
     166                stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
     167                pattynstokesflag=zeros(md.numberofelements,1);
     168                pattynstokesflag(find(commonelements))=1;
     169                gridonstokes=zeros(md.numberofgrids,1);
     170                gridonstokes(md.elements(find(stokesflag),:))=1;
     171                md.gridonstokes=gridonstokes;
     172
     173                %Create MacaAyealPattynApproximation where needed
     174                md.elements_type(find(pattynstokesflag))=PattynStokesApproximationEnum();
     175
     176                %Now recreate gridonmacayeal and gridonmacayealpattyn
     177                gridonpattynstokes(md.elements(find(pattynstokesflag),:))=1;
     178        elseif any(stokesflag) & any(macayealflag+hutterflag),
     179                error('type of coupling not supported yet');
     180        end
    150181end
    151182
     
    164195pos=find(gridonmacayealpattyn);
    165196md.vertices_type(pos)=MacAyealPattynApproximationEnum();
     197pos=find(gridonpattynstokes);
     198md.vertices_type(pos)=PattynStokesApproximationEnum();
    166199pos=find(gridonstokes);
    167200md.vertices_type(pos)=StokesApproximationEnum();
Note: See TracChangeset for help on using the changeset viewer.