fenics python perm

    from fenics import*
    from mshr import*
    from math import pi, sin, cos, tan, pi, fabs, sqrt
    import numpy


def gowniejsza_funkcja(alpha):


    # Angle of bypass
    # alpha = pi/4
    # Kinematic viscosity
    nu = 1.0
    # Radius of blood vessel
    r1 = 1.5
    # Radius of bypass
    r2 = 0.5
    # Length of blood vessel
    l1 = 10
    # Length of bypass
    l2 = 6
    # Left boundary of subdomain
    d1 = -4.5
    # Right boundary of subdomain
    d2 = 4.5

    def dimensions(alpha):
            # Check if dimensions of cylinders are correct
            condition0 = (alpha < pi/2)
            condition1 = (cos(alpha)*r2 < r1)
            condition2 = (tan(alpha)*l1/2 > r1)
            if (condition0) and (condition1) and (condition2):
                    print 'Dimensions of cylinders are correct'
            else:
                    print 'Dimensions of cylinders are incorrect'

            # Check if dimensions of subdomain are correct
            condition4 = (-l1/2 < d1) and (-r1/tan(alpha) - r2/sin(alpha) > d1)
            condition5 = (l1/2 > d2) and (-r1/tan(alpha) + r2/sin(alpha) < d2)
            if (condition4) and (condition5):
                    print 'Dimensions of subdomain are correct'
            else:
                    print 'Dimensions of subdomain are incorrect'
            return

    # Generate domain
    def domain(alpha):
        C1 = Rectangle(Point(-l1/2, -r1), Point(l1/2, r1))
        C2 = Rectangle(Point(-l2, -r2), Point(0.0, r2))
        C2_alpha = CSGRotation(C2, Point(0.0, 0.0), -alpha)
        return C1 + C2_alpha

    # Define left blood vessel boundary
    class Left(SubDomain):
        def inside(self, x, on_boundary):
            return ((x[0] + l1/2) < DOLFIN_EPS) and on_boundary
    # Define right blood vessel boundary
    class Right(SubDomain):
        def inside(self, x, on_boundary):
            return ((x[0] - l1/2) < DOLFIN_EPS) and on_boundary
    # Define upper bypass boundary
    class Bypass(SubDomain):
        def inside(self, x, on_boundary):
            return (((fabs(x[0]/tan(alpha) - x[1] + (cos(alpha)*l2)/tan(alpha) + \
            sin(alpha)*l2))/sqrt((1/tan(alpha))**2 + 1)) < DOLFIN_EPS) \
            and on_boundary
    # Define bypass inflow subdomain
    class Interior(SubDomain):
        def inside(self, x, on_boundary):
            return between(x[0], (d1, d2)) and between(x[1], (-r1, r1))


    # Initialize sub-domain instances
    left = Left()
    right = Right()
    bypass = Bypass()
    interior = Interior()

    # Generate mesh
    mesh = generate_mesh(domain(alpha), 32)

    # Mark domains
    domains = CellFunction('size_t', mesh)
    domains.set_all(0)

    # Mark boundaries
    boundaries = FacetFunction('size_t', mesh)
    boundaries.set_all(0)
    left.mark(boundaries, 1)
    right.mark(boundaries, 1)
    bypass.mark(boundaries, 1)

    # Define function spaces
    W2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)
    V2 = FunctionSpace(mesh, W2)
    W1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
    V1 = FunctionSpace(mesh, W1)
    W = W2 * W1
    V = FunctionSpace(mesh, W)

    # Define (u1, p1) and f1
    def u1():
        u0 = Expression(('2*x[0]*x[1]', '-x[1]*x[1]'), degree = 2)
        return interpolate(u0, V2)
    def p1():
        p0 = Expression('x[0] + x[1]', degree = 1)
        return interpolate(p0, V1)
    def f1(NS):
        f = Expression(('NS*(-sin(2*x[0])*cos(x[0])*cos(x[0]) + \
        2*x[1]*x[1]*cos(2*x[0])*sin(2*x[0])) + cos(x[0]) + \
        nu*(2*cos(2*x[0]))', 'NS*(sin(2*x[0])*sin(2*x[0])*x[1]) + \
        nu*(4*x[1]*sin(2*x[0])) + sin(x[1])'), degree = 2, nu = nu, NS = NS)
        return interpolate(f, V2)

    # Define (u2, p2) and f2
    def u2():
        u0 = Expression(('cos(x[0])*cos(x[0])', 'x[1]*sin(2*x[0])'), degree = 2)
        return interpolate(u0, V2)
    def p2():
        p0 = Expression('sin(x[0]) + cos(x[1])', degree = 1)
        return interpolate(p0, V1)
    def f2(NS):
        f = Expression(('NS*(4*x[0]*x[1]*x[1]) + 1', 'NS*(4*x[0]*x[0]*x[1] \
        + 2*x[1]*x[1]*x[1]) + 2*nu + 1'), degree = 2, nu = nu, NS = NS)
        return interpolate(f, V2)

    # Define Neumann boundaries
    def neumann_conditions(u_bc, p_bc):
            n = FacetNormal(mesh)
            I = Identity(mesh.geometry().dim())
            return (Constant(nu)*grad(u_bc) - p_bc*I)*n
    g = neumann_conditions(u1(), p1())

    # Define Dirichlet boundaries
    def dirichlet_conditions(u_bc, p_bc):
            b_u0 = DirichletBC(V.sub(0), u_bc, boundaries, 0)
            b_p0 = DirichletBC(V.sub(1), p_bc, boundaries, 0)
            return [b_u0, b_p0]
    bcs  = dirichlet_conditions(u1(), p1())

    # Define measures
    dx = Measure('dx', domain = mesh, subdomain_data = domains)
    ds = Measure('ds', domain = mesh, subdomain_data = boundaries)

    def navier_stokes(f, bc, NS):
            # Define variational forms for Stationary Navier-Stokes
            def A(u_h, v_h):
                    return Constant(nu)*inner(grad(u_h), grad(v_h))*dx(0)
            def B(p_h, v_h):
                    return p_h*div(v_h)*dx(0)
            def C(u_h, w_h, v_h):
                    return inner(grad(u_h)*w_h, v_h)*dx(0)
            def L(v_h):
                    return inner(f, v_h)*dx(0)
            def S(v_h):
                    return inner(g, v_h)*ds(1)

            # Define variational problem
            U = TrialFunction(V)
            (u, p) = split(U)
            (v, q) = TestFunctions(V)
            U_ = Function(V)
            (u_, p_) = split(U_)
            F = A(u, v) - B(p, v) + bool(NS)*C(u, u, v) - B(q, u) - L(v) - S(v)
            F = action(F, U_)
            J = derivative(F, U_, U)
            problem = NonlinearVariationalProblem(F, U_, bc, J)

            # Create Newton solver
            solver = NonlinearVariationalSolver(problem)
            solver.parameters.newton_solver['absolute_tolerance'] = 1E-8
            solver.parameters.newton_solver['relative_tolerance'] = 1E-7
            solver.parameters.newton_solver['maximum_iterations'] = 25
            solver.parameters.newton_solver['relaxation_parameter'] = 1.0
            solver.parameters.newton_solver['linear_solver'] = 'gmres'
            solver.parameters.newton_solver['preconditioner'] = 'ilu'
            solver.parameters.newton_solver.krylov_solver['absolute_tolerance'] = 1E-9
            solver.parameters.newton_solver.krylov_solver['relative_tolerance'] = 1E-7
            solver.parameters.newton_solver.krylov_solver['maximum_iterations'] = 1000
            solver.parameters.newton_solver.krylov_solver['monitor_convergence'] = True
            solver.parameters.newton_solver.krylov_solver['nonzero_initial_guess'] = False

            # Monitor solution
            PROGRESS = 16
            set_log_level(PROGRESS)
            info(solver.parameters.newton_solver, True)

            # Solve
            solver.solve()

            return U_

    # Call navier-stokes and get sub-functions
    U_s = navier_stokes(f1(0), bcs, 0)
    (u_s, p_s) = U_s.split(U_s)
    U_ns = navier_stokes(f1(1), bcs, 1)
    (u_ns, p_ns) = U_ns.split(U_ns)

    # Compute error parameters
    def errors(u0_V, u0, p0_V, p0):
            u_max = numpy.abs(u0.vector().array()).max()
            print 'Maximum norm of velocity: ', u_max
            E_u_max = numpy.abs(u0_V.vector().array() - u0.vector().array()).max()
            print 'Maximum error for velocity: ', E_u_max
            u_L2 = norm(u0, norm_type = 'L2')
            print 'L2 norm of velocity: ', u_L2
            E_u_L2 = errornorm(u0_V, u0, norm_type = 'L2', degree_rise = 1)
            print 'L2 norm of error for velocity: ', E_u_L2
            p_max = numpy.abs(p0.vector().array()).max()
            print 'Maximum error of pressure: ', p_max
            E_p_max = numpy.abs(p0_V.vector().array() - p0.vector().array()).max()
            print 'Maximum error for pressure: ', E_p_max
            u_L2 = norm(p0, norm_type = 'L2')
            print 'L2 norm of pressure: ', u_L2
            E_p_L2 = errornorm(p0_V, p0, norm_type = 'L2', degree_rise = 1)
            print 'L2 norm of error for pressure: ', E_p_L2
            div_u = assemble(abs(div(u0))*dx(0))
            print 'L2 norm of divergence of velocity : ', div_u
            return
    print 'Error parameters for Stokes: '
    errors(u1(), u_s, p1(), p_s)
    print 'Error parameters for Navier-Stokes: '
    errors(u1(), u_ns, p1(), p_ns)

    # Mark sudomain
    interior.mark(domains, 1)

    # Define cost functionals
    def functionals(g1, g2):
        # Viscious energy dissipation
        J_1 = nu/2*assemble(abs(inner(grad(g1), grad(g1)))*dx(1))
        print 'Viscious energy dissipation: ', J_1

        # Stokes-tracking type functional
        J_2 = assemble(abs(inner(g1 - g2, g1 - g2))*dx(1))
        print 'Stokes-tracking type functional: ', J_2

        # Vorticity
        rot = grad(g1)[1,0] - grad(g1)[0,1]
        J_3 = nu/2*assemble(abs(inner(rot, rot))*dx(1))
        print 'Vorticity: ', J_3

        # Galilean vortex
        J_4 = 1/2*assemble((det(grad(g1)) + abs(det(grad(g1))))*dx(1))
        print 'Galilean vortex: ', J_4
        return [J_1, J_2, J_3, J_4]

    print 'Evaluation of cost functionals: '
    J = functionals(u_ns, u_s)

    # Save solution in VTK format
    ufile_pvd = File('Solutions/Stokes/Velocity.pvd')
    ufile_pvd << u_s
    pfile_pvd = File('Solutions/Stokes/Pressure.pvd')
    pfile_pvd << p_s
    ufile_pvd = File('Solutions/Navier-Stokes/Velocity.pvd')
    ufile_pvd << u_ns
    pfile_pvd = File('Solutions/Navier-Stokes/Pressure.pvd')
    pfile_pvd << p_ns


gowniejsza_funkcja(pi/4)