source: anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py @ 4859

Last change on this file since 4859 was 4859, checked in by duncan, 16 years ago

check the last triangle fit/interp speed up. Upgrading the code to time runs

File size: 4.3 KB
Line 
1"""Least squares interpolation.
2
3   Implements a least-squares interpolation.
4
5   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
6   Geoscience Australia, 2004.
7
8DESIGN ISSUES
9* what variables should be global?
10- if there are no global vars functions can be moved around alot easier
11
12* The public interface
13__init__
14interpolate
15interpolate_block
16
17"""
18
19import time
20import os
21from warnings import warn
22
23from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
24     ArrayType, allclose, take, NewAxis, arange
25
26from anuga.caching.caching import cache
27from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
28
29from anuga.utilities.sparse import Sparse, Sparse_CSR
30from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
31from anuga.utilities.numerical_tools import ensure_numeric
32from anuga.utilities.polygon import in_and_outside_polygon
33from anuga.utilities.quad import build_quadtree
34
35from anuga.coordinate_transforms.geo_reference import Geo_reference
36from anuga.geospatial_data.geospatial_data import Geospatial_data, \
37     ensure_absolute
38from anuga.fit_interpolate.search_functions import set_last_triangle
39
40# tests fail if 2 is used
41MAX_VERTICES_PER_CELL = 13 # A value of 8 or lower can cause problems,
42                           # if a vert has 9 triangles.
43                           
44build_quadtree_time = 0
45
46def get_build_quadtree_time():
47    return build_quadtree_time
48
49class FitInterpolate:
50   
51    def __init__(self,
52                 vertex_coordinates=None,
53                 triangles=None,
54                 mesh=None,
55                 mesh_origin=None,
56                 verbose=False,
57                 max_vertices_per_cell=None):
58
59
60        """ Build interpolation matrix mapping from
61        function values at vertices to function values at data points
62
63        Pass in a mesh instance or vertex_coordinates and triangles
64        and optionally mesh_origin
65       
66        Inputs:
67
68          vertex_coordinates: List of coordinate pairs [xi, eta] of
69              points constituting a mesh (or an m x 2 Numeric array or
70              a geospatial object)
71              Points may appear multiple times
72              (e.g. if vertices have discontinuities)
73
74          triangles: List of 3-tuples (or a Numeric array) of
75              integers representing indices of all vertices in the mesh.
76
77        mesh: A mesh instance describing the mesh.
78
79          mesh_origin: A geo_reference object or 3-tuples consisting of
80              UTM zone, easting and northing.
81              If specified vertex coordinates are assumed to be
82              relative to their respective origins.
83
84          max_vertices_per_cell: Number of vertices in a quad tree cell
85          at which the cell is split into 4.
86
87          Note: Don't supply a vertex coords as a geospatial object and
88              a mesh origin, since geospatial has its own mesh origin.
89        """
90        global build_quadtree_time
91        if max_vertices_per_cell == None:
92            max_vertices_per_cell = MAX_VERTICES_PER_CELL
93        if mesh is None:
94            # Fixme (DSG) Throw errors if triangles or vertex_coordinates
95            # are None
96           
97            #Convert input to Numeric arrays
98            triangles = ensure_numeric(triangles, Int)
99            vertex_coordinates = ensure_absolute(vertex_coordinates,
100                                                 geo_reference = mesh_origin)
101
102            if verbose: print 'FitInterpolate: Building mesh'       
103            self.mesh = Mesh(vertex_coordinates, triangles)
104            self.mesh.check_integrity()
105        else:
106            self.mesh = mesh
107       
108        if verbose: print 'FitInterpolate: Building quad tree'
109        # This stores indices of vertices
110        t0 = time.time()
111        #print "self.mesh.get_extent(absolute=True)", \
112              #self.mesh.get_extent(absolute=True)
113        self.root = build_quadtree(self.mesh,
114                                   max_points_per_cell = max_vertices_per_cell)
115        #print "self.root",self.root.show()
116       
117        build_quadtree_time =  time.time()-t0
118        set_last_triangle()
119       
120    def __repr__(self):
121        return 'Interpolation object based on: ' + repr(self.mesh)
Note: See TracBrowser for help on using the repository browser.