Changeset 28028


Ignore:
Timestamp:
12/19/23 13:19:39 (15 months ago)
Author:
jdquinn
Message:

CHG: MATLAB -> Python

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

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/m/exp/contourlevelzero.py

    r27950 r28028  
    1 import os.path
    2 import numpy as np
    3 from collections import OrderedDict
    4 
    5 
    61def contourlevelzero(md,mask,level):
    7     """CONTOURLEVELZERO - figure out the zero level (or offset thereof,
    8     specified by the level value of a vectorial mask, and vectorialize it into
    9     an exp or shp compatible structure.
    10 
    11     Usage:
    12         contours=contourlevelzero(md,mask,level)
    13 
    14         See also: PLOT_CONTOUR
    15     """
    16 
    17     # Process data
    18     if md.mesh.dimension() == 3:
    19         x = md.mesh.x2d
    20         y = md.mesh.y2d
    21         z = md.mesh.z
    22         index = md.mesh.elements2d - 1
    23     else:
    24         x = md.mesh.x
    25         y = md.mesh.y
    26         index = md.mesh.elements - 1
    27         z = np.zeros((md.mesh.numberofvertices, 1))
    28 
    29     if len(mask) == 0:
    30         raise OSError("mask provided is empty")
    31 
    32     if md.mesh.dimension() == 3:
    33         if len(mask) != md.mesh.numberofvertices2d:
    34             raise OSError("mask provided should be specified at the vertices of the mesh")
    35     else:
    36         if len(mask) != md.mesh.numberofvertices:
    37             raise OSError("mask provided should be specified at the vertices of the mesh")
    38 
    39     # Initialization of some variables
    40     numberofelements = np.size(index, 0)
    41     elementslist = np.c_[0:numberofelements]
    42     c = []
    43     h = []
    44 
    45     # Get unique edges in mesh
    46     # 1: list of edges
    47     edges = np.vstack((np.vstack((index[:, (0, 1)], index[:, (1, 2)])), index[:, (2, 0)]))
    48     # 2: find unique edges
    49     [edges, J] = np.unique(np.sort(edges, 1), axis=0, return_inverse=True)
    50     # 3: unique edge numbers
    51     vec = J
    52     # 4: unique edges numbers in each triangle (2 triangles sharing the same
    53     # edge will have the same edge number)
    54     edges_tria = np.hstack((np.hstack((vec[elementslist], vec[elementslist + numberofelements])), vec[elementslist + 2 * numberofelements]))
    55 
    56     # Segments [nodes1 nodes2]
    57     Seg1 = index[:, (0, 1)]
    58     Seg2 = index[:, (1, 2)]
    59     Seg3 = index[:, (2, 0)]
    60 
    61     # Segment numbers [1;4;6;...]
    62     Seg1_num = edges_tria[:, 0]
    63     Seg2_num = edges_tria[:, 1]
    64     Seg3_num = edges_tria[:, 2]
    65 
    66     #value of data on each tips of the segments
    67     Data1 = mask[Seg1]
    68     Data2 = mask[Seg2]
    69     Data3 = mask[Seg3]
    70 
    71     # Get the ranges for each segment
    72     Range1 = np.sort(Data1, 1)
    73     Range2 = np.sort(Data2, 1)
    74     Range3 = np.sort(Data3, 1)
    75 
    76     # Find the segments that contain this value
    77     pos1 = (Range1[:, 0] < level) & (Range1[:, 1] >= level)
    78     pos2 = (Range2[:, 0] < level) & (Range2[:, 1] >= level)
    79     pos3 = (Range3[:, 0] < level) & (Range3[:, 1] >= level)
    80 
    81     # Get elements
    82     poselem12 = (pos1) & (pos2)
    83     poselem13 = (pos1) & (pos3)
    84     poselem23 = (pos2) & (pos3)
    85     poselem = np.where((poselem12) | (poselem13) | (poselem23))
    86     poselem = poselem[0]
    87     numelems = len(poselem)
    88 
    89     # If no element has been flagged, skip to the next level
    90     if numelems == 0:
    91         raise Exception('contourlevelzero warning message: no elements found with corresponding level value in mask')
    92         contours = []
    93         return contours
    94 
    95     # Go through the elements and build the coordinates for each segment (1 by element)
    96     x1 = np.zeros((numelems, 1))
    97     x2 = np.zeros((numelems, 1))
    98     y1 = np.zeros((numelems, 1))
    99     y2 = np.zeros((numelems, 1))
    100     z1 = np.zeros((numelems, 1))
    101     z2 = np.zeros((numelems, 1))
    102 
    103     edge_l = np.zeros((numelems, 2))
    104 
    105     for j in range(0, numelems):
    106         with np.errstate(divide='ignore', invalid='ignore'):
    107             weight1 = np.divide(level - Data1[poselem[j], 0],Data1[poselem[j], 1] - Data1[poselem[j], 0])
    108             weight2 = np.divide(level - Data2[poselem[j], 0],Data2[poselem[j], 1] - Data2[poselem[j], 0])
    109             weight3 = np.divide(level - Data3[poselem[j], 0],Data3[poselem[j], 1] - Data3[poselem[j], 0])
    110 
    111         if poselem12[poselem[j]] == True:
    112             x1[j] = x[Seg1[poselem[j], 0]] + weight1 * [x[Seg1[poselem[j], 1]] - x[Seg1[poselem[j], 0]]]
    113             x2[j] = x[Seg2[poselem[j], 0]] + weight2 * [x[Seg2[poselem[j], 1]] - x[Seg2[poselem[j], 0]]]
    114             y1[j] = y[Seg1[poselem[j], 0]] + weight1 * [y[Seg1[poselem[j], 1]] - y[Seg1[poselem[j], 0]]]
    115             y2[j] = y[Seg2[poselem[j], 0]] + weight2 * [y[Seg2[poselem[j], 1]] - y[Seg2[poselem[j], 0]]]
    116             z1[j] = z[Seg1[poselem[j], 0]] + weight1 * [z[Seg1[poselem[j], 1]] - z[Seg1[poselem[j], 0]]]
    117             z2[j] = z[Seg2[poselem[j], 0]] + weight2 * [z[Seg2[poselem[j], 1]] - z[Seg2[poselem[j], 0]]]
    118 
    119             edge_l[j, 0] = Seg1_num[poselem[j]]
    120             edge_l[j, 1] = Seg2_num[poselem[j]]
    121         elif poselem13[poselem[j]] == True:
    122             x1[j] = x[Seg1[poselem[j], 0]] + weight1 * [x[Seg1[poselem[j], 1]] - x[Seg1[poselem[j], 0]]]
    123             x2[j] = x[Seg3[poselem[j], 0]] + weight3 * [x[Seg3[poselem[j], 1]] - x[Seg3[poselem[j], 0]]]
    124             y1[j] = y[Seg1[poselem[j], 0]] + weight1 * [y[Seg1[poselem[j], 1]] - y[Seg1[poselem[j], 0]]]
    125             y2[j] = y[Seg3[poselem[j], 0]] + weight3 * [y[Seg3[poselem[j], 1]] - y[Seg3[poselem[j], 0]]]
    126             z1[j] = z[Seg1[poselem[j], 0]] + weight1 * [z[Seg1[poselem[j], 1]] - z[Seg1[poselem[j], 0]]]
    127             z2[j] = z[Seg3[poselem[j], 0]] + weight3 * [z[Seg3[poselem[j], 1]] - z[Seg3[poselem[j], 0]]]
    128 
    129             edge_l[j, 0] = Seg1_num[poselem[j]]
    130             edge_l[j, 1] = Seg3_num[poselem[j]]
    131         elif poselem23[poselem[j]] == True:
    132             x1[j] = x[Seg2[poselem[j], 0]] + weight2 * [x[Seg2[poselem[j], 1]] - x[Seg2[poselem[j], 0]]]
    133             x2[j] = x[Seg3[poselem[j], 0]] + weight3 * [x[Seg3[poselem[j], 1]] - x[Seg3[poselem[j], 0]]]
    134             y1[j] = y[Seg2[poselem[j], 0]] + weight2 * [y[Seg2[poselem[j], 1]] - y[Seg2[poselem[j], 0]]]
    135             y2[j] = y[Seg3[poselem[j], 0]] + weight3 * [y[Seg3[poselem[j], 1]] - y[Seg3[poselem[j], 0]]]
    136             z1[j] = z[Seg2[poselem[j], 0]] + weight2 * [z[Seg2[poselem[j], 1]] - z[Seg2[poselem[j], 0]]]
    137             z2[j] = z[Seg3[poselem[j], 0]] + weight3 * [z[Seg3[poselem[j], 1]] - z[Seg3[poselem[j], 0]]]
    138 
    139             edge_l[j, 0] = Seg2_num[poselem[j]]
    140             edge_l[j, 1] = Seg3_num[poselem[j]]
    141         # else:
    142         # Should never get here
    143 
    144     # Now that we have the segments, we must try to connect them...
    145 
    146     # Loop over the subcontours
    147     contours = []
    148 
    149     while len(edge_l) > 0:
    150         # Take the right edge of the second segment and connect it to the next segments if any
    151         e1 = edge_l[0, 0]
    152         e2 = edge_l[0, 1]
    153         xc = np.vstack((x1[0], x2[0]))
    154         yc = np.vstack((y1[0], y2[0]))
    155         zc = np.vstack((z1[0], z2[0]))
    156         # Erase the lines corresponding to this edge
    157         edge_l = np.delete(edge_l, 0, axis=0)
    158         x1 = np.delete(x1, 0, axis=0)
    159         x2 = np.delete(x2, 0, axis=0)
    160         y1 = np.delete(y1, 0, axis=0)
    161         y2 = np.delete(y2, 0, axis=0)
    162         z1 = np.delete(z1, 0, axis=0)
    163         z2 = np.delete(z2,0,axis=0)
    164         pos1 = np.where(edge_l == e1)
    165        
    166         while len(pos1[0]) > 0:
    167             if np.all(pos1[1] == 0):
    168                 xc = np.vstack((x2[pos1[0]], xc))
    169                 yc = np.vstack((y2[pos1[0]], yc))
    170                 zc = np.vstack((z2[pos1[0]], zc))
    171                 # Next edge:
    172                 e1 = edge_l[pos1[0], 1]
    173             else:
    174                 xc = np.vstack((x1[pos1[0]], xc))
    175                 yc = np.vstack((y1[pos1[0]], yc))
    176                 zc = np.vstack((z1[pos1[0]], zc))
    177                 # Next edge:
    178                 e1 = edge_l[pos1[0], 0]
    179 
    180             # Erase the lines of this
    181             edge_l = np.delete(edge_l, pos1[0], axis=0)
    182             x1 = np.delete(x1, pos1[0], axis=0)
    183             x2 = np.delete(x2, pos1[0], axis=0)
    184             y1 = np.delete(y1, pos1[0], axis=0)
    185             y2 = np.delete(y2, pos1[0], axis=0)
    186             z1 = np.delete(z1, pos1[0], axis=0)
    187             z2 = np.delete(z2, pos1[0], axis=0)
    188             # Next connection
    189             pos1 = np.where(edge_l == e1)
    190 
    191         # Same thing the other way (to the right)
    192         pos2 = np.where(edge_l == e2)
    193 
    194         while len(pos2[0]) > 0:
    195             if np.all(pos2[1] == 0):
    196                 xc = np.vstack((xc, x2[pos2[0]]))
    197                 yc = np.vstack((yc, y2[pos2[0]]))
    198                 zc = np.vstack((zc, z2[pos2[0]]))
    199                 # Next edge:
    200                 e2 = edge_l[pos2[0], 1]
    201             else:
    202                 xc = np.vstack((xc, x1[pos2[0]]))
    203                 yc = np.vstack((yc, y1[pos2[0]]))
    204                 zc = np.vstack((zc, z1[pos2[0]]))
    205                 # Next edge:
    206                 e2 = edge_l[pos2[0], 0]
    207 
    208             # Erase the lines of this
    209             edge_l = np.delete(edge_l, pos2[0], axis=0)
    210             x1 = np.delete(x1, pos2[0], axis=0)
    211             x2 = np.delete(x2, pos2[0], axis=0)
    212             y1 = np.delete(y1, pos2[0], axis=0)
    213             y2 = np.delete(y2, pos2[0], axis=0)
    214             z1 = np.delete(z1, pos2[0], axis=0)
    215             z2 = np.delete(z2, pos2[0], axis=0)
    216             # Next connection
    217             pos2 = np.where(edge_l == e2)
    218 
    219         # Save xc, yc contour:
    220         newcontour = OrderedDict()
    221         newcontour['nods'] = np.size(xc)
    222         newcontour['density'] = 1
    223         newcontour['closed'] = 0
    224         newcontour['x'] = np.ma.filled(xc.astype(float), np.nan)
    225         newcontour['y'] = np.ma.filled(yc.astype(float), np.nan)
    226         newcontour['z'] = np.ma.filled(zc.astype(float), np.nan)
    227         newcontour['name'] = ''
    228         contours.append(newcontour)
    229 
    230     return contours
     2    raise Exception('this function has been renamed: A = isoline(md, mask)')
  • issm/trunk-jpl/src/m/exp/isoline.m

    r27961 r28028  
    166166                edge_l(j,2)=Seg3_num(poselem(j));
    167167        else
    168                 %it shoud not go here
     168                %should never get here
    169169        end
    170170end
Note: See TracChangeset for help on using the changeset viewer.