source: anuga_work/debug/will_power_example/make_mesh.py @ 6842

""" This Script shows how exta points and segments can be added to a msh file.

# WARNING #
* this doesn't have tests assosiated with it.  For example that the
resolutions stay as described in create_mesh_from_regions is not tested
* Also, this code causes the mesh to be generated twice, so it is not
an efficient way of doing this.

If this is used a lot you may want to wrap the functionality into a function.
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

# Standard modules
import os
import time
import sys

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Dirichlet_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.pmesh.mesh import importMeshFromFile
from anuga.utilities.polygon import read_polygon, plot_polygons

from anuga.geospatial_data.geospatial_data import Geospatial_data



#------------------------------------------------------------------------------
# Filenames
#------------------------------------------------------------------------------

outname ='mesh'
meshname ='mesh.msh'

# Create a points file
ptfile = 'points_from_a_file.csv'
file = open(ptfile,"w")
file.write("x,y\n\
304200.0, 6185300.0\n\
304300.0, 6185400.0\n\
304400.0, 6185500.0\n")
file.close()

#------------------------------------------------------------------------------
# Create the triangular mesh based on overall clipping polygon with a tagged
# boundary and interior regions defined in project.py along with
# resolutions (maximal area of per triangle) for each polygon
#------------------------------------------------------------------------------

W=304180.0
S=6185270.0
E=307650.0
N=6189040.0

bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]



create_mesh_from_regions(bounding_polygon,
                         boundary_tags={'south': [0],
                                        'east': [1],
                                        'north': [2],
                                        'west': [3]},
                         maximum_triangle_area=10000, #background mesh size
                         filename=meshname,
                         use_cache=False,
                         verbose=True)

m = importMeshFromFile(meshname)

#m.add_vertices([(W +500,S +500),(W +1000,S +1000) ])

# load points from a points file
#m.add_vertices(Geospatial_data(ptfile))

p1 = (W +500,N - 500)
p2 = (W +1000,N - 1000)
p3 = (W +1000,N - 1500) 
#m.add_points_and_segments([p1,p2,p3], [[0,1],[1,2]])

# add points and segments, where the segments are automatically added.
# This will work on the next release of ANUGA.
# It doesn't work on the current sourceforge release.
p1 = (E -1000,N - 1000)
p2 = (E -1000,N - 1500) 
#m.add_points_and_segments([p1,p2])

# load points from a points file
m.add_points_and_segments(Geospatial_data(ptfile))

m.generate_mesh()
m.export_mesh_file('end.msh')

import sys; sys.exit() 
#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------

domain = Domain(meshname, use_cache=False, verbose=True)
domain.set_name(outname) 
print domain.statistics()

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
# Setting up initial conditions for a Flood model includes:
# Setting A Base topography over the model domain
# Setting the variation of Roughness over the terrain, which may be
# varied dependent on depth say
# Setting initial water levels at BOundaries and in Lakes and Water Bodies.


tide = 1.5 # 100yr tidal level as per WCC drainage design code
base_friction = 0.015

domain.set_quantity('friction', base_friction)
domain.set_quantity('stage', tide)
domain.set_quantity('elevation', 0)


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

print 'Available boundary tags', domain.get_boundary_tags()

Bd = Dirichlet_boundary([tide,0,0])

# Boundary conditions
domain.set_boundary({'west': Bd,
                     'south': Bd,
                     'north': Bd,
                     'east': Bd})
    

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

import time
t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 2): 
    print domain.timestepping_statistics()

print 'Finished'