| [4657] | 1 | """Create mesh and time boundary for the Okushiri island validation |
|---|
| 2 | """ |
|---|
| 3 | |
|---|
| 4 | |
|---|
| 5 | from Numeric import array, zeros, Float, allclose |
|---|
| 6 | |
|---|
| 7 | from anuga.pmesh.mesh import * |
|---|
| 8 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
|---|
| 9 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
|---|
| 10 | from anuga.geospatial_data import Geospatial_data |
|---|
| 11 | |
|---|
| 12 | import project |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | |
|---|
| 17 | #------------------------------------------------------------- |
|---|
| 18 | if __name__ == "__main__": |
|---|
| 19 | |
|---|
| 20 | |
|---|
| 21 | #-------------------------------------------------------------------------- |
|---|
| 22 | # Create the triangular mesh based on overall clipping polygon with a |
|---|
| 23 | # tagged |
|---|
| 24 | # boundary and interior regions defined in project.py along with |
|---|
| 25 | # resolutions (maximal area of per triangle) for each polygon |
|---|
| 26 | #-------------------------------------------------------------------------- |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | base_resolution = 100 # Use this to coarsen or refine entire mesh. |
|---|
| 30 | |
|---|
| 31 | # Basic geometry and bounding polygon |
|---|
| 32 | xleft = 0 |
|---|
| 33 | xright = 5.448 |
|---|
| 34 | ybottom = 0 |
|---|
| 35 | ytop = 3.402 |
|---|
| 36 | |
|---|
| 37 | point_sw = [xleft, ybottom] |
|---|
| 38 | point_se = [xright, ybottom] |
|---|
| 39 | point_nw = [xleft, ytop] |
|---|
| 40 | point_ne = [xright, ytop] |
|---|
| 41 | |
|---|
| 42 | bounding_polygon = [point_se, |
|---|
| 43 | point_ne, |
|---|
| 44 | point_nw, |
|---|
| 45 | point_sw] |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | # Localised refined area for gulleys |
|---|
| 49 | xl = 4.8 |
|---|
| 50 | xr = 5.3 |
|---|
| 51 | yb = 1.6 |
|---|
| 52 | yt = 2.3 |
|---|
| 53 | p0 = [xl, yb] |
|---|
| 54 | p1 = [xr, yb] |
|---|
| 55 | p2 = [xr, yt] |
|---|
| 56 | p3 = [xl, yt] |
|---|
| 57 | |
|---|
| 58 | gulleys = [p0, p1, p2, p3] |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | # Island area and drawdown region (original) |
|---|
| 62 | #island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5] |
|---|
| 63 | #island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3] |
|---|
| 64 | #island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3] |
|---|
| 65 | #island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3] |
|---|
| 66 | #island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3] |
|---|
| 67 | #island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2] |
|---|
| 68 | #island_6 = [xl-.01, yb] #OK |
|---|
| 69 | #island_7 = [xl-.01, yt] #OK |
|---|
| 70 | |
|---|
| 71 | |
|---|
| 72 | # Island area and drawdown region |
|---|
| 73 | island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5] |
|---|
| 74 | island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3] |
|---|
| 75 | island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3] |
|---|
| 76 | island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3] |
|---|
| 77 | island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3] |
|---|
| 78 | island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8] |
|---|
| 79 | island_6 = [xl-.01, yb] # Keep right edge just off the gulleys |
|---|
| 80 | island_7 = [xl-.01, yt] |
|---|
| 81 | |
|---|
| 82 | island = [island_0, island_1, island_2, |
|---|
| 83 | island_3, island_4, island_5, |
|---|
| 84 | island_6, island_7] |
|---|
| 85 | |
|---|
| 86 | |
|---|
| 87 | # Region spanning half right hand side of domain just inside boundary (org) |
|---|
| 88 | #rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02] |
|---|
| 89 | #rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02] |
|---|
| 90 | #rhs_se = [xright-0.02, ybottom+0.02] |
|---|
| 91 | #rhs_ne = [xright-0.02, ytop-0.02] |
|---|
| 92 | |
|---|
| 93 | # Region spanning half right hand side of domain just inside boundary |
|---|
| 94 | rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4] |
|---|
| 95 | rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5] |
|---|
| 96 | rhs_se = [xright-0.1, ybottom+0.2] |
|---|
| 97 | rhs_ne = [xright-0.1, ytop-0.2] |
|---|
| 98 | |
|---|
| 99 | rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw] |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | # Interior regions and creation of mesh |
|---|
| 103 | interior_regions = [[rhs_region, 0.0005], |
|---|
| 104 | [island, 0.0002*base_resolution], |
|---|
| 105 | [gulleys, 0.00002*base_resolution]] |
|---|
| 106 | |
|---|
| 107 | meshname = project.mesh_filename + '.msh' |
|---|
| 108 | m = create_mesh_from_regions(bounding_polygon, |
|---|
| 109 | boundary_tags={'wall': [0, 1, 3], |
|---|
| 110 | 'wave': [2]}, |
|---|
| 111 | maximum_triangle_area=0.1*base_resolution, |
|---|
| 112 | interior_regions=interior_regions, |
|---|
| 113 | filename=project.mesh_filename, |
|---|
| 114 | verbose=True) |
|---|
| 115 | |
|---|
