<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <h1>Table of Contents<span class="tocSkip"></span></h1>
<div class="toc"><ul class="toc-item"></ul></div>
Step1:  We test the system 
\begin{equation}\label{eq
Step2:  Since we have already tested the assembly, we focus here on the linear system. In particular
Step3:  To test extract_hanging_nodes and resolve_hanging_nodes we construct a simple mesh with hanging_nodes.
Step4:  For the assembly, we must define the bilinear form
\begin{equation}
a(u,v) = \int_\Omega \nabla u \cdot \nabla v dx =\int_\Omega u_x v_x + u_y v_y dx, \ \ \  \forall v \in H^1_0(\Omega)
\end{equation}
and the linear form
\begin{equation}
L(v) = \int_\Omega f v dx,  \qquad \ \ \  \forall v \in H_0^1(\Omega)
\end{equation}
where $\Omega = [0,1]^2$ and $H^1_0(\Omega) = {v\in H^1(\Omega)
Step5:  Let's visualize the meshes.
Step6:  It looks like the following dofs from mesh1 and mesh2 are equivalent 
mesh1 -> mesh2
- 0 -> 0
- 1 -> 1
- 4 -> 4
- 5 -> 5
- 8 -> 8
If we restrict to these, we should get the same matrix.
Step7:  Check that A1 and A2 coincide when restricting to the nodes
Step8:  Mark Dirichlet Regions on Meshes
Step9:  Now extract Dirichlet nodes
<ASSISTANT_TASK:>
Python Code:
# Set Path
import sys
sys.path.append('../../src/')
%autoreload 2
# Import Libraries
from fem import Function
from fem import QuadFE
from fem import DofHandler
from fem import Kernel
from fem import Basis
from fem import Form
from fem import Assembler
from fem import LinearSystem
from plot import Plot
from mesh import convert_to_array
from mesh import QuadMesh
from mesh import Mesh1D
import matplotlib.pyplot as plt
import scipy.sparse as sp
import numpy as np
% matplotlib inline
plt.rcParams['figure.figsize'] = [7, 7]
Explanation: <h1>Table of Contents<span class="tocSkip"></span></h1>
<div class="toc"><ul class="toc-item"></ul></div>
End of explanation
#
# Define the element
#
Q1 = QuadFE(2, 'Q1')
Explanation: We test the system 
\begin{equation}\label{eq:elliptic}
- u_{xx} - u_{yy} = 0 
\end{equation}
subject to Dirichlet conditions
\begin{align}
u(0,y) &= 0 \label{eq:dirichlet_at_xis0}\
u(1,y) &= 1 \label{eq:dirichlet_at_xis1}
\end{align}
Whose exact solution is 
\begin{equation}
u_e(x,y) = x. 
\end{equation}
We use a Galerkin approximation with $Q_1$ elements.
End of explanation
mesh1 = QuadMesh(resolution=(2,2))
Explanation: Since we have already tested the assembly, we focus here on the linear system. In particular:
Marking and extracting Dirichlet boundary conditions
Extracting hanging nodes
both by (i) eliminating the variables from the system (compressed=True) and (ii) by replacing affected equations with explicit Dirichlet data or interpolation formulae.  
We first test solving with Dirichlet conditions. To that end we define our first mesh.
End of explanation
mesh2 = QuadMesh(resolution=(2,2))
mesh2.cells.get_child(2).mark(1)
mesh2.cells.refine(refinement_flag=1)
Explanation: To test extract_hanging_nodes and resolve_hanging_nodes we construct a simple mesh with hanging_nodes.
End of explanation
#
# Weak form
#
# Kernel functions 
one = Function(1, 'constant')
zero = Function(0, 'constant')
# Basis functions 
u = Basis(Q1, 'u')
ux = Basis(Q1, 'ux')
uy = Basis(Q1, 'uy')
# Forms
ax = Form(kernel=Kernel(one), trial=ux, test=ux)
ay = Form(kernel=Kernel(one), trial=uy, test=uy)
L = Form(kernel=Kernel(zero), test=u)
# Assembler for mesh1
assembler1 = Assembler([ax, ay, L], mesh1)
assembler1.assemble()
# Assembler for mesh2
assembler2 = Assembler([ax,ay,L], mesh2)
assembler2.assemble()
Explanation: For the assembly, we must define the bilinear form
\begin{equation}
a(u,v) = \int_\Omega \nabla u \cdot \nabla v dx =\int_\Omega u_x v_x + u_y v_y dx, \ \ \  \forall v \in H^1_0(\Omega)
\end{equation}
and the linear form
\begin{equation}
L(v) = \int_\Omega f v dx,  \qquad \ \ \  \forall v \in H_0^1(\Omega)
\end{equation}
where $\Omega = [0,1]^2$ and $H^1_0(\Omega) = {v\in H^1(\Omega): v(0,\cdot) = v(1,\cdot)=0}$
End of explanation
# Get dofhandlers
dh1 = assembler1.dofhandlers['Q1']
dh2 = assembler2.dofhandlers['Q1']
# Plotting mesh 1
plot = Plot()
plot.mesh(mesh1, dofhandler=dh1, dofs=True)
# Plotting mesh 2
plot = Plot()
plot.mesh(mesh2, dofhandler=dh2, dofs=True)
Explanation: Let's visualize the meshes.
End of explanation
# Assembled matrices
# Mesh1 
# bilinear
rows = assembler1.af[0]['bilinear']['rows']
cols = assembler1.af[0]['bilinear']['cols']
vals = assembler1.af[0]['bilinear']['vals']
dofs = assembler1.af[0]['bilinear']['row_dofs']
A1 = sp.coo_matrix((vals, (rows, cols)))
A1 = A1.todense()
# linear
b1 = assembler1.af[0]['linear']['vals']
# number of dofs 
n = len(dofs)
# Print
print('Mesh 1')
print('A1 = \n', 6*A1)
print('b1 = \n', 6*b1)
print('n_dofs=', n)
print('='*60)
#
# Mesh2 
# 
# bilinear
rows = assembler2.af[0]['bilinear']['rows']
cols = assembler2.af[0]['bilinear']['cols']
vals = assembler2.af[0]['bilinear']['vals']
dofs = assembler2.af[0]['bilinear']['row_dofs']
A2 = sp.coo_matrix((vals, (rows, cols)))
A2 = A2.todense()
# linear
b2 = assembler1.af[0]['linear']['vals']
# number of dofs 
n = len(dofs)
# Print
print('Mesh 2')
print('A2 = \n', 6*A2)
print('b2 = \n', 6*b2)
print('n_dofs=', n)
Explanation: It looks like the following dofs from mesh1 and mesh2 are equivalent 
mesh1 -> mesh2
- 0 -> 0
- 1 -> 1
- 4 -> 4
- 5 -> 5
- 8 -> 8
If we restrict to these, we should get the same matrix.
End of explanation
print(A1[np.ix_([0,1,4,5,8],[0,1,4,5,8])] - A2[np.ix_([0,1,4,5,8],[0,1,4,5,8])])
# System for mesh1
system1 = LinearSystem(assembler1)
# Check that it's the same as before
assert np.allclose(A1, system1.A().todense())
Explanation: Check that A1 and A2 coincide when restricting to the nodes
End of explanation
# Mark Dirichlet Regions
f_left = lambda x,dummy: np.abs(x)<1e-9
f_right = lambda x,dummy: np.abs(x-1)<1e-9
# Mesh 1
mesh1.mark_region('left', f_left, on_boundary=True)
mesh1.mark_region('right', f_right, on_boundary=True)
# Mesh 2
mesh2.mark_region('left', f_left, on_boundary=True)
mesh2.mark_region('right', f_right, on_boundary=True)
#
# Check that we get the correct vertices back
#
for side in ['left', 'right']:
    # mesh1
    print('mesh1: ', side)
    for v in mesh1.get_region(side, entity_type='vertex', \
                          on_boundary=True, return_cells=False):
        print(v.coordinates())
 
    print('')
    
    # mesh2
    print('mesh2: ', side)
    for v in mesh2.get_region(side, entity_type='vertex', \
                              on_boundary=True, return_cells=False):
        print(v.coordinates())
        
        
    print('\n\n')
Explanation: Mark Dirichlet Regions on Meshes
End of explanation
#
# Extract Dirichlet conditions (uncompressed format)
#
system1a = LinearSystem(assembler1, compressed=False)
print('System matrix and vector before left Dirichlet nodes')
print('6A = \n', 6*system1a.A().todense())
print('6b = \n', 6*system1a.b() )
print('Extracting Dirichlet nodes on left')
system1a.extract_dirichlet_nodes('left', 0)
print('')
print('6A = \n', 6*system1a.A().todense())
print('6b = \n', 6*system1a.b() )
print('\n\n')
print('Extracting Dirichlet nodes on right')
system1a.extract_dirichlet_nodes('right',1)
print('')
print('6A = \n', 6*system1a.A().todense())
print('6b = \n', 6*system1a.b() )
#
# Extract Dirichlet conditions (compressed format)
#
system1b = LinearSystem(assembler1, compressed=True)
print('System matrix and vector before left Dirichlet nodes')
print('6A = \n', 6*system1b.A().todense())
print('6b = \n', 6*system1b.b() )
print('Extracting Dirichlet nodes on left')
system1b.extract_dirichlet_nodes('left', 0)
print('')
print('6A = \n', 6*system1b.A().todense())
print('6b = \n', 6*system1b.b() )
print('\n\n')
print('Extracting Dirichlet nodes on right')
system1b.extract_dirichlet_nodes('right',1)
print('')
print('6A = \n', 6*system1b.A().todense())
print('6b = \n', 6*system1b.b() )
#
# Check solutions 
# 
system1a.solve()
u1a = system1a.sol(as_function=True)
plot = Plot()
plot.wire(u1a)
# 
Explanation: Now extract Dirichlet nodes
End of explanation
<END_TASK>