| 1 | ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py |
|---|
| 2 | |
|---|
| 3 | #!/usr/bin/env python |
|---|
| 4 | # |
|---|
| 5 | |
|---|
| 6 | import unittest |
|---|
| 7 | |
|---|
| 8 | #from numpy.oldnumeric import allclose, array |
|---|
| 9 | from numpy import allclose, array |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | #from anuga.pyvolution.pmesh2domain import * |
|---|
| 13 | from pmesh2domain import * |
|---|
| 14 | |
|---|
| 15 | from anuga.shallow_water import Domain,\ |
|---|
| 16 | Reflective_boundary, Dirichlet_boundary,\ |
|---|
| 17 | Transmissive_boundary |
|---|
| 18 | |
|---|
| 19 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
|---|
| 20 | from anuga.pmesh.mesh import importMeshFromFile |
|---|
| 21 | |
|---|
| 22 | class Test_pmesh2domain(unittest.TestCase): |
|---|
| 23 | |
|---|
| 24 | def setUp(self): |
|---|
| 25 | pass |
|---|
| 26 | |
|---|
| 27 | def tearDown(self): |
|---|
| 28 | pass |
|---|
| 29 | |
|---|
| 30 | def test_pmesh2Domain(self): |
|---|
| 31 | import os |
|---|
| 32 | import tempfile |
|---|
| 33 | |
|---|
| 34 | fileName = tempfile.mktemp(".tsh") |
|---|
| 35 | file = open(fileName,"w") |
|---|
| 36 | file.write("4 3 # <vertex #> <x> <y> [attributes]\n \ |
|---|
| 37 | 0 0.0 0.0 0.0 0.0 0.01 \n \ |
|---|
| 38 | 1 1.0 0.0 10.0 10.0 0.02 \n \ |
|---|
| 39 | 2 0.0 1.0 0.0 10.0 0.03 \n \ |
|---|
| 40 | 3 0.5 0.25 8.0 12.0 0.04 \n \ |
|---|
| 41 | # Vert att title \n \ |
|---|
| 42 | elevation \n \ |
|---|
| 43 | stage \n \ |
|---|
| 44 | friction \n \ |
|---|
| 45 | 2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\ |
|---|
| 46 | 0 0 3 2 -1 -1 1 dsg\n\ |
|---|
| 47 | 1 0 1 3 -1 0 -1 ole nielsen\n\ |
|---|
| 48 | 4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\ |
|---|
| 49 | 0 1 0 2 \n\ |
|---|
| 50 | 1 0 2 3 \n\ |
|---|
| 51 | 2 2 3 \n\ |
|---|
| 52 | 3 3 1 1 \n\ |
|---|
| 53 | 3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \ |
|---|
| 54 | 0 216.0 -86.0 \n \ |
|---|
| 55 | 1 160.0 -167.0 \n \ |
|---|
| 56 | 2 114.0 -91.0 \n \ |
|---|
| 57 | 3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \ |
|---|
| 58 | 0 0 1 0 \n \ |
|---|
| 59 | 1 1 2 0 \n \ |
|---|
| 60 | 2 2 0 0 \n \ |
|---|
| 61 | 0 # <x> <y> ...Mesh Holes... \n \ |
|---|
| 62 | 0 # <x> <y> <attribute>...Mesh Regions... \n \ |
|---|
| 63 | 0 # <x> <y> <attribute>...Mesh Regions, area... \n\ |
|---|
| 64 | #Geo reference \n \ |
|---|
| 65 | 56 \n \ |
|---|
| 66 | 140 \n \ |
|---|
| 67 | 120 \n") |
|---|
| 68 | file.close() |
|---|
| 69 | |
|---|
| 70 | tags = {} |
|---|
| 71 | b1 = Dirichlet_boundary(conserved_quantities = array([0.0])) |
|---|
| 72 | b2 = Dirichlet_boundary(conserved_quantities = array([1.0])) |
|---|
| 73 | b3 = Dirichlet_boundary(conserved_quantities = array([2.0])) |
|---|
| 74 | tags["1"] = b1 |
|---|
| 75 | tags["2"] = b2 |
|---|
| 76 | tags["3"] = b3 |
|---|
| 77 | |
|---|
| 78 | domain = pmesh_to_domain_instance(fileName, Domain) |
|---|
| 79 | os.remove(fileName) |
|---|
| 80 | #print "domain.tagged_elements", domain.tagged_elements |
|---|
| 81 | ## check the quantities |
|---|
| 82 | #print domain.quantities['elevation'].vertex_values |
|---|
| 83 | answer = [[0., 8., 0.], |
|---|
| 84 | [0., 10., 8.]] |
|---|
| 85 | assert allclose(domain.quantities['elevation'].vertex_values, |
|---|
| 86 | answer) |
|---|
| 87 | |
|---|
| 88 | #print domain.quantities['stage'].vertex_values |
|---|
| 89 | answer = [[0., 12., 10.], |
|---|
| 90 | [0., 10., 12.]] |
|---|
| 91 | assert allclose(domain.quantities['stage'].vertex_values, |
|---|
| 92 | answer) |
|---|
| 93 | |
|---|
| 94 | #print domain.quantities['friction'].vertex_values |
|---|
| 95 | answer = [[0.01, 0.04, 0.03], |
|---|
| 96 | [0.01, 0.02, 0.04]] |
|---|
| 97 | assert allclose(domain.quantities['friction'].vertex_values, |
|---|
| 98 | answer) |
|---|
| 99 | |
|---|
| 100 | #print domain.quantities['friction'].vertex_values |
|---|
| 101 | assert allclose(domain.tagged_elements['dsg'][0],0) |
|---|
| 102 | assert allclose(domain.tagged_elements['ole nielsen'][0],1) |
|---|
| 103 | |
|---|
| 104 | self.failUnless( domain.boundary[(1, 0)] == '1', |
|---|
| 105 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 106 | self.failUnless( domain.boundary[(1, 2)] == '2', |
|---|
| 107 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 108 | self.failUnless( domain.boundary[(0, 1)] == '3', |
|---|
| 109 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 110 | self.failUnless( domain.boundary[(0, 0)] == 'exterior', |
|---|
| 111 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 112 | #print "domain.boundary",domain.boundary |
|---|
| 113 | self.failUnless( len(domain.boundary) == 4, |
|---|
| 114 | "test_pmesh2Domain Too many boundaries") |
|---|
| 115 | #FIXME change to use get_xllcorner |
|---|
| 116 | #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner |
|---|
| 117 | self.failUnless(domain.geo_reference.xllcorner == 140.0, |
|---|
| 118 | "bad geo_referece") |
|---|
| 119 | #************ |
|---|
| 120 | |
|---|
| 121 | def test_pmesh2Domain_instance(self): |
|---|
| 122 | import os |
|---|
| 123 | import tempfile |
|---|
| 124 | |
|---|
| 125 | fileName = tempfile.mktemp(".tsh") |
|---|
| 126 | file = open(fileName,"w") |
|---|
| 127 | file.write("4 3 # <vertex #> <x> <y> [attributes]\n \ |
|---|
| 128 | 0 0.0 0.0 0.0 0.0 0.01 \n \ |
|---|
| 129 | 1 1.0 0.0 10.0 10.0 0.02 \n \ |
|---|
| 130 | 2 0.0 1.0 0.0 10.0 0.03 \n \ |
|---|
| 131 | 3 0.5 0.25 8.0 12.0 0.04 \n \ |
|---|
| 132 | # Vert att title \n \ |
|---|
| 133 | elevation \n \ |
|---|
| 134 | stage \n \ |
|---|
| 135 | friction \n \ |
|---|
| 136 | 2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\ |
|---|
| 137 | 0 0 3 2 -1 -1 1 dsg\n\ |
|---|
| 138 | 1 0 1 3 -1 0 -1 ole nielsen\n\ |
|---|
| 139 | 4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\ |
|---|
| 140 | 0 1 0 2 \n\ |
|---|
| 141 | 1 0 2 3 \n\ |
|---|
| 142 | 2 2 3 \n\ |
|---|
| 143 | 3 3 1 1 \n\ |
|---|
| 144 | 3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \ |
|---|
| 145 | 0 216.0 -86.0 \n \ |
|---|
| 146 | 1 160.0 -167.0 \n \ |
|---|
| 147 | 2 114.0 -91.0 \n \ |
|---|
| 148 | 3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \ |
|---|
| 149 | 0 0 1 0 \n \ |
|---|
| 150 | 1 1 2 0 \n \ |
|---|
| 151 | 2 2 0 0 \n \ |
|---|
| 152 | 0 # <x> <y> ...Mesh Holes... \n \ |
|---|
| 153 | 0 # <x> <y> <attribute>...Mesh Regions... \n \ |
|---|
| 154 | 0 # <x> <y> <attribute>...Mesh Regions, area... \n\ |
|---|
| 155 | #Geo reference \n \ |
|---|
| 156 | 56 \n \ |
|---|
| 157 | 140 \n \ |
|---|
| 158 | 120 \n") |
|---|
| 159 | file.close() |
|---|
| 160 | |
|---|
| 161 | mesh_instance = importMeshFromFile(fileName) |
|---|
| 162 | |
|---|
| 163 | tags = {} |
|---|
| 164 | b1 = Dirichlet_boundary(conserved_quantities = array([0.0])) |
|---|
| 165 | b2 = Dirichlet_boundary(conserved_quantities = array([1.0])) |
|---|
| 166 | b3 = Dirichlet_boundary(conserved_quantities = array([2.0])) |
|---|
| 167 | tags["1"] = b1 |
|---|
| 168 | tags["2"] = b2 |
|---|
| 169 | tags["3"] = b3 |
|---|
| 170 | |
|---|
| 171 | domain = pmesh_instance_to_domain_instance(mesh_instance, Domain) |
|---|
| 172 | |
|---|
| 173 | os.remove(fileName) |
|---|
| 174 | #print "domain.tagged_elements", domain.tagged_elements |
|---|
| 175 | ## check the quantities |
|---|
| 176 | #print domain.quantities['elevation'].vertex_values |
|---|
| 177 | answer = [[0., 8., 0.], |
|---|
| 178 | [0., 10., 8.]] |
|---|
| 179 | assert allclose(domain.quantities['elevation'].vertex_values, |
|---|
| 180 | answer) |
|---|
| 181 | |
|---|
| 182 | #print domain.quantities['stage'].vertex_values |
|---|
| 183 | answer = [[0., 12., 10.], |
|---|
| 184 | [0., 10., 12.]] |
|---|
| 185 | assert allclose(domain.quantities['stage'].vertex_values, |
|---|
| 186 | answer) |
|---|
| 187 | |
|---|
| 188 | #print domain.quantities['friction'].vertex_values |
|---|
| 189 | answer = [[0.01, 0.04, 0.03], |
|---|
| 190 | [0.01, 0.02, 0.04]] |
|---|
| 191 | assert allclose(domain.quantities['friction'].vertex_values, |
|---|
| 192 | answer) |
|---|
| 193 | |
|---|
| 194 | #print domain.quantities['friction'].vertex_values |
|---|
| 195 | assert allclose(domain.tagged_elements['dsg'][0],0) |
|---|
| 196 | assert allclose(domain.tagged_elements['ole nielsen'][0],1) |
|---|
| 197 | |
|---|
| 198 | self.failUnless( domain.boundary[(1, 0)] == '1', |
|---|
| 199 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 200 | self.failUnless( domain.boundary[(1, 2)] == '2', |
|---|
| 201 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 202 | self.failUnless( domain.boundary[(0, 1)] == '3', |
|---|
| 203 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 204 | self.failUnless( domain.boundary[(0, 0)] == 'exterior', |
|---|
| 205 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 206 | #print "domain.boundary",domain.boundary |
|---|
| 207 | self.failUnless( len(domain.boundary) == 4, |
|---|
| 208 | "test_pmesh2Domain Too many boundaries") |
|---|
| 209 | #FIXME change to use get_xllcorner |
|---|
| 210 | #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner |
|---|
| 211 | self.failUnless(domain.geo_reference.xllcorner == 140.0, |
|---|
| 212 | "bad geo_referece") |
|---|
| 213 | |
|---|
| 214 | #*********** |
|---|
| 215 | def old_test_tags_to_boundaries (self): |
|---|
| 216 | meshDict = {} |
|---|
| 217 | p0 = [0.0, 0.0] |
|---|
| 218 | p1 = [1.0, 0.0] |
|---|
| 219 | p2 = [0.0, 1.0] |
|---|
| 220 | p3 = [0.646446609407, 0.353553390593] |
|---|
| 221 | meshDict['vertices'] = [p0,p1,p2,p3] |
|---|
| 222 | meshDict['vertex_attributes'] = [[0.0, 0.0,7.0],[10.0, 0.0,7.0],[0.0, 10.0,7.0],[6.46446609407, 3.53553390593,7.0]] |
|---|
| 223 | meshDict['triangles'] = [[0,3,2],[0,1,3]] |
|---|
| 224 | meshDict['triangle_tags'] = [6.6,6.6] |
|---|
| 225 | meshDict['triangle_neighbors'] = [[-1,-1,1],[-1,0,-1]] |
|---|
| 226 | meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]] |
|---|
| 227 | meshDict['segment_tags'] = [2,3,1,1] |
|---|
| 228 | |
|---|
| 229 | domain = Domain.pmesh_dictionary_to_domain(meshDict) |
|---|
| 230 | |
|---|
| 231 | #domain.set_tag_dict(tag_dict) |
|---|
| 232 | #Boundary tests |
|---|
| 233 | b1 = Dirichlet_boundary(conserved_quantities = array([0.0])) |
|---|
| 234 | b2 = Dirichlet_boundary(conserved_quantities = array([1.0])) |
|---|
| 235 | b3 = Dirichlet_boundary(conserved_quantities = array([1.0])) |
|---|
| 236 | #test adding a boundary |
|---|
| 237 | tags = {} |
|---|
| 238 | tags[1] = b1 |
|---|
| 239 | tags[2] = b2 |
|---|
| 240 | tags[3] = b3 |
|---|
| 241 | domain.set_boundary(tags) |
|---|
| 242 | inverted_id = Volume.instances[0].neighbours[0] |
|---|
| 243 | id = -(inverted_id+1) |
|---|
| 244 | boundary_value = Boundary_value.instances[id] |
|---|
| 245 | boundary_obj = boundary_value.boundary_object |
|---|
| 246 | |
|---|
| 247 | self.failUnless( boundary_obj == b1, |
|---|
| 248 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 249 | |
|---|
| 250 | inverted_id = Volume.instances[0].neighbours[1] |
|---|
| 251 | id = -(inverted_id+1) |
|---|
| 252 | boundary_value = Boundary_value.instances[id] |
|---|
| 253 | boundary_obj = boundary_value.boundary_object |
|---|
| 254 | self.failUnless( boundary_obj == b3, |
|---|
| 255 | "test_tags_to_boundaries failed. Single boundary wasn't added.") |
|---|
| 256 | |
|---|
| 257 | #------------------------------------------------------------- |
|---|
| 258 | if __name__ == "__main__": |
|---|
| 259 | suite = unittest.makeSuite(Test_pmesh2domain, 'test') |
|---|
| 260 | runner = unittest.TextTestRunner() |
|---|
| 261 | runner.run(suite) |
|---|
| 262 | |
|---|
| 263 | |
|---|
| 264 | |
|---|
| 265 | |
|---|
