removed growth eqns neumann

growth_eqns_neumann_moving_mesh.py

-import dolfin as df
-import numpy as np
-import progressbar
-
-df.set_log_level(df.LogLevel.ERROR)
-df.parameters['form_compiler']['optimize'] = True
-
-class Growth(object):
-    def __init__(self, parameters):
-        # read in parameters
-        for key in parameters:
-            setattr(self, key, parameters[key])
-
-        # create mesh, mixed element and function space
-        self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
-        conc_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
-        mixed_element = df.MixedElement([conc_element, conc_element, conc_element]) # u, v, rho
-
-        # define function space with this mixed element
-        self.function_space = df.FunctionSpace(self.mesh,mixed_element)
-
-        # define the reqd functions on this function space
-        self.f = df.Function(self.function_space)                 # f at current time
-        self.f0 = df.Function(self.function_space)                 # f at t =0
-
-
-    def advection(self, conc, v, tconc):
-        return (df.inner((v*conc).dx(0), tconc))
-
-    def reaction_rho(self, rho, trho):
-        return df.inner(self.turnover_rho * (rho - self.average_rho), trho)
-
-    def reaction_diffusion_rho(self, rho, trho):
-        return ( self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0))
-                + self.reaction_rho(rho, trho)
-                )
-
-    def setup_initial_conditions(self, u0, rho0):
-        if u0==None or rho0==None:
-            u0 = df.Expression('cos(x[0])', degree=1)
-            rho0 = df.Expression('cos(x[0])', degree=1)
-            u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
-            rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
-
-        # solve for v0
-        v0_FS = df.FunctionSpace(self.mesh, 'P', 1)
-        v0  = df.Function(v0_FS)
-        tv0 = df.TestFunction(v0_FS)
-        v0form = ( df.inner(v0, tv0)
-                + self.elasticity * df.inner(u0.dx(0), tv0.dx(0))
-                - self.lamda * df.inner(rho0.dx(0), tv0)
-                ) * df.dx
-        df.solve(v0form==0, v0)
-        df.assign(self.f0, [u0, v0, rho0])
-
-    def setup_weak_forms(self):
-        u0, v0, rho0 = df.split(self.f0)
-        u,  v, rho = df.split(self.f)
-        tu, tv, trho = df.TestFunctions(self.function_space)
-
-        uform = (df.inner((u-u0)/self.timestep, tu) 
-                - df.inner(v, tu)
-                )
-        
-        vform = (df.inner(v, tv)
-                + self.elasticity * df.inner(u.dx(0), tv.dx(0))
-                - self.lamda * df.inner(rho.dx(0), tv)
-                )
-        
-        rhoform = (df.inner((rho-rho0)/self.timestep,trho)
-                + self.advection(rho, v, trho)
-                + self.reaction_diffusion_rho(rho, trho)
-                )
-        
-        self.form = (uform + vform + rhoform) * df.dx
-
-
-    def solve(self, u0=None, rho0=None):
-        times = np.arange(0,self.maxtime,self.timestep)
-        x_array = np.zeros((len(times), self.mesh.num_vertices()))
-        u_array = np.zeros((len(times), self.mesh.num_vertices()))
-        v_array = np.zeros_like(u_array)
-        rho_array = np.zeros_like(u_array)
-
-        self.setup_initial_conditions(u0, rho0)
-        self.setup_weak_forms()
-
-        u, v, rho = self.f0.split(deepcopy=True)
-        u_array[0] = u.compute_vertex_values(self.mesh)
-        v_array[0] = v.compute_vertex_values(self.mesh) 
-        rho_array[0] = rho.compute_vertex_values(self.mesh)
-        x_array[0] = self.mesh.coordinates()[:,0]
-
-        for i in progressbar.progressbar(range(0, len(times))):
-            df.solve(self.form == 0, self.f)
-            u, v, rho = self.f.split(deepcopy=True)
-            u_array[i] = u.compute_vertex_values(self.mesh)
-            v_array[i] = v.compute_vertex_values(self.mesh)
-            rho_array[i] = rho.compute_vertex_values(self.mesh)
-            self.f0.assign(self.f)
-            
-            u0, v0, rho0 = self.f0.split()
-
-            df.ALE.move(self.mesh, u0)
-            
-            x_array[i] = self.mesh.coordinates()[:,0]
-            
-
-        return (x_array, u_array, v_array, rho_array)
-
-
-
-if __name__ == '__main__':
-    import dolfin as df
-
-    import json
-    import matplotlib.pyplot as plt
-    from matplotlib.widgets import Slider
-    
-    with open('growth_parameters.json') as jsonFile:
-        params = json.load(jsonFile)
-    assert(params['dimension']==1)
-
-    g = Growth(params)
-    (x, u, v, rho) = g.solve()
-
-    times = np.arange(0, g.maxtime, g.timestep)    
-    fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
-    axu.set_xlabel(r'$x$')
-    axu.set_ylabel(r'$u(x,t)$')
-    axv.set_ylabel(r'$v(x,t)$')
-    axrho.set_ylabel(r'$\rho(x,t)$')
-    axu.set_xlim(np.min(x), np.max(x))
-    axu.set_ylim(np.min(u), np.max(u))
-    axv.set_ylim(np.min(v), np.max(v))
-    axrho.set_ylim(np.min(rho), np.max(rho))
-
-    uplot, = axu.plot(x[0], u[0])
-    vplot, = axv.plot(x[0], v[0])
-    rhoplot, = axrho.plot(x[0], rho[0])
-
-    def update(value):
-        ti = np.abs(times-value).argmin()
-        uplot.set_xdata(x[ti])
-        uplot.set_ydata(u[ti])
-        
-        vplot.set_xdata(x[ti])
-        vplot.set_ydata(v[ti])
-        
-        rhoplot.set_xdata(x[ti])
-        rhoplot.set_ydata(rho[ti])
-
-        plt.draw()
-        
-    sax = plt.axes([0.1, 0.92, 0.7, 0.02])
-    slider = Slider(sax, r'$t/\tau$', min(times), max(times),
-            valinit=min(times), valfmt='%3.1f',
-            fc='#999999')
-    slider.drawon = False
-    slider.on_changed(update)
-    plt.show()
\ No newline at end of file