<jupyter_start><jupyter_text>Table of Contents<jupyter_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]<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>#
# Define the element
#
Q1 = QuadFE(2, 'Q1')<jupyter_output><empty_output><jupyter_text>Since we have already tested the assembly, we focus here on the linear system. In particular:- Marking and extracting Dirichlet boundary conditions- Extracting hanging nodesboth 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.<jupyter_code>mesh1 = QuadMesh(resolution=(2,2))<jupyter_output><empty_output><jupyter_text>To test ```extract_hanging_nodes``` and ```resolve_hanging_nodes``` we construct a simple mesh with hanging_nodes.<jupyter_code>mesh2 = QuadMesh(resolution=(2,2))
mesh2.cells.get_child(2).mark(1)
mesh2.cells.refine(refinement_flag=1)<jupyter_output><empty_output><jupyter_text>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\}$<jupyter_code>#
# 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()<jupyter_output><empty_output><jupyter_text>Let's visualize the meshes.<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>It looks like the following dofs from mesh1 and mesh2 are equivalent __mesh1 -> mesh2__- 0 -> 0- 1 -> 1- 4 -> 4- 5 -> 5- 8 -> 8If we restrict to these, we should get the same matrix.<jupyter_code># 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)<jupyter_output>Mesh 1
A1 = 
 [[ 4. -1. -2. -1.  0.  0.  0.  0.  0.]
 [-1.  8. -2. -2. -1. -2.  0.  0.  0.]
 [-2. -2. 16. -2. -2. -2. -2. -2. -2.]
 [-1. -2. -2.  8.  0.  0. -2. -1.  0.]
 [ 0. -1. -2.  0.  4. -1.  0.  0.  0.]
 [ 0. -2. -2.  0. -1.  8. -2.  0. -1.]
 [ 0.  0. -2. -2.  0. -2.  8. -1. -1.]
 [ 0.  0. -2. -1.  0.  0. -1.  4.  0.]
 [ 0.  0. -2.  0.  0. -1. -1.  0.  4.]]
b1 = 
 [0. 0. 0. 0. 0. 0. 0. 0. 0.]
n_dofs= 9
============================================================
Mesh 2
A2 = 
 [[ 4. -1. -2. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-1.  8. -2. -2. -1. -2.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-2. -2. 16. -1. -2. -2. -1.  0. -2. -1. -2.  0. -1.  0.]
 [-1. -2. -1.  8.  0.  0.  0.  0.  0. -1. -2. -1.  0.  0.]
 [ 0. -1. -2.  0.  4. -1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. -2. -2.  0. -1.  8. -2.  0. -1.  0.  0.  0.  0.  0.]
 [ 0.  0. -1.  0.  0. -2.  8.  0. -1.  0. -2.  0. -1. -1.]
 [ 0.  0.  0.  0.  0.  0.  0.  4.  0.  0. -2. -1.  0. -1.]
 [ 0.  0. -2.  0.  0. -1. -1.  0.  4.  0.[...]<jupyter_text>Check that ```A1``` and ```A2``` coincide when restricting to the nodes<jupyter_code>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())<jupyter_output><empty_output><jupyter_text>Mark Dirichlet Regions on Meshes<jupyter_code># 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')<jupyter_output>mesh1:  left
(0.0, 1.0)
(0.0, 0.0)
(0.0, 0.5)

mesh2:  left
(0.0, 1.0)
(0.0, 0.0)
(0.0, 0.75)
(0.0, 0.5)



mesh1:  right
(1.0, 1.0)
(1.0, 0.0)
(1.0, 0.5)

mesh2:  right
(1.0, 1.0)
(1.0, 0.0)
(1.0, 0.5)<jupyter_text>Now extract Dirichlet nodes<jupyter_code>#
# 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)
#<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19579484686436557, 'fraction_numerical': 0.060646572678658243, 'mean_word_length': 3.4747145187601958, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '638903', 'license': 'MIT', 'n_tokens_mistral': 3671, 'n_tokens_neox': 3484, 'n_words': 871, 'path': 'tests/test_fem/.ipynb_checkpoints/debug-checkpoint.ipynb', 'repo_name': 'hvanwyk/drifter', 'size': '8228'}