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