Nightly Commit

Started implementation of Magnetostatic solver
Work on Dirichlet boundary conditions (magnetic insulation)

.idea/dictionaries/jpries.xml

@@ -3,6 +3,7 @@
     <words>
       <w>discretization</w>
       <w>eigen</w>
+      <w>linearize</w>
       <w>supremum</w>
       <w>tangency</w>
       <w>verticies</w>

src/Physics/CMakeLists.txt

@@ -5,7 +5,7 @@ set(SOURCE_FILES
 
         #./src/PhysicsCommon.h
 
-        ./src/Condition.h ./src/Condition.cpp
+        ./src/Forcing.h ./src/Forcing.cpp
 
         ./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
 

src/Physics/include/Physics.hpp

 #ifndef OERSTED_PHYSICS_HPP
 #define OERSTED_PHYSICS_HPP
 
-
 //#include "../src/PhysicsCommon.h"
-//#include "../src/Physics.h"
 //#include "../src/Discretization.h"
 //#include "../src/Operator.h"
 
-#include "../src/Condition.h"
+#include "../src/Boundary.h"
+#include "../src/BoundaryCondition.h"
 #include "../src/FiniteElementMesh.h"
+#include "../src/Forcing.h"
 #include "../src/MatrixGroup.h"
 #include "../src/Node.h"
 #include "../src/Physics.h"

src/Physics/src/Boundary.cpp

+#include "Boundary.h"

src/Physics/src/Boundary.h

+#ifndef OERSTED_BOUNDARY_H
+#define OERSTED_BOUNDARY_H
+
+template<size_t Dimension>
+class Boundary {
+};
+
+template<>
+class Boundary<2> {
+public:
+    Boundary(std::vector<size_t> nodes) : Nodes{nodes} {};
+
+    std::vector<size_t> const &nodes() const { return Nodes; };
+
+    size_t const &node(size_t i) const { return Nodes[i]; };
+
+protected:
+    std::vector<size_t> Nodes;
+};
+
+#endif //OERSTED_BOUNDARY_H

src/Physics/src/BoundaryCondition.cpp

+#include "BoundaryCondition.h"

src/Physics/src/BoundaryCondition.h

+#ifndef OERSTED_BOUNDARYCONDITION_H
+#define OERSTED_BOUNDARYCONDITION_H
+
+#include <vector>
+#include <set>
+#include <algorithm>
+
+#include <Eigen/Sparse>
+
+#include "MatrixGroup.h"
+#include "QuadratureRule.h"
+#include "FiniteElementMesh.h"
+
+class BoundaryCondition {
+public:
+    BoundaryCondition(std::vector<size_t> boundaries) : Boundaries{boundaries} {};
+    virtual ~BoundaryCondition(){};
+
+    virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;
+    virtual void reduce(std::set<size_t,std::greater<size_t>> &index) const = 0;
+
+protected:
+    std::vector<size_t> Boundaries;
+};
+
+
+template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
+class ZeroDirichlet : public BoundaryCondition{};
+
+template<size_t ElementOrder, size_t QuadratureOrder>
+class ZeroDirichlet<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
+public:
+    ZeroDirichlet(std::vector<size_t> boundaries,
+                  FiniteElementMesh<2,ElementOrder> const &domain) : BoundaryCondition{boundaries} {
+
+        for(size_t i : boundaries) {
+            for(size_t j : domain.boundary(i).nodes()) {
+                Nodes.push_back(j);
+            }
+        }
+    };
+
+    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
+
+    void reduce(std::set<size_t,std::greater<size_t>> &index) const override {
+        for (size_t i = 0; i != Nodes.size(); ++i) {
+            index.insert(Nodes[i]);
+        }
+    };
+
+protected:
+    std::vector<size_t> Nodes;
+};
+#endif //OERSTED_BOUNDARYCONDITION_H

src/Physics/src/Condition.cpp

-//
-// Created by jpries on 1/2/17.
-//
-
-#include "Condition.h"

src/Physics/src/FiniteElementMesh.cpp

-//
-// Created by jpries on 12/20/16.
-//
-
 #include "FiniteElementMesh.h"

src/Physics/src/FiniteElementMesh.h

