Changeset 5372


Ignore:
Timestamp:
May 28, 2008, 2:47:30 PM (16 years ago)
Author:
jakeman
Message:

updated mux2 boundary interpolation methods

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r5361 r5372  
    1515
    1616from anuga.utilities.numerical_tools import ensure_numeric
    17 from Numeric import arange, choose, zeros, Float, array
     17from Numeric import arange, choose, zeros, Float, array, allclose, take, compress
    1818   
    1919from anuga.geospatial_data.geospatial_data import ensure_absolute
     
    3333                  time_thinning=1,
    3434                  verbose=False,
    35                   use_cache=False):
     35                  use_cache=False,
     36                  boundary_polygon=None):
    3637    """Read time history of spatial data from NetCDF file and return
    3738    a callable object.
     
    114115              'domain_starttime': domain_starttime,
    115116              'time_thinning': time_thinning,                   
    116               'verbose': verbose}
     117              'verbose': verbose,
     118              'boundary_polygon':boundary_polygon}
    117119
    118120
     
    130132                             dependencies=[filename],
    131133                             compression=False,                 
    132                              verbose=verbose)
     134                             verbose=verbose,
     135                             boundary_polygon=boundary_polygon)
    133136
    134137    else:
     
    166169                   domain_starttime=None,
    167170                   time_thinning=1,
    168                    verbose=False):
     171                   verbose=False,
     172                   boundary_polygon=None):
    169173    """Internal function
    170174   
     
    193197                                        domain_starttime,
    194198                                        time_thinning=time_thinning,
    195                                         verbose=verbose)
     199                                        verbose=verbose,
     200                                        boundary_polygon=boundary_polygon)
    196201    else:
    197202        raise 'Must be a NetCDF File'
     
    204209                             domain_starttime=None,                           
    205210                             time_thinning=1,                             
    206                              verbose=False):
     211                             verbose=False,
     212                             boundary_polygon=None):
    207213    """Read time history of spatial data from NetCDF sww file and
    208214    return a callable object f(t,x,y)
     
    282288        raise msg
    283289
     290    if filename[-3:] == 'sts' and boundary_polygon is None:
     291        #What if mux file only contains one point
     292        msg = 'Files of type sts require boundary polygon'       
     293        raise msg
     294
    284295    # Get first timestep
    285296    try:
     
    308319        y = reshape(y, (len(y),1))
    309320        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
     321
     322        if boundary_polygon is not None:
     323            #removes sts points that do not lie on boundary
     324            boundary_polygon=ensure_numeric(boundary_polygon)
     325            boundary_polygon[:,0] -= xllcorner
     326            boundary_polygon[:,1] -= yllcorner
     327            temp=[]
     328            boundary_id=[]
     329            gauge_id=[]
     330            for i in range(len(boundary_polygon)):
     331                for j in range(len(x)):
     332                    if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
     333                        #FIX ME:
     334                        #currently gauges lat and long is stored as float and
     335                        #then cast to double. This cuases slight repositioning
     336                        #of vertex_coordinates.
     337                        temp.append(boundary_polygon[i])
     338                        gauge_id.append(j)
     339                        boundary_id.append(i)
     340                        break
     341            gauge_neighbour_id=[]
     342            for i in range(len(boundary_id)-1):
     343                if boundary_id[i]+1==boundary_id[i+1]:
     344                    gauge_neighbour_id.append(i+1)
     345                else:
     346                    gauge_neighbour_id.append(-1)
     347            if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 and boundary_id[0]==0:
     348                gauge_neighbour_id.append(0)
     349            else:
     350                gauge_neighbour_id.append(-1)
     351            gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
     352            if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1:
     353                msg='incorrect number of segments'
     354                raise msg
     355            vertex_coordinates=ensure_numeric(temp)
     356            if len(vertex_coordinates)==0:
     357                msg = 'None of the sts gauges fall on the boundary'
     358                raise msg
     359        else:
     360            gauge_neighbour_id=None
    310361
    311362        if interpolation_points is not None:
     
    342393    for i, name in enumerate(quantity_names):
    343394        quantities[name] = fid.variables[name][:]
     395        if boundary_polygon is not None:
     396            #removes sts points that do not lie on boundary
     397            quantities[name] = take(quantities[name],gauge_id,1)
    344398    fid.close()
    345    
    346399
    347400    from anuga.fit_interpolate.interpolate import Interpolation_function
     
    362415                                  interpolation_points,
    363416                                  time_thinning=time_thinning,
    364                                   verbose=verbose), starttime
     417                                  verbose=verbose,gauge_neighbour_id=gauge_neighbour_id), starttime
    365418
    366419    # NOTE (Ole): Caching Interpolation function is too slow as
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5361 r5372  
    9393                             f,
    9494                             vertex_coordinates,
     95                             gauge_neighbour_id,
    9596                             point_coordinates=None,
    96                              start_blocking_len=500000,
    9797                             verbose=False):
    9898 
     
    112112        assert f.shape[0]>0,msg
    113113
    114         #print point_on_line_py(point_coordinates[3],[vertex_coordinates[0],vertex_coordinates[1]])
    115         #print vertex_coordinates[0],vertex_coordinates[1]
    116 
    117         #print point_coordinates
    118         #print vertex_coordinates
    119114        n=f.shape[0]
    120115        if n==1:
    121116            z=f*z
    122117        elif n>1:
    123             index=0#index of segment on which last point was found
    124             k=0#number of points on line
    125118            for i in range(len(point_coordinates)):
    126119                found = False
    127                 #find the segment the cetnroid lies on
    128                 #Start with the segment the last point was found on
    129                 #For n points there will be n-1 segments
    130                 for j in range(n-1):
    131                     #print 'searcing segment', index+j
    132                     #reached last segment look at first segment
    133                     if (index+j)==n-2:
    134                         index=0
    135                     if point_on_line_py(point_coordinates[i],[vertex_coordinates[(j)+index],vertex_coordinates[(j+1)+index]]):
    136                         found=True
    137                         x0=vertex_coordinates[j][0];y0=vertex_coordinates[j][1]
    138                         x1=vertex_coordinates[j+1][0];y1=vertex_coordinates[j+1][1]
    139                         x2=point_coordinates[i][0];y2=point_coordinates[i][1]
    140                        
    141                         segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
    142                         dist=sqrt((x2-x0)**2+(y2-y0)**2)
    143                         z[i]=(f[j+1]-f[j])/segment_len*dist+f[j]
    144                         #print 'element found on segment',j+index
    145                         index=j
    146                         break
    147                    
    148                
     120                for j in range(n):
     121                    if gauge_neighbour_id[j]>=0:
     122                        if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]):
     123                            found=True
     124                            x0=vertex_coordinates[j][0]
     125                            y0=vertex_coordinates[j][1]
     126                            x1=vertex_coordinates[gauge_neighbour_id[j]][0]
     127                            y1=vertex_coordinates[gauge_neighbour_id[j]][1]
     128                            x2=point_coordinates[i][0]
     129                            y2=point_coordinates[i][1]
     130                           
     131                            segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
     132                            dist=sqrt((x2-x0)**2+(y2-y0)**2)
     133                            z[i]=(f[gauge_neighbour_id[j]]-f[j])/segment_len*dist+f[j]
     134                            #print 'element found on segment'
     135                            break
     136                                 
    149137                if not found:
    150138                    z[i]=0.0
    151139                    #print 'point not on urs boundary'
    152                 else:
    153                     k+=1
    154         #print 'number of midpoints on urs boundary',k
    155         #print z
    156140        return z
    157141
     
    550534                 interpolation_points=None,
    551535                 time_thinning=1,
    552                  verbose=False):
     536                 verbose=False,
     537                 gauge_neighbour_id=None):
    553538        """Initialise object and build spatial interpolation if required
    554539
     
    749734                                                      verbose=False) # Don't clutter
    750735                    elif triangles is None and vertex_coordinates is not None:
    751                         result=interpol.interpolate_polyline(Q,vertex_coordinates,point_coordinates=self.interpolation_points)
     736                        result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points)
    752737
    753738                    #assert len(result), len(interpolation_points)
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5369 r5372  
    64416441        domain_fbound = pmesh_to_domain_instance(meshname, Domain)
    64426442        domain_fbound.set_quantity('stage', tide)
    6443         Bf = File_boundary(sts_file+'.sts', domain_fbound)
     6443        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
    64446444        Br = Reflective_boundary(domain_fbound)
    64456445        Bd=Dirichlet_boundary([2.0,220,-220])
     
    64666466            i+=1
    64676467
    6468         assert temp_fbound-temp_drchlt<epsilon
     6468        assert allclose(temp_fbound,temp_drchlt)
     6469        os.remove(sts_file+'.sts')
     6470        os.remove(meshname)
     6471
     6472    def test_file_boundary_sts2(self):
     6473        # mux2 file has points not included in boundary
     6474        # mux2 gauges are not stored with the same order as they are
     6475        # found in bounding_polygon
     6476        from anuga.shallow_water import Domain
     6477        from anuga.shallow_water import Reflective_boundary
     6478        from anuga.shallow_water import Dirichlet_boundary
     6479        from anuga.shallow_water import File_boundary
     6480        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     6481        from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
     6482        bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]]
     6483        tide=0.
     6484        time_step_count = 5
     6485        time_step = 2
     6486        lat_long_points=bounding_polygon[0:2]
     6487        lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1])
     6488        lat_long_points.insert(0,[6.0,97.01])
     6489        n=len(lat_long_points)
     6490        first_tstep=ones(n,Int)
     6491        last_tstep=(time_step_count)*ones(n,Int)
     6492        gauge_depth=20*ones(n,Float)
     6493        ha=2*ones((n,time_step_count),Float)
     6494        ua=10*ones((n,time_step_count),Float)
     6495        va=-10*ones((n,time_step_count),Float)
     6496        base_name, files = self.write_mux2(lat_long_points,
     6497                                   time_step_count, time_step,
     6498                                   first_tstep, last_tstep,
     6499                                   depth=gauge_depth,
     6500                                   ha=ha,
     6501                                   ua=ua,
     6502                                   va=va)
     6503
     6504        sts_file=base_name
     6505        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
     6506        self.delete_mux(files)
     6507
     6508        #print 'start create mesh from regions'
     6509        for i in range(len(bounding_polygon)):
     6510            zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1])
     6511        extent_res=1000000
     6512        meshname = 'urs_test_mesh' + '.tsh'
     6513        interior_regions=None
     6514        boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]}
     6515        #have to change boundary tags from last example because now bounding
     6516        #polygon starts in different place.
     6517        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
     6518                         maximum_triangle_area=extent_res,filename=meshname,
     6519                         interior_regions=interior_regions,verbose=False)
     6520       
     6521        domain_fbound = pmesh_to_domain_instance(meshname, Domain)
     6522        domain_fbound.set_quantity('stage', tide)
     6523        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
     6524        Br = Reflective_boundary(domain_fbound)
     6525        Bd=Dirichlet_boundary([2.0,220,-220])
     6526        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
     6527        finaltime=time_step*(time_step_count-1)
     6528        yieldstep=time_step
     6529        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
     6530        i=0
     6531        for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
     6532                                      skip_initial_step = False):
     6533            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
     6534            i+=1
     6535       
     6536        domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
     6537        domain_drchlt.set_quantity('stage', tide)
     6538        Br = Reflective_boundary(domain_drchlt)
     6539        Bd=Dirichlet_boundary([2.0,220,-220])
     6540        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
     6541        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
     6542        i=0
     6543        for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     6544                                      skip_initial_step = False):
     6545            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
     6546            i+=1
     6547
     6548        assert allclose(temp_fbound,temp_drchlt)
    64696549        os.remove(sts_file+'.sts')
    64706550        os.remove(meshname)
  • anuga_core/source/anuga/utilities/polygon.py

    r5321 r5372  
    1616from anuga.utilities.numerical_tools import ensure_numeric
    1717from anuga.geospatial_data.geospatial_data import ensure_absolute
     18
     19def point_on_line_py(point, line):
     20    from Numeric import fabs
     21    point = ensure_numeric(point)
     22    line = ensure_numeric(line)
     23
     24
     25    x=point[0];y=point[1]
     26    x0=line[0,0];y0=line[0,1]
     27    x1=line[1,0];y1=line[1,1]
     28    #from pylab import plot,show
     29    #plot(line[:,0],line[:,1])
     30    #plot([x],[y],'o')
     31    #show()
     32    a0 = x - x0;
     33    a1 = y - y0;
     34
     35    a_normal0 = a1;
     36    a_normal1 = -a0;
     37
     38    b0 = x1 - x0;
     39    b1 = y1 - y0;
     40
     41    #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ):
     42    eps=200
     43    #print 'normal',a_normal0*b0 + a_normal1*b1
     44    if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ):
     45        #for some reason (perhaps catastrophic cancellation) urs_boundary.py
     46        #example only works for eps=2
     47        #point is somewhere on the infinite extension of the line
     48        # FIXME (Ole): Perhaps add a tolerance here instead of 0.0
     49
     50        len_a = sqrt(a0*a0 + a1*a1);
     51        len_b = sqrt(b0*b0 + b1*b1);
     52
     53        if (a0*b0 + a1*b1 >= 0 and len_a <= len_b):
     54            return 1
     55        else:
     56            return 0
     57    else:
     58        return 0
     59
    1860
    1961def point_on_line(point, line):
Note: See TracChangeset for help on using the changeset viewer.