From 1804fa37608e340a9599db2aacb9f1412b421245 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 9 Jan 2019 09:22:54 +0100
Subject: [PATCH] Changes necessary for n-phase model

---
 phasefield/analytical.py | 11 +++++++----
 phasefield/kerneleqs.py  |  5 +++--
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index c38469dd..1318b341 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -303,10 +303,13 @@ def extract_gamma(free_energy, order_parameters):
     return result
 
 
-def pressure_tensor_bulk_component(free_energy, order_parameters):
+def pressure_tensor_bulk_component(free_energy, order_parameters, bulk_chemical_potential):
     """Diagonal component of pressure tensor in bulk"""
     bulk_free_energy, _ = separate_into_bulk_and_interface(free_energy)
-    mu_bulk = chemical_potentials_from_free_energy(bulk_free_energy, order_parameters)
+    if bulk_chemical_potential is None:
+        mu_bulk = chemical_potentials_from_free_energy(bulk_free_energy, order_parameters)
+    else:
+        mu_bulk = bulk_chemical_potential
     return sum(c_i * mu_i for c_i, mu_i in zip(order_parameters, mu_bulk)) - bulk_free_energy
 
 
@@ -341,14 +344,14 @@ def pressure_tensor_interface_component_new(free_energy, order_parameters, dim,
     return result
 
 
-def pressure_tensor_from_free_energy(free_energy, order_parameters, dim,
+def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, bulk_chemical_potential=None,
                                      include_bulk=True, include_interface=True):
     op = order_parameters
 
     def get_entry(i, j):
         p_if = pressure_tensor_interface_component(free_energy, op, dim, i, j) if include_interface else 0
         if include_bulk:
-            p_b = pressure_tensor_bulk_component(free_energy, op) if i == j else 0
+            p_b = pressure_tensor_bulk_component(free_energy, op, bulk_chemical_potential) if i == j else 0
         else:
             p_b = 0
         return sp.expand(p_if + p_b)
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index 75f421d9..87bc9cf0 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -29,10 +29,11 @@ def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1, discretization
 
 
 def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field,
-                           dx=1, discretization='standard'):
+                           dx=1, discretization='standard', bulk_chemical_potential=None):
     dim = phi_field.spatial_dimensions
 
-    p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim)
+    p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim,
+                                         bulk_chemical_potential=bulk_chemical_potential)
 
     p = p.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
     index_map = symmetric_tensor_linearization(dim)
-- 
GitLab