Changeset 8699

Timestamp:

Feb 15, 2013, 10:28:56 AM (12 years ago)

Author:

steve

Message:

Some fixes picked up by other users

Location:

trunk/anuga_core

Files:

• source/anuga/__init__.py (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• source/anuga/fit_interpolate/fit.py (modified) (1 diff)
• source/anuga/pmesh/mesh_interface.py (modified) (1 diff)
• validation_tests/Tests/Case_studies/Okushiri/create_okushiri.py (modified) (1 diff)

trunk/anuga_core/source/anuga/__init__.py

@@ -117 +117 @@
 from anuga.file_conversion.ferret2sww import ferret2sww     
 from anuga.file_conversion.dem2dem import dem2dem
-from anuga.file_conversion.sww2png import sww2png
+from anuga.file_conversion.sww2array import sww2array
 
 

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -740 +740 @@
         V = self.get_vertex_coordinates()
 
-        # Check each triangle
-        for i in xrange(N):
-
-            x0, y0 = V[3*i, :]
-            x1, y1 = V[3*i+1, :]
-            x2, y2 = V[3*i+2, :]
-
-            # Check that area hasn't been compromised
-            area = self.areas[i]
-            ref = -((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
-            msg += 'Wrong area: %f  %f'\
-                  %(area, ref)
-            assert abs((area - ref)/area) < epsilon, msg
-
-            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
-            msg += ' is degenerate:  area == %f' % self.areas[i]
-            assert area > 0.0, msg
-
-            # Check that points are arranged in counter clock-wise order
-            v0 = [x1-x0, y1-y0]
-            v1 = [x2-x1, y2-y1]
-            v2 = [x0-x2, y0-y2]
-            a0 = anglediff(v1, v0)
-            a1 = anglediff(v2, v1)
-            a2 = anglediff(v0, v2)
-
-            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
-            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
-            assert a0 < pi and a1 < pi and a2 < pi, msg
-
-            # Check that normals are orthogonal to edge vectors
-            # Note that normal[k] lies opposite vertex k
-
-            normal0 = self.normals[i, 0:2]
-            normal1 = self.normals[i, 2:4]
-            normal2 = self.normals[i, 4:6]
-
-            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
-
-                # Normalise
-                l_u = num.sqrt(u[0]*u[0] + u[1]*u[1])
-                l_v = num.sqrt(v[0]*v[0] + v[1]*v[1])
-
-                msg = 'Normal vector in triangle %d does not have unit length' %i
-                assert num.allclose(l_v, 1), msg
-
-                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
-
-                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
-                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
-                msg += ' Inner product is %e.' %x
-                assert x < epsilon, msg
+#        # Check each triangle
+#        for i in xrange(0):
+#
+#            x0, y0 = V[3*i, :]
+#            x1, y1 = V[3*i+1, :]
+#            x2, y2 = V[3*i+2, :]
+#
+#            # Check that area hasn't been compromised
+#            area = self.areas[i]
+#            ref = -((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+#            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
+#            msg += 'Wrong area: %f  %f'\
+#                  %(area, ref)
+#            assert abs((area - ref)/area) < epsilon, msg
+#
+#            msg = 'Triangle %i (%f,%f), (%f,%f), (%f, %f)' % (i, x0,y0,x1,y1,x2,y2)
+#            msg += ' is degenerate:  area == %f' % self.areas[i]
+#            assert area > 0.0, msg
+#
+#            # Check that points are arranged in counter clock-wise order
+#            v0 = [x1-x0, y1-y0]
+#            v1 = [x2-x1, y2-y1]
+#            v2 = [x0-x2, y0-y2]
+#            a0 = anglediff(v1, v0)
+#            a1 = anglediff(v2, v1)
+#            a2 = anglediff(v0, v2)
+#
+#            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
+#            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
+#            assert a0 < pi and a1 < pi and a2 < pi, msg
+#
+#            # Check that normals are orthogonal to edge vectors
+#            # Note that normal[k] lies opposite vertex k
+#
+#            normal0 = self.normals[i, 0:2]
+#            normal1 = self.normals[i, 2:4]
+#            normal2 = self.normals[i, 4:6]
+#
+#            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
+#
+#                # Normalise
+#                l_u = num.sqrt(u[0]*u[0] + u[1]*u[1])
+#                l_v = num.sqrt(v[0]*v[0] + v[1]*v[1])
+#
+#                msg = 'Normal vector in triangle %d does not have unit length' %i
+#                assert num.allclose(l_v, 1), msg
+#
+#                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
+#
+#                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
+#                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
+#                msg += ' Inner product is %e.' %x
+#                assert x < epsilon, msg
+
+
+        # let's try numpy constructs
+
+        x0 = V[0::3, 0]
+        y0 = V[0::3, 1]
+        x1 = V[1::3, 0]
+        y1 = V[1::3, 1]
+        x2 = V[2::3, 0]
+        y2 = V[2::3, 1]
+
+        #print 'check areas'
+        area = self.areas
+
+        ref = -((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+
+
+
+        assert num.sum(num.abs((area - ref)/area)) < epsilon, 'Error in areas'
+
+        assert num.all(area > 0.0), 'A negative area'
+
+
+        tx0 = x2 - x1
+        ty0 = y2 - y1
+        tx1 = x0 - x2
+        ty1 = y0 - y2
+        tx2 = x1 - x0
+        ty2 = y1 - y0
+
+        nx0 = self.normals[:,0]
+        ny0 = self.normals[:,1]
+        nx1 = self.normals[:,2]
+        ny1 = self.normals[:,3]
+        nx2 = self.normals[:,4]
+        ny2 = self.normals[:,5]
+
+        #print nx0.shape
+        #print tx0.shape
+        #print ny0.shape
+        #print ty0.shape
+
+
+        #print 'check edge |_ to normals'
+        assert num.all(tx0*nx0 + ty0*ny0 < epsilon), 'Normal not perpendicular to edge'
+        assert num.all(tx1*nx1 + ty1*ny1 < epsilon), 'Normal not perpendicular to edge'
+        assert num.all(tx2*nx2 + ty2*ny2 < epsilon), 'Normal not perpendicular to edge'
+
+
+        #print 'check normals are unit length'
+        assert num.all(num.abs(nx0**2 + ny0**2 - 1) < epsilon), 'Normal are not normalised'
+        assert num.all(num.abs(nx1**2 + ny1**2 - 1) < epsilon), 'Normal are not normalised'
+        assert num.all(num.abs(nx2**2 + ny2**2 - 1) < epsilon), 'Normal are not normalised'
+
+
+
+
+        #print 'check neighbours'
+
+        # 0 neighbours
+        neighs = self.neighbours
+        nids = neighs[:,0]
+        eids = self.neighbour_edges[:,0]
+
+        ids = num.arange(len(nids))
+        nnid = num.where(nids>-1, 3*nids+eids, 0)
+        
+        #assert num.all(num.where(nids>-1, neighs[nnid]-ids, 0)==0)
+
+
+
 
 

trunk/anuga_core/source/anuga/fit_interpolate/fit.py

@@ -580 +580 @@
 
         # This is very slow!
-        #mesh.check_integrity()
+        mesh.check_integrity()
 
     interp = Fit(mesh=mesh,

trunk/anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -358 +358 @@
                         verbose=verbose)
         m.export_mesh_file(filename)
+
+        return m
         

trunk/anuga_core/validation_tests/Tests/Case_studies/Okushiri/create_okushiri.py

@@ -170 +170 @@
                                  use_cache=False,
                                  verbose=True)
+
+
     
     if elevation_in_mesh is True: