| 1 | function md2=modelextract(md1,area)
|
|---|
| 2 | %modelextract - extract a model according to an Argus contour or flag list
|
|---|
| 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:
|
|---|
| 14 | % md2=modelextract(md1,area);
|
|---|
| 15 | %
|
|---|
| 16 | % Examples:
|
|---|
| 17 | % md2=modelextract(md,'Domain.exp');
|
|---|
| 18 | % md2=modelextract(md,md.elementoniceshelf);
|
|---|
| 19 | %
|
|---|
| 20 | % See also: EXTRUDE, COLLAPSE
|
|---|
| 21 |
|
|---|
| 22 | %some checks
|
|---|
| 23 | if ((nargin~=2) | (nargout~=1)),
|
|---|
| 24 | help modelextract
|
|---|
| 25 | error('modelextract error message: bad usage');
|
|---|
| 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
|
|---|
| 36 | flag_elem=FlagElements(md1,area);
|
|---|
| 37 |
|
|---|
| 38 | %kick out all elements with 3 dirichlets
|
|---|
| 39 | spc_elem=find(~flag_elem);
|
|---|
| 40 | spc_node=sort(unique(md1.elements(spc_elem,:)));
|
|---|
| 41 | flag=ones(md1.numberofnodes,1);
|
|---|
| 42 | flag(spc_node)=0;
|
|---|
| 43 | pos=find(sum(flag(md1.elements),2)==0);
|
|---|
| 44 | flag_elem(pos)=0;
|
|---|
| 45 |
|
|---|
| 46 | %extracted elements and nodes lists
|
|---|
| 47 | pos_elem=find(flag_elem);
|
|---|
| 48 | pos_node=sort(unique(md1.elements(pos_elem,:)));
|
|---|
| 49 |
|
|---|
| 50 | %keep track of some fields
|
|---|
| 51 | numberofnodes1=md1.numberofnodes;
|
|---|
| 52 | numberofelements1=md1.numberofelements;
|
|---|
| 53 | numberofnodes2=length(pos_node);
|
|---|
| 54 | numberofelements2=length(pos_elem);
|
|---|
| 55 | flag_node=zeros(numberofnodes1,1);
|
|---|
| 56 | flag_node(pos_node)=1;
|
|---|
| 57 |
|
|---|
| 58 | %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
|
|---|
| 59 | Pelem=zeros(numberofelements1,1);
|
|---|
| 60 | Pelem(pos_elem)=[1:numberofelements2]';
|
|---|
| 61 | Pnode=zeros(numberofnodes1,1);
|
|---|
| 62 | Pnode(pos_node)=[1:numberofnodes2]';
|
|---|
| 63 |
|
|---|
| 64 | %renumber the elements (some node won't exist anymore)
|
|---|
| 65 | elements_1=md1.elements;
|
|---|
| 66 | elements_2=elements_1(pos_elem,:);
|
|---|
| 67 | elements_2(:,1)=Pnode(elements_2(:,1));
|
|---|
| 68 | elements_2(:,2)=Pnode(elements_2(:,2));
|
|---|
| 69 | elements_2(:,3)=Pnode(elements_2(:,3));
|
|---|
| 70 | if md1.dim==3,
|
|---|
| 71 | elements_2(:,4)=Pnode(elements_2(:,4));
|
|---|
| 72 | elements_2(:,5)=Pnode(elements_2(:,5));
|
|---|
| 73 | elements_2(:,6)=Pnode(elements_2(:,6));
|
|---|
| 74 | end
|
|---|
| 75 |
|
|---|
| 76 | %OK, now create the new model !
|
|---|
| 77 |
|
|---|
| 78 | %take every fields from model
|
|---|
| 79 | md2=md1;
|
|---|
| 80 |
|
|---|
| 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
|
|---|
| 88 | field=md1.(model_fields(i));
|
|---|
| 89 | fieldsize=size(field);
|
|---|
| 90 |
|
|---|
| 91 | %size = number of nodes * n
|
|---|
| 92 | if fieldsize(1)==numberofnodes1
|
|---|
| 93 | md2.(model_fields(i))=field(pos_node,:);
|
|---|
| 94 |
|
|---|
| 95 | %size = number of elements * n
|
|---|
| 96 | elseif fieldsize(1)==numberofelements1
|
|---|
| 97 | md2.(model_fields(i))=field(pos_elem,:);
|
|---|
| 98 | end
|
|---|
| 99 | end
|
|---|
| 100 |
|
|---|
| 101 | %modify some specific fields
|
|---|
| 102 |
|
|---|
| 103 | %Mesh
|
|---|
| 104 | md2.numberofelements=numberofelements2;
|
|---|
| 105 | md2.numberofnodes=numberofnodes2;
|
|---|
| 106 | md2.elements=elements_2;
|
|---|
| 107 |
|
|---|
| 108 | %uppernodes lowernodes
|
|---|
| 109 | if md1.dim==3
|
|---|
| 110 | md2.uppernodes=md1.uppernodes(pos_node);
|
|---|
| 111 | pos=find(~isnan(md2.uppernodes));
|
|---|
| 112 | md2.uppernodes(pos)=Pnode(md2.uppernodes(pos));
|
|---|
| 113 |
|
|---|
| 114 | md2.lowernodes=md1.lowernodes(pos_node);
|
|---|
| 115 | pos=find(~isnan(md2.lowernodes));
|
|---|
| 116 | md2.lowernodes(pos)=Pnode(md2.lowernodes(pos));
|
|---|
| 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));
|
|---|
| 125 | end
|
|---|
| 126 |
|
|---|
| 127 | %Initial 2d mesh
|
|---|
| 128 | if md1.dim==3
|
|---|
| 129 | flag_elem_2d=flag_elem(1:md1.numberofelements2d);
|
|---|
| 130 | pos_elem_2d=find(flag_elem_2d);
|
|---|
| 131 | flag_node_2d=flag_node(1:md1.numberofnodes2d);
|
|---|
| 132 | pos_node_2d=find(flag_node_2d);
|
|---|
| 133 |
|
|---|
| 134 | md2.numberofelements2d=length(pos_elem_2d);
|
|---|
| 135 | md2.numberofnodes2d=length(pos_node_2d);
|
|---|
| 136 | md2.elements2d=md1.elements2d(pos_elem_2d,:);
|
|---|
| 137 | md2.elements2d(:,1)=Pnode(md2.elements2d(:,1));
|
|---|
| 138 | md2.elements2d(:,2)=Pnode(md2.elements2d(:,2));
|
|---|
| 139 | md2.elements2d(:,3)=Pnode(md2.elements2d(:,3));
|
|---|
| 140 |
|
|---|
| 141 | if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d); end;
|
|---|
| 142 | md2.x2d=md1.x(pos_node_2d);
|
|---|
| 143 | md2.y2d=md1.y(pos_node_2d);
|
|---|
| 144 | md2.z2d=md1.z(pos_node_2d);
|
|---|
| 145 | end
|
|---|
| 146 |
|
|---|
| 147 | %Edges
|
|---|
| 148 | if size(md2.edges,2)>1, %do not use ~isnan because there are some NaNs...
|
|---|
| 149 | %renumber first two columns
|
|---|
| 150 | pos=find(~isnan(md2.edges(:,4)));
|
|---|
| 151 | md2.edges(: ,1)=Pnode(md2.edges(:,1));
|
|---|
| 152 | md2.edges(: ,2)=Pnode(md2.edges(:,2));
|
|---|
| 153 | md2.edges(: ,3)=Pelem(md2.edges(:,3));
|
|---|
| 154 | md2.edges(pos,4)=Pelem(md2.edges(pos,4));
|
|---|
| 155 | %remove edges when the 2 vertices are not in the domain.
|
|---|
| 156 | md2.edges=md2.edges(find(md2.edges(:,1) & md2.edges(:,2)),:);
|
|---|
| 157 | %Replace all zeros by NaN in the last two columns;
|
|---|
| 158 | pos=find(md2.edges(:,3)==0);
|
|---|
| 159 | md2.edges(pos,3)=NaN;
|
|---|
| 160 | pos=find(md2.edges(:,4)==0);
|
|---|
| 161 | md2.edges(pos,4)=NaN;
|
|---|
| 162 | %Invert NaN of the third column with last column (Also invert first two columns!!)
|
|---|
| 163 | pos=find(isnan(md2.edges(:,3)));
|
|---|
| 164 | md2.edges(pos,3)=md2.edges(pos,4);
|
|---|
| 165 | md2.edges(pos,4)=NaN;
|
|---|
| 166 | values=md2.edges(pos,2);
|
|---|
| 167 | md2.edges(pos,2)=md2.edges(pos,1);
|
|---|
| 168 | md2.edges(pos,1)=values;
|
|---|
| 169 | %Finally remove edges that do not belong to any element
|
|---|
| 170 | pos=find(isnan(md2.edges(:,3)) & isnan(md2.edges(:,4)));
|
|---|
| 171 | md2.edges(pos,:)=[];
|
|---|
| 172 | end
|
|---|
| 173 |
|
|---|
| 174 | %Penalties
|
|---|
| 175 | if ~isnan(md2.penalties),
|
|---|
| 176 | for i=1:size(md1.penalties,1);
|
|---|
| 177 | md2.penalties(i,:)=Pnode(md1.penalties(i,:));
|
|---|
| 178 | end
|
|---|
| 179 | md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
|
|---|
| 180 | end
|
|---|
| 181 |
|
|---|
| 182 | %recreate segments
|
|---|
| 183 | if md1.dim==2
|
|---|
| 184 | md2.nodeconnectivity=NodeConnectivity(md2.elements,md2.numberofnodes);
|
|---|
| 185 | md2.elementconnectivity=ElementConnectivity(md2.elements,md2.nodeconnectivity);
|
|---|
| 186 | md2.segments=contourenvelope(md2);
|
|---|
| 187 | md2.nodeonboundary=zeros(numberofnodes2,1); md2.nodeonboundary(md2.segments(:,1:2))=1;
|
|---|
| 188 | end
|
|---|
| 189 |
|
|---|
| 190 | %Boundary conditions: Dirichlets on new boundary
|
|---|
| 191 | %Catch the elements that have not been extracted
|
|---|
| 192 | orphans_elem=find(~flag_elem);
|
|---|
| 193 | orphans_node=unique(md1.elements(orphans_elem,:))';
|
|---|
| 194 | %Figure out which node are on the boundary between md2 and md1
|
|---|
| 195 | nodestoflag1=intersect(orphans_node,pos_node);
|
|---|
| 196 | nodestoflag2=Pnode(nodestoflag1);
|
|---|
| 197 | if ~isnan(md1.spcvelocity),
|
|---|
| 198 | md2.spcvelocity(nodestoflag2,1:3)=1;
|
|---|
| 199 | if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
|
|---|
| 200 | md2.spcvelocity(nodestoflag2,4)=md2.vx_obs(nodestoflag2);
|
|---|
| 201 | md2.spcvelocity(nodestoflag2,5)=md2.vy_obs(nodestoflag2);
|
|---|
| 202 | else
|
|---|
| 203 | md2.spcvelocity(nodestoflag2,4:5)=zeros(length(nodestoflag2),2);
|
|---|
| 204 | disp(' ')
|
|---|
| 205 | disp('!! modelextract warning: spc values should be checked !!')
|
|---|
| 206 | disp(' ')
|
|---|
| 207 | end
|
|---|
| 208 | end
|
|---|
| 209 | if ~isnan(md1.spctemperature),
|
|---|
| 210 | md2.spctemperature(nodestoflag2,1)=1;
|
|---|
| 211 | if ~isnan(md1.observed_temperature)
|
|---|
| 212 | md2.spctemperature(nodestoflag2,2)=md2.observed_temperature(nodestoflag2);
|
|---|
| 213 | else
|
|---|
| 214 | md2.spctemperature(nodestoflag2,2)=zeros(length(nodestoflag2),2);
|
|---|
| 215 | disp(' ')
|
|---|
| 216 | disp('!! modelextract warning: spc values should be checked !!')
|
|---|
| 217 | disp(' ')
|
|---|
| 218 | end
|
|---|
| 219 | end
|
|---|
| 220 |
|
|---|
| 221 | %Diagnostic
|
|---|
| 222 | if ~isnan(md2.pressureload)
|
|---|
| 223 | md2.pressureload(:,1)=Pnode(md1.pressureload(:,1));
|
|---|
| 224 | md2.pressureload(:,2)=Pnode(md1.pressureload(:,2));
|
|---|
| 225 | md2.pressureload(:,end-1)=Pelem(md1.pressureload(:,end-1));
|
|---|
| 226 | if md1.dim==3
|
|---|
| 227 | md2.pressureload(:,3)=Pnode(md1.pressureload(:,3));
|
|---|
| 228 | md2.pressureload(:,4)=Pnode(md1.pressureload(:,4));
|
|---|
| 229 | end
|
|---|
| 230 | md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2) & md2.pressureload(:,end)),:);
|
|---|
| 231 | end
|
|---|
| 232 |
|
|---|
| 233 | %Results fields
|
|---|
| 234 | if isstruct(md1.results),
|
|---|
| 235 | md2.results=struct();
|
|---|
| 236 | solutionfields=fields(md1.results);
|
|---|
| 237 | for i=1:length(solutionfields),
|
|---|
| 238 | %get subfields
|
|---|
| 239 | solutionsubfields=fields(md1.results.(solutionfields(i)));
|
|---|
| 240 | for j=1:length(solutionsubfields),
|
|---|
| 241 | field=md1.results.(solutionfields(i)).(solutionsubfields(j));
|
|---|
| 242 | if length(field)==numberofnodes1,
|
|---|
| 243 | md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_node);
|
|---|
| 244 | elseif length(field)==numberofelements1,
|
|---|
| 245 | md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_elem);
|
|---|
| 246 | else
|
|---|
| 247 | md2.results.(solutionfields(i)).(solutionsubfields(j))=field;
|
|---|
| 248 | end
|
|---|
| 249 | end
|
|---|
| 250 | end
|
|---|
| 251 | end
|
|---|
| 252 |
|
|---|
| 253 | %reinitialize output parameters
|
|---|
| 254 | md2.viscousheating=NaN;
|
|---|
| 255 | md2.pressure_elem=NaN;
|
|---|
| 256 | md2.stress=NaN;
|
|---|
| 257 | md2.stress_surface=NaN;
|
|---|
| 258 | md2.stress_bed=NaN;
|
|---|
| 259 | md2.deviatoricstress=NaN;
|
|---|
| 260 | md2.strainrate=NaN;
|
|---|
| 261 |
|
|---|
| 262 | %Keep track of pos_node and pos_elem
|
|---|
| 263 | md2.extractednodes=pos_node;
|
|---|
| 264 | md2.extractedelements=pos_elem;
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
