Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 13856) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 13857) @@ -38,5 +38,5 @@ pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0] if not numpy.size(pos): - print("SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually.") + print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually." md.diagnostic.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) @@ -50,9 +50,9 @@ #Dirichlet Values if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices: - print(" boundary conditions for diagnostic model: spc set as observed velocities") + print " boundary conditions for diagnostic model: spc set as observed velocities" md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos] md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos] else: - print(" boundary conditions for diagnostic model: spc set as zero") + print " boundary conditions for diagnostic model: spc set as zero" md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2)) @@ -83,14 +83,14 @@ if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1): md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) - print(" no surfaceforcings.precipitation specified: values set as zero") + print " no surfaceforcings.precipitation specified: values set as zero" if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0): md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1)) - print(" no surfaceforcings.mass_balance specified: values set as zero") + print " no surfaceforcings.mass_balance specified: values set as zero" if numpy.all(numpy.isnan(md.basalforcings.melting_rate)): md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) - print(" no basalforcings.melting_rate specified: values set as zero") + print " no basalforcings.melting_rate specified: values set as zero" if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) - print(" no balancethickness.thickening_rate specified: values set as zero") + print " no balancethickness.thickening_rate specified: values set as zero" md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) @@ -104,7 +104,7 @@ if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) - md.basalforcings.geothermalflux[numpy.nonzero(md.mask.vertexongroundedice)]=50.*10.**-3; #50mW/m2 + md.basalforcings.geothermalflux[numpy.nonzero(md.mask.vertexongroundedice)]=50.*10.**-3 #50mW/m2 else: - print(" no thermal boundary conditions created: no observed temperature found"); + print " no thermal boundary conditions created: no observed temperature found" return md Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 13856) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 13857) @@ -252,5 +252,4 @@ % an empty string '' will be considered as an empty domain % a string 'all' will be considered as the entire domain - % add an argument 0 if you do not want the elements to be checked (faster) % % Usage: @@ -272,11 +271,4 @@ end - %get check option - if (nargin==3 & varargin{1}==0), - checkoutline=0; - else - checkoutline=1; - end - %get elements that are inside area flag_elem=FlagElements(md1,area); @@ -311,5 +303,5 @@ Pnode(pos_node)=[1:numberofvertices2]'; - %renumber the elements (some node won't exist anymore) + %renumber the elements (some nodes won't exist anymore) elements_1=md1.mesh.elements; elements_2=elements_1(pos_elem,:); @@ -323,7 +315,7 @@ end - %OK, now create the new model ! - - %take every fields from model + %OK, now create the new model! + + %take every field from model md2=md1; @@ -347,5 +339,5 @@ elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); @@ -358,5 +350,5 @@ elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i})=field(pos_elem,:); @@ -413,11 +405,11 @@ %renumber first two columns pos=find(md2.mesh.edges(:,4)~=-1); - md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); - md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); + md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); + md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); %remove edges when the 2 vertices are not in the domain. md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); - %Replace all zeros by -1 in the last two columns; + %Replace all zeros by -1 in the last two columns pos=find(md2.mesh.edges(:,3)==0); md2.mesh.edges(pos,3)=-1; Index: /issm/trunk-jpl/src/m/classes/model/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.py (revision 13856) +++ /issm/trunk-jpl/src/m/classes/model/model.py (revision 13857) @@ -1,4 +1,5 @@ #module imports {{{ import numpy +import copy from mesh import mesh from mask import mask @@ -39,4 +40,8 @@ from iluasmoptions import * from project3d import * +from FlagElements import * +from NodeConnectivity import * +from ElementConnectivity import * +from contourenvelope import * #}}} @@ -170,293 +175,258 @@ # }}} - """ - function md2 = extract(md,area) % {{{ - %extract - extract a model according to an Argus contour or flag list - % - % This routine extracts a submodel from a bigger model with respect to a given contour - % md must be followed by the corresponding exp file or flags list - % It can either be a domain file (argus type, .exp extension), or an array of element flags. - % If user wants every element outside the domain to be - % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); - % an empty string '' will be considered as an empty domain - % a string 'all' will be considered as the entire domain - % add an argument 0 if you do not want the elements to be checked (faster) - % - % Usage: - % md2=extract(md,area); - % - % Examples: - % md2=extract(md,'Domain.exp'); - % md2=extract(md,md.mask.elementonfloatingice); - % - % See also: EXTRUDE, COLLAPSE - - %copy model - md1=md; - - %some checks - if ((nargin~=2) | (nargout~=1)), - help extract - error('extract error message: bad usage'); - end - - %get check option - if (nargin==3 & varargin{1}==0), - checkoutline=0; - else - checkoutline=1; - end - - %get elements that are inside area - flag_elem=FlagElements(md1,area); - if ~any(flag_elem), - error('extracted model is empty'); - end - - %kick out all elements with 3 dirichlets - spc_elem=find(~flag_elem); - spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); - flag=ones(md1.mesh.numberofvertices,1); - flag(spc_node)=0; - pos=find(sum(flag(md1.mesh.elements),2)==0); - flag_elem(pos)=0; - - %extracted elements and nodes lists - pos_elem=find(flag_elem); - pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); - - %keep track of some fields - numberofvertices1=md1.mesh.numberofvertices; - numberofelements1=md1.mesh.numberofelements; - numberofvertices2=length(pos_node); - numberofelements2=length(pos_elem); - flag_node=zeros(numberofvertices1,1); - flag_node(pos_node)=1; - - %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) - Pelem=zeros(numberofelements1,1); - Pelem(pos_elem)=[1:numberofelements2]'; - Pnode=zeros(numberofvertices1,1); - Pnode(pos_node)=[1:numberofvertices2]'; - - %renumber the elements (some node won't exist anymore) - elements_1=md1.mesh.elements; - elements_2=elements_1(pos_elem,:); - elements_2(:,1)=Pnode(elements_2(:,1)); - elements_2(:,2)=Pnode(elements_2(:,2)); - elements_2(:,3)=Pnode(elements_2(:,3)); - if md1.mesh.dimension==3, - elements_2(:,4)=Pnode(elements_2(:,4)); - elements_2(:,5)=Pnode(elements_2(:,5)); - elements_2(:,6)=Pnode(elements_2(:,6)); - end - - %OK, now create the new model ! - - %take every fields from model - md2=md1; - - %automatically modify fields - - %loop over model fields - model_fields=fields(md1); - for i=1:length(model_fields), - %get field - field=md1.(model_fields{i}); - fieldsize=size(field); - if isobject(field), %recursive call - object_fields=fields(md1.(model_fields{i})); - for j=1:length(object_fields), - %get field - field=md1.(model_fields{i}).(object_fields{j}); - fieldsize=size(field); - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); - end - end - else - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i})=field(pos_elem,:); - end - end - end - - %modify some specific fields - - %Mesh - md2.mesh.numberofelements=numberofelements2; - md2.mesh.numberofvertices=numberofvertices2; - md2.mesh.elements=elements_2; - - %mesh.uppervertex mesh.lowervertex - if md1.mesh.dimension==3 - md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); - pos=find(~isnan(md2.mesh.uppervertex)); - md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); - - md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); - pos=find(~isnan(md2.mesh.lowervertex)); - md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); - - md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); - pos=find(~isnan(md2.mesh.upperelements)); - md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); - - md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); - pos=find(~isnan(md2.mesh.lowerelements)); - md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); - end - - %Initial 2d mesh - if md1.mesh.dimension==3 - flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); - pos_elem_2d=find(flag_elem_2d); - flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); - pos_node_2d=find(flag_node_2d); - - md2.mesh.numberofelements2d=length(pos_elem_2d); - md2.mesh.numberofvertices2d=length(pos_node_2d); - md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); - md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); - md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); - md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); - - md2.mesh.x2d=md1.mesh.x(pos_node_2d); - md2.mesh.y2d=md1.mesh.y(pos_node_2d); - end - - %Edges - if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... - %renumber first two columns - pos=find(md2.mesh.edges(:,4)~=-1); - md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); - md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); - md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); - md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); - %remove edges when the 2 vertices are not in the domain. - md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); - %Replace all zeros by -1 in the last two columns; - pos=find(md2.mesh.edges(:,3)==0); - md2.mesh.edges(pos,3)=-1; - pos=find(md2.mesh.edges(:,4)==0); - md2.mesh.edges(pos,4)=-1; - %Invert -1 on the third column with last column (Also invert first two columns!!) - pos=find(md2.mesh.edges(:,3)==-1); - md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); - md2.mesh.edges(pos,4)=-1; - values=md2.mesh.edges(pos,2); - md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); - md2.mesh.edges(pos,1)=values; - %Finally remove edges that do not belong to any element - pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); - md2.mesh.edges(pos,:)=[]; - end - - %Penalties - if ~isnan(md2.diagnostic.vertex_pairing), - for i=1:size(md1.diagnostic.vertex_pairing,1); - md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); - end - md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); - end - if ~isnan(md2.prognostic.vertex_pairing), - for i=1:size(md1.prognostic.vertex_pairing,1); - md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); - end - md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); - end - - %recreate segments - if md1.mesh.dimension==2 - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - else - %First do the connectivity for the contourenvelope in 2d - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); - %Then do it for 3d as usual - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - end - - %Boundary conditions: Dirichlets on new boundary - %Catch the elements that have not been extracted - orphans_elem=find(~flag_elem); - orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; - %Figure out which node are on the boundary between md2 and md1 - nodestoflag1=intersect(orphans_node,pos_node); - nodestoflag2=Pnode(nodestoflag1); - if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, - if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 - md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); - md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); - else - md2.diagnostic.spcvx(nodestoflag2)=NaN; - md2.diagnostic.spcvy(nodestoflag2)=NaN; - disp(' ') - disp('!! extract warning: spc values should be checked !!') - disp(' ') - end - %put 0 for vz - md2.diagnostic.spcvz(nodestoflag2)=0; - end - if ~isnan(md1.thermal.spctemperature), - md2.thermal.spctemperature(nodestoflag2,1)=1; - end - - %Diagnostic - if ~isnan(md2.diagnostic.icefront) - md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); - md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); - md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); - if md1.mesh.dimension==3 - md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); - md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); - end - md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); - end - - %Results fields - if isstruct(md1.results), - md2.results=struct(); - solutionfields=fields(md1.results); - for i=1:length(solutionfields), - %get subfields - solutionsubfields=fields(md1.results.(solutionfields{i})); - for j=1:length(solutionsubfields), - field=md1.results.(solutionfields{i}).(solutionsubfields{j}); - if length(field)==numberofvertices1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); - elseif length(field)==numberofelements1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); - else - md2.results.(solutionfields{i}).(solutionsubfields{j})=field; - end - end - end - end - - %Keep track of pos_node and pos_elem - md2.mesh.extractedvertices=pos_node; - md2.mesh.extractedelements=pos_elem; - end % }}} - """ + def extract(md,area): # {{{ + """ + extract - extract a model according to an Argus contour or flag list + + This routine extracts a submodel from a bigger model with respect to a given contour + md must be followed by the corresponding exp file or flags list + It can either be a domain file (argus type, .exp extension), or an array of element flags. + If user wants every element outside the domain to be + extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + an empty string '' will be considered as an empty domain + a string 'all' will be considered as the entire domain + + Usage: + md2=extract(md,area); + + Examples: + md2=extract(md,'Domain.exp'); + md2=extract(md,md.mask.elementonfloatingice); + + See also: EXTRUDE, COLLAPSE + """ + + #copy model + md1=copy.deepcopy(md) + + #get elements that are inside area + flag_elem=FlagElements(md1,area) + if not numpy.any(flag_elem): + raise RuntimeError("extracted model is empty") + + #kick out all elements with 3 dirichlets + spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] + spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1 + flag=numpy.ones(md1.mesh.numberofvertices) + flag[spc_node]=0 + pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0] + flag_elem[pos]=0 + + #extracted elements and nodes lists + pos_elem=numpy.nonzero(flag_elem)[0] + pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1 + + #keep track of some fields + numberofvertices1=md1.mesh.numberofvertices + numberofelements1=md1.mesh.numberofelements + numberofvertices2=numpy.size(pos_node) + numberofelements2=numpy.size(pos_elem) + flag_node=numpy.zeros(numberofvertices1) + flag_node[pos_node]=1 + + #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) + Pelem=numpy.zeros(numberofelements1) + Pelem[pos_elem]=numpy.arange(1,numberofelements2+1) + Pnode=numpy.zeros(numberofvertices1) + Pnode[pos_node]=numpy.arange(1,numberofvertices2+1) + + #renumber the elements (some node won't exist anymore) + elements_1=copy.deepcopy(md1.mesh.elements) + elements_2=elements_1[pos_elem,:] + elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1] + elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1] + elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1] + if md1.mesh.dimension==3: + elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1] + elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1] + elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1] + + #OK, now create the new model! + + #take every field from model + md2=copy.deepcopy(md1) + + #automatically modify fields + + #loop over model fields + model_fields=vars(md1) + for fieldi in model_fields: + #get field + field=getattr(md1,fieldi) + fieldsize=numpy.shape(field) + if hasattr(field,'__dict__') and not ismember(fieldi,['results'])[0]: #recursive call + object_fields=vars(field) + for fieldj in object_fields: + #get field + field=getattr(getattr(md1,fieldi),fieldj) + fieldsize=numpy.shape(field) + if len(fieldsize): + #size = number of nodes * n + if fieldsize[0]==numberofvertices1: + setattr(getattr(md2,fieldi),fieldj,field[pos_node,:]) + elif fieldsize[0]==numberofvertices1+1: + setattr(getattr(md2,fieldi),fieldj,numpy.vstack((field[pos_node,:],field[-1,:]))) + #size = number of elements * n + elif fieldsize[0]==numberofelements1: + setattr(getattr(md2,fieldi),fieldj,field[pos_elem,:]) + else: + if len(fieldsize): + #size = number of nodes * n + if fieldsize[0]==numberofvertices1: + setattr(md2,fieldi,field[pos_node,:]) + elif fieldsize[0]==numberofvertices1+1: + setattr(md2,fieldi,numpy.hstack((field[pos_node,:],field[-1,:]))) + #size = number of elements * n + elif fieldsize[0]==numberofelements1: + setattr(md2,fieldi,field[pos_elem,:]) + + #modify some specific fields + + #Mesh + md2.mesh.numberofelements=numberofelements2 + md2.mesh.numberofvertices=numberofvertices2 + md2.mesh.elements=elements_2 + + #mesh.uppervertex mesh.lowervertex + if md1.mesh.dimension==3: + md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node] + pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.uppervertex)))[0] + md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1] + + md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node] + pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowervertex)))[0] + md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1] + + md2.mesh.upperelements=md1.mesh.upperelements[pos_elem] + pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.upperelements)))[0] + md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1] + + md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem] + pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowerelements)))[0] + md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1] + + #Initial 2d mesh + if md1.mesh.dimension==3: + flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)] + pos_elem_2d=numpy.nonzero(flag_elem_2d)[0] + flag_node_2d=flag_node[numpy.arange(0,md1.mesh.numberofvertices2d)] + pos_node_2d=numpy.nonzero(flag_node_2d)[0] + + md2.mesh.numberofelements2d=numpy.size(pos_elem_2d) + md2.mesh.numberofvertices2d=numpy.size(pos_node_2d) + md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:] + md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1] + md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1] + md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1] + + md2.mesh.x2d=md1.mesh.x[pos_node_2d] + md2.mesh.y2d=md1.mesh.y[pos_node_2d] + + #Edges + if len(numpy.shape(md2.mesh.edges))>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs... + #renumber first two columns + pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0] + md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0].astype(int)-1] + md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1] + md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1] + md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1] + #remove edges when the 2 vertices are not in the domain. + md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:] + #Replace all zeros by -1 in the last two columns + pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0] + md2.mesh.edges[pos,2]=-1 + pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0] + md2.mesh.edges[pos,3]=-1 + #Invert -1 on the third column with last column (Also invert first two columns!!) + pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0] + md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3] + md2.mesh.edges[pos,3]=-1 + values=md2.mesh.edges[pos,1] + md2.mesh.edges[pos,1]=md2.mesh.edges[pos,0] + md2.mesh.edges[pos,0]=values + #Finally remove edges that do not belong to any element + pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0] + md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0) + + #Penalties + if numpy.any(numpy.logical_not(numpy.isnan(md2.diagnostic.vertex_pairing))): + for i in xrange(numpy.size(md1.diagnostic.vertex_pairing,axis=0)): + md2.diagnostic.vertex_pairing[i,:]=Pnode[md1.diagnostic.vertex_pairing[i,:]] + md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing[numpy.nonzero(md2.diagnostic.vertex_pairing[:,0])[0],:] + if numpy.any(numpy.logical_not(numpy.isnan(md2.prognostic.vertex_pairing))): + for i in xrange(numpy.size(md1.prognostic.vertex_pairing,axis=0)): + md2.prognostic.vertex_pairing[i,:]=Pnode[md1.prognostic.vertex_pairing[i,:]] + md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing[numpy.nonzero(md2.prognostic.vertex_pairing[:,0])[0],:] + + #recreate segments + if md1.mesh.dimension==2: + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity) + md2.mesh.segments=contourenvelope(md2) + md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2) + md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1 + else: + #First do the connectivity for the contourenvelope in 2d + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity) + md2.mesh.segments=contourenvelope(md2) + md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers) + md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1 + md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) + #Then do it for 3d as usual + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity) + + #Boundary conditions: Dirichlets on new boundary + #Catch the elements that have not been extracted + orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] + orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1 + #Figure out which node are on the boundary between md2 and md1 + nodestoflag1=numpy.intersect1d(orphans_node,pos_node) + nodestoflag2=Pnode[nodestoflag1].astype(int)-1 + if numpy.size(md1.diagnostic.spcvx)>1 and numpy.size(md1.diagnostic.spcvy)>2 and numpy.size(md1.diagnostic.spcvz)>2: + if numpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1: + md2.diagnostic.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] + md2.diagnostic.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2] + else: + md2.diagnostic.spcvx[nodestoflag2]=float('NaN') + md2.diagnostic.spcvy[nodestoflag2]=float('NaN') + print "\n!! extract warning: spc values should be checked !!\n\n" + #put 0 for vz + md2.diagnostic.spcvz[nodestoflag2]=0 + if numpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))): + md2.thermal.spctemperature[nodestoflag2,0]=1 + + #Diagnostic + if numpy.any(numpy.logical_not(numpy.isnan(md2.diagnostic.icefront))): + md2.diagnostic.icefront[:,0]=Pnode[md1.diagnostic.icefront[:,0].astype(int)-1] + md2.diagnostic.icefront[:,1]=Pnode[md1.diagnostic.icefront[:,1].astype(int)-1] + md2.diagnostic.icefront[:,-2]=Pelem[md1.diagnostic.icefront[:,-2].astype(int)-1] + if md1.mesh.dimension==3: + md2.diagnostic.icefront[:,2]=Pnode[md1.diagnostic.icefront[:,2].astype(int)-1] + md2.diagnostic.icefront[:,3]=Pnode[md1.diagnostic.icefront[:,3].astype(int)-1] + md2.diagnostic.icefront=md2.diagnostic.icefront[numpy.nonzero(numpy.logical_and(numpy.logical_and(md2.diagnostic.icefront[:,0],md2.diagnostic.icefront[:,1]),md2.diagnostic.icefront[:,-1]))[0],:] + + #Results fields + if md1.results: + md2.results=OrderedDict() + for solutionfield in md1.results.iterkeys(): + #get time step + for i in m1.results[solutionfield].iterkeys(): + #get subfields + for solutionsubfield,field in m1.results[solutionfield][i].iteritems(): + if numpy.size(field)==numberofvertices1: + md2.results[solutionfield][i][solutionsubfield]=field[pos_node] + elif numpy.size(field)==numberofelements1: + md2.results[solutionfield][i][solutionsubfield]=field[pos_elem] + else: + md2.results[solutionfield][i][solutionsubfield]=field + + #Keep track of pos_node and pos_elem + md2.mesh.extractedvertices=pos_node.astype(float)+1 + md2.mesh.extractedelements=pos_elem.astype(float)+1 + + return md2 + # }}} def extrude(md,*args): # {{{ Index: /issm/trunk-jpl/src/m/classes/verbose.py =================================================================== --- /issm/trunk-jpl/src/m/classes/verbose.py (revision 13856) +++ /issm/trunk-jpl/src/m/classes/verbose.py (revision 13857) @@ -61,5 +61,5 @@ #Cast to logicals listproperties=vars(self) - for [fieldname,fieldvalue] in listproperties.iteritems(): + for fieldname,fieldvalue in listproperties.iteritems(): if isinstance(fieldvalue,bool) or isinstance(fieldvalue,(int,long,float)): setattr(self,fieldname,bool(fieldvalue)) Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13856) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13857) @@ -123,5 +123,5 @@ elif numpy.any(numpy.logical_not(flags)): #We LOTS of work to do - print("Rift tip outside of or on the domain has been detected and is being processed...") + print "Rift tip outside of or on the domain has been detected and is being processed..." #check that only one point is outside (for now) @@ -166,5 +166,5 @@ if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0): - print("moving tip-domain intersection point") + print "moving tip-domain intersection point" #Get position of the closer point @@ -350,5 +350,5 @@ raise RuntimeError("bamg.py/processgeometry is not complete.") #Deal with edges - print("Checking Edge crossing...") + print "Checking Edge crossing..." i=0 while (i1: + raise RuntimeError("contourenvelope error message: bad usage") + + if len(args)==1: + flags=args[0] + + if isinstance(flags,(str,unicode)): + file=flags + if not os.path.exists(file): + raise IOError("contourenvelope error message: file '%s' not found" % file) + isfile=1 + elif isinstance(flags,(bool,int,long,float)): + #do nothing for now + isfile=0 + else: + raise TypeError("contourenvelope error message: second argument should be a file or an elements flag") + + #Now, build the connectivity tables for this mesh. + #Computing connectivity + if numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d: + [md.mesh.vertexconnectivity]=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices) + if numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d: + [md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity) + + #get nodes inside profile + mesh.elementconnectivity=copy.deepcopy(md.mesh.elementconnectivity) + if md.mesh.dimension==2: + mesh.elements=copy.deepcopy(md.mesh.elements) + mesh.x=copy.deepcopy(md.mesh.x) + mesh.y=copy.deepcopy(md.mesh.y) + mesh.numberofvertices=copy.deepcopy(md.mesh.numberofvertices) + mesh.numberofelements=copy.deepcopy(md.mesh.numberofelements) + else: + mesh.elements=copy.deepcopy(md.mesh.elements2d) + mesh.x=copy.deepcopy(md.mesh.x2d) + mesh.y=copy.deepcopy(md.mesh.y2d) + mesh.numberofvertices=copy.deepcopy(md.mesh.numberofvertices2d) + mesh.numberofelements=copy.deepcopy(md.mesh.numberofelements2d) + + if len(args)==1: + + if isfile: + #get flag list of elements and nodes inside the contour + nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1) + elemin=(numpy.sum(nodein(mesh.elements),axis=1)==numpy.size(mesh.elements,axis=1)) + #modify element connectivity + elemout=numpy.nonzero(numpy.logical_not(elemin))[0] + mesh.elementconnectivity[elemout,:]=0 + mesh.elementconnectivity[numpy.nonzero(ismember(mesh.elementconnectivity,elemout+1))]=0 + else: + #get flag list of elements and nodes inside the contour + nodein=numpy.zeros(mesh.numberofvertices) + elemin=numpy.zeros(mesh.numberofelements) + + pos=numpy.nonzero(flags) + elemin[pos]=1 + nodein[mesh.elements[pos,:].astype(int)-1]=1 + + #modify element connectivity + elemout=numpy.nonzero(numpy.logical_not(elemin))[0] + mesh.elementconnectivity[elemout,:]=0 + mesh.elementconnectivity[numpy.nonzero(ismember(mesh.elementconnectivity,elemout+1))]=0 + + #Find element on boundary + #First: find elements on the boundary of the domain + flag=copy.deepcopy(mesh.elementconnectivity) + if len(args)==1: + flag[numpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]] + elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0).astype(float) + + #Find segments on boundary + pos=numpy.nonzero(elementonboundary)[0] + num_segments=numpy.size(pos) + segments=numpy.zeros((num_segments*3,3)) + count=0 + + for el1 in pos: + els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1 + if numpy.size(els2)>1: + flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:]) + nods1=mesh.elements[el1,:] + nods1=numpy.delete(nods1,numpy.nonzero(nods1==flag)) + segments[count,:]=[nods1[0],nods1[1],el1+1] + + ord1=numpy.nonzero(nods1[0]==mesh.elements[el1,:])[0][0] + ord2=numpy.nonzero(nods1[1]==mesh.elements[el1,:])[0][0] + + #swap segment nodes if necessary + if ( (ord1==0 and ord2==1) or (ord1==1 and ord2==2) or (ord1==2 and ord2==0) ): + temp=segments[count,0] + segments[count,0]=segments[count,1] + segments[count,1]=temp + segments[count,0:2]=numpy.flipud(segments[count,0:2]) + count+=1 + else: + nods1=mesh.elements[el1,:] + flag=numpy.setdiff1d(nods1,mesh.elements[els2,:]) + for j in xrange(0,3): + nods=numpy.delete(nods1,j) + if numpy.any(ismember(flag,nods)): + segments[count,:]=[nods[0],nods[1],el1+1] + ord1=numpy.nonzero(nods[0]==mesh.elements[el1,:])[0][0] + ord2=numpy.nonzero(nods[1]==mesh.elements[el1,:])[0][0] + if ( (ord1==0 and ord2==1) or (ord1==1 and ord2==2) or (ord1==2 and ord2==0) ): + temp=segments[count,0] + segments[count,0]=segments[count,1] + segments[count,1]=temp + segments[count,0:2]=numpy.flipud(segments[count,0:2]) + count+=1 + segments=segments[0:count,:] + + return segments +