from firedrake import *
from firedrake.utils import cached_property

from abc import ABCMeta, abstractproperty


class Problem(object):

    __metaclass__ = ABCMeta

    def __init__(self, N=None, degree=None):
        super(Problem, self).__init__()

        self.degree = degree
        self.N = N

    @property
    def comm(self):
        return self.mesh.comm

    @cached_property
    def mesh(self):
        return UnitCubeMesh(self.N, self.N, self.N)

    @abstractproperty
    def name(self):
        pass

    @abstractproperty
    def function_space(self):
        pass

    @cached_property
    def u(self):
        return Function(self.function_space, name="Solution")

    @abstractproperty
    def a(self):
        pass

    @abstractproperty
    def L(self):
        pass

    @abstractproperty
    def bcs(self):
        pass

    @cached_property
    def analytic_solution(self):
        x, y, z = SpatialCoordinate(self.mesh)
        return exp(sin(pi*x)*sin(pi*y)*sin(pi*z))

    def solver(self, parameters=None):

        problem = LinearVariationalProblem(self.a, self.L, self.u,
                                           bcs=self.bcs)
        solver = LinearVariationalSolver(problem, solver_parameters=parameters)

        return solver

    @abstractproperty
    def output(self):
        pass

    @abstractproperty
    def err(self):
        pass

    @abstractproperty
    def true_err(self):
        pass

    @abstractproperty
    def sol(self):
        pass

@@ -201,16 +201,71 @@
 ree=None
+, dimension=None,%0A                 quadrilaterals=False
 ):%0A     
@@ -347,16 +347,81 @@
 lf.N = N
+%0A        self.dim = dimension%0A        self.quads = quadrilaterals
 %0A%0A    @p
@@ -512,32 +512,389 @@
 def mesh(self):%0A
+        if self.dim == 2:%0A            return UnitSquareMesh(self.N, self.N, quadrilateral=self.quads)%0A        else:%0A            assert self.dim == 3%0A            if self.quads:%0A                base = UnitSquareMesh(self.N, self.N, quadrilateral=self.quads)%0A                return ExtrudedMesh(base, self.N, layer_height=1.0/self.N)%0A            else:%0A        
         return U
@@ -1049,39 +1049,40 @@
      pass%0A%0A    @
-cached_
+abstract
 property%0A    def
@@ -1103,61 +1103,12 @@
     
-return Function(self.function_space, name=%22Solution%22)
+pass
 %0A%0A  
@@ -1333,14 +1333,8 @@
    x
-, y, z
  = S
@@ -1361,16 +1361,46 @@
 f.mesh)%0A
+        if self.dim == 2:%0A    
         
@@ -1418,32 +1418,148 @@
 sin(
+3*
 pi*x
+%5B0%5D
 )*sin(
+3*
 pi*
-y
+x%5B1%5D))%0A        else:%0A            assert self.dim == 3%0A            return exp(sin(3*pi*x%5B0%5D
 )*sin(
+3*
 pi*
-z
+x%5B1%5D)*sin(3*pi*x%5B2%5D
 ))%0A%0A

Support 2 and 3 dim prbs and quads