[1495] | 1 | function 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] | 23 | if ((nargin~=2) | (nargout~=1)),
|
---|
[3589] | 24 | help modelextract
|
---|
| 25 | error('modelextract error message: bad usage');
|
---|
[1348] | 26 | end
|
---|
| 27 |
|
---|
| 28 | %get check option
|
---|
| 29 | if (nargin==3 & varargin{1}==0),
|
---|
| 30 | checkoutline=0;
|
---|
| 31 | else
|
---|
| 32 | checkoutline=1;
|
---|
| 33 | end
|
---|
| 34 |
|
---|
| 35 | %first
|
---|
[2988] | 36 | flag_elem=FlagElements(md1,area);
|
---|
[1348] | 37 |
|
---|
[1361] | 38 | %kick out all elements with 3 dirichlets
|
---|
| 39 | spc_elem=find(~flag_elem);
|
---|
[8298] | 40 | spc_node=sort(unique(md1.elements(spc_elem,:)));
|
---|
| 41 | flag=ones(md1.numberofnodes,1);
|
---|
| 42 | flag(spc_node)=0;
|
---|
[1361] | 43 | pos=find(sum(flag(md1.elements),2)==0);
|
---|
| 44 | flag_elem(pos)=0;
|
---|
| 45 |
|
---|
| 46 | %extracted elements and nodes lists
|
---|
[1348] | 47 | pos_elem=find(flag_elem);
|
---|
[8298] | 48 | pos_node=sort(unique(md1.elements(pos_elem,:)));
|
---|
[1348] | 49 |
|
---|
| 50 | %keep track of some fields
|
---|
[8298] | 51 | numberofnodes1=md1.numberofnodes;
|
---|
[1348] | 52 | numberofelements1=md1.numberofelements;
|
---|
[8298] | 53 | numberofnodes2=length(pos_node);
|
---|
[1348] | 54 | numberofelements2=length(pos_elem);
|
---|
[8298] | 55 | flag_node=zeros(numberofnodes1,1);
|
---|
| 56 | flag_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] | 59 | Pelem=zeros(numberofelements1,1);
|
---|
| 60 | Pelem(pos_elem)=[1:numberofelements2]';
|
---|
[8298] | 61 | Pnode=zeros(numberofnodes1,1);
|
---|
| 62 | Pnode(pos_node)=[1:numberofnodes2]';
|
---|
[1348] | 63 |
|
---|
[8298] | 64 | %renumber the elements (some node won't exist anymore)
|
---|
[1348] | 65 | elements_1=md1.elements;
|
---|
| 66 | elements_2=elements_1(pos_elem,:);
|
---|
[8298] | 67 | elements_2(:,1)=Pnode(elements_2(:,1));
|
---|
| 68 | elements_2(:,2)=Pnode(elements_2(:,2));
|
---|
| 69 | elements_2(:,3)=Pnode(elements_2(:,3));
|
---|
[4348] | 70 | if 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] | 74 | end
|
---|
| 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,:);
|
---|
[9627] | 94 | elseif (fieldsize(1)==numberofnodes1+1)
|
---|
| 95 | md2.(model_fields(i))=[field(pos_node,:); field(end,:)];
|
---|
[1361] | 96 | %size = number of elements * n
|
---|
| 97 | elseif fieldsize(1)==numberofelements1
|
---|
[9436] | 98 | md2.(model_fields{i})=field(pos_elem,:);
|
---|
[1361] | 99 | end
|
---|
| 100 | end
|
---|
| 101 |
|
---|
[1348] | 102 | %modify some specific fields
|
---|
| 103 |
|
---|
| 104 | %Mesh
|
---|
| 105 | md2.numberofelements=numberofelements2;
|
---|
[8298] | 106 | md2.numberofnodes=numberofnodes2;
|
---|
[1348] | 107 | md2.elements=elements_2;
|
---|
| 108 |
|
---|
[1368] | 109 | %uppernodes lowernodes
|
---|
[4348] | 110 | if md1.dim==3
|
---|
[8298] | 111 | md2.uppernodes=md1.uppernodes(pos_node);
|
---|
| 112 | pos=find(~isnan(md2.uppernodes));
|
---|
| 113 | md2.uppernodes(pos)=Pnode(md2.uppernodes(pos));
|
---|
[1368] | 114 |
|
---|
[8298] | 115 | md2.lowernodes=md1.lowernodes(pos_node);
|
---|
| 116 | pos=find(~isnan(md2.lowernodes));
|
---|
| 117 | md2.lowernodes(pos)=Pnode(md2.lowernodes(pos));
|
---|
[4690] | 118 |
|
---|
| 119 | md2.upperelements=md1.upperelements(pos_elem);
|
---|
| 120 | pos=find(~isnan(md2.upperelements));
|
---|
| 121 | md2.upperelements(pos)=Pelem(md2.upperelements(pos));
|
---|
| 122 |
|
---|
| 123 | md2.lowerelements=md1.lowerelements(pos_elem);
|
---|
| 124 | pos=find(~isnan(md2.lowerelements));
|
---|
| 125 | md2.lowerelements(pos)=Pelem(md2.lowerelements(pos));
|
---|
[1368] | 126 | end
|
---|
| 127 |
|
---|
[1348] | 128 | %Initial 2d mesh
|
---|
[4348] | 129 | if md1.dim==3
|
---|
[1348] | 130 | flag_elem_2d=flag_elem(1:md1.numberofelements2d);
|
---|
| 131 | pos_elem_2d=find(flag_elem_2d);
|
---|
[8298] | 132 | flag_node_2d=flag_node(1:md1.numberofnodes2d);
|
---|
| 133 | pos_node_2d=find(flag_node_2d);
|
---|
[1348] | 134 |
|
---|
| 135 | md2.numberofelements2d=length(pos_elem_2d);
|
---|
[8298] | 136 | md2.numberofnodes2d=length(pos_node_2d);
|
---|
[1348] | 137 | md2.elements2d=md1.elements2d(pos_elem_2d,:);
|
---|
[8298] | 138 | md2.elements2d(:,1)=Pnode(md2.elements2d(:,1));
|
---|
| 139 | md2.elements2d(:,2)=Pnode(md2.elements2d(:,2));
|
---|
| 140 | md2.elements2d(:,3)=Pnode(md2.elements2d(:,3));
|
---|
[1348] | 141 |
|
---|
[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
|
---|
| 247 | md2.extractednodes=pos_node;
|
---|
[1348] | 248 | md2.extractedelements=pos_elem;
|
---|