From 14981b075db4adc8bc980c15b5f5e33c5fbeb99f Mon Sep 17 00:00:00 2001
From: KulaginVladimir <kulagin.vladimir.l@yandex.ru>
Date: Fri, 16 Feb 2024 23:33:38 +0300
Subject: [PATCH 1/2] upd docs

---
 docs/source/api/temperature.rst                 |  4 ++++
 docs/source/theory.rst                          | 15 ++++++++++++---
 docs/source/userguide/boundary_conditions.rst   |  2 ++
 festim/boundary_conditions/dirichlets/dc_imp.py |  2 +-
 4 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/docs/source/api/temperature.rst b/docs/source/api/temperature.rst
index 5143338c1..bf5e3fc2c 100644
--- a/docs/source/api/temperature.rst
+++ b/docs/source/api/temperature.rst
@@ -7,6 +7,10 @@ Temperature
     :members:
     :show-inheritance:
 
+.. autoclass:: TemperatureFromXDMF
+    :members:
+    :show-inheritance:
+
 .. autoclass:: HeatTransferProblem
     :members:
     :show-inheritance:
\ No newline at end of file
diff --git a/docs/source/theory.rst b/docs/source/theory.rst
index 2f64c3ad7..5a738528d 100644
--- a/docs/source/theory.rst
+++ b/docs/source/theory.rst
@@ -52,7 +52,7 @@ Continuity of local partial pressure :math:`P` at interfaces between materials h
 .. math::
     :label: eq_Sievert   
 
-    P = \frac{c_\mathrm{m}^2}{K_S^2}
+    P = \left(\frac{c_\mathrm{m}}{K_S}\right)^2
 
 where :math:`K_S` is the material solubility (or Sivert's constant).
 
@@ -210,10 +210,19 @@ Assuming second order recombination, :math:`\varphi_{\mathrm{recomb}}` can also
 By substituting Equation :eq:`eq_flux_balance_approx2` into :eq:`eq_flux_balance_approx1` one can obtain:
 
 .. math::
-    :label: eq_DirichletBC_triangle_full
+    :label: eq_DirichletBC_triangle_recomb
     
     c_\mathrm{m} = \frac{\varphi_{\mathrm{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\mathrm{imp}}}{K_r}}
 
+Similarly, dissociation can be accounted for:
+
+.. math::
+    :label: eq_DirichletBC_triangle_full
+    
+    c_\mathrm{m} = \frac{\varphi_{\mathrm{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\mathrm{imp}}+K_d P}{K_r}}
+
+where :math:`K_d` is the dissociation coefficient. 
+
 When recombination is fast (i.e. :math:`K_r\rightarrow\infty`), Equation :eq:`eq_DirichletBC_triangle_full` can be reduced to:
 
 .. math::
@@ -252,7 +261,7 @@ Recombination and dissociation fluxes can also be applied:
     J \cdot \mathrm{\textbf{n}} = -D\nabla c_\mathrm{m} \cdot \mathrm{\textbf{n}}
     = K_d P - K_r c_\mathrm{m}^{\{1,2\}} ~\text{on}~\delta\Omega
 
-where :math:`K_d` is the dissociation coefficient and :math:`K_r` is the recombination coefficient. In Equation :eq:`eq_NeumannBC_DisRec`, the exponent of :math:`c_\mathrm{m}` is either 1 or 2 depending on the reaction order. 
+In Equation :eq:`eq_NeumannBC_DisRec`, the exponent of :math:`c_\mathrm{m}` is either 1 or 2 depending on the reaction order. 
 These boundary conditions are Robin boundary conditions since the gradient is imposed as a function of the solution. 
 
 Finally, convective heat fluxes can be applied to boundaries:
diff --git a/docs/source/userguide/boundary_conditions.rst b/docs/source/userguide/boundary_conditions.rst
index cde7f9942..355148a0d 100644
--- a/docs/source/userguide/boundary_conditions.rst
+++ b/docs/source/userguide/boundary_conditions.rst
@@ -141,6 +141,8 @@ Plasma implantation can be approximated by a Dirichlet boundary condition with t
     # non-instantaneous recombination
     my_bc = ImplantationDirichlet(surfaces=3, phi=1e10 + t, R_p=1e-9, D_0=1, E_D=0.1, Kr_0=2, E_Kr=0.2)
 
+    # non-instantaneous recombination and dissociation
+    my_bc = ImplantationDirichlet(surfaces=3, phi=1e10 + t, R_p=1e-9, D_0=1, E_D=0.1, Kr_0=2, E_Kr=0.2, Kd_0=3, E_Kd=0.3, P=4)
 
 -----------------
 Heat transfer BCs
diff --git a/festim/boundary_conditions/dirichlets/dc_imp.py b/festim/boundary_conditions/dirichlets/dc_imp.py
index b82e85315..f53b43641 100644
--- a/festim/boundary_conditions/dirichlets/dc_imp.py
+++ b/festim/boundary_conditions/dirichlets/dc_imp.py
@@ -23,7 +23,7 @@ class ImplantationDirichlet(DirichletBC):
     The details of the approximation can be found in
     https://www.nature.com/articles/s41598-020-74844-w
 
-    c = phi*R_p/D + (phi+Kd*P/Kr)**0.5
+    c = phi*R_p/D + ((phi+Kd*P)/Kr)**0.5
 
     Args:
         surfaces (list or int): the surfaces of the BC

From 092032902a25f33ba40296f12b54cff5496a6204 Mon Sep 17 00:00:00 2001
From: KulaginVladimir <kulagin.vladimir.l@yandex.ru>
Date: Sat, 17 Feb 2024 17:02:07 +0300
Subject: [PATCH 2/2] upd theory

---
 docs/source/theory.rst | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/source/theory.rst b/docs/source/theory.rst
index 5a738528d..251a9f646 100644
--- a/docs/source/theory.rst
+++ b/docs/source/theory.rst
@@ -94,8 +94,9 @@ In FESTIM, the conservation of chemical potential is obtained by a change of var
 
     \theta = 
     \begin{cases}
-    \frac{c_\mathrm{m}^2}{K_S^2} & \text{in Sievert materials} \\
-    \frac{c_\mathrm{m}}{K_H}     & \text{in Henry materials}
+    \left(\dfrac{c_\mathrm{m}}{K_S}\right)^2 & \text{in Sievert materials} \\
+    \\
+    \dfrac{c_\mathrm{m}}{K_H}     & \text{in Henry materials}
     \end{cases}
 
 The variable :math:`\theta` is continuous at interfaces.