From 2a4518d487dc6200f94c65b149e37aceb9e909e4 Mon Sep 17 00:00:00 2001
From: JanHab <Jan.Habscheid@web.de>
Date: Fri, 16 Aug 2024 17:31:33 +0200
Subject: [PATCH] Added ParameterAnalysis + InstructiveExample
---
...nsionalEquilibrium-AnInstructiveExample.py | 129 ++++++++++++++++++
examples/ParameterAnalysis/Lambda.py | 102 ++++++++++++++
examples/ParameterAnalysis/ScalingPhi.py | 72 ++++++++++
examples/ParameterAnalysis/a.py | 68 +++++++++
src/Eq04.py | 44 +++++-
src/__pycache__/Eq04.cpython-312.pyc | Bin 17910 -> 20310 bytes
src/oneD_refined_mesh.py | 50 -------
7 files changed, 412 insertions(+), 53 deletions(-)
create mode 100644 examples/OneDimensionalEquilibrium-AnInstructiveExample.py
create mode 100644 examples/ParameterAnalysis/Lambda.py
create mode 100644 examples/ParameterAnalysis/ScalingPhi.py
create mode 100644 examples/ParameterAnalysis/a.py
delete mode 100644 src/oneD_refined_mesh.py
diff --git a/examples/OneDimensionalEquilibrium-AnInstructiveExample.py b/examples/OneDimensionalEquilibrium-AnInstructiveExample.py
new file mode 100644
index 0000000..aebf81f
--- /dev/null
+++ b/examples/OneDimensionalEquilibrium-AnInstructiveExample.py
@@ -0,0 +1,129 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+import numpy as np
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+
+# Define the parameters and boundary conditions
+Lambda2 = 8.553e-6
+z_A = -1
+z_C = 1
+y_L = 1/3
+
+# Define the equilibrium concentrations
+def y_C(phi, phi_prime):
+ return y_L * np.exp(-z_C * phi - 1/2 * Lambda2 * phi_prime**2)
+
+def y_A(phi, phi_prime):
+ return y_L * np.exp(-z_A * phi - 1/2 * Lambda2 * phi_prime**2)
+
+# Define the ODEs as a system of first-order ODEs
+def func(PHI, x):
+ phi, phi_prime = PHI
+ f = -1/Lambda2 * (z_C * y_C(phi, phi_prime) + z_A * y_A(phi, phi_prime))
+ return [phi_prime, f]
+
+# Initial conditions
+Y0 = [0, 1e-9]
+x_end = np.array([-0.95, -1.0, -1.05, -1.1])*1e-1
+
+# Storage for the results
+y_C_end, y_A_end, y_S_end = [], [], []
+phi_end, x_vec = [], []
+
+# Solve the ODEs for different values of x_end
+for x_ in x_end:
+ x = np.linspace(x_, 0, 1024)
+
+ Y = odeint(func, Y0, x)
+
+ Y_mirror = np.flip(Y, axis=0)
+ phi_end.append(Y_mirror[:, 0])
+ y_C_end.append(y_C(Y_mirror[:, 0], Y_mirror[:, 1]))
+ y_A_end.append(y_A(Y_mirror[:, 0], Y_mirror[:, 1]))
+ y_S_end.append(1-y_C_end[-1]-y_A_end[-1])
+ x_vec.append(x)
+
+# Visualize the results
+# create 2x1 subplots
+fig, axs = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(14-14/3, 12), constrained_layout=True)
+
+# clear subplots
+for ax in axs:
+ ax.remove()
+
+# add subfigure per subplot
+gridspec = axs[0].get_subplotspec().get_gridspec()
+subfigs = [fig.add_subfigure(gs) for gs in gridspec]
+
+for row, subfig in enumerate(subfigs):
+ subfig.suptitle(f'$L_x = $ {abs(round(x_end[row], 4))}', fontsize=22)
+
+ # create 1x2 subplots per subfig
+ axs = subfig.subplots(nrows=1, ncols=2)
+
+ for col, ax in enumerate(axs):
+ labelsize = 20
+ ax.set_ylim(-0.05, 1.05)
+ match col:
+ case 0:
+ ax.plot(x_vec[row], phi_end[row])
+ ax.set_ylim(-0.5, phi_end[-1][0]+0.5)
+ ax.set_ylabel('$\\varphi$', fontsize=labelsize)
+ case 1:
+ ax.plot(x_vec[row], y_C_end[row], label='$y_C$')
+ ax.plot(x_vec[row], y_A_end[row], label='$y_A$')
+ ax.plot(x_vec[row], y_S_end[row], label='$y_S$')
+ ax.legend()
+ ax.set_ylabel('$y_\\alpha$', fontsize=labelsize)
+ ax.set_xlim(min(x_end), 0)
+ ax.set_xlabel('x [-]', fontsize=labelsize)
+ ax.grid(True)
+ ax.tick_params(axis='x', labelrotation=45)
+fig.show()
+
+
+
+
+# Streamplot
+phi_range = np.arange(-0.1, 0.1, 0.01) # Range for x-axis
+dphi_range = np.linspace(-5e-3, 5e-3, 10) # Range for y-axis
+
+# Create a meshgrid with different sizes for x and y axes
+phi_values, dphi_values = np.meshgrid(phi_range, dphi_range)
+initial_conditions = np.array([phi_values.ravel(), dphi_values.ravel()]).T
+
+# Points at which the solution is requested
+x_domain = np.linspace(0, 1.0e-1, 1024)
+
+# Prepare storage for the results
+results = np.zeros((len(initial_conditions), len(x_domain), 2))
+
+# Numerical integration for each initial condition
+for i, Y0 in enumerate(initial_conditions):
+ Y = odeint(func, Y0, x_domain)
+ results[i, :, :] = Y
+
+# Reshape results for plotting
+Y_reshaped = results[:, :, 0].reshape(len(dphi_range), len(phi_range), len(x_domain))
+
+# Select a specific slice for plotting, e.g., the last time point -> last point of phi
+Y_slice = Y_reshaped[:, :, -1]
+
+# Compute derivatives for the streamplot
+U, V = np.gradient(Y_slice, dphi_range, phi_range)
+
+# Plotting
+plt.subplots(tight_layout=True)
+# lw = 5*Y_slice/Y_slice.max()
+strm = plt.streamplot(phi_values, dphi_values, U, V, color=U, linewidth=2, cmap='autumn', density=10)
+plt.colorbar(strm.lines)
+plt.xlim(-0.01, 0.01)
+plt.ylim(dphi_values.min(), dphi_values.max())
+plt.xlabel('$\\varphi$')
+plt.xticks(rotation='vertical')
+plt.ylabel("$\phi$")
+plt.show()
\ No newline at end of file
diff --git a/examples/ParameterAnalysis/Lambda.py b/examples/ParameterAnalysis/Lambda.py
new file mode 100644
index 0000000..b5fc5f3
--- /dev/null
+++ b/examples/ParameterAnalysis/Lambda.py
@@ -0,0 +1,102 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+a2 = 7.5412e-4
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+# Define left values of electric potential and values for lambda^2
+phi_left_vec = [-1.0, 1.0, 4.0, 8.0]
+Lambda2_vec = [8.553e-5, 8.553e-6, 8.553e-7]
+y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
+
+# Solve the system for different values of phi_left
+for phi_left_ in phi_left_vec:
+ y_A_bcs, y_C_bcs, y_S_bcs, phi_bcs, p_bcs, x_bcs = [], [], [], [], [], []
+ for l2 in Lambda2_vec:
+ y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, l2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
+ y_A_bcs.append(y_A_)
+ y_C_bcs.append(y_C_)
+ y_S_bcs.append(1 - y_A_ - y_C_)
+ phi_bcs.append(phi_)
+ p_bcs.append(p_)
+ x_bcs.append(x_)
+ y_A.append(y_A_bcs)
+ y_C.append(y_C_bcs)
+ y_S.append(y_S_bcs)
+ phi.append(phi_bcs)
+ p.append(p_bcs)
+ x.append(x_bcs)
+
+# Visualize the results
+fig, axs = plt.subplots(2, 2)
+xlim = 0.08
+markers = ['--', '-', ':']
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+for i in range(len(phi_left_vec)):
+ # Plot phi
+
+ match i:
+ case 0:
+ i_ = 0
+ j_ = 0
+ case 1:
+ i_ = 0
+ j_ = 1
+ case 2:
+ i_ = 1
+ j_ = 0
+ case 3:
+ i_ = 1
+ j_ = 1
+ axs[i_,j_].set_title(f'$\delta \\varphi$ = {phi_left_vec[i]}')#, fontsize=labelsize)
+ for j in range(len(Lambda2_vec)):
+ clr = colors[j]
+ axs[i_,j_].plot(x[i][j], y_A[i][j], markers[0], color=clr)#, lw=lw)
+ axs[i_,j_].plot(x[i][j], y_C[i][j], markers[1], color=clr)#, lw=lw)
+ axs[i_,j_].plot(x[i][j], y_S[i][j], markers[2], color=clr)#, lw=lw)
+ axs[i_,j_].grid()
+ axs[i_,j_].set_xlim(0,xlim)
+ axs[i_,j_].set_xlabel('$x$ [-]')#, fontsize=labelsize)
+ axs[i_,j_].set_ylabel('$y_\\alpha$ [-]')#, fontsize=labelsize)
+ axs[i_,j_].tick_params(axis='both')#, labelsize=labelsize)
+
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle='--', label='$y_A$')
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle='-', label='$y_C$')
+dummy, = axs[i_,j_].plot(2, y_A_R, color='grey', linestyle=':', label='$y_S$')
+for index, Lambda2_cur in enumerate(Lambda2_vec):
+ dummy, = axs[i_,j_].plot(10, y_A_R, color=colors[index], linestyle='-', label=f'$\lambda^2$ = {Lambda2_cur}')
+
+lines_labels = [ax.get_legend_handles_labels() for ax in fig.axes]
+lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]
+
+# Finally, the legend (that maybe you'll customize differently)
+lgnd = fig.legend(lines, labels, bbox_to_anchor=(0.88,1.15), ncol=3)
+fig.tight_layout()
+fig.show()
\ No newline at end of file
diff --git a/examples/ParameterAnalysis/ScalingPhi.py b/examples/ParameterAnalysis/ScalingPhi.py
new file mode 100644
index 0000000..95da151
--- /dev/null
+++ b/examples/ParameterAnalysis/ScalingPhi.py
@@ -0,0 +1,72 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+phi_left_vector = [1,2,4,8]
+
+phi, x = [], []
+# Solve the system for different values of phi_left
+for phi_L in phi_left_vector:
+ y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_L, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', rtol=rtol, max_iter=max_iter)
+ phi.append(phi_)
+ x.append(x_)
+
+# Normalize the electric potentials
+phi_scaled = []
+for i, phi_L in enumerate(phi_left_vector):
+ phi_scaled.append(phi[i]/phi_L)
+
+# Visualize the results
+# Non-normalized electric potential
+for i, phi_L in enumerate(phi_left_vector):
+ plt.plot(x[i], phi[i], label='$\delta \\varphi$ = ' + str(phi_L))
+plt.legend()
+plt.xlabel('$x$ [-]')
+plt.ylabel('$\\varphi$ [-]')
+plt.ylim(0, max(phi_left_vector) + 0.1)
+plt.xlim(0,0.05)
+plt.grid()
+plt.tight_layout()
+plt.show()
+
+# Normalized electric potential
+for i, phi_L in enumerate(phi_left_vector):
+ plt.plot(x[i], phi_scaled[i], label='$\delta \\varphi$ = ' + str(phi_L))
+plt.legend()
+plt.xlabel('$x$ [-]')
+plt.ylabel('$\\varphi/\delta \\varphi$ [-]')
+plt.ylim(0, 1.05)
+plt.xlim(0,0.05)
+plt.grid()
+plt.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/examples/ParameterAnalysis/a.py b/examples/ParameterAnalysis/a.py
new file mode 100644
index 0000000..042e3b4
--- /dev/null
+++ b/examples/ParameterAnalysis/a.py
@@ -0,0 +1,68 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+sys.path.insert(0, src_path)
+
+from Eq04 import solve_System_4eq
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Define Parameter and buondary conditions
+phi_right = 0.0
+p_right = 0.0
+y_A_R, y_C_R = 1/3, 1/3
+z_A, z_C = -1.0, 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+number_cells = 1024*4
+refinement_style = 'hard_log'
+rtol = 1e-8
+max_iter = 10_000
+
+phi_left_vec = [0.1, 1.0, 4.0, 8.0]
+a2_vec = [7.5412e-3, 7.5412e-4, 7.5412e-5]
+y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
+for phi_left_ in phi_left_vec:
+ y_A_bcs, y_C_bcs, y_S_bcs, phi_bcs, p_bcs, x_bcs = [], [], [], [], [], []
+ for a2 in a2_vec:
+ y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=0.05, x0=0, x1=1, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
+ y_A_bcs.append(y_A_)
+ y_C_bcs.append(y_C_)
+ y_S_bcs.append(1 - y_A_ - y_C_)
+ phi_bcs.append(phi_)
+ p_bcs.append(p_)
+ x_bcs.append(x_)
+ y_A.append(y_A_bcs)
+ y_C.append(y_C_bcs)
+ y_S.append(y_S_bcs)
+ phi.append(phi_bcs)
+ p.append(p_bcs)
+ x.append(x_bcs)
+
+p = np.array(p)
+
+plt.figure()
+for j in range(len(a2_vec)):
+ plt.plot(x[-1][j], p[-1][j], label='$a^2$ = {}'.format(a2_vec[j]))
+plt.grid()
+plt.ylim(0.0, np.max(p[-1][0]))
+plt.xlim(0,0.08)
+plt.xlabel('$x$ [-]')
+plt.ylabel('$p$ [-]')
+plt.tick_params(axis='both')
+plt.legend(loc='best')
+plt.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/src/Eq04.py b/src/Eq04.py
index bed5462..37c292f 100644
--- a/src/Eq04.py
+++ b/src/Eq04.py
@@ -8,12 +8,48 @@ from mpi4py import MPI
from dolfinx import mesh, fem, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
-from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Cell
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
from basix.ufl import element, mixed_element
-from petsc4py import PETSc
import matplotlib.pyplot as plt
-from oneD_refined_mesh import create_refined_mesh
+# Define mesh
+def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
+ '''
+ Creates a one-dimensional mesh with a refined region at the left boundary
+
+ Parameters
+ ----------
+ refinement_style : str
+ How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
+ number_cells : int
+ Number of cells in the mesh
+
+ Returns
+ -------
+ Mesh
+ One-dimensional mesh, ready for use in FEniCSx
+ '''
+ if refinement_style == 'log':
+ coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
+ elif refinement_style == 'hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
+ coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ elif refinement_style == 'hard_hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
+ coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ num_vertices = len(coordinates_np)
+ num_cells = num_vertices - 1
+ cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
+ gdim = 1
+ shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
+ degree = 1
+ domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
+ coordinates_np_ = []
+ [coordinates_np_.append([coord]) for coord in coordinates_np]
+ msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
+ return msh
def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
'''
@@ -26,6 +62,8 @@ def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
Parameters
----------
phi_left : float
diff --git a/src/__pycache__/Eq04.cpython-312.pyc b/src/__pycache__/Eq04.cpython-312.pyc
index d115d2ed6ed42485613169199be3432cb7a86dc5..e94b79d647eeb82bdf59ecca6b780f96647b96e5 100644
GIT binary patch
delta 4847
zcmbVPeN0=|6@Tyf0~@e`Kp@xz9*_XWfcZ{Z(f}a|NgyE!u%uox$ML;@Q-5)vNpNar
zXUnv7tJY?BTFJUDiB_pZVN%j`TJ}d(_QxjFwH1!L7Cl9p)J>hbbuDS6b*(CG=RVsQ
z5<0be;+}iYyXXAQx#ygFc>dQ9`FDQ78~<o9WODF${jZI@)_l>pfxmLG`g}z26njl}
zGtVJCE<u@iBg(?-N3<yWyxLxhwxe8JhVo#ngOLxT0LFS4g)oY687?`S!EtB<aK$i~
zz6oYC7EuYz8)1~fVAvpwL|B*OvXL58ab9CD#}yG&`9uxJJ<F5V(I$H(+HBthV>7C<
zSHY-8)%F@xW49oSy%yElt;lMxLv{9gRBvxU4R#yKu{R1FF2|L46W)w#umxA)YFvx0
zxDMCj2E-!+G9nvp#J1VS>+4tmIPRE+<LpiLW_yd>)@pA90P`tNVrwPEZMfySka)v$
zij<GuG!)P!K1XFf6`)Oe6J1hiU`^=Ftb7_$8(|$(*54G2y!ppMdF0EWn@}riL+!I$
zu7fx1Nt#&Ok!G|Nq#yZe(#U|D&O)+r^Vv)jWam*(?RZFa<yj_IoaKL>D`(qvex<u`
zLd%z(u=W2A-V16ac%}4p^>xpFER9(NIRksJ7yE*;!RwyH$dO`m2_A9#lr4u840*7e
z?VIqP#Khsk9*-nz2@XyWpDUS2#d$qY6^#3LfIi2CImELnBo*{rb7><IkSfe0)mf0Z
z&cxsW4g!kQXR|nNN_dv{sX$4iC`~b?-ApN%%p}Kwy7nw^sbD3*P{yq8y82;2q-U~I
zDx_Z<O4(iwwU%EAlg#TWb&#i@YB!dt!s-zL3A350R7?rJbKrTF;(2&&b6>~v?4;+}
zYeRii&kYP_rRT@j@I1)UraTcL6^$uPSd+-jls1^GOpZIHbh?yuVqUSEP5F%c_1YB1
zsvWq;5mi{5%zm<HM$+~&kZMYYj8l4;GlMxP8^OHP3gv`h59QA0rBa=;GnF;UO=X9L
zo!m19moUn?gyRK}D?>Q@RDMdE&D4;QQDhk9Q+gMLrwi2}J!RlI`m2o6?BIr!_$d=R
zN+ctk9nMObn=*!tgCNc;DrY#$I_D%rCR2o`hhS$AOQKWs`*1yS!-bOEexK7Ldf}Fd
zXWYRt*h}0f1oJ3tiq4=I9K)gqj|9b&{s|v)l2E1schE_kUL2(NN{Z-rDhgYbm(0!{
z{~1L`L6ye*6CNa<#3|qgaUc)`WblJv`p=3=db4aB4Q<*OVRwBQWMB;!>j1-vLgI*@
zh!YYP-9B;mE}y$+a8gktT8Io`S^Iz$!w0NlWZjsPAcx04I{M}t?D?>Bw}qE=*oOjs
zw=ej&pIu_lZ#yk~L{>W|-I8ojG?vxAxOVT({*!^u(b*C^>iF}=ow5$}Fgfk?$hv;#
zC~^8mv4tl^U=5Hp(wH-VEk;@J1!NuAl>$x|mQ|o5XSw`77o@-kLsogPPtJDvJriCZ
z^e(3>@`S8~kOCUX1uVjf3f=?);45K`1N4;)I#^%{>rHZI&%plujw1ty`*&J$2rDE4
zRmI6hm)}p2n;DWE;7HDDzUjdi&CoX^6JEz@OoDC~mSjD&(PR(RDXo;Ds7Ik_h{z0R
z5gsKN!|eAu-M$C$9I{5SD672En8Yp$D~uJ-=YE^l7qq!-!xAQv?Ud6uM(=O#h$M>&
zlWwTXkmQzZ4Grv1A}-sm@y6DMKuFF{h9Oz0iEwNM8!T+n`;hx6w<L194byFRGQ_*4
zlFQah)+^4LyqW&ndDjlaO&zzLF;nMUQ>=LB4~&cY;%WOGef3@Q#>*X%OC49nXYyuS
zX2Q37Vz$n>wQDXfX5BHr>uvMlg~6eid1&DqM`ETU3-;qN)A0*Ai>8e)u6xOR)jV@(
zX8e`>8^`7<V{N<QO}ppDV@-SF<$LGLW95DGhhpM^g+oIzaVTaQx{!0X*nGM4Qt3?D
zZQDZ4zJ)`F-!2|}ig^;5cHGg|DXx5VPaH9mgLzW+V)IKKS373LZ{)qwaw9y~6Wh`o
zZ`m`S7i-xUuk4%ej#ciTACFZWTo^nWt2i1n9ev1?h~mlx<K6tiXTKfKuZs3B=9Wy4
z$8yE!-n%(@&qSVzyePb+zp9_H#KqQEPX5My(>=$}Ju-jjcOCKefrYYz@ti}^J&Re{
z(VdIgh0$G$24i&hT|?eas;6tO;F&`=(A(=e;)ZR}-HU}Aq63S~UGJ#7XVym@?`q63
zjrp<Z;pyH#X*R{x-FGvM?`E1}nWmq&T;6_Z`^?rGj=7O|=Kg5UVwN$w@3T)doMs0v
zu{-|<vAB1WniSDrYmQR0w)vJriwpP<Z~fl1Sr{l6_l_v-=m0(wgoY;hJy7K$aswU3
zQk@v`Pl#S;NDLAferPdh!i5zLf<4YjCu>$>AkpkFx%HQlb2`3=o@gFyxF#G>+RZhA
zl!F^`hFfxZ{1S!_7JicKVb#dGkr*XgV4~Zqs`xN<Rh=#!0G1pCViCw;w%!Y5S(biU
zWvz*@*F%g9Gx7)!)?H9A6u_hmxDS*@rvgWYtPC&Zgnneg<8gT0Pk<W<@DA8bgTGbX
z!aq;ns43Rh#B;0*Ix87w7zv+_F*8ZEZL8FNTGPhA0BYs@mxw>MO8kVSiNCB^sgCE=
zF6e3>v;yubC+}LNR#n@^U!})tAF`4X?_MRoRNI#IOL#*@7#XIG*2fFBvt=g}V?9$8
zqY!=if^~D%>+D@2kbg(H#e$;g%8PrisIIuK>SG0#g=}yTwp?m@u%Uo{VeQ|^`uR1U
zu-;|KAU>FgABe09jJah#Kz6b>yBJ{>CHNr(#+bQ8wq?~mM>oA^E92jyvvqCSRp|ND
zSl_;lFl7nIN_-MzZzmW@TTqU($^wnnS7*hT)EP#G>C5$7s}w*l@Q^|UVF#?2GY@N-
z_nkFU-f9wj@6)=5Qpjar!@dsotAi+d+Vnva2v4%DbV%84TX9GdpbHJ<{NL#(4K;w}
zXsn~XwiCVSgXvb_r4Pp3T5&ohe!2uQi(evu<TU-rR>6Nn^BVUn*iV05TqSA`<AwmY
zT<Ez6z@=1bvH>{Mf01%0S^bkSAj?ss#^xqL#nbL)7A+7gNeiNtf=n|?7n-ZH)bJmO
zoMZej6<R{;NiAC^I>d;_?+muLlCAW)mKfhh-K`|T?t!I`uxm=#r6cSrrGrko7|RNu
zD1qx?E2TErzxqT8?q;iAMp)NT+DYPTU>{@q80lwZkdY&dv;&cKJ-tniO7>$))*exQ
zRLT4inTKCYOoWwtvK;ADn|1xS7=MD1lZ;?S-1OD9mvg_%mH|diGZLbac0d0l{kVOo
zmUTWkS82MA!H^WhUPmh)C-d-{xga?pLI3#z{$cOD<=t{gr~MuSI7PF&Ol`yr>(Tn2
zl7GC{Zt5i82M%t^m!48NU>|jKoo{4!S818#0)u$;>Bjed6xy|&Jjx&*lQSF+<aasX
zBL%*40wGdP^Shs?uXc+Ol@uhTl4-N=30b40OIDxs`#p-wvd-b~!B?-t@nE{Tn39sD
z<=e7mx)k@*30WQJjKrloOi%A<tz;K6(fnBqKi6&Z`q6|3Zzq2RmWLs=(>Hh6>}5;Z
zO?=Kw{t^e%?G|O4YrW67d;NS!#c!U*OB|bLCKB`QBtO4C!7mIwo?JLbk_(>%3)<aN
z7-5P_9M>lt6cT&D3=8%&;}$$=W>T0+GefGQs%0IIYCO$6sp?WM@ARts(~MJ%(o9$#
fsZ4VNnjy`y64JOe_mxvEP2T#f=QCc-I_!S|8gl3j
delta 2696
zcmb7FeQXp}5Z~RqYma^cZTUFbyDR;?Yw6v!ch{C;d(u*b7Wt}tERdz#du<Q)?#kOM
z&;qu95dMJ}>Prl$A(D^?Dbd9KBY%O3F-9aQiBYmDF(m$BOiYc5QA5<3-P&zI65V8e
zZ|1$3d2eRky#4a3`0E9+;9h=yjsV|G{LX+_KUuI!ygpfTZhVQlI_Z*KqCi9<scVv@
za+$&BsP3dk_5j~b3aFP90xbnv2DBWg2&fpSgnDWDShgUL6(FqyT1D5A)iAnf1z7`Q
zDNq@Zo0by~kQexJ>K&*hYtK2x<w{ybDvs3(!X=S?Mk?iMQYF^_)skx22jnL;aviCa
z>xoZZNBnXFsgnbwUJjCVa)>m@jU-=gk_1{wt7#3brFFEP`lz3-qYX4bgERzz<&r!S
zpiQ*#R&Z)LCP5GmIs`!u%efIbDmTk5;Gm2jwig%E78;$BOdV0c>7yC$wN;2^e8yI8
zD?56MdnFvlc$YDr;W4`hR>9Z~#$o=BRAP!Q5F@08L`dtc<|)_@Khq~_*U}b}OT%MO
z6&gV#>p`{cX;pK3*znc4jh;Tnd&aRA90Ox${*N)7w__}5T*^POPsKW^s`?~7o2oQP
zQyD!!nK(s>GAH^(-LYd+-|jdtMR7*^GgM)8Adv!<B-Ms~5c#gceffE*k>nA|lsGuk
z_|3wmWlK#hbg5{WQI$@AT{uuIzA?SAvP$BAE^$wPU476gmhiVjyZMDcj#rmb!+Lf?
zi>qmk!hyt7tQxYXVP(7%f5`lPUT;A|LqKEkz^3EDaKrE@za6L+_wermC;Sbd#R336
zi8Uf$12AXHnt<0GyAp##nK2LlDB!9y`2DCCL5L!(0KjlbW^|b9&ffk()}I=r{BH0-
z9lF0Mu^<fkl7SIbRn){WIHkGQDwv+%i+6{b#1TFgDi@;=T%X%GT01(QQ(wJEy`V87
zp5$GP&sBFVQh(SO5=VJX(=+ODJR|V3Me0|ZBI0TO#q-tOi`0k0k=(N|!(Kv&@*Cm(
zWowWvMKwEs*E=y`CR9pwkCkFJZK>|-ye;C%!ZgP1{934(XCl2DN<h3RvO*Yq*=z?4
zY$t$j8y-sN;xH@6nF<8lLW36tyrX~EfIo@Uh&Om;G_u%h5kDV|Zs<d!R{$3JWqWXT
zFM{RSIQL@DzeQ_ur%=nNf0X;1Tik}j?dVW4J$Qh(N;=V8QFRw*Fr?omgoni?m-;6=
zFE5*HzPxj?Z6?ocT%owli2a-9Q0!$;WQM~mD{|BIAzZSOng!xNQRB)k05h-7RQJaq
zXdWJmtSGp-^TzS1z6V92nY<8G=Bw6vey8Q|Rx1_G1shqZz(5w{!^pG~=3q<Av{_5G
z4`yC(c1tgQ%I~jt^LJYVaCyIKz4(9K8|i-5{Q<+hN4y8#x0!9Ra=tCtR#>%ieZDOu
zeJk+WZEjPs+*0gbs0che0x-4~?;wA_&6)cM=i&%aKDJ@BRL|+$zHOVbcl)m1F0*&D
z2h?<bCK_hN-1))-(a-f4*f_S0uZLk5Gwe-<9cNXDVMj6S6J{jSi=1)8)*@$I;7a7G
z5b(`19v=ojIDz>P{0MajVT2X{-ImmbbZ6)0kWzt9t)8{lcye`duP%-}Q%U|)dtFHk
zjW;56A#6tI;SC)hETzce#W4!W4)EEIw0MXgj`jJCw{NMIR!>l6_o$YkNhM5=vs<9N
z97t23T7Ny65N2c3k78GxS%c`Mm!I9@Y8`y@3Vydc2G@0#{`ffR>R<_2g`RtEq(+c_
z{`r=3L2MP{O=U-c8(e=e`1rNaP3>$Yx`^r7ibB$HMbRDb>lhwoK|Zndtz73Y&1ms3
z$VGf*Pv!Wsbc%K@`r)u*u-6?%9d(<QVY)2=e=YlwbXql{(w&Nu>Q7QddAg!j9p}66
zUn^=vXVlDP#54bLEqGxX`A0qB3T$xm=YeVd!?_`uCL=0sXZJu9fwXeIsyDF7^{wFf
zUI_hcKWP)IF4I{7#~UN&_~D%VX^AP%^zENx`UmEilm-mH*1KXH4QB<RN7^BoYaqNL
z$(E3iR7*G|owbBM+dkX8jcPk?38!r_`~1pgd#@$*+leKdv5!|+(soCmW8R235{_St
QsAe<M`?vir+VNlTH!RPK`2YX_
diff --git a/src/oneD_refined_mesh.py b/src/oneD_refined_mesh.py
deleted file mode 100644
index d34bca8..0000000
--- a/src/oneD_refined_mesh.py
+++ /dev/null
@@ -1,50 +0,0 @@
-'''
-Jan Habscheid
-Jan.Habscheid@rwth-aachen.de
-'''
-
-import numpy as np
-from mpi4py import MPI
-from dolfinx import mesh
-from basix.ufl import element
-from ufl import Mesh
-
-
-# Define mesh
-def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
- '''
- Creates a one-dimensional mesh with a refined region at the left boundary
-
- Parameters
- ----------
- refinement_style : str
- How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
- number_cells : int
- Number of cells in the mesh
-
- Returns
- -------
- Mesh
- One-dimensional mesh, ready for use in FEniCSx
- '''
- if refinement_style == 'log':
- coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
- elif refinement_style == 'hard_log':
- coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
- coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
- coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
- elif refinement_style == 'hard_hard_log':
- coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
- coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
- coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
- num_vertices = len(coordinates_np)
- num_cells = num_vertices - 1
- cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
- gdim = 1
- shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
- degree = 1
- domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
- coordinates_np_ = []
- [coordinates_np_.append([coord]) for coord in coordinates_np]
- msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
- return msh
\ No newline at end of file
--
GitLab