@@ -17,7 +17,7 @@ class FiniteElementMesh<2,Order> {
 public:
     FiniteElementMesh() {};
 
-    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<Region<2>> r) : Nodes(nodes), Triangles(tris), Regions(r) {};
+    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<Region<2>> r, std::vector<Boundary<2>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
 
     std::vector<XY> const &nodes() const { return Nodes; };
 
@@ -25,12 +25,16 @@ public:
 
     std::vector<Region<2>> const &regions() const { return Regions; };
 
+    std::vector<Boundary<2>> const &boundaries() const { return Boundaries; };
+
     XY const &node(size_t i) const { return Nodes[i]; };
 
     Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
 
     Region<2> const &region(size_t i) const { return Regions[i]; };
 
+    Boundary<2> const &boundary(size_t i) const { return Boundaries[i]; };
+
     template<size_t Q>
     DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> determinant() const {
         DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Triangles.size());
@@ -69,6 +73,7 @@ protected:
     std::vector<Triangle<Order>> Triangles;
 
     std::vector<Region<2>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
+    std::vector<Boundary<2>> Boundaries;
 };
 
 #endif //OERSTED_FINITEELEMENTMESH_H

src/Physics/src/Forcing.cpp

+#include "Forcing.h"

src/Physics/src/Condition.h → src/Physics/src/Forcing.h

-#ifndef OERSTED_CONDITION_H
-#define OERSTED_CONDITION_H
+#ifndef OERSTED_FORCING_H
+#define OERSTED_FORCING_H
 
 #include <functional>
 #include <vector>
 
 #include "Eigen"
-
 #include "Region.h"
-#include "Physics.h"
+#include "MatrixGroup.h"
+#include "QuadratureRule.h"
+#include "FiniteElementMesh.h"
 
