# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# L2 error squared estimation
# ----------------------------
#
#
# ### Bilinear quad

from sympy import *
from sympy.integrals.intpoly import polytope_integrate
from sympy.abc import x, y

# +
points = [ Point2D(-1, -1), Point2D(2, -2), Point2D(4, 1), Point2D(-2, 3)]

def phi_alpha_beta(alpha, beta, x, y):
    return (1 + alpha * x) * (1 + beta * y) / 4

# Define basis functions phi(x, y)
def phi_local(x, y):
    return [
        phi_alpha_beta(-1, -1, x, y),
        phi_alpha_beta(1, -1, x, y),
        phi_alpha_beta(1, 1, x, y),
        phi_alpha_beta(-1, 1, x, y)
    ]

# Define transformation from reference element T: K_hat -> K,
# with K being the element defined by quad.
def T(x, y):
    p = phi_local(x, y)
    return points[0] * p[0] + points[1] * p[1] + points[2] * p[2] + points[3] * p[3]
    
def u(x, y):
    return 5 * x * y + 3 * x - 2 * y - 5

def u_local(xi, eta):
    (x, y) = T(xi, eta)
    return u(x, y)

u_h_weights = [u(p[0], p[1]) for p in points]

def u_h_local(xi, eta):
    p = phi_local(xi, eta)
    u = u_h_weights
    return u[0] * p[0] + u[1] * p[1] + u[2] * p[2] + u[3] * p[3]


# +
det_J_K = Matrix(T(x, y)).jacobian(Matrix([x, y])).det()

integrand = expand(det_J_K * (u_h_local(x, y) - u_local(x, y))**2)

# Note: It may be necessary to expand the polynomial for use with polytope_integrate
#integral = polytope_integrate(reference_quad, 1)
# Note: polytope_integrate did not seem to work so well. Since we anyway integrate in the reference domain,
# which is a simple square, we can just integrate normally with simple limits
integral = integrate(integrand, (x, -1, 1), (y, -1, 1))
integral
# -

expand(u_h_local(x, y))

expand(u_local(x, y))





StarcoderdataJupyterScriptsDedupFiltered