diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H
index e58af394a7a..d855f74f417 100644
--- a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H
+++ b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H
@@ -136,6 +136,15 @@ public:
         std::array<amrex::Real, 3> beta
     ) const;
 
+    /** Compute the potential `phi` by solving the Poisson equation with the
+       simulation specific boundary conditions and boundary values, then add the
+       E field due to that `phi` to `Efield_fp`.
+    * \param[in] Efield Efield updated to include potential gradient from boundary condition
+    */
+    void AddBoundaryField (
+        ablastr::fields::MultiLevelVectorField& Efield
+    );    
+
     /** Maximum levels for the electrostatic solver grid */
     int num_levels;
 
diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
index 0b1dca675be..946703ef550 100644
--- a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
+++ b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
@@ -541,3 +541,39 @@ void ElectrostaticSolver::computeB (
         }
     }
 }
+
+void ElectrostaticSolver::AddBoundaryField (ablastr::fields::MultiLevelVectorField& Efield_fp)
+{
+    WARPX_PROFILE("RelativisticExplicitES::AddBoundaryField");
+
+    auto & warpx = WarpX::GetInstance();
+
+    // Allocate fields for charge and potential
+    Vector<std::unique_ptr<MultiFab>> rho(num_levels);
+    Vector<std::unique_ptr<MultiFab>> phi(num_levels);
+    // Use number of guard cells used for local deposition of rho
+    const amrex::IntVect ng = warpx.get_ng_depos_rho();
+    for (int lev = 0; lev < num_levels; lev++) {
+        BoxArray nba = warpx.boxArray(lev);
+        nba.surroundingNodes();
+        rho[lev] = std::make_unique<amrex::MultiFab>(nba, warpx.DistributionMap(lev), 1, ng);
+        rho[lev]->setVal(0.);
+        phi[lev] = std::make_unique<amrex::MultiFab>(nba, warpx.DistributionMap(lev), 1, 1);
+        phi[lev]->setVal(0.);
+    }
+
+    // Set the boundary potentials appropriately
+    setPhiBC( amrex::GetVecOfPtrs(phi), warpx.gett_new(0));
+
+    // beta is zero for boundaries
+    const std::array<Real, 3> beta = {0._rt};
+
+    // Compute the potential phi, by solving the Poisson equation
+    computePhi( amrex::GetVecOfPtrs(rho), amrex::GetVecOfPtrs(phi),
+                beta, self_fields_required_precision,
+                self_fields_absolute_tolerance, self_fields_max_iters,
+                self_fields_verbosity );
+
+    // Compute the corresponding electric field, from the potential phi.
+    computeE( Efield_fp, amrex::GetVecOfPtrs(phi), beta );
+}
diff --git a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H
index cf831a7ab10..0ee7d3c35bd 100644
--- a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H
+++ b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H
@@ -65,14 +65,6 @@ public:
         ablastr::fields::MultiLevelVectorField& Bfield_fp
     );
 
-    /** Compute the potential `phi` by solving the Poisson equation with the
-       simulation specific boundary conditions and boundary values, then add the
-       E field due to that `phi` to `Efield_fp`.
-    * \param[in] Efield Efield updated to include potential gradient from boundary condition
-    */
-    void AddBoundaryField (
-        ablastr::fields::MultiLevelVectorField& Efield
-    );
 };
 
 #endif // WARPX_RELATIVISTICEXPLICITES_H_
diff --git a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp
index 69647da1702..dff9f729f62 100644
--- a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp
+++ b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp
@@ -136,40 +136,4 @@ void RelativisticExplicitES::AddSpaceChargeField (
     computeE( Efield_fp, amrex::GetVecOfPtrs(phi), beta );
     computeB( Bfield_fp, amrex::GetVecOfPtrs(phi), beta );
 
-}
-
-void RelativisticExplicitES::AddBoundaryField (ablastr::fields::MultiLevelVectorField& Efield_fp)
-{
-    WARPX_PROFILE("RelativisticExplicitES::AddBoundaryField");
-
-    auto & warpx = WarpX::GetInstance();
-
-    // Allocate fields for charge and potential
-    Vector<std::unique_ptr<MultiFab>> rho(num_levels);
-    Vector<std::unique_ptr<MultiFab>> phi(num_levels);
-    // Use number of guard cells used for local deposition of rho
-    const amrex::IntVect ng = warpx.get_ng_depos_rho();
-    for (int lev = 0; lev < num_levels; lev++) {
-        BoxArray nba = warpx.boxArray(lev);
-        nba.surroundingNodes();
-        rho[lev] = std::make_unique<amrex::MultiFab>(nba, warpx.DistributionMap(lev), 1, ng);
-        rho[lev]->setVal(0.);
-        phi[lev] = std::make_unique<amrex::MultiFab>(nba, warpx.DistributionMap(lev), 1, 1);
-        phi[lev]->setVal(0.);
-    }
-
-    // Set the boundary potentials appropriately
-    setPhiBC( amrex::GetVecOfPtrs(phi), warpx.gett_new(0));
-
-    // beta is zero for boundaries
-    const std::array<Real, 3> beta = {0._rt};
-
-    // Compute the potential phi, by solving the Poisson equation
-    computePhi( amrex::GetVecOfPtrs(rho), amrex::GetVecOfPtrs(phi),
-                beta, self_fields_required_precision,
-                self_fields_absolute_tolerance, self_fields_max_iters,
-                self_fields_verbosity );
-
-    // Compute the corresponding electric field, from the potential phi.
-    computeE( Efield_fp, amrex::GetVecOfPtrs(phi), beta );
-}
+}
\ No newline at end of file
diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp
index 01ab2d3e48f..d71efd246d7 100644
--- a/Source/Python/WarpX.cpp
+++ b/Source/Python/WarpX.cpp
@@ -24,6 +24,7 @@
 #       include <FieldSolver/SpectralSolver/SpectralSolver.H>
 #   endif // RZ ifdef
 #endif // use PSATD ifdef
+#include <FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H>
 #include <FieldSolver/WarpX_FDTD.H>
 #include <Filter/NCIGodfreyFilter.H>
 #include <Initialization/ExternalField.H>
@@ -274,6 +275,13 @@ The physical fields in WarpX have the following naming:
             [] (WarpX& wx) { wx.Synchronize(); },
             "Synchronize particle velocities and positions."
         )
+        .def("add_boundary_electrostatic_field",
+            [] (WarpX& wx) { 
+                auto Efield_fp = wx.m_fields.get_mr_levels_alldirs("Efield_fp", wx.maxLevel());
+                wx.GetElectrostaticSolver().AddBoundaryField( Efield_fp );
+            },
+            "Compute the electric field due to the potential specified on the domain boundaries and embedded boundaries."
+        )
     ;
 
     py::class_<warpx::Config>(m, "Config")