Changeset 13436


Ignore:
Timestamp:
09/25/12 11:41:00 (13 years ago)
Author:
jschierm
Message:

NEW: Working python bamg et al.

Location:
issm/trunk-jpl/src/m
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/bamggeom.py

    r13270 r13436  
     1import numpy
     2
    13class bamggeom(object):
    24        """
    … …  
    810
    911        def __init__(self,*args):    # {{{
    10                 self.Vertices=[]
    11                 self.Edges=[]
    12                 self.TangentAtEdges=[]
    13                 self.Corners=[]
    14                 self.RequiredVertices=[]
    15                 self.RequiredEdges=[]
    16                 self.CrackedEdges=[]
    17                 self.SubDomains=[]
     12                self.Vertices=numpy.empty((0,3))
     13                self.Edges=numpy.empty((0,3))
     14                self.TangentAtEdges=numpy.empty((0,0))
     15                self.Corners=numpy.empty((0,1))
     16                self.RequiredVertices=numpy.empty((0,1))
     17                self.RequiredEdges=numpy.empty((0,1))
     18                self.CrackedEdges=numpy.empty((0,0))
     19                self.SubDomains=numpy.empty((0,4))
    1820
    1921                if not len(args):

  • issm/trunk-jpl/src/m/classes/bamgmesh.py

    r13270 r13436  
     1import numpy
     2
    13class bamgmesh(object):
    24        """
    … …  
    810
    911        def __init__(self,*args):    # {{{
    10                 self.Vertices=[]
    11                 self.Edges=[]
    12                 self.Triangles=[]
    13                 self.Quadrilaterals=[]
    14                 self.IssmEdges=[]
    15                 self.IssmSegments=[]
    16                 self.VerticesOnGeomVertex=[]
    17                 self.VerticesOnGeomEdge=[]
    18                 self.EdgesOnGeomEdge=[]
    19                 self.SubDomains=[]
    20                 self.SubDomainsFromGeom=[]
    21                 self.ElementConnectivity=[]
    22                 self.NodalConnectivity=[]
    23                 self.NodalElementConnectivity=[]
    24                 self.CrackedVertices=[]
    25                 self.CrackedEdges=[]
     12                self.Vertices=numpy.empty((0,3))
     13                self.Edges=numpy.empty((0,3))
     14                self.Triangles=numpy.empty((0,0))
     15                self.Quadrilaterals=numpy.empty((0,0))
     16                self.IssmEdges=numpy.empty((0,0))
     17                self.IssmSegments=numpy.empty((0,0))
     18                self.VerticesOnGeomVertex=numpy.empty((0,0))
     19                self.VerticesOnGeomEdge=numpy.empty((0,0))
     20                self.EdgesOnGeomEdge=numpy.empty((0,0))
     21                self.SubDomains=numpy.empty((0,4))
     22                self.SubDomainsFromGeom=numpy.empty((0,0))
     23                self.ElementConnectivity=numpy.empty((0,0))
     24                self.NodalConnectivity=numpy.empty((0,0))
     25                self.NodalElementConnectivity=numpy.empty((0,0))
     26                self.CrackedVertices=numpy.empty((0,0))
     27                self.CrackedEdges=numpy.empty((0,0))
    2628
    2729                if not len(args):

  • issm/trunk-jpl/src/m/mesh/bamg.m

    r13353 r13436  
    5555%initialize the structures required as input of Bamg
    5656bamg_options=struct();
    57 bamg_geometry=bamggeom;
    58 bamg_mesh=bamgmesh;
     57bamg_geometry=bamggeom();
     58bamg_mesh=bamgmesh();
    5959
    6060% Bamg Geometry parameters {{{
    … …  
    6363        %Check that file exists
    6464        domainfile=getfieldvalue(options,'domain');
    65         if ~exist(domainfile,'file') error(['bamg error message: file ' domainfile ' not found ']); end
     65        if ~exist(domainfile,'file') error(['bamg error message: file ' domainfile ' not found']); end
    6666        domain=expread(domainfile);
    6767
    … …  
    7979                        flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0);
    8080                        if any(~flags),
    81                                 error('bamg error message: All holes should be stricly inside the principal domain');
     81                                error('bamg error message: All holes should be strictly inside the principal domain');
    8282                        end
    8383                end
    … …  
    8686                nods=domain(i).nods-1; %the domain are closed 1=end;
    8787                bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
    88                 bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1*ones(nods,1)]];
     88                bamg_geometry.Edges   =[bamg_geometry.Edges;    [transpose(count+1:count+nods) transpose([count+2:count+nods count+1])  1.*ones(nods,1)]];
    8989                if i>1, bamg_geometry.SubDomains=[bamg_geometry.SubDomains; 2 count+1 1 1]; end
    9090
    … …  
    105105                for i=1:length(rift),
    106106
    107                         %detect wether all points of the rift are inside the domain
     107                        %detect whether all points of the rift are inside the domain
    108108                        flags=ContourToNodes(rift(i).x,rift(i).y,domain(1),0);
    109109                        if ~flags,
    110                                 error('one Rift has all his points outside of the domain outline'),
     110                                error('one rift has all its points outside of the domain outline'),
    111111
    112112                        elseif any(~flags),
    … …  
    114114                                disp('Rift tip outside of or on the domain has been detected and is being processed...');
    115115
    116                                 %check that only one point is outsie (for now)
     116                                %check that only one point is outside (for now)
    117117                                if sum(~flags)~=1,
    118118                                        error('bamg error message: only one point outside of the domain is supported yet');
    … …  
    160160
    161161                                                        %This point is only in Vertices (number pos).
    162                                                         %OK, no we can add our own rift
     162                                                        %OK, now we can add our own rift
    163163                                                        nods=rift(i).nods-1;
    164164                                                        bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
    … …  
    183183                                                                bamg_geometry.Edges(pos+1:end,:)];
    184184
    185                                                         %OK, no we can add our own rift
     185                                                        %OK, now we can add our own rift
    186186                                                        nods=rift(i).nods-1;
    187187                                                        bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
    … …  
    210210                track=getfieldvalue(options,'tracks');
    211211                if all(ischar(track)),
    212                         A=expread(track); 
     212                        A=expread(track);
    213213                        track=[A.x A.y];
    214214                else
    215215                        track=double(track); %for some reason, it is of class "single"
    216216                end
    217                 if(size(track,2)==2), track=[track 3*ones(size(track,1),1)]; end
     217                if(size(track,2)==2), track=[track 3.*ones(size(track,1),1)]; end
    218218
    219219                %only keep those inside
    … …  
    224224                nods=size(track,1);
    225225                bamg_geometry.Vertices=[bamg_geometry.Vertices; track];
    226                 bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose([count+2:count+nods])  3*ones(nods-1,1)]];
     226                bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose([count+2:count+nods])  3.*ones(nods-1,1)]];
    227227
    228228                %update counter
    … …  
    235235                %recover RequiredVertices
    236236                requiredvertices=double(getfieldvalue(options,'RequiredVertices')); %for some reason, it is of class "single"
    237                 if(size(requiredvertices,2)==2), requiredvertices=[requiredvertices 4*ones(size(requiredvertices,1),1)]; end
     237                if(size(requiredvertices,2)==2), requiredvertices=[requiredvertices 4.*ones(size(requiredvertices,1),1)]; end
    238238       
    239239                %only keep those inside
    … …  
    276276% Bamg Options {{{
    277277bamg_options.Crack=getfieldvalue(options,'Crack',0);
    278 bamg_options.anisomax=getfieldvalue(options,'anisomax',10^30);
    279 bamg_options.coeff=getfieldvalue(options,'coeff',1);
    280 bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
     278bamg_options.anisomax=getfieldvalue(options,'anisomax',10.^30);
     279bamg_options.coeff=getfieldvalue(options,'coeff',1.);
     280bamg_options.cutoff=getfieldvalue(options,'cutoff',10.^-5);
    281281bamg_options.err=getfieldvalue(options,'err',0.01);
    282282bamg_options.errg=getfieldvalue(options,'errg',0.1);
    … …  
    284284bamg_options.gradation=getfieldvalue(options,'gradation',1.5);
    285285bamg_options.Hessiantype=getfieldvalue(options,'Hessiantype',0);
    286 bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
    287 bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
     286bamg_options.hmin=getfieldvalue(options,'hmin',10.^-100);
     287bamg_options.hmax=getfieldvalue(options,'hmax',10.^100);
    288288bamg_options.hminVertices=getfieldvalue(options,'hminVertices',[]);
    289289bamg_options.hmaxVertices=getfieldvalue(options,'hmaxVertices',[]);
    290290bamg_options.hVertices=getfieldvalue(options,'hVertices',[]);
    291291bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1);
    292 bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10);
     292bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10.);
    293293bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6);
    294 bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10);
     294bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10.);
    295295bamg_options.metric=getfieldvalue(options,'metric',[]);
    296296bamg_options.Metrictype=getfieldvalue(options,'Metrictype',0);
    … …  
    298298bamg_options.nbsmooth=getfieldvalue(options,'nbsmooth',3);
    299299bamg_options.omega=getfieldvalue(options,'omega',1.8);
    300 bamg_options.power=getfieldvalue(options,'power',1);
     300bamg_options.power=getfieldvalue(options,'power',1.);
    301301bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1);
    302302bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0);
    … …  
    334334
    335335%Check for orphan
    336 reshape(md.mesh.elements,3*md.mesh.numberofelements,1);
    337 (~ismember(1:md.mesh.numberofvertices,sort(unique(reshape(md.mesh.elements,3*md.mesh.numberofelements,1)))));
    338336if any(~ismember(1:md.mesh.numberofvertices,sort(unique(reshape(md.mesh.elements,3*md.mesh.numberofelements,1)))))
    339337        error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)');

  • issm/trunk-jpl/src/m/mesh/meshconvert.py

    r13279 r13436  
    5252        md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
    5353        md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
    54         md.mesh.vertexonboundary[md.mesh.segments[:,0:1].astype(int)-1]=1
     54        md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
    5555
    5656        return md

  • issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py

    r13170 r13436  
    5353
    5454        else:
    55                 b=numpy.empty_like(a)
    56                 for i,item in enumerate(a.flat):
    57                         b.flat[i]=item in s
     55                if not isinstance(s,numpy.ndarray):
     56                        b=numpy.empty_like(a)
     57                        for i,item in enumerate(a.flat):
     58                                b.flat[i]=item in s
     59                else:
     60                        b=numpy.in1d(a.flat,s.flat).reshape(a.shape)
    5861
    5962        return b
    6063
     64def det(a):
     65        import numpy
     66
     67        if   a.shape==(1,):
     68                return a[0]
     69        elif a.shape==(1,1):
     70                return a[0,0]
     71        elif a.shape==(2,2):
     72                return a[0,0]*a[1,1]-a[0,1]*a[1,0]
     73        else:
     74                raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
     75
Note: See TracChangeset for help on using the changeset viewer.