-class Condition {
+class Forcing {
 public:
-    virtual ~Condition() {};
+    virtual ~Forcing() {};
 
     virtual Eigen::VectorXd operator()(double const t) const = 0;
 };
 
 template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
-class HomogeneousForcing : public Condition {
+class HomogeneousForcing : public Forcing {
 };
 
 template<size_t ElementOrder, size_t QuadratureOrder>
-class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Condition {
+class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Forcing {
 public:
-    HomogeneousForcing(Physics<2, ElementOrder, QuadratureOrder> &op, std::vector<size_t> reg, std::function<double(double)> f) : Operator{op}, Regions{reg}, Function{f} {
-        Indicator.resize(op.domain().triangles().size());  // TODO: num_elements (once quadrilateral elements added)
+    HomogeneousForcing(std::function<double(double)> f,
+                       std::vector<size_t> reg,
+                       SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &basis,
+                       DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &weights,
+                       FiniteElementMesh<2,ElementOrder> const &domain) : Regions{reg}, Basis(basis), Weights(weights), Function{f} {
+
+        Indicator.resize(domain.triangles().size());  // TODO: num_elements (once quadrilateral elements added)
         Indicator.setZero();
 
         for (size_t const &i : reg) {
-            for (size_t const &j : op.domain().region(i).triangles()) {
+            for (size_t const &j : domain.region(i).triangles()) {
                 Indicator(j) = 1.0;
             }
-    }
+        }
     };
 
     Eigen::VectorXd operator()(double const t) const override {
-        Eigen::VectorXd output = (Operator.basis()[0] * (Operator.weights()(0) * Indicator)); // TODO: ?SparseVector?
+        Eigen::VectorXd output = (Basis[0] * (Weights(0) * Indicator)); // TODO: ?SparseVector?
         for (size_t i = 1; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
-            output += (Operator.basis()[i] * (Operator.weights()(i) * Indicator));
+            output += (Basis[i] * (Weights(i) * Indicator));
         }
 
         output *= Function(t);
@@ -47,11 +53,12 @@ public:
     auto const &indicator() const { return Indicator; };
 
 protected:
-    Physics<2, ElementOrder, QuadratureOrder> &Operator;
+    std::function<double(double)> Function;
 
     std::vector<size_t> Regions;
 
-    std::function<double(double)> Function;
+    SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Basis;
+    DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Weights;
 
     Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
 };

src/Physics/src/MatrixGroup.h

@@ -24,6 +24,12 @@ public:
 
     auto &operator[](size_t const q) const { return Matrices[q]; };
 
+    void transform(Eigen::SparseMatrix<double> &A) {
+        for (size_t i = 0; i != Q; ++i) {
+            Matrices[i] = A * Matrices[i];
+        }
+    }
+
 protected:
     std::array<Eigen::SparseMatrix<double>, Q> Matrices;
 };
@@ -49,14 +55,19 @@ class DerivativeMatrixGroup {
 public:
     DerivativeMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) : Dx{rows, cols, nnz}, Dy{rows, cols, nnz} {};
 
-    auto &dx(size_t const q) { return Dx[q]; };
+    auto const &dx(size_t const q) const { return Dx[q]; };
 
-    auto &dy(size_t const q) { return Dy[q]; };
+    auto const &dy(size_t const q) const { return Dy[q]; };
 
     auto &dx(size_t const q, size_t const i, size_t const j) { return Dx(q, i, j); };
 
     auto &dy(size_t const q, size_t const i, size_t const j) { return Dy(q, i, j); };
 
+    void transform(Eigen::SparseMatrix<double> &A) {
+        Dx.transform(A);
+        Dy.transform(A);
+    }
+
 protected:
     SparseMatrixGroup<Q> Dx;
     SparseMatrixGroup<Q> Dy;

src/Physics/src/Physics.h

@@ -3,10 +3,14 @@
 
 #include <memory>
 
+#include "Eigen/SparseCholesky"
+#include "Eigen/IterativeLinearSolvers"
+
 #include "Node.h"
 #include "Triangle.h"
 #include "QuadratureRule.h"
 #include "FiniteElementMesh.h"
+#include "Forcing.h"
 
 enum class Variable {
     MagneticVectorPotential = 1,
@@ -27,8 +31,21 @@ enum class Variable {
     Psi = 26 + 23
 };
 
+class PhysicsInterface {
+public:
+    double time() { return Time; };
+
+    void time(double t) { Time = t; };
+
+protected:
+    double Time;
+
+    std::vector<std::unique_ptr<Forcing>> ForcingCondtions;
+    std::vector<std::unique_ptr<BoundaryCondition>> BoundaryConditions;
+};
+
 template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
-class Physics {
+class PhysicsData {
 public:
     static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
 
@@ -40,12 +57,12 @@ public:
     virtual void apply_conditions() = 0;
      */
 
-    Physics(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
+    PhysicsData(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
                                                      ElementWeights{d.triangles().size()},
                                                      Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
                                                      Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
 
-    FiniteElementMesh<Dimension, ElementOrder> &domain() const { return Domain; };
+    FiniteElementMesh<Dimension, ElementOrder> const &domain() const { return Domain; };
 
     DiagonalMatrixGroup<QuadraturePoints> const &weights() { return ElementWeights; };
 
@@ -62,29 +79,35 @@ protected:
 };
 
 template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, Variable V>
-class Magnetostatic : public Physics<Dimension, ElementOrder, QuadratureOrder> {
+class Magnetostatic : public PhysicsData<Dimension, ElementOrder, QuadratureOrder>, public PhysicsInterface {
 };
 
 template<size_t ElementOrder, size_t QuadratureOrder>
-class Magnetostatic<2, ElementOrder, QuadratureOrder, Variable::MagneticVectorPotential> : public Physics<2, ElementOrder, QuadratureOrder> {
+class Magnetostatic<2, ElementOrder, QuadratureOrder, Variable::MagneticVectorPotential> : public PhysicsData<2, ElementOrder, QuadratureOrder>, public PhysicsInterface {
 public:
-    using Physics<2, ElementOrder, QuadratureOrder>::QuadraturePoints;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::QuadraturePoints;
 
-    using Physics<2, ElementOrder, QuadratureOrder>::domain;
-    using Physics<2, ElementOrder, QuadratureOrder>::derivatives;
-    using Physics<2, ElementOrder, QuadratureOrder>::weights;
-    using Physics<2, ElementOrder, QuadratureOrder>::basis;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::domain;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::derivatives;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::weights;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::basis;
 
-    Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : Physics<2, ElementOrder, QuadratureOrder>{d} {};
+    Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.nodes().size()} {};
 
-    Eigen::VectorXd init_residual() { return Eigen::VectorXd(Domain.nodes().size()); };
+    Eigen::SparseMatrix<double> init_matrix() { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); };
 
-    Eigen::SparseMatrix<double> init_jacobian() { return Eigen::SparseMatrix<double>(Domain.nodes().size(), Domain.nodes().size()); };
+    Eigen::VectorXd init_vector() { return Eigen::VectorXd(Unknowns); };
 
-    void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) {
-        J.setZero();
-        r.setZero(); // TODO: Unecessary for r?
+    void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f) {
+        // Calculate residual
+        r = -f;
+        for (size_t i = 0; i!= TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
+            r += Derivatives.dx(i) * (ElementWeights(i) * (Derivatives.dx(i).transpose() * v))
+                 + Derivatives.dy(i) * (ElementWeights(i) * (Derivatives.dy(i).transpose() * v));
+        }
 
+        // Calculate jacobian
+        J.setZero();
         for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
             J += (Derivatives.dx(i) * ElementWeights(i) * Derivatives.dx(i).transpose())
                  + (Derivatives.dy(i) * ElementWeights(i) * Derivatives.dy(i).transpose());
@@ -101,13 +124,70 @@ public:
         }
     };
 
-    void apply_conditions() {};
+    void calculate_forcing(Eigen::VectorXd &f) {
+        f.setZero();
+        for(size_t i = 0; i!=ForcingCondtions.size();++i) {
+            f += ForcingCondtions[i]->operator()(Time);
+        }
+    }
+
+    void add_current_density(std::function<double(double)> func, std::vector<size_t> regions) {
+        ForcingCondtions.push_back(std::make_unique<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights, Domain));
+    }
+
+    void remove_current_density(size_t i) {
+        ForcingCondtions.erase(ForcingCondtions.begin() + i);
+    }
+
+    void add_magnetic_insulation(std::vector<size_t> boundaries) {
+        BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries, Domain));
+    }
+
+    void apply_conditions() {
+        std::vector<Eigen::Triplet<double>> triplets;
+        triplets.reserve(Domain.nodes().size());
+
+        for(size_t i = 0; i!=Domain.nodes().size();++i) {
+            triplets.push_back(Eigen::Triplet<double>(i,i,1.0));
+        }
+
+        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
+            BoundaryConditions[i]->apply(triplets);
+        }
+
+        std::set<size_t, std::greater<size_t>> index;
+        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
+            BoundaryConditions[i]->reduce(index);
+        }
+
+        size_t rows{Domain.nodes().size()};
+        size_t cols{Domain.nodes().size()};
+        for (auto &i : index) {
+            triplets.erase(triplets.begin()+i);
+            --rows;
+        }
+
+        for (size_t i = 0; i != triplets.size(); ++i) {
+            triplets[i] = Eigen::Triplet<double>(i, triplets[i].col(), triplets[i].value());
+        }
+
+        Eigen::SparseMatrix<double> bc_matrix(rows,cols);
+        bc_matrix.setFromTriplets(triplets.begin(),triplets.end());
+
+        Basis.transform(bc_matrix);
+        Derivatives.transform(bc_matrix);
+
+        Unknowns = rows;
+    };
 
 protected:
-    using Physics<2, ElementOrder, QuadratureOrder>::Domain;
-    using Physics<2, ElementOrder, QuadratureOrder>::Derivatives;
-    using Physics<2, ElementOrder, QuadratureOrder>::ElementWeights;
-    using Physics<2, ElementOrder, QuadratureOrder>::Basis;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::Domain;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights;
+    using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis;
+
+private:
+    size_t Unknowns;
 };
 
 #endif //OERSTED_PHYSICS_H

src/Physics/src/Region.h

@@ -2,7 +2,6 @@
 #define OERSTED_REGION_H
 
 #include <vector>
-#include <cstddef>
 
 template<size_t Dimension>
 class Region {

test/Physics/test_FiniteElementMesh.cpp

@@ -15,12 +15,17 @@ public:
 
         reg.push_back(Region<2>(std::vector<size_t>{0}));
 
-        femesh = FiniteElementMesh<2,1>(node, tri, reg);
+        bnd.push_back(Boundary<2>(std::vector<size_t>{0,1}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{0,2}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{1,2}));
+
+        femesh = FiniteElementMesh<2,1>(node, tri, reg, bnd);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
     std::vector<Region<2>> reg;
+    std::vector<Boundary<2>> bnd;
 
     FiniteElementMesh<2,1> femesh;
 };
@@ -66,23 +71,27 @@ TEST_F(MasterTriangle, magnetostatic) {
     static constexpr size_t Q = 1;
 
     Magnetostatic<2, 1, Q, Variable::A> ms(femesh);
+    ms.time(1.0);
 
-    Magnetostatic<2, 1, Q, Variable::A> *ms_ptr{&ms};
+    auto ff = [](double t) { return 1.0 * t; };
 
-    Physics<2,1,Q> *ph_ptr = dynamic_cast<Physics<2,1,Q>*>(ms_ptr);
+    ms.add_current_density(ff, std::vector<size_t>{0});
 
     ms.build_matrices();
 
-    auto J = ms.init_jacobian();
-    auto r = ms.init_residual();
+    auto J = ms.init_matrix();
+    auto v = ms.init_vector();
+    auto r = ms.init_vector();
+    auto f = ms.init_vector();
 
-    ms.residual(J, r);
+    ms.calculate_forcing(f);
+    ms.linearize(J, v, r, f);
 
-    auto ff = [](double t) { return 1.0 * t; };
+    //auto ff = [](double t) { return 1.0 * t; };
 
-    HomogeneousForcing<2,1,Q> hf{ms, std::vector<size_t>{0}, ff};
+    //HomogeneousForcing<2,1,Q> hf{ff, std::vector<size_t>{0}, ms.basis(), ms.weights(), ms.domain()};
 
-    std::cout << hf(1.0) << std::endl;
+    std::cout << f << std::endl;
 }
 
 class MasterTriangleRotationReflection : public ::testing::Test {
@@ -104,12 +113,14 @@ public:
 
         reg.push_back(Region<2>(std::vector<size_t>{0,1,2,3})); // Triangles 4,5,6 are duplicates
 
-        femesh = FiniteElementMesh<2, 1>(node, tri, reg);
+        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
     std::vector<Region<2>> reg;
+    std::vector<Boundary<2>> bnd;
+
 
     FiniteElementMesh<2, 1> femesh;
 };
@@ -231,7 +242,13 @@ public:
         reg.push_back(Region<2>(std::vector<size_t>{0}));
         reg.push_back(Region<2>(std::vector<size_t>{1}));
 
-        femesh = FiniteElementMesh<2, 1>(node, tri, reg);
+        bnd.push_back(Boundary<2>(std::vector<size_t>{0,1}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{0,2}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{1,2}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{1,3}));
+        bnd.push_back(Boundary<2>(std::vector<size_t>{2,3}));
+
+        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
 
         vx_m << 1.0, 0.0, 1.0, 0.0;
         vx_p << 0.0, 1.0, 0.0, 1.0;
@@ -242,7 +259,7 @@ public:
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
     std::vector<Region<2>> reg;
-
+    std::vector<Boundary<2>> bnd;
 
     Eigen::Vector4d vx_m;
     Eigen::Vector4d vx_p;
@@ -338,33 +355,222 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     //      TODO: Applied in reverse order, Motion->BoundaryCondition->Forcing,
     //      TODO: Or, TimeVarying->TimeInvariant->Forcing?
 
-
     Magnetostatic<2, 1, 2, Variable::A> ms{femesh};
-
+    ms.time(1.0);
     ms.build_matrices();
 
     auto ff = [](double t) { return 1.0 * t; };
+    auto f = ms.init_vector();
+
+    ms.add_current_density(ff, std::vector<size_t>{0});
+    ms.calculate_forcing(f);
+    EXPECT_NEAR(f(0),1.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(1),1.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(2),1.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(3),0.0/6.0,FLT_EPSILON);
+    ms.remove_current_density(0);
+
+    ms.add_current_density(ff, std::vector<size_t>{1});
+    ms.calculate_forcing(f);
+    EXPECT_NEAR(f(0),0.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(1),1.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(2),1.0/6.0,FLT_EPSILON);
+    EXPECT_NEAR(f(3),1.0/6.0,FLT_EPSILON);
+    ms.remove_current_density(0);
+
+    ms.add_current_density(ff, std::vector<size_t>{0,1});