Updating values in SolveStep

src/electrophysiology/solvers/LinearImplicitSolver.cpp

@@ -123,7 +123,6 @@ void LinearImplicitSolver::SolvePDE(IElphyAssemble &A) {
 
 }
 
-
 void LinearImplicitSolver::updateValues(Vectors &values) {
   values[0] = (*gating)[0];
   if (M.Type() == TENSION) {
@@ -214,7 +213,7 @@ void ImplicitSolver::Solve(IElphyAssemble &A, Vector &V) {
 
 void ImplicitSolver::SolveStep(IElphyAssemble &A, Vectors &values) {
   SolveGating(A.Time(), A.StepSize());
-  SolveConcentration(A.Time(), A.StepSize());
+  SolveConcentration(A.Time(), A.StepSize(), values);
   A.updateVCW(*gating);
   vout(10).StartBlock("Solve PDE");
   bool conv = SolvePDE(A, values);
@@ -285,46 +284,46 @@ void SemiImplicitSolver::SolvePDE(IElphyAssemble &A) {
 }
 
 void LinearImplicitSolver::SolveStep(IElphyAssemble &A, Vectors &values) {
-  StaggeredStep(A);
+  StaggeredStep(A, values);
   updateMaxMin();
   updateValues(values);
 }
 
 void LinearImplicitSolver::selectOrder(std::string o) {
   if (o == "Vcw") {
-    StaggeredStep = [this](IElphyAssemble &A) {
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
       SolvePDE(A);
-      SolveConcentration(A.Time(), A.StepSize());
+      SolveConcentration(A.Time(), A.StepSize(), values);
       SolveGating(A.Time(), A.StepSize());
     };
   } else if (o == "Vwc") {
-    StaggeredStep = [this](IElphyAssemble &A) {
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
       SolvePDE(A);
       SolveGating(A.Time(), A.StepSize());
-      SolveConcentration(A.Time(), A.StepSize());
+      SolveConcentration(A.Time(), A.StepSize(), values);
     };
   } else if (o == "cVw") {
-    StaggeredStep = [this](IElphyAssemble &A) {
-      SolveConcentration(A.Time(), A.StepSize());
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
+      SolveConcentration(A.Time(), A.StepSize(), values);
       SolvePDE(A);
       SolveGating(A.Time(), A.StepSize());
     };
   } else if (o == "cwV") {
-    StaggeredStep = [this](IElphyAssemble &A) {
-      SolveConcentration(A.Time(), A.StepSize());
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
+      SolveConcentration(A.Time(), A.StepSize(), values);
       SolveGating(A.Time(), A.StepSize());
       SolvePDE(A);
     };
   } else if (o == "wVc") {
-    StaggeredStep = [this](IElphyAssemble &A) {
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
       SolveGating(A.Time(), A.StepSize());
       SolvePDE(A);
-      SolveConcentration(A.Time(), A.StepSize());
+      SolveConcentration(A.Time(), A.StepSize(), values);
     };
   } else {//wcV
-    StaggeredStep = [this](IElphyAssemble &A) {
+    StaggeredStep = [this](IElphyAssemble &A, const Vectors &values) {
       SolveGating(A.Time(), A.StepSize());
-      SolveConcentration(A.Time(), A.StepSize());
+      SolveConcentration(A.Time(), A.StepSize(), values);
       SolvePDE(A);
     };
   }
@@ -345,8 +344,8 @@ void LinearImplicitSolver::selectIonScheme(const string &solvingIonScheme) {
   }
 }
 
-void LinearImplicitSolver::SolveConcentration(double t, double dt) {
-  std::vector<double> vcw(M.Size());
+void LinearImplicitSolver::SolveConcentration(double t, double dt, const Vectors &values) {
+  std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
   std::vector<vector<double>> G(ionSteps);
   for (int i = 0; i < ionSteps; ++i) {
     G[i] = std::vector<double>(M.Size());
@@ -356,6 +355,9 @@ void LinearImplicitSolver::SolveConcentration(double t, double dt) {
       for (int i = 0; i < M.Size(); ++i) {
         vcw[i] = (*gating)[i](r, j);
       }
+      if (M.Type() == TENSION) {
+        vcw[M.Size()] = values[2](r, j);
+      }
       SolveConcentrationStep(M, t, M.TimeScale() * dt, vcw, G);
       for (int i = 1; i < M.GatingIndex(); i++) {
         (*gating)[i](r, j) = vcw[i];

src/electrophysiology/solvers/LinearImplicitSolver.hpp

@@ -78,7 +78,7 @@ public:
 
   virtual void Solve(IElphyAssemble &A, Vector &V);
 
-  void SolveConcentration(double t, double dt);
+  void SolveConcentration(double t, double dt, const Vectors &values);
 
   void SolveGating(double t, double dt);
 
@@ -86,7 +86,7 @@ public:
 
   void SolveStep(IElphyAssemble &A, Vectors &values) override;
 
-  std::function<void(IElphyAssemble &A)> StaggeredStep;
+  std::function<void(IElphyAssemble &A, const Vectors &values)> StaggeredStep;
 
 
   void selectOrder(std::string o);