import meshplex
import meshzoo
import numpy as np
from scipy.sparse import linalg
import pyfvm
from pyfvm.form_language import Boundary, dS, dV, integrate, n_dot_grad
class Poisson:
    def apply(self, u):
        return integrate(lambda x: -n_dot_grad(u(x)), dS) - integrate(lambda x: 1.0, dV)
    def dirichlet(self, u):
        return [(lambda x: u(x) - 0.0, Boundary())]
# Create mesh using meshzoo
vertices, cells = meshzoo.rectangle_tri(
    np.linspace(0.0, 2.0, 401), np.linspace(0.0, 1.0, 201)
)
mesh = meshplex.Mesh(vertices, cells)
matrix, rhs = pyfvm.discretize_linear(Poisson(), mesh)
u = linalg.spsolve(matrix, rhs)
mesh.write("out.vtk", point_data={"u": u})

mesh = meshplex.read("pacman.e")

import pyamg
ml = pyamg.smoothed_aggregation_solver(matrix)
u = ml.solve(rhs, tol=1e-10)

import pyfvm
from pyfvm.form_language import *
import meshzoo
import numpy as np
from sympy import exp
import meshplex
class Bratu:
    def apply(self, u):
        return integrate(lambda x: -n_dot_grad(u(x)), dS) - integrate(
            lambda x: 2.0 * exp(u(x)), dV
        )
    def dirichlet(self, u):
        return [(u, Boundary())]
vertices, cells = meshzoo.rectangle_tri(
    np.linspace(0.0, 2.0, 101), np.linspace(0.0, 1.0, 51)
)
mesh = meshplex.Mesh(vertices, cells)
f, jacobian = pyfvm.discretize(Bratu(), mesh)
def jacobian_solver(u0, rhs):
    from scipy.sparse import linalg
    jac = jacobian.get_linear_operator(u0)
    return linalg.spsolve(jac, rhs)
u0 = np.zeros(len(vertices))
u = pyfvm.newton(f.eval, jacobian_solver, u0)
mesh.write("out.vtk", point_data={"u": u})

import scipy.optimize
u = scipy.optimize.newton_krylov(f.eval, u0)