"""
Script for playing with and plotting solutions assuming these are
precomputed.
"""

# Copyright (C) 2013 Marie E. Rognes

from dolfin import *

# -----------------------------------------------------------------------------
# Create and initialize application parameters
application_parameters = Parameters("Discretization")
application_parameters.add("family", "RT")
application_parameters.add("order", 1)
application_parameters.add("exact_case", 2)
application_parameters.add("meshid", "sphere_ico6")
application_parameters.parse()
info(application_parameters, True)

family = application_parameters["family"]
order = application_parameters["order"]
meshid = application_parameters["meshid"]
filename = "results/w_%s_%d_%s.xml.gz" % (family, order, meshid)

# Define approximation spaces and basis functions
mesh = Mesh("meshes/%s.xml.gz" % meshid)
global_normal = Expression(("x[0]", "x[1]", "x[2]"))
mesh.init_cell_orientations(global_normal)

V = FunctionSpace(mesh, family, order)
Q = FunctionSpace(mesh, "DG", order-1)
R = FunctionSpace(mesh, "R", 0)
W = MixedFunctionSpace((V, Q, R))

w = Function(W, filename)

(sigma, u, r) = w.split(deepcopy=True)

p = plot(u)
p.zoom(1.7)
p.write_png("results/u_%s_%d_%s" % (family, order, meshid))

S = VectorFunctionSpace(mesh, "CG", 1)
s = plot(project(sigma, S))
s.zoom(1.7)
s.write_png("results/sigma_%s_%d_%s" % (family, order, meshid))

interactive()
