#!/usr/bin/env python
# coding: utf-8

# You can actually do a lot of interesting things -- not related to hyperbolic PDEs -- with the building blocks of PyClaw.  Let's see how.
# 
# To run this notebook, you'll first need to [install PyClaw](http://www.clawpack.org/pyclaw/started.html) if you haven't already.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from clawpack import pyclaw
import numpy as np
import matplotlib.pyplot as plt


# The `pyclaw.Geometry` module contains things of general use when dealing with structured grids in 1, 2, or 3 dimensions.  If you place the cursor at the end of the next cell and hit `<tab>`, you'll see what is included:
# 
# - Dimension
# - Patch
# - Grid
# - Domain

# In[ ]:


pyclaw.geometry.


# # `pyclaw.Dimension`

# The building block for all PyClaw geometry is the Dimension object.  The docstring explains most of what it can do.

# In[2]:


print pyclaw.Dimension.__doc__


# Let's instantiate a `Dimension` object and see what we can do with it.

# In[3]:


x = pyclaw.Dimension(0.,1.,10) # Dimension with 10 intervals in [0,1]


# A Dimension object is essentially an equipartitioning of an interval.  The four arguments used to initialize it are, in order:
# 
# - A name
# - The left end of the interval (`lower`)
# - The right end of the interval (`upper`)
# - The number of partitions, or cells (`n`)
# 
# Printing the Dimension object gives us essentially this information back:

# In[4]:


print x
x.centers


# Notice that the printed statement also includes `delta`, the width of a single partition.  The Dimension knows a lot more.  For instance, it can tell us the locations of the cell centers:

# In[5]:


print x.centers


# and of the cell edges (often referred to as *interfaces* in finite volume terminology):

# In[6]:


print x.nodes


# ## Modifying a Dimension
# The Dimension is an interactive object: if we change one of its properties, the others are updated automatically.  Here, we change the number of cells from 10 to 8 and the cell centers are automatically respaced.

# In[7]:


x.num_cells=8
print x
x.centers


# Similarly, we can change the boundary locations:

# In[8]:


x.lower=-0.5
print x
x.centers


# Using the cell and edge coordinates, we can easily plot a function that is defined piecewise over these intervals:

# In[9]:


q = np.exp(x.centers)
print q
for i in range(x.num_cells):
    plt.plot([x.nodes[i],x.edges[i+1]],[q[i],q[i]],'b',lw=2)
    if i<x.num_cells-1:
        plt.plot([x.nodes[i+1],x.edges[i+1]],[q[i],q[i+1]],'b',lw=2)


# # `pyclaw.Grid`

# We can take one or more Dimension objects and form a `Grid` with them.  The Grid is just a partitioning of the tensor product of the Dimensions, where the partitioning is given by the tensor product of the Dimension intervals.  If that sounds complicated, just know that the grid can give us multidimensional coordinates, similar to the MATLAB or Numpy command `meshgrid`.

# One-dimensional grids are not very interesting, so let's create a two-dimensional grid:

# In[10]:


x = pyclaw.Dimension(-1.,1.,10,name='x')
y = pyclaw.Dimension(-1.,1.,10,name='y')
print x
print y


# In[11]:


grid = pyclaw.geometry.Grid((x,y))
print grid


# The Grid is also fully interactive, so we can change its properties (or those of its constituent Dimensions) and all other properties automatically update.

# In[12]:


x.num_cells = 12
x.upper = 3.
print grid


# We can get the coordinates of all the cell centers:

# In[13]:


X, Y = grid.c_centers
print Y


# We can plot the grid, optionally marking the locations of nodes and centers:

# In[14]:


grid.plot(mark_nodes=True,mark_centers=True);


# The plot above is simply drawing all the lines between adjacent nodes of the grid.

# In numerical PDEs, we often use *ghost cells* around the edge of the grid to implement boundary conditions.  The `grid.plot` function knows about those too:

# In[15]:


grid.plot(num_ghost=2);


# # Mappings

# PyClaw grids need not be rectangular.  We can define a grid over any region that results from a mapping of a rectangle.  For instance, here is a mapping that takes a square to a circle (from [the SIAM Review paper of Calhoun, Helzel, & LeVeque](http://faculty.washington.edu/rjl/pubs/circles/)):

# In[16]:


def square2circle(xc,yc,r1=1.0):
    d = np.maximum(np.abs(xc),np.abs(yc))
    r = np.sqrt(xc**2 + yc**2)
    r = np.maximum(r, 1.e-10)
    xp = r1 * d * xc/r
    yp = r1 * d * yc/r
    return [xp, yp]


# The `grid.mapc2p` function is the identity mapping by default.  Let's create a new grid and associate with it the circle mapping.

# In[17]:


x.upper = 1.

grid.mapc2p = square2circle
grid.plot(num_ghost=1);