Index: /issm/trunk-jpl/src/m/classes/bamggeom.py =================================================================== --- /issm/trunk-jpl/src/m/classes/bamggeom.py (revision 13435) +++ /issm/trunk-jpl/src/m/classes/bamggeom.py (revision 13436) @@ -1,2 +1,4 @@ +import numpy + class bamggeom(object): """ @@ -8,12 +10,12 @@ def __init__(self,*args): # {{{ - self.Vertices=[] - self.Edges=[] - self.TangentAtEdges=[] - self.Corners=[] - self.RequiredVertices=[] - self.RequiredEdges=[] - self.CrackedEdges=[] - self.SubDomains=[] + self.Vertices=numpy.empty((0,3)) + self.Edges=numpy.empty((0,3)) + self.TangentAtEdges=numpy.empty((0,0)) + self.Corners=numpy.empty((0,1)) + self.RequiredVertices=numpy.empty((0,1)) + self.RequiredEdges=numpy.empty((0,1)) + self.CrackedEdges=numpy.empty((0,0)) + self.SubDomains=numpy.empty((0,4)) if not len(args): Index: /issm/trunk-jpl/src/m/classes/bamgmesh.py =================================================================== --- /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 13435) +++ /issm/trunk-jpl/src/m/classes/bamgmesh.py (revision 13436) @@ -1,2 +1,4 @@ +import numpy + class bamgmesh(object): """ @@ -8,20 +10,20 @@ def __init__(self,*args): # {{{ - self.Vertices=[] - self.Edges=[] - self.Triangles=[] - self.Quadrilaterals=[] - self.IssmEdges=[] - self.IssmSegments=[] - self.VerticesOnGeomVertex=[] - self.VerticesOnGeomEdge=[] - self.EdgesOnGeomEdge=[] - self.SubDomains=[] - self.SubDomainsFromGeom=[] - self.ElementConnectivity=[] - self.NodalConnectivity=[] - self.NodalElementConnectivity=[] - self.CrackedVertices=[] - self.CrackedEdges=[] + self.Vertices=numpy.empty((0,3)) + self.Edges=numpy.empty((0,3)) + self.Triangles=numpy.empty((0,0)) + self.Quadrilaterals=numpy.empty((0,0)) + self.IssmEdges=numpy.empty((0,0)) + self.IssmSegments=numpy.empty((0,0)) + self.VerticesOnGeomVertex=numpy.empty((0,0)) + self.VerticesOnGeomEdge=numpy.empty((0,0)) + self.EdgesOnGeomEdge=numpy.empty((0,0)) + self.SubDomains=numpy.empty((0,4)) + self.SubDomainsFromGeom=numpy.empty((0,0)) + self.ElementConnectivity=numpy.empty((0,0)) + self.NodalConnectivity=numpy.empty((0,0)) + self.NodalElementConnectivity=numpy.empty((0,0)) + self.CrackedVertices=numpy.empty((0,0)) + self.CrackedEdges=numpy.empty((0,0)) if not len(args): Index: /issm/trunk-jpl/src/m/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.m (revision 13435) +++ /issm/trunk-jpl/src/m/mesh/bamg.m (revision 13436) @@ -55,6 +55,6 @@ %initialize the structures required as input of Bamg bamg_options=struct(); -bamg_geometry=bamggeom; -bamg_mesh=bamgmesh; +bamg_geometry=bamggeom(); +bamg_mesh=bamgmesh(); % Bamg Geometry parameters {{{ @@ -63,5 +63,5 @@ %Check that file exists domainfile=getfieldvalue(options,'domain'); - if ~exist(domainfile,'file') error(['bamg error message: file ' domainfile ' not found ']); end + if ~exist(domainfile,'file') error(['bamg error message: file ' domainfile ' not found']); end domain=expread(domainfile); @@ -79,5 +79,5 @@ flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0); if any(~flags), - error('bamg error message: All holes should be stricly inside the principal domain'); + error('bamg error message: All holes should be strictly inside the principal domain'); end end @@ -86,5 +86,5 @@ nods=domain(i).nods-1; %the domain are closed 1=end; bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]]; - bamg_geometry.Edges =[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose([count+2:count+nods count+1]) 1*ones(nods,1)]]; + bamg_geometry.Edges =[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose([count+2:count+nods count+1]) 1.*ones(nods,1)]]; if i>1, bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 1]; end @@ -105,8 +105,8 @@ for i=1:length(rift), - %detect wether all points of the rift are inside the domain + %detect whether all points of the rift are inside the domain flags=ContourToNodes(rift(i).x,rift(i).y,domain(1),0); if ~flags, - error('one Rift has all his points outside of the domain outline'), + error('one rift has all its points outside of the domain outline'), elseif any(~flags), @@ -114,5 +114,5 @@ disp('Rift tip outside of or on the domain has been detected and is being processed...'); - %check that only one point is outsie (for now) + %check that only one point is outside (for now) if sum(~flags)~=1, error('bamg error message: only one point outside of the domain is supported yet'); @@ -160,5 +160,5 @@ %This point is only in Vertices (number pos). - %OK, no we can add our own rift + %OK, now we can add our own rift nods=rift(i).nods-1; bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; @@ -183,5 +183,5 @@ bamg_geometry.Edges(pos+1:end,:)]; - %OK, no we can add our own rift + %OK, now we can add our own rift nods=rift(i).nods-1; bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; @@ -210,10 +210,10 @@ track=getfieldvalue(options,'tracks'); if all(ischar(track)), - A=expread(track); + A=expread(track); track=[A.x A.y]; else track=double(track); %for some reason, it is of class "single" end - if(size(track,2)==2), track=[track 3*ones(size(track,1),1)]; end + if(size(track,2)==2), track=[track 3.*ones(size(track,1),1)]; end %only keep those inside @@ -224,5 +224,5 @@ nods=size(track,1); bamg_geometry.Vertices=[bamg_geometry.Vertices; track]; - bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) 3*ones(nods-1,1)]]; + bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) 3.*ones(nods-1,1)]]; %update counter @@ -235,5 +235,5 @@ %recover RequiredVertices requiredvertices=double(getfieldvalue(options,'RequiredVertices')); %for some reason, it is of class "single" - if(size(requiredvertices,2)==2), requiredvertices=[requiredvertices 4*ones(size(requiredvertices,1),1)]; end + if(size(requiredvertices,2)==2), requiredvertices=[requiredvertices 4.*ones(size(requiredvertices,1),1)]; end %only keep those inside @@ -276,7 +276,7 @@ % Bamg Options {{{ bamg_options.Crack=getfieldvalue(options,'Crack',0); -bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30); -bamg_options.coeff=getfieldvalue(options,'coeff',1); -bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5); +bamg_options.anisomax=getfieldvalue(options,'anisomax',10.^30); +bamg_options.coeff=getfieldvalue(options,'coeff',1.); +bamg_options.cutoff=getfieldvalue(options,'cutoff',10.^-5); bamg_options.err=getfieldvalue(options,'err',0.01); bamg_options.errg=getfieldvalue(options,'errg',0.1); @@ -284,13 +284,13 @@ bamg_options.gradation=getfieldvalue(options,'gradation',1.5); bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0); -bamg_options.hmin=getfieldvalue(options,'hmin',10^-100); -bamg_options.hmax=getfieldvalue(options,'hmax',10^100); +bamg_options.hmin=getfieldvalue(options,'hmin',10.^-100); +bamg_options.hmax=getfieldvalue(options,'hmax',10.^100); bamg_options.hminVertices=getfieldvalue(options,'hminVertices',[]); bamg_options.hmaxVertices=getfieldvalue(options,'hmaxVertices',[]); bamg_options.hVertices=getfieldvalue(options,'hVertices',[]); bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1); -bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10); +bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10.); bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6); -bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10); +bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10.); bamg_options.metric=getfieldvalue(options,'metric',[]); bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0); @@ -298,5 +298,5 @@ bamg_options.nbsmooth=getfieldvalue(options,'nbsmooth',3); bamg_options.omega=getfieldvalue(options,'omega',1.8); -bamg_options.power=getfieldvalue(options,'power',1); +bamg_options.power=getfieldvalue(options,'power',1.); bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1); bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0); @@ -334,6 +334,4 @@ %Check for orphan -reshape(md.mesh.elements,3*md.mesh.numberofelements,1); -(~ismember(1:md.mesh.numberofvertices,sort(unique(reshape(md.mesh.elements,3*md.mesh.numberofelements,1))))); if any(~ismember(1:md.mesh.numberofvertices,sort(unique(reshape(md.mesh.elements,3*md.mesh.numberofelements,1))))) error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)'); Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13436) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13436) @@ -0,0 +1,486 @@ +import os.path +import numpy +from collections import OrderedDict +from pairoptions import * +from bamggeom import * +from bamgmesh import * +from expread import * +from expwrite import * +from MatlabFuncs import * +from BamgMesher import * + +def bamg(md,*args): + """ + BAMG - mesh generation + + Available options (for more details see ISSM website http://issm.jpl.nasa.gov/): + + - domain : followed by an ARGUS file that prescribes the domain outline + - hmin : minimum edge length (default is 10^-100) + - hmax : maximum edge length (default is 10^100) + - hVertices : imposed edge length for each vertex (geometry or mesh) + - hminVertices : minimum edge length for each vertex (mesh) + - hmaxVertices : maximum edge length for each vertex (mesh) + + - anisomax : maximum ratio between the smallest and largest edges (default is 10^30) + - coeff : coefficient applied to the metric (2-> twice as many elements, default is 1) + - cutoff : scalar used to compute the metric when metric type 2 or 3 are applied + - err : error used to generate the metric from a field + - errg : geometric error (default is 0.1) + - field : field of the model that will be used to compute the metric + to apply several fields, use one column per field + - gradation : maximum ratio between two adjacent edges + - Hessiantype : 0 -> use double P2 projection (default) + 1 -> use Green formula + - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) + - MaxCornerAngle : maximum angle of corners in degree (default is 10) + - maxnbv : maximum number of vertices used to allocate memory (default is 10^6) + - maxsubdiv : maximum subdivision of exisiting elements (default is 10) + - metric : matrix (numberofnodes x 3) used as a metric + - Metrictype : 1 -> absolute error c/(err coeff^2) * Abs(H) (default) + 2 -> relative error c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s)) + 3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin) + - nbjacoby : correction used by Hessiantype=1 (default is 1) + - nbsmooth : number of metric smoothing procedure (default is 3) + - omega : relaxation parameter of the smoothing procedure (default is 1.8) + - power : power applied to the metric (default is 1) + - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) + - geometricalmetric : take the geometry into account to generate the metric (default is 0) + - verbose : level of verbosity (default is 1) + + - rifts : followed by an ARGUS file that prescribes the rifts + - toltip : tolerance to move tip on an existing point of the domain outline + - tracks : followed by an ARGUS file that prescribes the tracks that the mesh will stick to + - RequiredVertices : mesh vertices that are required. [x,y,ref]; ref is optional + - tol : if the distance between 2 points of the domain outline is less than tol, they + will be merged + + Examples: + md=bamg(md,'domain','DomainOutline.exp','hmax',3000); + md=bamg(md,'field',[md.inversion.vel_obs md.geometry.thickness],'hmax',20000,'hmin',1000); + md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1); + """ + + #process options + options=pairoptions(*args) +# options=deleteduplicates(options,1); + + #initialize the structures required as input of Bamg + bamg_options=OrderedDict() + bamg_geometry=bamggeom() + bamg_mesh=bamgmesh() + + # Bamg Geometry parameters {{{ + if options.exist('domain'): + + #Check that file exists + domainfile=options.getfieldvalue('domain') + if not os.path.exists(domainfile): + raise IOError("bamg error message: file '%s' not found" % domainfile) + domain=expread(domainfile) + + #Build geometry + count=0 + for i,domaini in enumerate(domain): + + #Check that the domain is closed + if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]): + raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed") + + #Checks that all holes are INSIDE the principle domain outline + if i: + flags=ContourToNodes(domaini['x'],domaini['y'],domain[0],0) + if numpy.any(numpy.logical_not(flags)): + raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") + + #Add all points to bamg_geometry + nods=domaini['nods']-1 #the domain are closed 0=end + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((domaini['x'][0:nods].reshape(-1,1),domaini['y'][0:nods].reshape(-1,1),numpy.ones((nods,1)))))) + bamg_geometry.Edges =numpy.vstack((bamg_geometry.Edges, numpy.hstack((numpy.arange(count+1,count+nods+1).reshape(-1,1),numpy.hstack((numpy.arange(count+2,count+nods+1),count+1)).reshape(-1,1),1.*numpy.ones((nods,1)))))) + if i: + bamg_geometry.SubDomains=numpy.vstack((bamg_geometry.SubDomains,[2,count+1,1,1])) + + #update counter + count=+nods + + #take care of rifts + if options.exist('rifts'): + + #Check that file exists + riftfile=options.getfieldvalue('rifts') + if not os.path.exists(riftfile): + raise IOError("bamg error message: file '%s' not found" % riftfile) + rift=expread(riftfile) + + for i,rifti in enumerate(rift): + + #detect whether all points of the rift are inside the domain + flags=ContourToNodes(rifti['x'],rifti['y'],domain[0],0) + if numpy.all(numpy.logical_not(flags)): + raise RuntimeError("one rift has all its points outside of the domain outline") + + elif numpy.any(numpy.logical_not(flags)): + raise RuntimeError("bamg.m for rifts is not complete.") + #We LOTS of work to do + 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) + if sum(~flags)~=1, + error('bamg error message: only one point outside of the domain is supported yet'); + end + + %Move tip outside to the first position + if flags(1)==0, + %OK, first point is outside (do nothing), + elseif (flags(end)==0), + rift(i).x=flipud(rift(i).x); + rift(i).y=flipud(rift(i).y); + else + error('bamg error message: only a rift tip can be outside of the domain'); + end + + %Get cordinate of intersection point + x1=rift(i).x(1); y1=rift(i).y(1); + x2=rift(i).x(2); y2=rift(i).y(2); + for j=1:length(domain[0]['x'])-1; + if SegIntersect([x1 y1; x2 y2],[domain[0]['x'](j) domain[0]['y'](j); domain[0]['x'](j+1) domain[0]['y'](j+1)]), + + %Get position of the two nodes of the edge in domain + i1=j; + i2=mod(j+1,domain[0]['nods']); + + %rift is crossing edge [i1 i2] of the domain + %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) + x3=domain[0]['x'](i1); y3=domain[0]['y'](i1); + x4=domain[0]['x'](i2); y4=domain[0]['y'](i2); + x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); + y=det([det([x1 y1; x2 y2]) y1-y2;det([x3 y3; x4 y4]) y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); + + segdis= sqrt((x4-x3)**2+(y4-y3)**2); + tipdis=[sqrt((x-x3)**2+(y-y3)**2) sqrt((x-x4)**2+(y-y4)**2)]; + + if (min(tipdis)/segdis) < getfieldvalue(options,'toltip',0), + disp('moving tip-domain intersection point'); + + %Get position of the closest point + if tipdis(1)>tipdis(2), + pos=i2; + else + pos=i1; + end + + %This point is only in Vertices (number pos). + %OK, now we can add our own rift + nods=rift(i).nods-1; + bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; + bamg_geometry.Edges=[bamg_geometry.Edges;... + pos count+1 (1+i);... + [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; + count=count+nods; + + break; + + else + %Add intersection point to Vertices + bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1]; + count=count+1; + + %Decompose the crossing edge in 2 subedges + pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2); + if isempty(pos) error('bamg error message: a problem occured...'); end + bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);... + bamg_geometry.Edges(pos,1) count bamg_geometry.Edges(pos,3);... + count bamg_geometry.Edges(pos,2) bamg_geometry.Edges(pos,3);... + bamg_geometry.Edges(pos+1:end,:)]; + + %OK, now we can add our own rift + nods=rift(i).nods-1; + bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]]; + bamg_geometry.Edges=[bamg_geometry.Edges;... + count count+1 2 ;... + [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; + count=count+nods; + + break; + end + end + end + """ + else: + nods=rifti['nods']-1 + bamg_geometry.Vertices=numpy.vstack(bamg_geometry.Vertices, numpy.hstack(rifti['x'][:],rifti['y'][:],numpy.ones((nods+1,1)))) + bamg_geometry.Edges =numpy.vstack(bamg_geometry.Edges, numpy.hstack(numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),i*numpy.ones((nods,1)))) + count=+nods+1 + + #Deal with tracks + if options.exist('tracks'): + + #read tracks + track=options.getfieldvalue('tracks') + if all(isinstance(track,(str,unicode))): + A=expread(track) + track=numpy.hstack((A.x.reshape(-1,1),A.y.reshape(-1,1))) + else: + track=float(track) #for some reason, it is of class "single" + if numpy.size(track,axis=1)==2: + track=numpy.hstack((track,3.*numpy.ones((size(track,axis=0),1)))) + + #only keep those inside + flags=ContourToNodes(track[:,0],track[:,1],domainfile,0) + track=track[numpy.nonzero(flags),:] + + #Add all points to bamg_geometry + nods=numpy.size(track,axis=0) + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,track)) + bamg_geometry.Edges =numpy.vstack((bamg_geometry.Edges,numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),3.*numpy.ones((nods-1,1)))))) + + #update counter + count+=nods + + #Deal with vertices that need to be kept by mesher + if options.exist('RequiredVertices'): + + #recover RequiredVertices + requiredvertices=float(options.getfieldvalue('RequiredVertices')) #for some reason, it is of class "single" + if numpy.size(requiredvertices,axis=1)==2: + requiredvertices=numpy.hstack((requiredvertices,4.*numpy.ones((numpy.size(requiredvertices,axis=0),1)))) + + #only keep those inside + flags=ContourToNodes(requiredvertices[:,0],domain[0],0) + requiredvertices=requiredvertices[numpy.nonzero(flags),:] + + #Add all points to bamg_geometry + nods=numpy.size(requiredvertices,axis=0) + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,requiredvertices)) + + #update counter + count+=nods + + #process geom + #bamg_geometry=processgeometry(bamg_geometry,options.getfieldvalue('tol',float(nan)),domain[0]) + + elif isinstance(md.private.bamg,dict) and 'geometry' in md.private.bamg: + bamg_geometry=bamggeom(md.private.bamg['geometry']) + else: + #do nothing... + pass + #}}} + # Bamg Mesh parameters {{{ + if options.exist('domain') and md.mesh.numberofvertices and md.mesh.dimension==2: + + if isinstance(md.private.bamg,dict) and 'mesh' in md.private.bamg: + bamg_mesh=bamgmesh(md.private.bamg['mesh']) + else: + bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1)))) + bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1)))) + + if isinstance(md.rifts.riftstruct,dict): + raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg") + #}}} + # Bamg Options {{{ + bamg_options['Crack']=options.getfieldvalue('Crack',0) + bamg_options['anisomax']=options.getfieldvalue('anisomax',10.**30) + bamg_options['coeff']=options.getfieldvalue('coeff',1.) + bamg_options['cutoff']=options.getfieldvalue('cutoff',10.**-5) + bamg_options['err']=options.getfieldvalue('err',numpy.array([[0.01]])) + bamg_options['errg']=options.getfieldvalue('errg',0.1) + bamg_options['field']=options.getfieldvalue('field',numpy.empty((0,1))) + bamg_options['gradation']=options.getfieldvalue('gradation',1.5) + bamg_options['Hessiantype']=options.getfieldvalue('Hessiantype',0) + bamg_options['hmin']=options.getfieldvalue('hmin',10.**-100) + bamg_options['hmax']=options.getfieldvalue('hmax',10.**100) + bamg_options['hminVertices']=options.getfieldvalue('hminVertices',numpy.empty((0,1))) + bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',numpy.empty((0,1))) + bamg_options['hVertices']=options.getfieldvalue('hVertices',numpy.empty((0,1))) + bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1) + bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.) + bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6) + bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.) + bamg_options['metric']=options.getfieldvalue('metric',numpy.empty((0,1))) + bamg_options['Metrictype']=options.getfieldvalue('Metrictype',0) + bamg_options['nbjacobi']=options.getfieldvalue('nbjacobi',1) + bamg_options['nbsmooth']=options.getfieldvalue('nbsmooth',3) + bamg_options['omega']=options.getfieldvalue('omega',1.8) + bamg_options['power']=options.getfieldvalue('power',1.) + bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1) + bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0) + bamg_options['verbose']=options.getfieldvalue('verbose',1) + #}}} + + #call Bamg + bamgmesh_out,bamggeom_out=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options) + + # plug results onto model + md.private.bamg=OrderedDict() + md.private.bamg['mesh']=bamgmesh(bamgmesh_out) + md.private.bamg['geometry']=bamggeom(bamggeom_out) + md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() + md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy() + md.mesh.edges=bamgmesh_out['IssmEdges'].copy() + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy() + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy() + + #Fill in rest of fields: + md.mesh.dimension=2 + md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=numpy.size(md.mesh.x); + md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) + md.mesh.z=numpy.zeros(md.mesh.numberofvertices) + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) + md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) + md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity + md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0 + + #Check for orphan + if numpy.any(numpy.logical_not(numpy.in1d(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))): + raise RuntimeError("Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)") + + return md + +def processgeometry(geom,tol,outline): # {{{ + + raise RuntimeError("bamg.m/processgeometry is not complete.") + #Deal with edges + print("Checking Edge crossing...") + i=0 + while (ii) + geom.Edges[posedges]=geom.Edges[posedges]-1 + + #update counter + i-=1 + + if num: + print("WARNING: %d points outside the domain outline have been removed" % num) + + """ + %Check point spacing + if ~isnan(tol), + disp('Checking point spacing...'); + i=0; + while (ij); + geom.Edges(posedges)=geom.Edges(posedges)-1; + + %update counter + j=j-1; + + end + end + end + end + %remove empty edges + geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[]; + """ + return geom +# }}} + Index: /issm/trunk-jpl/src/m/mesh/meshconvert.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 13435) +++ /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 13436) @@ -52,5 +52,5 @@ md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary[md.mesh.segments[:,0:1].astype(int)-1]=1 + md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 return md Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 13435) +++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 13436) @@ -53,8 +53,23 @@ else: - b=numpy.empty_like(a) - for i,item in enumerate(a.flat): - b.flat[i]=item in s + if not isinstance(s,numpy.ndarray): + b=numpy.empty_like(a) + for i,item in enumerate(a.flat): + b.flat[i]=item in s + else: + b=numpy.in1d(a.flat,s.flat).reshape(a.shape) return b +def det(a): + import numpy + + if a.shape==(1,): + return a[0] + elif a.shape==(1,1): + return a[0,0] + elif a.shape==(2,2): + return a[0,0]*a[1,1]-a[0,1]*a[1,0] + else: + raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape)) +