source: issm/trunk/src/m/model/modelextract.m@ 9451

Last change on this file since 9451 was 9451, checked in by Mathieu Morlighem, 14 years ago

removed all other model fields to be deleted

File size: 7.7 KB
RevLine 
[1495]1function md2=modelextract(md1,area)
2%modelextract - extract a model according to an Argus contour or flag list
[1348]3%
4% This routine extracts a submodel from a bigger model with respect to a given contour
5% md must be followed by the corresponding exp file or flags list
6% It can either be a domain file (argus type, .exp extension), or an array of element flags.
7% If user wants every element outside the domain to be
8% modelextract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp');
9% an empty string '' will be considered as an empty domain
10% a string 'all' will be considered as the entire domain
11% add an argument 0 if you do not want the elements to be checked (faster)
12%
13% Usage:
[1495]14% md2=modelextract(md1,area);
[1348]15%
16% Examples:
[1495]17% md2=modelextract(md,'Domain.exp');
18% md2=modelextract(md,md.elementoniceshelf);
[1348]19%
20% See also: EXTRUDE, COLLAPSE
21
22%some checks
[1368]23if ((nargin~=2) | (nargout~=1)),
[3589]24 help modelextract
25 error('modelextract error message: bad usage');
[1348]26end
27
28%get check option
29if (nargin==3 & varargin{1}==0),
30 checkoutline=0;
31else
32 checkoutline=1;
33end
34
35%first
[2988]36flag_elem=FlagElements(md1,area);
[1348]37
[1361]38%kick out all elements with 3 dirichlets
39spc_elem=find(~flag_elem);
[8298]40spc_node=sort(unique(md1.elements(spc_elem,:)));
41flag=ones(md1.numberofnodes,1);
42flag(spc_node)=0;
[1361]43pos=find(sum(flag(md1.elements),2)==0);
44flag_elem(pos)=0;
45
46%extracted elements and nodes lists
[1348]47pos_elem=find(flag_elem);
[8298]48pos_node=sort(unique(md1.elements(pos_elem,:)));
[1348]49
50%keep track of some fields
[8298]51numberofnodes1=md1.numberofnodes;
[1348]52numberofelements1=md1.numberofelements;
[8298]53numberofnodes2=length(pos_node);
[1348]54numberofelements2=length(pos_elem);
[8298]55flag_node=zeros(numberofnodes1,1);
56flag_node(pos_node)=1;
[1348]57
[8298]58%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
[1348]59Pelem=zeros(numberofelements1,1);
60Pelem(pos_elem)=[1:numberofelements2]';
[8298]61Pnode=zeros(numberofnodes1,1);
62Pnode(pos_node)=[1:numberofnodes2]';
[1348]63
[8298]64%renumber the elements (some node won't exist anymore)
[1348]65elements_1=md1.elements;
66elements_2=elements_1(pos_elem,:);
[8298]67elements_2(:,1)=Pnode(elements_2(:,1));
68elements_2(:,2)=Pnode(elements_2(:,2));
69elements_2(:,3)=Pnode(elements_2(:,3));
[4348]70if md1.dim==3,
[8298]71 elements_2(:,4)=Pnode(elements_2(:,4));
72 elements_2(:,5)=Pnode(elements_2(:,5));
73 elements_2(:,6)=Pnode(elements_2(:,6));
[1348]74end
75
76%OK, now create the new model !
77
78 %take every fields from model
79 md2=md1;
80
[1361]81 %automatically modify fields
82
83 %loop over model fields
84 model_fields=fields(md1);
85 for i=1:length(model_fields),
86
87 %get field
[9436]88 field=md1.(model_fields{i});
[1361]89 fieldsize=size(field);
90
[8298]91 %size = number of nodes * n
92 if fieldsize(1)==numberofnodes1
[9436]93 md2.(model_fields{i})=field(pos_node,:);
[1361]94
95 %size = number of elements * n
96 elseif fieldsize(1)==numberofelements1
[9436]97 md2.(model_fields{i})=field(pos_elem,:);
[1361]98 end
99 end
100
[1348]101 %modify some specific fields
102
103 %Mesh
104 md2.numberofelements=numberofelements2;
[8298]105 md2.numberofnodes=numberofnodes2;
[1348]106 md2.elements=elements_2;
107
[1368]108 %uppernodes lowernodes
[4348]109 if md1.dim==3
[8298]110 md2.uppernodes=md1.uppernodes(pos_node);
111 pos=find(~isnan(md2.uppernodes));
112 md2.uppernodes(pos)=Pnode(md2.uppernodes(pos));
[1368]113
[8298]114 md2.lowernodes=md1.lowernodes(pos_node);
115 pos=find(~isnan(md2.lowernodes));
116 md2.lowernodes(pos)=Pnode(md2.lowernodes(pos));
[4690]117
118 md2.upperelements=md1.upperelements(pos_elem);
119 pos=find(~isnan(md2.upperelements));
120 md2.upperelements(pos)=Pelem(md2.upperelements(pos));
121
122 md2.lowerelements=md1.lowerelements(pos_elem);
123 pos=find(~isnan(md2.lowerelements));
124 md2.lowerelements(pos)=Pelem(md2.lowerelements(pos));
[1368]125 end
126
[1348]127 %Initial 2d mesh
[4348]128 if md1.dim==3
[1348]129 flag_elem_2d=flag_elem(1:md1.numberofelements2d);
130 pos_elem_2d=find(flag_elem_2d);
[8298]131 flag_node_2d=flag_node(1:md1.numberofnodes2d);
132 pos_node_2d=find(flag_node_2d);
[1348]133
134 md2.numberofelements2d=length(pos_elem_2d);
[8298]135 md2.numberofnodes2d=length(pos_node_2d);
[1348]136 md2.elements2d=md1.elements2d(pos_elem_2d,:);
[8298]137 md2.elements2d(:,1)=Pnode(md2.elements2d(:,1));
138 md2.elements2d(:,2)=Pnode(md2.elements2d(:,2));
139 md2.elements2d(:,3)=Pnode(md2.elements2d(:,3));
[1348]140
[5433]141 if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d); end;
[8298]142 md2.x2d=md1.x(pos_node_2d);
143 md2.y2d=md1.y(pos_node_2d);
[1348]144 end
145
[3583]146 %Edges
147 if size(md2.edges,2)>1, %do not use ~isnan because there are some NaNs...
148 %renumber first two columns
[5442]149 pos=find(~isnan(md2.edges(:,4)));
[8298]150 md2.edges(: ,1)=Pnode(md2.edges(:,1));
151 md2.edges(: ,2)=Pnode(md2.edges(:,2));
[5442]152 md2.edges(: ,3)=Pelem(md2.edges(:,3));
153 md2.edges(pos,4)=Pelem(md2.edges(pos,4));
154 %remove edges when the 2 vertices are not in the domain.
155 md2.edges=md2.edges(find(md2.edges(:,1) & md2.edges(:,2)),:);
156 %Replace all zeros by NaN in the last two columns;
157 pos=find(md2.edges(:,3)==0);
158 md2.edges(pos,3)=NaN;
159 pos=find(md2.edges(:,4)==0);
160 md2.edges(pos,4)=NaN;
[5451]161 %Invert NaN of the third column with last column (Also invert first two columns!!)
[5442]162 pos=find(isnan(md2.edges(:,3)));
163 md2.edges(pos,3)=md2.edges(pos,4);
164 md2.edges(pos,4)=NaN;
[5451]165 values=md2.edges(pos,2);
166 md2.edges(pos,2)=md2.edges(pos,1);
167 md2.edges(pos,1)=values;
[5442]168 %Finally remove edges that do not belong to any element
169 pos=find(isnan(md2.edges(:,3)) & isnan(md2.edges(:,4)));
170 md2.edges(pos,:)=[];
[3583]171 end
172
[1348]173 %Penalties
174 if ~isnan(md2.penalties),
175 for i=1:size(md1.penalties,1);
[8298]176 md2.penalties(i,:)=Pnode(md1.penalties(i,:));
[1348]177 end
178 md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
179 end
180
181 %recreate segments
[4348]182 if md1.dim==2
[8298]183 md2.nodeconnectivity=NodeConnectivity(md2.elements,md2.numberofnodes);
[2093]184 md2.elementconnectivity=ElementConnectivity(md2.elements,md2.nodeconnectivity);
[2530]185 md2.segments=contourenvelope(md2);
[8298]186 md2.nodeonboundary=zeros(numberofnodes2,1); md2.nodeonboundary(md2.segments(:,1:2))=1;
[1368]187 end
[1348]188
[1361]189 %Boundary conditions: Dirichlets on new boundary
[1348]190 %Catch the elements that have not been extracted
191 orphans_elem=find(~flag_elem);
[8298]192 orphans_node=unique(md1.elements(orphans_elem,:))';
193 %Figure out which node are on the boundary between md2 and md1
194 nodestoflag1=intersect(orphans_node,pos_node);
195 nodestoflag2=Pnode(nodestoflag1);
[8830]196 if numel(md1.spcvx)>1 & numel(md1.spcvy)>2 & numel(md1.spcvz)>2,
197 if numel(md1.vx_obs)>1 & numel(md1.vy_obs)>1
[8823]198 md2.spcvx(nodestoflag2)=md2.vx_obs(nodestoflag2);
199 md2.spcvy(nodestoflag2)=md2.vy_obs(nodestoflag2);
[1361]200 else
[8823]201 md2.spcvx(nodestoflag2)=NaN;
202 md2.spcvy(nodestoflag2)=NaN;
[1361]203 disp(' ')
[1755]204 disp('!! modelextract warning: spc values should be checked !!')
[1361]205 disp(' ')
206 end
[8833]207 %put 0 for vz
208 md2.spcvz(nodestoflag2)=0;
[1348]209 end
[1755]210 if ~isnan(md1.spctemperature),
[8298]211 md2.spctemperature(nodestoflag2,1)=1;
[1361]212 end
[1348]213
214 %Diagnostic
[1755]215 if ~isnan(md2.pressureload)
[8298]216 md2.pressureload(:,1)=Pnode(md1.pressureload(:,1));
217 md2.pressureload(:,2)=Pnode(md1.pressureload(:,2));
[3099]218 md2.pressureload(:,end-1)=Pelem(md1.pressureload(:,end-1));
[4348]219 if md1.dim==3
[8298]220 md2.pressureload(:,3)=Pnode(md1.pressureload(:,3));
221 md2.pressureload(:,4)=Pnode(md1.pressureload(:,4));
[1348]222 end
[2274]223 md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2) & md2.pressureload(:,end)),:);
[1348]224 end
225
226 %Results fields
227 if isstruct(md1.results),
[2530]228 md2.results=struct();
[1348]229 solutionfields=fields(md1.results);
230 for i=1:length(solutionfields),
231 %get subfields
232 solutionsubfields=fields(md1.results.(solutionfields(i)));
233 for j=1:length(solutionsubfields),
234 field=md1.results.(solutionfields(i)).(solutionsubfields(j));
[8298]235 if length(field)==numberofnodes1,
236 md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_node);
[1348]237 elseif length(field)==numberofelements1,
238 md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_elem);
[2530]239 else
240 md2.results.(solutionfields(i)).(solutionsubfields(j))=field;
[1348]241 end
242 end
243 end
244 end
245
[8298]246%Keep track of pos_node and pos_elem
247md2.extractednodes=pos_node;
[1348]248md2.extractedelements=pos_elem;
Note: See TracBrowser for help on using the repository browser.