#705 fix sign error

CHANGELOG.md

@@ -21,6 +21,7 @@
 
 ## Bug fixes
 
+-   Correct a sign error in Dirichlet boundary conditions in the Finite Element Method ([#706](https://github.com/pybamm-team/PyBaMM/pull/706))
 -   Pass the correct dimensional temperature to open circuit potential ([#702](https://github.com/pybamm-team/PyBaMM/pull/702))
 -   Adds missing temperature dependence in electrolyte and interface submodels ([#698](https://github.com/pybamm-team/PyBaMM/pull/698))
 -   Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))

pybamm/spatial_methods/scikit_finite_element.py

@@ -121,7 +121,7 @@ def unit_bc_load_form(v, dv, w):
             # set Dirichlet value at facets corresponding to tab
             neg_bc_load = np.zeros(mesh.npts)
             neg_bc_load[mesh.negative_tab_dofs] = 1
-            boundary_load = boundary_load - neg_bc_value * pybamm.Vector(neg_bc_load)
+            boundary_load = boundary_load + neg_bc_value * pybamm.Vector(neg_bc_load)
         else:
             raise ValueError(
                 "boundary condition must be Dirichlet or Neumann, not '{}'".format(
@@ -138,7 +138,7 @@ def unit_bc_load_form(v, dv, w):
             # set Dirichlet value at facets corresponding to tab
             pos_bc_load = np.zeros(mesh.npts)
             pos_bc_load[mesh.positive_tab_dofs] = 1
-            boundary_load = boundary_load - pos_bc_value * pybamm.Vector(pos_bc_load)
+            boundary_load = boundary_load + pos_bc_value * pybamm.Vector(pos_bc_load)
         else:
             raise ValueError(
                 "boundary condition must be Dirichlet or Neumann, not '{}'".format(

tests/unit/test_spatial_methods/test_scikit_finite_element.py

@@ -461,6 +461,43 @@ def test_pure_neumann_poisson(self):
         u_exact = z ** 2 / 2 - 1 / 6
         np.testing.assert_array_almost_equal(solution.y[:-1], u_exact, decimal=1)
 
+    def test_dirichlet_bcs(self):
+        # manufactured solution u = a*z^2 + b*z + c
+        model = pybamm.BaseModel()
+        a = 3
+        b = 4
+        c = 5
+        u = pybamm.Variable("variable", domain="current collector")
+        model.algebraic = {u: -pybamm.laplacian(u) + pybamm.source(2 * a, u)}
+        # set boundary conditions ("negative tab" = bottom of unit square,
+        # "positive tab" = top of unit square, elsewhere normal derivative is zero)
+        model.boundary_conditions = {
+            u: {
+                "negative tab": (pybamm.Scalar(c), "Dirichlet"),
+                "positive tab": (pybamm.Scalar(a + b + c), "Dirichlet"),
+            }
+        }
+        # bad initial guess (on purpose)
+        model.initial_conditions = {u: pybamm.Scalar(1)}
+        model.variables = {"u": u}
+        # create discretisation
+        mesh = get_unit_2p1D_mesh_for_testing(ypts=8, zpts=32)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "current collector": pybamm.ScikitFiniteElement,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        # solve model
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model)
+
+        # indepdent of y, so just check values for one y
+        z = mesh["current collector"][0].edges["z"][:, np.newaxis]
+        u_exact = a * z ** 2 + b * z + c
+        np.testing.assert_array_almost_equal(solution.y[0 : len(z)], u_exact)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")