L2 error squared estimation
----------------------------


### Bilinear quad


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


```python
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]
```


```python
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
```




$\displaystyle \frac{9955}{12}$




```python
expand(u_h_local(x, y))
```




$\displaystyle \frac{43 x y}{2} + \frac{29 x}{2} - \frac{3 y}{2} - \frac{19}{2}$




```python
expand(u_local(x, y))
```




$\displaystyle - \frac{15 x^{2} y^{2}}{16} - \frac{45 x^{2} y}{8} - \frac{135 x^{2}}{16} + \frac{25 x y^{2}}{4} + \frac{43 x y}{2} + \frac{33 x}{4} + \frac{35 y^{2}}{16} + \frac{33 y}{8} - \frac{37}{16}$




```python

```


```python

```

bea9d8022650f9cfba8e632e49ba455f80c860e9
ipynb
Jupyter Notebook
notebooks/unit_tests_analytic_solutions.ipynb
InteractiveComputerGraphics/higher_order_embedded_fem
868fbc25f93cae32aa3caaa41a60987d4192cf1b
MIT
2021-10-19T17:11:52.000Z
2021-12-26T10:20:53.000Z
notebooks/unit_tests_analytic_solutions.ipynb
InteractiveComputerGraphics/higher_order_embedded_fem
868fbc25f93cae32aa3caaa41a60987d4192cf1b
MIT
notebooks/unit_tests_analytic_solutions.ipynb
InteractiveComputerGraphics/higher_order_embedded_fem
868fbc25f93cae32aa3caaa41a60987d4192cf1b
MIT
2021-10-20T16:13:05.000Z
2022-03-16T01:50:35.000Z
Qwen/Qwen-72B
1. YES
2. YES