diff --git a/mpp b/mpp
index 3c5bdc99ad054f3accb9a024179325ded12e9192..cfd2f66c53e10c51415ddc429e56a002ad7b97b0 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit 3c5bdc99ad054f3accb9a024179325ded12e9192
+Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0
diff --git a/src/CardiacAssemble.hpp b/src/CardiacAssemble.hpp
index 57f794183142ab578b41dd3275114f89320cbab9..df092a4e24cf319260000a56b140f9d6b132d68f 100644
--- a/src/CardiacAssemble.hpp
+++ b/src/CardiacAssemble.hpp
@@ -20,7 +20,7 @@ protected:
void setVtkSteps(const std::string& configString){
Config::Get(configString, vtksteps);
- vtksteps = max(vtksteps, 1);
+ vtksteps = std::max(vtksteps, 1);
}
public:
explicit ICardiacAssemble(const std::string& plotEntry= "PlotVTK", const std::string& stepEntry = "PlottingSteps") {
diff --git a/src/CardiacData.cpp b/src/CardiacData.cpp
index 7bfc2a845a9fe6bb7a979f5b620d34e1b3e6e842..15ad6925580824b4c74bf58c872d12c69ca81a67 100644
--- a/src/CardiacData.cpp
+++ b/src/CardiacData.cpp
@@ -54,7 +54,7 @@ void emplace_chamber(CHAMBER c, std::vector<int> & chambers){
chambers.emplace_back(c + 100 + i);
}
}
-std::vector<int> getIncludedChambers(const string &chamberName) {
+std::vector<int> getIncludedChambers(const std::string &chamberName) {
std::vector<int> chambers{};
if (chamberName == "leftventricle") {
emplace_chamber(LEFT_VENTRICLE, chambers);
diff --git a/src/CardiacData.hpp b/src/CardiacData.hpp
index 1bbaca31d046f7ba9e098128a18244c75e9016d4..10b97a16ec621b38c22457ed670348d9f19f6528 100644
--- a/src/CardiacData.hpp
+++ b/src/CardiacData.hpp
@@ -28,8 +28,7 @@ enum CHAMBER {
PERICARD = 11000
};
-std::vector<int> getIncludedChambers(const string &chamberName = "");
-
+std::vector<int> getIncludedChambers(const std::string &chamberName = "");
static bool inChamber(int domain, int chamber) {
switch (chamber) {
@@ -168,10 +167,10 @@ static Tensor DeformationGradient(const ElemenT &E, const Point &z, const Vector
double Volume(CELLTYPE type, const std::vector<Point> &corners);
-static string stepAsString(int step) {
+static std::string stepAsString(int step) {
char buffer[256];
sprintf(buffer, "%04d", step);
- string str(buffer);
- return string("-") + str;
+ std::string str(buffer);
+ return std::string("-") + str;
}
#endif //CARDIACDATA_HPP
diff --git a/src/CardiacInterpolation.cpp b/src/CardiacInterpolation.cpp
index 7256d6645bac1ed8e7ae3123472f23969ebddb0a..70de09f7d2eaab32223330f42b6baf7babc6f797 100644
--- a/src/CardiacInterpolation.cpp
+++ b/src/CardiacInterpolation.cpp
@@ -46,10 +46,9 @@ void fillLagrangeData(Vector &vec, row r, DataContainer d) {
}*/
}
-
-void
-fillLagrangeCellData(const vector<Point> &corners, vector<DataContainer> data, vector<Point> &cells,
- Vector &vec) {
+void fillLagrangeCellData(const std::vector<Point> &corners,
+ std::vector<DataContainer> data,
+ std::vector<Point> &cells, Vector &vec) {
auto min_index = std::distance(data.begin(),
std::min_element(data.begin(), data.end(), compare_excitation));
auto max_index = std::distance(data.begin(),
@@ -399,7 +398,8 @@ void smoothMeshData(const Mesh &mesh, Vector &vec) {
}
if (max_exc_time >= 0.0 && max_amplitude > 0.0) {
for (int i = 0; i < c.Corners(); ++i) {
- double amplitude = min(mesh.find_vertex(c.Corner(i)).GetData()[2], max_amplitude);
+ double amplitude =
+ std::min(mesh.find_vertex(c.Corner(i)).GetData()[2], max_amplitude);
DataContainer exc_data({max_exc_time, max_duration, amplitude});
row cornerRow = vec.find_row(c.Corner(i));
mesh.find_vertex(c.Corner(i)).SetData(exc_data);
diff --git a/src/CardiacProblem.hpp b/src/CardiacProblem.hpp
index 4202228ffb9221b9e952f25903318c3e802c77d2..4c16fe75136b9bd057123dee7d8c5db4b87382fd 100644
--- a/src/CardiacProblem.hpp
+++ b/src/CardiacProblem.hpp
@@ -11,8 +11,8 @@ static constexpr int PRESSURE_BC_RA = 233;
static constexpr int ROBIN_BC_EPI = 330;
static constexpr int ROBIN_BC_BASE = 331;
-
-static std::pair<int, bool> isIn(const Point &p, const vector<Point> &pointList) {
+static std::pair<int, bool> isIn(const Point &p,
+ const std::vector<Point> &pointList) {
for (int i = 0; i < pointList.size(); ++i) {
if (p.isNear(pointList[i], 1e-4)) {
return {i, true};
diff --git a/src/ExternalCurrent.cpp b/src/ExternalCurrent.cpp
index 7ddd0ca0caea198f5e2acc0e952fa0b26e999c3b..361f829ef214beb3af99960acf2e1dd0c23fa6e9 100644
--- a/src/ExternalCurrent.cpp
+++ b/src/ExternalCurrent.cpp
@@ -1,4 +1,7 @@
#include "ExternalCurrent.hpp"
+
+#include <numbers>
+
#include "GlobalDefinitions.hpp"
std::function<double(double,double,double,double)> excitationFunc=IextStepfct;
@@ -17,7 +20,7 @@ double IextTanh(double t, double amp, double exc, double dur) {
double IextArcTan(double t, double amp, double exc, double dur) {
double start_t = 4 * (t - exc);
double end_t = 4 * (t - (exc + dur));
- return (atan(start_t) - atan(end_t)) / Pi * amp;
+ return (atan(start_t) - atan(end_t)) / std::numbers::pi * amp;
}
void SelectExternalCurrentFunction(const std::string &externalCurrentFunction) {
diff --git a/src/MainElphy.cpp b/src/MainElphy.cpp
index 0b630e35fcfecc371d77bdb647bbe251ae61b6b3..26bde2f68a0e7c8c09df38f9a4e10fd787373bf4 100644
--- a/src/MainElphy.cpp
+++ b/src/MainElphy.cpp
@@ -32,7 +32,7 @@ int main(int argc, char **argv) {
Date Start;
- string model = "";
+ std::string model = "";
Config::Get("Model", model);
MainCardMech *cmMain;
diff --git a/src/cellmodels/MBeelerReuter.cpp b/src/cellmodels/MBeelerReuter.cpp
index 358d522357969a22d849dd5359919ef47b116b5b..3a4a2badf344de6df1aa5ff019ee314e0366d8e3 100644
--- a/src/cellmodels/MBeelerReuter.cpp
+++ b/src/cellmodels/MBeelerReuter.cpp
@@ -45,7 +45,7 @@ void MBeelerReuter::Initialize(Vectors &VW) {
BR_x_1(VW) = 0.00565;//x_1
}
-void MBeelerReuter::Initialize(vector<double> &vw) {
+void MBeelerReuter::Initialize(std::vector<double> &vw) {
BR_Pot(vw) = initialPotential(); //V
BR_Ca(vw) = 0.0000002; // Ca
BR_d(vw) = 0.002980; //d
@@ -97,7 +97,7 @@ double MBeelerReuter::JacobiEntry(int i, int j, const std::vector<double> &VW) c
}
}
-double MBeelerReuter::dIiondV(const vector<double> &VW) const {
+double MBeelerReuter::dIiondV(const std::vector<double> &VW) const {
double V = BR_Pot(VW);
double dI_NaToV = g_Na * BR_m(VW) * BR_m(VW) * BR_m(VW) * BR_h(VW) * BR_j(VW) + g_NaC;
double dI_sToV = g_s * BR_d(VW) * BR_f(VW);
@@ -163,7 +163,7 @@ double MBeelerReuter::dGw(int i, int j, const std::vector<double> &VW) const {
}
}
-vector<double> MBeelerReuter::Tau(const vector<double> &var) const {
+std::vector<double> MBeelerReuter::Tau(const std::vector<double> &var) const {
auto V = BR_Pot(var);
return {
BR_tau_d(V),
@@ -174,7 +174,8 @@ vector<double> MBeelerReuter::Tau(const vector<double> &var) const {
BR_tau_x_1(V)};
}
-vector<double> MBeelerReuter::YInfty(const vector<double> &var) const {
+std::vector<double> MBeelerReuter::YInfty(
+ const std::vector<double> &var) const {
auto V = BR_Pot(var);
return {
BR_yinfty_d(V),
@@ -186,8 +187,9 @@ vector<double> MBeelerReuter::YInfty(const vector<double> &var) const {
};
}
-void
-MBeelerReuter::ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const {
+void MBeelerReuter::ExpIntegratorUpdate(double dt,
+ const std::vector<double> &var,
+ std::vector<double> &g) const {
if (!useExpInt) return;
double V = BR_Pot(var);
@@ -200,7 +202,7 @@ MBeelerReuter::ExpIntegratorUpdate(double dt, const vector<double> &var, vector<
}
void
-MBeelerReuter::IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+MBeelerReuter::IIonUpdate(const std::vector<double> &var, std::vector<double> &g, double iExt) const {
auto V = BR_Pot(var);
@@ -219,14 +221,16 @@ MBeelerReuter::IIonUpdate(const vector<double> &var, vector<double> &g, double i
}
}
-void MBeelerReuter::ConcentrationUpdate(const vector<double> &var, vector<double> &g,
+void MBeelerReuter::ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g,
double iExt) const {
- double I_s = g_s * BR_d(var) * BR_f(var) * (BR_Pot(var) - (-82.3 - 13.0287 * log(BR_Ca(var))));
+ double I_s = g_s * BR_d(var) * BR_f(var) *
+ (BR_Pot(var) - (-82.3 - 13.0287 * log(BR_Ca(var))));
BR_Ca(g) = -0.0000001 * I_s + 0.07 * (0.0000001 - BR_Ca(var));
}
-
-void MBeelerReuter::GatingUpdate(const std::vector<double> &var, vector<double> &g) const {
+void MBeelerReuter::GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const {
if (useExpInt) return;
auto V = BR_Pot(var);
@@ -238,7 +242,8 @@ void MBeelerReuter::GatingUpdate(const std::vector<double> &var, vector<double>
BR_x_1(g) = (BR_yinfty_x_1(V) - BR_x_1(var)) / BR_tau_x_1(V);
}
-void MBeelerReuter::UpdateValues(double dt, vector<double> &var, const vector<double> &g) const {
+void MBeelerReuter::UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const {
BR_Pot(var) += dt * BR_Pot(g);
BR_Ca(var) += dt * BR_Ca(g);
if (useExpInt) {
diff --git a/src/cellmodels/MBeelerReuter.hpp b/src/cellmodels/MBeelerReuter.hpp
index 3373d4c67c79dc1759734fcf457e7cca158e04dd..f33230ffe4e1031306d7045cd0b9e20d7c716708 100644
--- a/src/cellmodels/MBeelerReuter.hpp
+++ b/src/cellmodels/MBeelerReuter.hpp
@@ -52,7 +52,7 @@ class MBeelerReuter : public MElphyModel {
[[nodiscard]] double dBeta(int i, double V) const;
- [[nodiscard]] double dIiondV(const vector<double> &VW) const;
+ [[nodiscard]] double dIiondV(const std::vector<double> &VW) const;
[[nodiscard]] double dIion(int j, const std::vector<double> &VW) const;
@@ -78,20 +78,24 @@ public:
void Initialize(std::vector<double> &vw) override;
- void IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void IIonUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override;
- void
- ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const override;
- void GatingUpdate(const vector<double> &var, vector<double> &g) const override;
+ void GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const override;
- void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override;
+ void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const override;
- void ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const override;
+ void ExpIntegratorUpdate(double dt, const std::vector<double> &var,
+ std::vector<double> &g) const override;
- vector<double> Tau(const vector<double> &var) const override;
+ std::vector<double> Tau(const std::vector<double> &var) const override;
- vector<double> YInfty(const vector<double> &var) const override;
+ std::vector<double> YInfty(const std::vector<double> &var) const override;
double GetIion(const std::vector<double> &var) const override;
diff --git a/src/cellmodels/MCellModel.cpp b/src/cellmodels/MCellModel.cpp
index a6b398dc0ffd7c119b25ecedae0d8d4dc2dd96ec..a3f75cf7db6283aa2738c1a8e37dcc08cb77a7fa 100644
--- a/src/cellmodels/MCellModel.cpp
+++ b/src/cellmodels/MCellModel.cpp
@@ -7,18 +7,18 @@
std::unique_ptr<MElphyModel> GetElphyModel() {
- string cmClass = "MElphyModel";
+ std::string cmClass = "MElphyModel";
Config::Get("ElphyModelClass", cmClass);
if(contains(cmClass,"IBT"))
return GetElphyModel("IBTCellModel");
- string cmName = "FitzHughNagumo";
+ std::string cmName = "FitzHughNagumo";
Config::Get("ElphyModelName", cmName);
return GetElphyModel(cmName);
}
-std::unique_ptr<MElphyModel> GetElphyModel(const string &name) {
+std::unique_ptr<MElphyModel> GetElphyModel(const std::string &name) {
if (name == "FitzHughNagumo") {
return std::make_unique<MFitzHughNagumo>();
}
@@ -28,16 +28,15 @@ std::unique_ptr<MElphyModel> GetElphyModel(const string &name) {
if (name == "TenTusscher") {
return std::make_unique<MTenTusscher>();
}
- if(contains(name, "IBT")){
+ if (contains(name, "IBT")) {
return nullptr;
}
THROW(name + " cellmodel is not implemented.")
-
}
-
-std::unique_ptr<MTensionModel> GetTensionModel(const string &name, std::unique_ptr<MElphyModel> &&M) {
+std::unique_ptr<MTensionModel> GetTensionModel(
+ const std::string &name, std::unique_ptr<MElphyModel> &&M) {
if (name == "Rossi") {
return std::make_unique<MRossi>(std::move(M));
}
@@ -46,7 +45,7 @@ std::unique_ptr<MTensionModel> GetTensionModel(const string &name, std::unique_p
}
std::unique_ptr<MTensionModel> GetTensionModel(std::unique_ptr<MElphyModel> &&M) {
- string cmName = "Rossi";
+ std::string cmName = "Rossi";
Config::Get("TensionModelName", cmName);
return GetTensionModel(cmName, std::move(M));
}
diff --git a/src/cellmodels/MCellModel.hpp b/src/cellmodels/MCellModel.hpp
index 4f4ba7f8d9bd33984fec9f33558788b900a9ec37..242ff30c76f18f759b77deb7c296ad268b5cb4bf 100644
--- a/src/cellmodels/MCellModel.hpp
+++ b/src/cellmodels/MCellModel.hpp
@@ -78,8 +78,8 @@ public:
* @param var current values at time t
* @param[out] g stores the next value var(t+dt)
*/
- virtual void
- ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const {};
+ virtual void ExpIntegratorUpdate(double dt, const std::vector<double> &var,
+ std::vector<double> &g) const {};
/**
* Calculates the change rates of concentration-type variables and stores the results in g.
@@ -88,7 +88,8 @@ public:
* @param[out] g change rates
* @param iExt external Current
*/
- virtual void ConcentrationUpdate(const std::vector<double> &var, vector<double> &g,
+ virtual void ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g,
double iExt) const {};
/**
@@ -97,7 +98,8 @@ public:
* @param var current values at time t
* @param[out] g change rates
*/
- virtual void GatingUpdate(const std::vector<double> &var, vector<double> &g) const {};
+ virtual void GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const {};
/**
* Calculates the change rates of all variables and stores the results in g.
@@ -105,8 +107,8 @@ public:
* @param[out] g change rates
* @param iExt external current
*/
- virtual void
- RHSUpdate(const std::vector<double> &var, vector<double> &g, double iExt) const = 0;
+ virtual void RHSUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const = 0;
/**
* Updates the values with the given timestep and change rates
@@ -114,7 +116,8 @@ public:
* @param[out] var values which are updated
* @param g change rates (time derivatives of the correspondig ODE, depending on the method)
*/
- virtual void UpdateValues(double dt, std::vector<double> &var, const vector<double> &g) const = 0;
+ virtual void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const = 0;
/**
* Returns a value to initialize external components
@@ -131,7 +134,7 @@ public:
* for the given cell model.
* @param var current values to calculate tau
*/
- virtual std::vector<double> Tau(const vector<double> &var) const {
+ virtual std::vector<double> Tau(const std::vector<double> &var) const {
THROW("Not implemented in base class")
}
@@ -140,7 +143,7 @@ public:
* variables for the given cell model.
* @param var current values to calculate y_infty
*/
- virtual std::vector<double> YInfty(const vector<double> &var) const {
+ virtual std::vector<double> YInfty(const std::vector<double> &var) const {
THROW("Not implemented in base class")
}
@@ -176,10 +179,11 @@ public:
return initialPotential();
};
- virtual void
- IIonUpdate(const std::vector<double> &var, vector<double> &g, double iExt) const {};
+ virtual void IIonUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const {};
- void RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const override {
+ void RHSUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override {
IIonUpdate(var, g, iExt);
ConcentrationUpdate(var, g, iExt);
GatingUpdate(var, g);
@@ -231,43 +235,47 @@ public:
elphy->UseExpInt(expInt);
}
- vector<double> Tau(const vector<double> &var) const override {
+ std::vector<double> Tau(const std::vector<double> &var) const override {
return elphy->Tau(var);
}
- vector<double> YInfty(const vector<double> &var) const override {
+ std::vector<double> YInfty(const std::vector<double> &var) const override {
return elphy->YInfty(var);
}
- double JacobiEntry(int i, int j, const vector<double> &VW) const override {
+ double JacobiEntry(int i, int j,
+ const std::vector<double> &VW) const override {
return elphy->JacobiEntry(i, j, VW);
};
- void AlgebraicUpdate(double dt, vector<double> &var) const override {
+ void AlgebraicUpdate(double dt, std::vector<double> &var) const override {
elphy->AlgebraicUpdate(dt, var);
};
- void ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const override {
+ void ExpIntegratorUpdate(double dt, const std::vector<double> &var,
+ std::vector<double> &g) const override {
elphy->ExpIntegratorUpdate(dt, var, g);
};
- void
- ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override {
+ void ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const override {
elphy->ConcentrationUpdate(var, g, iExt);
TensionUpdate(var, g, iExt);
};
- void GatingUpdate(const vector<double> &var, vector<double> &g) const override {
+ void GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const override {
elphy->GatingUpdate(var, g);
};
-
- void RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const override {
+ void RHSUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override {
elphy->RHSUpdate(var, g, iExt);
TensionUpdate(var, g, iExt);
};
- virtual void TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const = 0;
+ virtual void TensionUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const = 0;
};
static std::vector<int> ConcentrationIndices(const MCellModel &M) {
@@ -283,9 +291,10 @@ static std::vector<int> ConcentrationIndices(const MCellModel &M) {
std::unique_ptr<MElphyModel> GetElphyModel();
-std::unique_ptr<MElphyModel> GetElphyModel(const string &name);
+std::unique_ptr<MElphyModel> GetElphyModel(const std::string &name);
-std::unique_ptr<MTensionModel> GetTensionModel(const string &name, std::unique_ptr<MElphyModel> &&);
+std::unique_ptr<MTensionModel> GetTensionModel(const std::string &name,
+ std::unique_ptr<MElphyModel> &&);
std::unique_ptr<MTensionModel> GetTensionModel(std::unique_ptr<MElphyModel> &&);
diff --git a/src/cellmodels/MFitzHughNagumo.cpp b/src/cellmodels/MFitzHughNagumo.cpp
index 1e9899342059989c9a382c74845d6bc2c1a638c8..bbf3f921cfb1358c4e7db36ee655674aac75644b 100644
--- a/src/cellmodels/MFitzHughNagumo.cpp
+++ b/src/cellmodels/MFitzHughNagumo.cpp
@@ -7,32 +7,35 @@ void MFitzHughNagumo::Initialize(Vectors &var) {
var[1] = 0.0;
}
-void MFitzHughNagumo::Initialize(vector<double> &var) {
+void MFitzHughNagumo::Initialize(std::vector<double> &var) {
var[0] = initialPotential();
var[1] = 0.0;
}
-vector<double> MFitzHughNagumo::Tau(const vector<double> &var) const {
+std::vector<double> MFitzHughNagumo::Tau(const std::vector<double> &var) const {
return {};
}
-vector<double> MFitzHughNagumo::YInfty(const vector<double> &var) const {
+std::vector<double> MFitzHughNagumo::YInfty(
+ const std::vector<double> &var) const {
return {};
}
-void
-MFitzHughNagumo::IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+void MFitzHughNagumo::IIonUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const {
double mVolt = var[0];
g[0] = c1 * mVolt * (mVolt - a) * (1 - mVolt) - c2 * var[1] + iExt;
}
-void MFitzHughNagumo::ConcentrationUpdate(const vector<double> &var, vector<double> &g,
+void MFitzHughNagumo::ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g,
double iExt) const {
double mVolt = var[0];
g[1] = b * (mVolt - c3 * var[1]);
}
-void MFitzHughNagumo::UpdateValues(double dt, vector<double> &var, const vector<double> &g) const {
+void MFitzHughNagumo::UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const {
var[0] += dt * g[0];
var[1] += dt * g[1];
}
diff --git a/src/cellmodels/MFitzHughNagumo.hpp b/src/cellmodels/MFitzHughNagumo.hpp
index 3d02e5e6392ff74ae678e63fd9de714e339ca032..a38db43a259de2843b3f4360f131544f9c6db893 100644
--- a/src/cellmodels/MFitzHughNagumo.hpp
+++ b/src/cellmodels/MFitzHughNagumo.hpp
@@ -30,16 +30,18 @@ public:
void Initialize(std::vector<double> &vw) override;
- void IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void IIonUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override;
- void
- ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const override;
- void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override;
+ void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const override;
- vector<double> Tau(const vector<double> &var) const override;
+ std::vector<double> Tau(const std::vector<double> &var) const override;
- vector<double> YInfty(const vector<double> &var) const override;
+ std::vector<double> YInfty(const std::vector<double> &var) const override;
};
#endif //MFITZHUGHNAGUMO_HPP
diff --git a/src/cellmodels/MRossi.cpp b/src/cellmodels/MRossi.cpp
index 43d95401e6df50bbe4cd77a8a85fc451ab1714cb..bb910f390631ff35366f82c66577b01da5ffb78a 100644
--- a/src/cellmodels/MRossi.cpp
+++ b/src/cellmodels/MRossi.cpp
@@ -41,7 +41,7 @@ void MRossi::Initialize(Vectors &vectors) {
restingCA = elphyValues[elphy->CalciumIndex()] * elphy->CAScale();
}
-void MRossi::Initialize(vector<double> &vw) {
+void MRossi::Initialize(std::vector<double> &vw) {
elphy->Initialize(vw);
R_gamma(vw) = 0.0;
R_force(vw) = 0.0;
@@ -49,7 +49,8 @@ void MRossi::Initialize(vector<double> &vw) {
restingCA = vw[elphy->CalciumIndex()] * elphy->CAScale();
}
-void MRossi::TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+void MRossi::TensionUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const {
double gamma = R_gamma(var);
double ca = var[elphy->CalciumIndex()] * elphy->CAScale();
double iota = R_iota(var);
@@ -73,7 +74,8 @@ void MRossi::TensionUpdate(const vector<double> &var, vector<double> &g, double
R_force(g) = activeForce;
}
-void MRossi::UpdateValues(double dt, vector<double> &var, const vector<double> &g) const {
+void MRossi::UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const {
elphy->UpdateValues(dt, var, g);
R_gamma(var) += dt * R_gamma(g);
R_force(var) += dt * R_force(g);
diff --git a/src/cellmodels/MRossi.hpp b/src/cellmodels/MRossi.hpp
index 1012ed2b7b53e3c1bcfc2d92bca4aa730866f195..42cbc800629e1ddfb744c95f433d55aef01a3b17 100644
--- a/src/cellmodels/MRossi.hpp
+++ b/src/cellmodels/MRossi.hpp
@@ -58,11 +58,11 @@ public:
* 1. calcium
* 2. I4 (macroscopic second material invariant)
*/
- void TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
-
- void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override;
-
+ void TensionUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override;
+ void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const override;
};
#endif //MROSSI_HPP
diff --git a/src/cellmodels/MTenTusscher.cpp b/src/cellmodels/MTenTusscher.cpp
index 92cc64fbb7187a85f1b57ddf443e4da5fdb41efd..97f5809ddc30db3515dd076b276de2172ef54a47 100644
--- a/src/cellmodels/MTenTusscher.cpp
+++ b/src/cellmodels/MTenTusscher.cpp
@@ -23,7 +23,7 @@ void MTenTusscher::Initialize(Vectors &VW) {
}
-void MTenTusscher::Initialize(vector<double> &vw) {
+void MTenTusscher::Initialize(std::vector<double> &vw) {
TT_Pot(vw) = initialPotential();
TT_Ca(vw) = 0.000126;//0.00007//0.000153;
TT_Ca_SS(vw) = 0.00036;//0.00007//0.00042;
@@ -43,7 +43,6 @@ void MTenTusscher::Initialize(vector<double> &vw) {
TT_f_2(vw) = 0.9755; // 1.0//0.9526;
TT_f_Ca(vw) = 0.9953; // 1.0// 0.9942;
TT_R_quer(vw) = 0.9073; // 1.0//0.8978;
-
}
double tauFormula_h(double V) {
@@ -94,42 +93,51 @@ double tauFormula_j(double V) {
#define TT_yinfty_f_2(V) (0.67 / (1. + exp((V + 35.) / 7)) + 0.33)
#define TT_yinfty_f_Ca(CaSS) (0.6 / (1 + (CaSS / 0.05) * (CaSS / 0.05)) + 0.4)
-
-void
-MTenTusscher::ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const {
- if (!useExpInt) return;
+void MTenTusscher::ExpIntegratorUpdate(double dt,
+ const std::vector<double> &var,
+ std::vector<double> &g) const {
+ if (!useExpInt)
+ return;
// This should not be here!
double kcasr =
- max_sr - ((max_sr - min_sr) / (1.00000 + (EC / TT_Ca_SR(var)) * (EC / TT_Ca_SR(var)))); // x
- double k_2 = k_2prime * kcasr; // x
- TT_R_quer(g) =
- TT_R_quer(var) + dt * (k_4 * (1 - TT_R_quer(var)) - k_2 * TT_Ca_SS(var) * TT_R_quer(var));
-
+ max_sr - ((max_sr - min_sr) /
+ (1.00000 + (EC / TT_Ca_SR(var)) * (EC / TT_Ca_SR(var)))); // x
+ double k_2 = k_2prime * kcasr; // x
+ TT_R_quer(g) = TT_R_quer(var) + dt * (k_4 * (1 - TT_R_quer(var)) -
+ k_2 * TT_Ca_SS(var) * TT_R_quer(var));
double V = TT_Pot(var);
- TT_m(g) = (TT_yinfty_m(V) - (TT_yinfty_m(V) - TT_m(var)) * exp(-dt / TT_tau_m(V)));
- TT_h(g) = (TT_yinfty_h(V) - (TT_yinfty_h(V) - TT_h(var)) * exp(-dt / TT_tau_h(V)));
- TT_j(g) = (TT_yinfty_j(V) - (TT_yinfty_j(V) - TT_j(var)) * exp(-dt / TT_tau_j(V)));
+ TT_m(g) =
+ (TT_yinfty_m(V) - (TT_yinfty_m(V) - TT_m(var)) * exp(-dt / TT_tau_m(V)));
+ TT_h(g) =
+ (TT_yinfty_h(V) - (TT_yinfty_h(V) - TT_h(var)) * exp(-dt / TT_tau_h(V)));
+ TT_j(g) =
+ (TT_yinfty_j(V) - (TT_yinfty_j(V) - TT_j(var)) * exp(-dt / TT_tau_j(V)));
TT_x_r1(g) = (TT_yinfty_x_r1(V) -
- (TT_yinfty_x_r1(V) - TT_x_r1(var)) * exp(-dt / TT_tau_x_r1(V)));
+ (TT_yinfty_x_r1(V) - TT_x_r1(var)) * exp(-dt / TT_tau_x_r1(V)));
TT_x_r2(g) = (TT_yinfty_x_r2(V) -
- (TT_yinfty_x_r2(V) - TT_x_r2(var)) * exp(-dt / TT_tau_x_r2(V)));
- TT_x_s(g) = (TT_yinfty_x_s(V) - (TT_yinfty_x_s(V) - TT_x_s(var)) * exp(-dt / TT_tau_x_s(V)));
- TT_s(g) = (TT_yinfty_s(V) - (TT_yinfty_s(V) - TT_s(var)) * exp(-dt / TT_tau_s(V)));
- TT_r(g) = (TT_yinfty_r(V) - (TT_yinfty_r(V) - TT_r(var)) * exp(-dt / TT_tau_r(V)));
- TT_d(g) =(TT_yinfty_d(V) - (TT_yinfty_d(V) - TT_d(var)) * exp(-dt / TT_tau_d(V)));
- TT_f(g) = (TT_yinfty_f(V) -
- (TT_yinfty_f(V) - TT_f(var)) * exp(-dt / TT_tau_f(V)));
- TT_f_2(g) = (TT_yinfty_f_2(V) - (TT_yinfty_f_2(V) - TT_f_2(var)) * exp(-dt / TT_tau_f_2(V)));
+ (TT_yinfty_x_r2(V) - TT_x_r2(var)) * exp(-dt / TT_tau_x_r2(V)));
+ TT_x_s(g) = (TT_yinfty_x_s(V) -
+ (TT_yinfty_x_s(V) - TT_x_s(var)) * exp(-dt / TT_tau_x_s(V)));
+ TT_s(g) =
+ (TT_yinfty_s(V) - (TT_yinfty_s(V) - TT_s(var)) * exp(-dt / TT_tau_s(V)));
+ TT_r(g) =
+ (TT_yinfty_r(V) - (TT_yinfty_r(V) - TT_r(var)) * exp(-dt / TT_tau_r(V)));
+ TT_d(g) =
+ (TT_yinfty_d(V) - (TT_yinfty_d(V) - TT_d(var)) * exp(-dt / TT_tau_d(V)));
+ TT_f(g) =
+ (TT_yinfty_f(V) - (TT_yinfty_f(V) - TT_f(var)) * exp(-dt / TT_tau_f(V)));
+ TT_f_2(g) = (TT_yinfty_f_2(V) -
+ (TT_yinfty_f_2(V) - TT_f_2(var)) * exp(-dt / TT_tau_f_2(V)));
auto caSS = TT_Ca_SS(var);
- TT_f_Ca(g) =
- TT_yinfty_f_Ca(caSS) - (TT_yinfty_f_Ca(caSS) - TT_f_Ca(var)) * exp(-dt / TT_tau_f_Ca(caSS));
+ TT_f_Ca(g) = TT_yinfty_f_Ca(caSS) - (TT_yinfty_f_Ca(caSS) - TT_f_Ca(var)) *
+ exp(-dt / TT_tau_f_Ca(caSS));
}
-
-std::vector<double> MTenTusscher::calculateCurrents(const vector<double> &var) const {
+std::vector<double> MTenTusscher::calculateCurrents(
+ const std::vector<double> &var) const {
double V = TT_Pot(var);
// reverse potentials
@@ -180,8 +188,7 @@ double MTenTusscher::oQuer(const std::vector<double> &var) const {
(k_3 + k_1 * TT_Ca_SS(var) * TT_Ca_SS(var));
}
-
-vector<double> MTenTusscher::Tau(const vector<double> &var) const {
+std::vector<double> MTenTusscher::Tau(const std::vector<double> &var) const {
auto V = TT_Pot(var);
return {
TT_tau_m(V),
@@ -197,10 +204,9 @@ vector<double> MTenTusscher::Tau(const vector<double> &var) const {
TT_tau_f_2(V) /*,
TT_tau_f_Ca(CaSS)*/
};
-
}
-vector<double> MTenTusscher::YInfty(const vector<double> &var) const {
+std::vector<double> MTenTusscher::YInfty(const std::vector<double> &var) const {
auto V = TT_Pot(var);
return {
TT_yinfty_m(V),
@@ -217,7 +223,8 @@ vector<double> MTenTusscher::YInfty(const vector<double> &var) const {
TT_yinfty_f_Ca(CaSS)*/ };
}
-void MTenTusscher::IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+void MTenTusscher::IIonUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const {
auto currents = calculateCurrents(var);
TT_Pot(g)=- TT_I_Na(currents)
- TT_I_CaL(currents)
@@ -236,8 +243,9 @@ void MTenTusscher::IIonUpdate(const vector<double> &var, vector<double> &g, doub
}
}
-void
-MTenTusscher::ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+void MTenTusscher::ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g,
+ double iExt) const {
auto ion = calculateCurrents(var);
double I_leak = V_leak * (TT_Ca_SR(var) - TT_Ca(var)); // x
double I_up = V_maxup / (1.00000 + ((K_up * K_up) / (TT_Ca(var) * TT_Ca(var)))); //x
@@ -261,8 +269,10 @@ MTenTusscher::ConcentrationUpdate(const vector<double> &var, vector<double> &g,
(F * V_c) * Capacitance; // x
}
-void MTenTusscher::GatingUpdate(const vector<double> &var, vector<double> &g) const {
- if (useExpInt) return;
+void MTenTusscher::GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const {
+ if (useExpInt)
+ return;
double V = TT_Pot(var);
TT_m(g) = (TT_yinfty_m(V) - TT_m(var)) / TT_tau_m(V);
@@ -276,28 +286,28 @@ void MTenTusscher::GatingUpdate(const vector<double> &var, vector<double> &g) co
TT_d(g) = (TT_yinfty_d(V) - TT_d(var)) / TT_tau_d(V);
TT_f(g) = (TT_yinfty_f(V) - TT_f(var)) / TT_tau_f(V);
TT_f_2(g) = (TT_yinfty_f_2(V) - TT_f_2(var)) / TT_tau_f_2(V);
- TT_f_Ca(g) = (TT_yinfty_f_Ca(TT_Ca_SS(var)) - TT_f_Ca(var)) / TT_tau_f_Ca(TT_Ca_SS(var));
-
+ TT_f_Ca(g) = (TT_yinfty_f_Ca(TT_Ca_SS(var)) - TT_f_Ca(var)) /
+ TT_tau_f_Ca(TT_Ca_SS(var));
double kcasr =
- max_sr - ((max_sr - min_sr) / (1.00000 + (EC / TT_Ca_SR(var)) * (EC / TT_Ca_SR(var)))); // x
- double k_2 = k_2prime * kcasr; // x
- TT_R_quer(g) =
- k_4 * (1 - TT_R_quer(var)) - k_2 * TT_Ca_SS(var) * TT_R_quer(var);//von tentusscher // x
+ max_sr - ((max_sr - min_sr) /
+ (1.00000 + (EC / TT_Ca_SR(var)) * (EC / TT_Ca_SR(var)))); // x
+ double k_2 = k_2prime * kcasr; // x
+ TT_R_quer(g) = k_4 * (1 - TT_R_quer(var)) -
+ k_2 * TT_Ca_SS(var) * TT_R_quer(var); // von tentusscher // x
/* Vorher:
- * double alpha_quer = k_4;
- * double beta_quer = k_2 * Ca_SS(var);
- * double y_infty_R = alpha_quer / (alpha_quer + beta_quer);
- * double TT_tau_R = 1 / (alpha_quer + beta_quer);
- * return y_infty_R + (R_quer - y_infty_R) * exp (-dt* TimeScale() * 1/TT_tau_R);
- */
+ * double alpha_quer = k_4;
+ * double beta_quer = k_2 * Ca_SS(var);
+ * double y_infty_R = alpha_quer / (alpha_quer + beta_quer);
+ * double TT_tau_R = 1 / (alpha_quer + beta_quer);
+ * return y_infty_R + (R_quer - y_infty_R) * exp (-dt* TimeScale() *
+ * 1/TT_tau_R);
+ */
}
-
-void MTenTusscher::UpdateValues(double dt, vector<double> &var, const vector<double> &g) const {
-
-
+void MTenTusscher::UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const {
if (useExpInt) {
TT_R_quer(var) = TT_R_quer(g);
} else {
@@ -359,7 +369,7 @@ void MTenTusscher::UpdateValues(double dt, vector<double> &var, const vector<dou
TT_Pot(var) += dt * TT_Pot(g);
}
-double MTenTusscher::GetIion(const vector<double> &var) const {
+double MTenTusscher::GetIion(const std::vector<double> &var) const {
auto ionCurrents = calculateCurrents(var);
return std::accumulate(ionCurrents.begin(), ionCurrents.end(), 0.0);
}
diff --git a/src/cellmodels/MTenTusscher.hpp b/src/cellmodels/MTenTusscher.hpp
index 02f1ed56ce52a8ef148b2b33d80b6e2374db4abf..709cf6aac9ff47b9714ee49ec46d7fabddbd34df 100644
--- a/src/cellmodels/MTenTusscher.hpp
+++ b/src/cellmodels/MTenTusscher.hpp
@@ -97,9 +97,9 @@ class MTenTusscher : public MElphyModel {
double initialPotential() const override { return -85.23; } //-86.2//-85.423
- double oQuer(const vector<double> &var) const;
+ double oQuer(const std::vector<double> &var) const;
- std::vector<double> calculateCurrents(const vector<double> &var) const;
+ std::vector<double> calculateCurrents(const std::vector<double> &var) const;
public:
explicit MTenTusscher() : MElphyModel(19, 1000.0, 1e3) {
@@ -118,21 +118,24 @@ public:
void Initialize(std::vector<double> &vw) override;
- void IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void IIonUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double iExt) const override;
- void
- ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+ void ConcentrationUpdate(const std::vector<double> &var,
+ std::vector<double> &g, double iExt) const override;
- void GatingUpdate(const vector<double> &var, vector<double> &g) const override;
+ void GatingUpdate(const std::vector<double> &var,
+ std::vector<double> &g) const override;
- void ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const override;
+ void ExpIntegratorUpdate(double dt, const std::vector<double> &var,
+ std::vector<double> &g) const override;
+ std::vector<double> Tau(const std::vector<double> &var) const override;
- vector<double> Tau(const vector<double> &var) const override;
+ std::vector<double> YInfty(const std::vector<double> &var) const override;
- vector<double> YInfty(const vector<double> &var) const override;
-
- void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const;
+ void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const;
double GetIion(const std::vector<double> &var) const override;
};
diff --git a/src/cellmodels/solvers/MCellSolver.cpp b/src/cellmodels/solvers/MCellSolver.cpp
index b0e25041ebbd71e5d668742f4511c08390bdb899..5189d32d7561ef731b1098af665a5e52346d6caf 100644
--- a/src/cellmodels/solvers/MCellSolver.cpp
+++ b/src/cellmodels/solvers/MCellSolver.cpp
@@ -1,13 +1,15 @@
#include "MCellSolver.hpp"
-void ExplicitEulerStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void ExplicitEulerStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
cellModel.RHSUpdate(VW, g[0], cellModel.IExt(t));
cellModel.UpdateValues(dt, VW, g[0]);
}
-void RungeKutte2Step(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void RungeKutte2Step(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
std::vector<double> k(VW.size());
std::uninitialized_copy(VW.begin(), VW.end(), k.begin());
@@ -17,8 +19,8 @@ void RungeKutte2Step(MCellModel &cellModel, double t, double dt, vector<double>
cellModel.RHSUpdate(k, g[1], cellModel.IExt(t+0.5*dt));
cellModel.UpdateValues(dt, VW, g[1]);
}
-void HeunStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g){
+void HeunStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW, std::vector<std::vector<double>> &g) {
std::vector<double> k(VW.size());
std::uninitialized_copy(VW.begin(), VW.end(), k.begin());
@@ -32,8 +34,9 @@ void HeunStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
cellModel.UpdateValues(0.5*dt, VW, g[0]);
}
-void RungeKutta4Step(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void RungeKutta4Step(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
std::vector<double> k(VW.size());
std::uninitialized_copy(VW.begin(), VW.end(), k.begin());
diff --git a/src/cellmodels/solvers/MElphySolver.cpp b/src/cellmodels/solvers/MElphySolver.cpp
index c0f30225e48ecdcddd63e77854915bd4e303ff36..145d4d5d74b62a607ae81e1c954c936ae91ec39b 100644
--- a/src/cellmodels/solvers/MElphySolver.cpp
+++ b/src/cellmodels/solvers/MElphySolver.cpp
@@ -98,7 +98,7 @@ void MElphySolver::ComputeAndPrintEigenvalues(const std::vector<double> &VW) {
}
}
-void MElphySolver::fillTimeIndependant(vector<double> vector, Vectors &vectors, row row, int j) {
+void MElphySolver::fillTimeIndependant(std::vector<double> vector, Vectors &vectors, row row, int j) {
cellModel.SetExternal(Excitation(excData, row, j), Duration(excData, row, j),
Amplitude(excData, row, j));
@@ -118,7 +118,7 @@ void MElphySolver::Solve(TimeSeries µTS, Vector &u) {
u = V[0];
}
-void MElphySolver::Solve(TimeSeries µTS, vector<double> &u) {
+void MElphySolver::Solve(TimeSeries µTS, std::vector<double> &u) {
u.resize(microTS.Steps() + 1);
std::vector<double> vw(size);
@@ -172,7 +172,7 @@ void MElphySolver::PrintMaxMinGating() const {
vout(1) << endl;
}
-void MElphySolver::selectScheme(const string &solvingScheme) {
+void MElphySolver::selectScheme(const std::string &solvingScheme) {
if (getSplitting() == STEIN) {
steps = 2;
SolveStep = HeunStep;
@@ -214,7 +214,8 @@ void MElphySolver::SetExternal(std::vector<double> values) {
cellModel.SetExternal(values[0], values[1], values[2]);
}
-void MElphySolver::Solve(TimeSeries µTS, vector<std::vector<double>> &u) {
+void MElphySolver::Solve(TimeSeries µTS,
+ std::vector<std::vector<double>> &u) {
u.resize(microTS.Steps() + 1);
std::vector<double> vw(cellModel.Size() + (cellModel.Type() == TENSION));
@@ -224,7 +225,8 @@ void MElphySolver::Solve(TimeSeries µTS, vector<std::vector<double>> &u) {
}
cellModel.Initialize(vw);
if (cellModel.Type() == TENSION) {
- vw[cellModel.Size()] = (1.0 + vw[cellModel.ElphySize()]) * (1.0 + vw[cellModel.ElphySize()]);
+ vw[cellModel.Size()] =
+ (1.0 + vw[cellModel.ElphySize()]) * (1.0 + vw[cellModel.ElphySize()]);
}
u[0] = vw;
@@ -239,9 +241,9 @@ void MElphySolver::Solve(TimeSeries µTS, vector<std::vector<double>> &u) {
}
}
-
-void ExponentialIntegratorStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void ExponentialIntegratorStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
int position = cellModel.GatingIndex();
cellModel.ExpIntegratorUpdate(dt, VW, g[0]);
@@ -253,12 +255,11 @@ void ExponentialIntegratorStep(MCellModel &cellModel, double t, double dt, vecto
for (int i =1;i<position;i++){
VW[i] += dt * g[0][i];
}*/
-
-
}
-void EIWithNewtonStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void EIWithNewtonStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
double mVolt = VW[0];
int position = cellModel.GatingIndex();
@@ -272,8 +273,9 @@ void EIWithNewtonStep(MCellModel &cellModel, double t, double dt, vector<double>
}
}
-void EigenvalueStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g) {
+void EigenvalueStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g) {
double mVolt = VW[0];
double volt = mVolt / 1000.;
@@ -291,7 +293,6 @@ void EigenvalueStep(MCellModel &cellModel, double t, double dt, vector<double> &
}
}
-
/*void MEigenvalueSolver::ComputeAndPrintEigenvalues(std::vector<double> &VW) {
int n = cellModel.GatingDim()+1;
std::vector<std::vector<double>> J(n);
diff --git a/src/cellmodels/solvers/MElphySolver.hpp b/src/cellmodels/solvers/MElphySolver.hpp
index 4da76d8633362bba650924c9888ea3a587ceadad..38120d3ff3dadfbbc941f03d4063d2d7704ef524 100644
--- a/src/cellmodels/solvers/MElphySolver.hpp
+++ b/src/cellmodels/solvers/MElphySolver.hpp
@@ -10,8 +10,9 @@ void ExponentialIntegratorStep(MCellModel &cellModel, double t, double dt, std::
void EigenvalueStep(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
std::vector<std::vector<double>> &g);
-void EIWithNewtonStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
- vector<std::vector<double>> &g);
+void EIWithNewtonStep(MCellModel &cellModel, double t, double dt,
+ std::vector<double> &VW,
+ std::vector<std::vector<double>> &g);
class MElphySolver : public MCellSolver {
MCellModel &cellModel;
@@ -23,22 +24,23 @@ class MElphySolver : public MCellSolver {
- void selectScheme(const string &solvingScheme);
+ void selectScheme(const std::string &solvingScheme);
std::function<void(MCellModel &, double, double, std::vector<double> &,
std::vector<std::vector<double>> &)> SolveStep;
- void fillTimeIndependant(vector<double> vector, Vectors &vectors, row row, int j = 0);
+ void fillTimeIndependant(std::vector<double> vector, Vectors &vectors, row row, int j = 0);
public :
explicit MElphySolver(MCellModel &M, const Vector &data) : MCellSolver(M.Size()),
cellModel(M), excData(data) {
- string solvingScheme{"ExplicitEuler"};
+ std::string solvingScheme{"ExplicitEuler"};
Config::Get("CellModelScheme", solvingScheme);
selectScheme(solvingScheme);
}
- explicit MElphySolver(MCellModel &M, const Vector &data, const string &solvingScheme)
+ explicit MElphySolver(MCellModel &M, const Vector &data,
+ const std::string &solvingScheme)
: MCellSolver(M.Size() + 1), cellModel(M), excData(data) {
selectScheme(solvingScheme);
}
diff --git a/src/cellmodels/solvers/MTensionSolver.cpp b/src/cellmodels/solvers/MTensionSolver.cpp
index ebeb383d870fa98002dfb8e3b627fd79013ad75a..37f43d3218a0565335073edadde0f0b9591ee30f 100644
--- a/src/cellmodels/solvers/MTensionSolver.cpp
+++ b/src/cellmodels/solvers/MTensionSolver.cpp
@@ -56,9 +56,7 @@ void MTensionSolver::Solve(TimeSeries µTS, Vectors &vecs) {
}
}
-void MTensionSolver::Solve(TimeSeries µTS, vector<double> &u) {
-
-}
+void MTensionSolver::Solve(TimeSeries µTS, std::vector<double> &u) {}
void MTensionSolver::initialize(Vector &u) {
gating = std::make_unique<Vectors>(cellModel.Size(), u);
@@ -67,7 +65,7 @@ void MTensionSolver::initialize(Vector &u) {
init = true;
}
-void MTensionSolver::selectScheme(const string &solvingScheme) {
+void MTensionSolver::selectScheme(const std::string &solvingScheme) {
if (solvingScheme == "ExplicitEuler") {
steps = 1;
SolveStep = ExplicitEulerStep;
@@ -80,10 +78,10 @@ void MTensionSolver::selectScheme(const string &solvingScheme) {
} else {
THROW(("Scheme \"" + solvingScheme + "\" not implemented").c_str())
}
-
}
-void MTensionSolver::PlotEntries(const vector<double> &entries, row row, int j, int step) const {
+void MTensionSolver::PlotEntries(const std::vector<double> &entries, row row,
+ int j, int step) const {
if (row() == Origin) {
mout.PrintInfo("CellModelValues", 1,
PrintInfoEntry("Calcium", entries[cellModel.ElphyModel().CalciumIndex()]),
diff --git a/src/cellmodels/solvers/MTensionSolver.hpp b/src/cellmodels/solvers/MTensionSolver.hpp
index a43764e496d4181b8389759b36dc40c6e2bd190a..97f9ed450b44256fa32f59838685ed94ce701d13 100644
--- a/src/cellmodels/solvers/MTensionSolver.hpp
+++ b/src/cellmodels/solvers/MTensionSolver.hpp
@@ -11,7 +11,7 @@ class MTensionSolver : public MCellSolver {
void initialize(Vector &u);
- void selectScheme(const string &solvingScheme);
+ void selectScheme(const std::string &solvingScheme);
std::function<void(MTensionModel &, double, double, std::vector<double> &,
std::vector<std::vector<double>> &)> SolveStep;
@@ -22,13 +22,14 @@ class MTensionSolver : public MCellSolver {
public :
explicit MTensionSolver(MTensionModel &M, const Vector &data) : MCellSolver(M.Size() + 1),
cellModel(M), excData(data) {
- string solvingScheme{"ExplicitEuler"};
+ std::string solvingScheme{"ExplicitEuler"};
Config::Get("CellModelScheme", solvingScheme);
selectScheme(solvingScheme);
Config::Get("IextInPDE", iextInPDE);
}
- explicit MTensionSolver(MTensionModel &M, const Vector &data, const string &solvingScheme)
+ explicit MTensionSolver(MTensionModel &M, const Vector &data,
+ const std::string &solvingScheme)
: MCellSolver(M.Size() + 1), cellModel(M), excData(data) {
selectScheme(solvingScheme);
Config::Get("IextInPDE", iextInPDE);
@@ -42,7 +43,8 @@ public :
void Solve(TimeSeries µTS, std::vector<double> &u) override;
- void PlotEntries(const vector<double> &vector, row row, int j, int step) const;
+ void PlotEntries(const std::vector<double> &vector, row row, int j,
+ int step) const;
void PlotResult(const Vectors &vectors, double beginT, double endT) const;
diff --git a/src/coupled/MainCoupled.cpp b/src/coupled/MainCoupled.cpp
index 8215a43c44eea7e85dc13e811ce45ab70d7cd243..29e9becc5082fb869115b855382c1e3f71922708 100644
--- a/src/coupled/MainCoupled.cpp
+++ b/src/coupled/MainCoupled.cpp
@@ -18,7 +18,7 @@ void MainCoupled::Initialize() {
}
Vector &MainCoupled::Run() {
- string eModelName = "";
+ std::string eModelName = "";
Config::Get("ElphyModel", eModelName);
// Load ElphySolver from src/electrophysiology/solvers/ElphySolver.cpp
@@ -129,7 +129,7 @@ std::vector<double> MainCoupled::Evaluate(const Vector &solution) {
std::string elphyModelClass = "IBTElphyModel";
Config::Get("ElphyModelClass", elphyModelClass);
- string cmName = "FitzHughNagumo";
+ std::string cmName = "FitzHughNagumo";
if (elphyModelClass == "IBTElphyModel") {
Config::Get("IBTVentricleModel", cmName);
} else {
@@ -210,7 +210,7 @@ std::vector<double> MainCoupled::Evaluate(const Vector &solution) {
SetDisplacement(true, exactSolution);
solution.GetMesh().PrintInfo();
- string m = "Passive";
+ std::string m = "Passive";
Config::Get("MechDiscretization", m);
mout.PrintInfo("Assemble", verbose,
diff --git a/src/coupled/MainCoupled.hpp b/src/coupled/MainCoupled.hpp
index 42f5e25da6ef5594c36bd55a746039c571251d2e..0c3b66e60d2426fd83c42e668613cfa8f77d5afe 100644
--- a/src/coupled/MainCoupled.hpp
+++ b/src/coupled/MainCoupled.hpp
@@ -61,7 +61,8 @@ public:
void SetDisplacement(bool exact, Vector &u) const;
- double EvaluateQuantity(const Vector &solution, const string &quantity) const;
+ double EvaluateQuantity(const Vector &solution,
+ const std::string &quantity) const;
};
#endif //MAINCOUPLED_H
diff --git a/src/coupled/problem/CoupledBenchmarkProblems.hpp b/src/coupled/problem/CoupledBenchmarkProblems.hpp
index b016bb0d4b5bb3e823c6fc7007493e92926dc3b0..c0908dbad6397103835ca08943e8da6f6febd9ad 100644
--- a/src/coupled/problem/CoupledBenchmarkProblems.hpp
+++ b/src/coupled/problem/CoupledBenchmarkProblems.hpp
@@ -9,8 +9,6 @@ public:
defaultMeshName("BenchmarkBeam3DTet");
}
- string Name() const override {
- return "Coupled Beam Problem";
- }
+ std::string Name() const override { return "Coupled Beam Problem"; }
};
#endif //COUPLEDBENCHMARKPROBLEMS_HPP
diff --git a/src/coupled/problem/CoupledProblem.cpp b/src/coupled/problem/CoupledProblem.cpp
index bfd47c218d11908cccc9f2b301a7781ae061e223..168f40fe435e45e4ab6279963374f938bcd2bf20 100644
--- a/src/coupled/problem/CoupledProblem.cpp
+++ b/src/coupled/problem/CoupledProblem.cpp
@@ -11,7 +11,8 @@ std::unique_ptr<CoupledProblem> GetCoupledProblem() {
return GetCoupledProblem(problemName);
}
-std::unique_ptr<CoupledProblem> undeformedCoupledProblem(const string &problemName) {
+std::unique_ptr<CoupledProblem> undeformedCoupledProblem(
+ const std::string &problemName) {
if (contains(problemName, "ContractingBeam")) {
return std::make_unique<CoupledBeamProblem>();
}
@@ -33,7 +34,8 @@ std::unique_ptr<CoupledProblem> undeformedCoupledProblem(const string &problemNa
return std::make_unique<CompositeCoupledProblem>(std::move(elphyProblem), std::move(mechProblem));
}
-std::unique_ptr<CoupledProblem> GetCoupledProblem(const string &problemName) {
+std::unique_ptr<CoupledProblem> GetCoupledProblem(
+ const std::string &problemName) {
if (contains(problemName, "Deformed")) {
return std::make_unique<DeformedCoupledProblem>(undeformedCoupledProblem(problemName));
}
@@ -48,7 +50,7 @@ std::unique_ptr<CoupledProblem> GetCoupledProblem(const string &problemName) {
#define mechMeshes this->ElasticityProblem::meshes
#define mechGeo this->ElasticityProblem::cGeo
-void CoupledProblem::defaultMeshName(const string &name) {
+void CoupledProblem::defaultMeshName(const std::string &name) {
if (!ElasticityProblem::meshFromConfig()) {
elphyMeshName = name;
mechMeshName = name;
diff --git a/src/coupled/problem/CoupledProblem.hpp b/src/coupled/problem/CoupledProblem.hpp
index e12614e1d67fac1c70c5c6362241dca73424d868..9549e1fdc0393165afcead85f68bcb8c81b6d710 100644
--- a/src/coupled/problem/CoupledProblem.hpp
+++ b/src/coupled/problem/CoupledProblem.hpp
@@ -94,35 +94,41 @@ public:
virtual std::string Name() const override = 0;
virtual void InitializeEvaluationPoints() override;
- string Evaluate(const Vector &solution) const override {
+ std::string Evaluate(const Vector &solution) const override {
return EvaluateElphy(solution) + "\n" + EvaluateMech(solution);
}
- vector<double> EvaluationResults(const Vector &solution) const override {
+ std::vector<double> EvaluationResults(const Vector &solution) const override {
auto elphyResults = EvaluationResultsElphy(solution);
auto mechResults = EvaluationResultsMech(solution);
return mechResults;
}
- vector<std::string> EvaluationMetrics() const override {
+ std::vector<std::string> EvaluationMetrics() const override {
auto elphyMetrics = EvaluationMetricsElphy();
auto mechMetrics = EvaluationMetricsMech();
return mechMetrics;
}
- virtual string EvaluateElphy(const Vector &solution) const { return ""; }
+ virtual std::string EvaluateElphy(const Vector &solution) const { return ""; }
- virtual string EvaluateMech(const Vector &solution) const { return ""; }
+ virtual std::string EvaluateMech(const Vector &solution) const { return ""; }
- virtual vector<double> EvaluationResultsElphy(const Vector &solution) const { return {}; }
+ virtual std::vector<double> EvaluationResultsElphy(
+ const Vector &solution) const {
+ return {};
+ }
- virtual vector<double> EvaluationResultsMech(const Vector &solution) const { return {}; }
+ virtual std::vector<double> EvaluationResultsMech(
+ const Vector &solution) const {
+ return {};
+ }
- virtual vector<std::string> EvaluationMetricsElphy() const { return {}; }
+ virtual std::vector<std::string> EvaluationMetricsElphy() const { return {}; }
- virtual vector<std::string> EvaluationMetricsMech() const { return {}; }
+ virtual std::vector<std::string> EvaluationMetricsMech() const { return {}; }
virtual ~CoupledProblem() = default;
};
@@ -144,13 +150,11 @@ public:
}
}
- string Name() const override {
+ std::string Name() const override {
return "Composite Coupled Problem [" + elphyProblem->Name() + " | " + mechProblem->Name() + "]";
}
- const string &Model() const override {
- return mechProblem->Model();
- }
+ const std::string &Model() const override { return mechProblem->Model(); }
bool IsStatic() const override {
return mechProblem->IsStatic();
@@ -241,15 +245,15 @@ public:
return elphyProblem->Conductivity(c);
}
- vector<double> EvaluationResults(const Vector &solution) const override {
+ std::vector<double> EvaluationResults(const Vector &solution) const override {
return elphyProblem->EvaluationResults(solution);
}
- virtual string EvaluateElphy(const Vector &solution) const {
+ virtual std::string EvaluateElphy(const Vector &solution) const {
return elphyProblem->Evaluate(solution);
}
- virtual string EvaluateMech(const Vector &solution) const {
+ virtual std::string EvaluateMech(const Vector &solution) const {
return mechProblem->Evaluate(solution);
}
@@ -261,12 +265,13 @@ public:
return elphyProblem->IonCurrent(time, x);
}
- void chamberSpecialization(const string &chamberName) const override {
+ void chamberSpecialization(const std::string &chamberName) const override {
elphyProblem->chamberSpecialization(chamberName);
}
};
-std::unique_ptr<CoupledProblem> GetCoupledProblem(const string &problemName);
+std::unique_ptr<CoupledProblem> GetCoupledProblem(
+ const std::string &problemName);
std::unique_ptr<CoupledProblem> GetCoupledProblem();
diff --git a/src/coupled/problem/DeformedProblem.hpp b/src/coupled/problem/DeformedProblem.hpp
index f08c98e0090c21e7f1a920ec121ed9df462b2bd7..e633bb5bf03d383a29c1e4e76d5ed21afa235a39 100644
--- a/src/coupled/problem/DeformedProblem.hpp
+++ b/src/coupled/problem/DeformedProblem.hpp
@@ -20,49 +20,49 @@ public:
InitializeEvaluationPoints();
Config::Get("deformFibersInConductivity", withDeformFibersInConductivity);
}
- string Name() const override {
+ std::string Name() const override {
return originalProblem->Name() + " (Deformed)";
}
- string Evaluate(const Vector &solution) const override {
+ std::string Evaluate(const Vector &solution) const override {
return originalProblem->Evaluate(solution);
}
- vector<double> EvaluationResults(const Vector &solution) const override {
+ std::vector<double> EvaluationResults(const Vector &solution) const override {
return originalProblem->EvaluationResults(solution);
}
- vector<std::string> EvaluationMetrics() const override {
+ std::vector<std::string> EvaluationMetrics() const override {
return originalProblem->EvaluationMetrics();
}
- string EvaluateElphy(const Vector &solution) const override {
+ std::string EvaluateElphy(const Vector &solution) const override {
return originalProblem->EvaluateElphy(solution);
}
- string EvaluateMech(const Vector &solution) const override {
+ std::string EvaluateMech(const Vector &solution) const override {
return originalProblem->EvaluateMech(solution);
}
- vector<double> EvaluationResultsElphy(const Vector &solution) const override {
+ std::vector<double> EvaluationResultsElphy(
+ const Vector &solution) const override {
return originalProblem->EvaluationResultsElphy(solution);
}
- vector<double> EvaluationResultsMech(const Vector &solution) const override {
+ std::vector<double> EvaluationResultsMech(
+ const Vector &solution) const override {
return originalProblem->EvaluationResultsMech(solution);
}
- vector<std::string> EvaluationMetricsElphy() const override {
+ std::vector<std::string> EvaluationMetricsElphy() const override {
return originalProblem->EvaluationMetricsElphy();
}
- vector<std::string> EvaluationMetricsMech() const override {
+ std::vector<std::string> EvaluationMetricsMech() const override {
return originalProblem->EvaluationMetricsMech();
}
- const string &Model() const override {
- return originalProblem->Model();
- }
+ const std::string &Model() const override { return originalProblem->Model(); }
bool IsStatic() const override {
return originalProblem->IsStatic();
@@ -211,7 +211,7 @@ public:
return det(F) * FInv * D * transpose(FInv);
}
- void chamberSpecialization(const string &chamberName) const override {
+ void chamberSpecialization(const std::string &chamberName) const override {
originalProblem->chamberSpecialization(chamberName);
}
};
diff --git a/src/coupled/problem/PhysiologicalProblems.cpp b/src/coupled/problem/PhysiologicalProblems.cpp
index e5864c9b144ff0ae79a9e4145bb7058bc04dceae..12fe100e71cc90af27da436a30cb9f52740873d9 100644
--- a/src/coupled/problem/PhysiologicalProblems.cpp
+++ b/src/coupled/problem/PhysiologicalProblems.cpp
@@ -5,8 +5,8 @@
#include "VentriclePVLoops.hpp"
-
-vector<double> ContractingVentricleProblem::EvaluationResults(const Vector &solution) const {
+std::vector<double> ContractingVentricleProblem::EvaluationResults(
+ const Vector &solution) const {
auto apexPot = findPotential(solution, {ventricleApex});
mout.PrintInfo("LeftVentricle", coupledVerbose(),
diff --git a/src/coupled/problem/PhysiologicalProblems.hpp b/src/coupled/problem/PhysiologicalProblems.hpp
index 1445056e5ff6abf25fb5aa0e24a08a7706a3e591..4c5cb4ac39b4954fdd56efa41fd30fad9a93097e 100644
--- a/src/coupled/problem/PhysiologicalProblems.hpp
+++ b/src/coupled/problem/PhysiologicalProblems.hpp
@@ -29,21 +29,18 @@ public:
ElphyProblem::meshesName = ventricleProblem.GetMeshName();
}
- ContractingVentricleProblem(const string &matName, const vector<double> &matPar) :
- ventricleProblem(matName, matPar) {
+ ContractingVentricleProblem(const std::string &matName,
+ const std::vector<double> &matPar)
+ : ventricleProblem(matName, matPar) {
Config::Get("TractionK", k);
Config::Get("TractionC", c);
ElasticityProblem::meshesName = ventricleProblem.GetMeshName();
ElphyProblem::meshesName = ventricleProblem.GetMeshName();
}
- string Name() const override {
- return "Contracting Ventricle";
- }
+ std::string Name() const override { return "Contracting Ventricle"; }
- const string &Model() const override {
- return ventricleProblem.Model();
- }
+ const std::string &Model() const override { return ventricleProblem.Model(); }
bool IsStatic() const override {
return false;
@@ -77,7 +74,7 @@ public:
double TractionViscosity(double t, const Point &x, int sd) const override { return c; };
- vector<double> EvaluationResults(const Vector &solution) const;
+ std::vector<double> EvaluationResults(const Vector &solution) const;
};
class BiVentricleProblem : public CoupledProblem {
@@ -126,7 +123,7 @@ public:
Config::Get("pressureES_RV", pressureES[1]);
}
- string Name() const override {
+ std::string Name() const override {
return "BiVentricular Problem on " + ElasticityProblem::GetMeshName();
}
diff --git a/src/coupled/solvers/CellModelSolver.hpp b/src/coupled/solvers/CellModelSolver.hpp
index fd15ae63309cb04716b274d9c375c7ead3b53cd5..1dc4c9359a44ebd14d56c61e4af48121ef881194 100644
--- a/src/coupled/solvers/CellModelSolver.hpp
+++ b/src/coupled/solvers/CellModelSolver.hpp
@@ -11,7 +11,7 @@ class CellModelSolver : public CardiacSolver {
int plottingsteps;
void setPlottingSteps(const std::string& configString = "PlottingSteps"){
Config::Get(configString, plottingsteps);
- plottingsteps = min(plottingsteps, 1);
+ plottingsteps = std::min(plottingsteps, 1);
}
public:
explicit CellModelSolver(IElphyAssemble &A, ElphySolver &eSolver) : elphyA(A), elphySolver(eSolver) {
diff --git a/src/coupled/transfers/MultiPartTransfer.cpp b/src/coupled/transfers/MultiPartTransfer.cpp
index 2a3a7e29fd24d99024f873c8cf5963fe3856e1c4..c81ce50b25372bf081e128ebb3d010ea25679e96 100644
--- a/src/coupled/transfers/MultiPartTransfer.cpp
+++ b/src/coupled/transfers/MultiPartTransfer.cpp
@@ -55,7 +55,7 @@ void MultiPartLagrangeTransfer::constructRowList(
recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
for (int l = 1; l <= dlevel; ++l) {
for (const auto &parent: recursiveCells[l - 1]) {
- const vector<Rule> &R = parent->Refine(childrenVerteces);
+ const std::vector<Rule> &R = parent->Refine(childrenVerteces);
for (int i = 0; i < R.size(); ++i) {
R[i](childrenVerteces, childVerteces);
recursiveCells[l].emplace_back(
@@ -66,7 +66,7 @@ void MultiPartLagrangeTransfer::constructRowList(
}
for (const auto &child: recursiveCells[dlevel]) {
LagrangeNodalPoints(*child, coarseNodalPoints, 1);
- vector<int> coarseIds = mpp::nodalIds(coarse, coarseNodalPoints);
+ std::vector<int> coarseIds = mpp::nodalIds(coarse, coarseNodalPoints);
mpp::linearTransfer(child->ReferenceType(), coarseNodalPoints, coarseIds, lvl, fine, rowList);
}
}
diff --git a/src/elasticity/MainElasticity.cpp b/src/elasticity/MainElasticity.cpp
index 920d99f6836e12eec2d612130749be713bd4d549..bad17c051dce7f649ebef055b39b3ee0fd7a640b 100644
--- a/src/elasticity/MainElasticity.cpp
+++ b/src/elasticity/MainElasticity.cpp
@@ -110,7 +110,7 @@ Vector &MainElasticity::Run() {
return *displacement;
}
-void MainElasticity::generateMesh(const string &model) const {
+void MainElasticity::generateMesh(const std::string &model) const {
// search for included chambers (src/CardiacData.cpp)
std::vector<int> chambers = getIncludedChambers(chamber);
auto mCreator = MeshesCreator().
diff --git a/src/elasticity/MainElasticity.hpp b/src/elasticity/MainElasticity.hpp
index 32d8f189bc37e8e87291290cdafd6425208a2f38..e17c4e066a3a072c940c6e1e1de5b3caccf1ec08 100644
--- a/src/elasticity/MainElasticity.hpp
+++ b/src/elasticity/MainElasticity.hpp
@@ -28,8 +28,7 @@ class MainElasticity : public MainCardMech {
friend class MultiElasticity;
friend class MultiDiscElasticity;
- void generateMesh(const string &model) const;
-
+ void generateMesh(const std::string &model) const;
int discDegree{1};
@@ -92,7 +91,8 @@ public:
/// Creates vector- and scalar-valued discretizations and generates a displacement vector.
void generateStartVectors();
/// Evaluates some error measures for the solution.
- double EvaluateQuantity(const Vector &solution, const string &quantity) const;
+ double EvaluateQuantity(const Vector &solution,
+ const std::string &quantity) const;
std::array<double, 4> GetVolume(const Vector &solution) const;
diff --git a/src/elasticity/MultiDiscElasticity.cpp b/src/elasticity/MultiDiscElasticity.cpp
index 24603cbfca46754cdf2857fb2a7cd39711651e89..924e01238b61ceea4aa77d5b64b421ec7a91544c 100644
--- a/src/elasticity/MultiDiscElasticity.cpp
+++ b/src/elasticity/MultiDiscElasticity.cpp
@@ -86,8 +86,7 @@ Vector &MultiDiscElasticity::Run() {
return solutions[discs.size() - 1]();
}
-vector<double> MultiDiscElasticity::Evaluate(const Vector &solution) {
+std::vector<double> MultiDiscElasticity::Evaluate(const Vector &solution) {
mout << endl << "=============================================================" << endl << endl;
return {};
}
-
diff --git a/src/elasticity/MultiDiscElasticity.hpp b/src/elasticity/MultiDiscElasticity.hpp
index b12cc4ecfcb2955a125f92e667e942c1e032f70c..18836d80f9c90dc1dd96ce77fd05bcae35c4149e 100644
--- a/src/elasticity/MultiDiscElasticity.hpp
+++ b/src/elasticity/MultiDiscElasticity.hpp
@@ -25,9 +25,9 @@ class MultiDiscElasticity : public MainCardMech {
bool useDifferentMeshes{false};
bool useDifferentDegs{false};
std::vector<VectorStruct> solutions{};
- std::vector<string> discs{};
- std::vector<string> meshes;
- std::vector<string> meshesList;
+ std::vector<std::string> discs{};
+ std::vector<std::string> meshes;
+ std::vector<std::string> meshesList;
std::vector<int> degrees;
std::vector<int> degreeList;
std::vector<int> DGPenalties;
@@ -40,7 +40,6 @@ public:
Vector &Run() override;
- vector<double> Evaluate(const Vector &solution) override;
-
+ std::vector<double> Evaluate(const Vector &solution) override;
};
diff --git a/src/elasticity/MultiElasticity.cpp b/src/elasticity/MultiElasticity.cpp
index b1d1645c1ccf8973780c6f607bd859bc3d1b3b6e..000dacde0aedbf8bbb60fea3a4b9b4c6a9a7780c 100644
--- a/src/elasticity/MultiElasticity.cpp
+++ b/src/elasticity/MultiElasticity.cpp
@@ -31,11 +31,10 @@ Vector &MultiElasticity::Run() {
return solutions[referenceLevel]();
}
-vector<double> MultiElasticity::Evaluate(const Vector &solution) {
+std::vector<double> MultiElasticity::Evaluate(const Vector &solution) {
for (int level = initialLevel; level < referenceLevel; ++level) {
main.EvaluateReferece(solutions[level](), solution);
}
mout << endl << "=============================================================" << endl << endl;
return {};
}
-
diff --git a/src/elasticity/MultiElasticity.hpp b/src/elasticity/MultiElasticity.hpp
index 88407e999867a13251670cc00f6758fe279c1668..5a692278dee7b598e7dde663a26bdc6517ad76f9 100644
--- a/src/elasticity/MultiElasticity.hpp
+++ b/src/elasticity/MultiElasticity.hpp
@@ -24,8 +24,7 @@ public:
Vector &Run() override;
- vector<double> Evaluate(const Vector &solution) override;
-
+ std::vector<double> Evaluate(const Vector &solution) override;
};
#endif //MULTIELASTICITY_HPP
diff --git a/src/elasticity/assemble/DGElasticity.cpp b/src/elasticity/assemble/DGElasticity.cpp
index ab9428a22cb218db3ce73129999623e2629a7cd7..7c01614f2b9a105a0edf41bab9c44caeee7e24ef 100644
--- a/src/elasticity/assemble/DGElasticity.cpp
+++ b/src/elasticity/assemble/DGElasticity.cpp
@@ -125,8 +125,8 @@ void DGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
}
}
for (int f = 0; f < c.Faces(); ++f) {
- Scalar s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
- U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
DGVectorFieldFaceElement FE(U, *c, f);
if (E.Bnd(f) >= PRESSURE_BC_LV && E.Bnd(f) <= PRESSURE_BC_RA) {// Neumann-Randbedingung
@@ -270,8 +270,8 @@ void DGElasticity::DirichletBoundary(const cell &c, int face, int bc, const Vect
DGVectorFieldFaceElement FE(U, *c, face);
DGSizedRowBndValues r_c(r, *c, E.shape_size());
const auto &cellMat = eProblem.GetMaterial(c);
- Scalar s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
- U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
Tensor Q = OrientationMatrix(c.GetData());
for (int q = 0; q < FE.nQ(); ++q) {
@@ -302,8 +302,8 @@ void DGElasticity::FixationBoundary(const cell &c, int face, int bc, const Vecto
DGVectorFieldFaceElement FE(U, *c, face);
DGSizedRowBndValues r_c(r, *c, E.shape_size());
const auto &cellMat = eProblem.GetMaterial(c);
- Scalar s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
- U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
Tensor Q = OrientationMatrix(c.GetData());
for (int q = 0; q < FE.nQ(); ++q) {
@@ -431,8 +431,8 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
}
}
for (int face = 0; face < c.Faces(); ++face) {
- Scalar s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
- U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
DGVectorFieldFaceElement FE(U, *c, face);
if (E.Bnd(face) == DIRICHLET_BC || E.Bnd(face) == DIRICHLET_FIX_BC) {
for (int q = 0; q < FE.nQ(); ++q) {
@@ -580,7 +580,7 @@ void DGElasticity::Project(const Vector &fine, Vector &u) const {
recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
for (int l = 1; l <= dlevel; ++l) {
for (const auto &parent: recursiveCells[l - 1]) {
- const vector<Rule> &R = parent->Refine(childrenVerteces);
+ const std::vector<Rule> &R = parent->Refine(childrenVerteces);
for (int i = 0; i < R.size(); ++i) {
R[i](childrenVerteces, childVerteces);
recursiveCells[l].emplace_back(
@@ -627,7 +627,7 @@ void DGElasticity::Interpolate(const Vector &coarse, Vector &u) const {
recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
for (int l = 1; l <= dlevel; ++l) {
for (const auto &parent: recursiveCells[l - 1]) {
- const vector<Rule> &R = parent->Refine(childrenVerteces);
+ const std::vector<Rule> &R = parent->Refine(childrenVerteces);
for (int i = 0; i < R.size(); ++i) {
R[i](childrenVerteces, childVerteces);
recursiveCells[l].emplace_back(
@@ -966,7 +966,7 @@ double DGElasticity::MaxStress(const Vector &U) const {
cellS += w * Frobenius(P, P);
}
- maxS = max(maxS, cellS);
+ maxS = std::max(maxS, cellS);
}
return maxS;
@@ -1103,8 +1103,10 @@ void DGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
for (int face = 0; face < c.Faces(); ++face) {
DGVectorFieldFaceElement FE(U, *c, face);
int bc = u_c.bc(face);
- bnd(c(), 0) = max(bnd(c(), 0), (double) bc);
- bnd(c(), 1) = max(bnd(c(), 1), bc == 199 ? 1.0 : 0.0); // Fixation Boundary bc == 199
+ bnd(c(), 0) = std::max(bnd(c(), 0), (double)bc);
+ bnd(c(), 1) =
+ std::max(bnd(c(), 1),
+ bc == 199 ? 1.0 : 0.0); // Fixation Boundary bc == 199
}
}
plot.AddData("Boundary", bnd, 0);
@@ -1125,12 +1127,10 @@ Vector makePlottableDGVector(const Vector &U) {
return plotU;
}
-
void DGElasticity::PlotDisplacement(const Vector &U, int step,
- const string &varname) const {
+ const std::string &varname) const {
// Plotting the deformation in src/CardiacAssemble.hpp
PlotDeformed(U, {{"Displacement", makePlottableDGVector(U)}}, step, varname);
-
}
void DGElasticity::SetDisplacement(Vector &U) const {
@@ -1269,8 +1269,8 @@ std::pair<double, double> DGElasticity::detF(const Vector &u) const {
for (int q = 0; q < E.nQ(); ++q) {
Tensor DU = E.VectorGradient(q, u);
double J = det(One + DU);
- Jmin = min(Jmin, J);
- Jmax = max(Jmax, J);
+ Jmin = std::min(Jmin, J);
+ Jmax = std::max(Jmax, J);
}
}
Jmin = PPM->Min(Jmin);
diff --git a/src/elasticity/assemble/DGElasticity.hpp b/src/elasticity/assemble/DGElasticity.hpp
index 010a62e350c6efea650988639a5979707aef3445..9d501c10ff83d97964bf889f876ed6745cdaf463 100644
--- a/src/elasticity/assemble/DGElasticity.hpp
+++ b/src/elasticity/assemble/DGElasticity.hpp
@@ -90,7 +90,8 @@ public:
void PlotBoundary(const Vector &U, const Vector &scal) const override;
- void PlotDisplacement(const Vector &u, int step, const string &varname) const override;
+ void PlotDisplacement(const Vector &u, int step,
+ const std::string &varname) const override;
//std::array<double, 4> ChamberVolume(const Vector &U) const override;
void FixationBoundary(const cell &c, int face, int bc, const Vector &U, Vector &r) const;
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index 7b835427f0d6bd76535268443522fd2de1f12ac1..c4ec4b1086d5e82265d7ab564916186ea4c3fae1 100644
--- a/src/elasticity/assemble/EGElasticity.cpp
+++ b/src/elasticity/assemble/EGElasticity.cpp
@@ -36,11 +36,9 @@ void EGElasticity::Energy(const cell &c, const Vector &U, double &energy) const
// energy += w * rho * (v * v);
// energy += w * cellMat.TotalEnergy(F, Q);
-
- //energy = min(energy,det(F));
+ // energy = std::min(energy,det(F));
- energy = max(energy,cellMat.q(E, Q[0]));
-
+ energy = std::max(energy, cellMat.q(E, Q[0]));
}
return;
@@ -128,7 +126,8 @@ void EGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
w * (cellMat.TotalDerivative(Fiso, F, Q, Dtheta) - f * theta);
}
for (int f = 0; f < c.Faces(); ++f) {
- Scalar s = cellMat.Isochoric().GetMaxParameter()* (pow(degree, 2) * penalty) /U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
MixedEGVectorFieldFaceElement FE(U, *c, f);
if (E.Bnd(f) >= PRESSURE_BC_LV && E.Bnd(f) <= PRESSURE_BC_RA) {
for (int q = 0; q < FE.nQ(); ++q) {
@@ -373,7 +372,8 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
w * (cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta, Dtheta));
}
for (int face = 0; face < c.Faces(); ++face) {
- Scalar s = cellMat.Isochoric().GetMaxParameter()* (pow(degree, 2) * penalty) /U.GetMesh().MaxMeshWidth();
+ double s = cellMat.Isochoric().GetMaxParameter() *
+ (pow(degree, 2) * penalty) / U.GetMesh().MaxMeshWidth();
MixedEGVectorFieldFaceElement FE(U, *c, face);
if (E.Bnd(face) == DIRICHLET_BC || E.Bnd(face) == DIRICHLET_FIX_BC) {
for (int q = 0; q < FE.nQ(); ++q) {
@@ -852,7 +852,7 @@ double EGElasticity::MaxStress(const Vector &U) const {
cellS += w * Frobenius(P, P);
}
- maxS = max(maxS, cellS);
+ maxS = std::max(maxS, cellS);
}
return maxS;
@@ -980,7 +980,7 @@ void EGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
for (int i = 0; i < c.Faces(); ++i) {
for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
int id = U.IdOfNodalPointOnFace(*c, i, j);
- u_c(id, 0) = max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
+ u_c(id, 0) = std::max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
}
}
}
@@ -1044,12 +1044,12 @@ void EGElasticity::Project(const Vector &fine, Vector &u) const {
c.Subdomain(), c.AsVector(), false));
for (int l = 1; l <= dlevel; ++l) {
for (const auto &parent: recursiveCells[l - 1]) {
- const vector<Rule> &R = parent->Refine(childrenVerteces);
- for (int i = 0; i < R.size(); ++i) {
- R[i](childrenVerteces, childVerteces);
- recursiveCells[l].emplace_back(
- CreateCell(R[i].type(), c.Subdomain(),
- childVerteces, false));
+ const std::vector<Rule> &R = parent->Refine(childrenVerteces);
+ for (int i = 0; i < R.size(); ++i) {
+ R[i](childrenVerteces, childVerteces);
+ recursiveCells[l].emplace_back(
+ CreateCell(R[i].type(), c.Subdomain(), childVerteces,
+ false));
}
}
}
@@ -1234,8 +1234,8 @@ std::pair<double, double> EGElasticity::detF(const Vector &u) const {
for (int q = 0; q < E.nQ(); ++q) {
Tensor DU = E.VectorGradient(q, u);
double J = det(One + DU);
- Jmin = min(Jmin, J);
- Jmax = max(Jmax, J);
+ Jmin = std::min(Jmin, J);
+ Jmax = std::max(Jmax, J);
}
}
Jmin = PPM->Min(Jmin);
diff --git a/src/elasticity/assemble/IElasticity.cpp b/src/elasticity/assemble/IElasticity.cpp
index fd163f5b1cb256f9799094269b53179a96f7d992..cf29f0878031dab1fcfac5cf4546e8f1a21e026d 100644
--- a/src/elasticity/assemble/IElasticity.cpp
+++ b/src/elasticity/assemble/IElasticity.cpp
@@ -66,10 +66,9 @@ void IElasticity::FinishTimeStep(const Vector &u) {
}
void IElasticity::PlotDisplacement(const Vector &U, int step,
- const string &varname) const {
+ const std::string &varname) const {
// Plotting the deformation in src/CardiacAssemble.hpp
PlotDeformed(U, {{"Displacement", U}}, step, varname);
-
}
void IElasticity::Initialize(Vectors &vecs) const {
@@ -97,10 +96,13 @@ void IElasticity::CalculateAndPrintMinMaxGammaF() const{
vout(0) << "maximal value of gamma_f is " << max_gamma << endl ;
}
-void IElasticity::PointEvaluationGamma(const std::string &varname, const vector<Point> &point) {
+void IElasticity::PointEvaluationGamma(const std::string &varname,
+ const std::vector<Point> &point) {
IElasticity::PrintPointEvaluation(*gamma, varname, point);
};
-void IElasticity::PrintPointEvaluation(const Vector &var,const std::string &varname, const vector<Point> &point) {
+void IElasticity::PrintPointEvaluation(const Vector &var,
+ const std::string &varname,
+ const std::vector<Point> &point) {
std::vector<VectorField> values(point.size());
for (row r = var.rows(); r != var.rows_end(); ++r) {
auto [i, isNear] = isIn(r(), point);
diff --git a/src/elasticity/assemble/IElasticity.hpp b/src/elasticity/assemble/IElasticity.hpp
index 413875218fb9ccd697b01850a76000c2611a9fb5..c8fb19966fc7850c7afe8a112d1fa9c7d8a016fa 100644
--- a/src/elasticity/assemble/IElasticity.hpp
+++ b/src/elasticity/assemble/IElasticity.hpp
@@ -144,7 +144,8 @@ public:
virtual void GetInvariant(const Vector &u, Vector &iota4, Vector &iota) const { return GetInvariant(u,iota4); }
- virtual void PlotDisplacement(const Vector &u, int step, const string &varname) const;
+ virtual void PlotDisplacement(const Vector &u, int step,
+ const std::string &varname) const;
/**
* Creates a new vector with the discretization of target (coarse vector) and projects fine
@@ -283,9 +284,11 @@ public:
void SetInitialValue(Vectors &values);
- void PointEvaluationGamma(const std::string &varname, const vector<Point> &point);
+ void PointEvaluationGamma(const std::string &varname,
+ const std::vector<Point> &point);
- void PrintPointEvaluation(const Vector &var,const std::string &varname, const vector<Point> &point);
+ void PrintPointEvaluation(const Vector &var, const std::string &varname,
+ const std::vector<Point> &point);
/// Every time step will be plotted if we are in the dynamic case.
void PlotIteration(const Vector &u) const override;
diff --git a/src/elasticity/assemble/LagrangeElasticity.cpp b/src/elasticity/assemble/LagrangeElasticity.cpp
index 9c262f1e462d5fa821972d231c981bfe55a0e252..ecd738354e1d1b1c7529180b96b9c4683c9b39af 100644
--- a/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/src/elasticity/assemble/LagrangeElasticity.cpp
@@ -918,7 +918,7 @@ double LagrangeElasticity::MaxStress(const Vector &U) const {
cellS += w * Frobenius(P, P);
}
- maxS = max(maxS, cellS);
+ maxS = std::max(maxS, cellS);
}
return maxS;
@@ -1124,7 +1124,7 @@ void LagrangeElasticity::PlotBoundary(const Vector &U, const Vector &scal) const
for (int i = 0; i < c.Faces(); ++i) {
for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
int id = U.IdOfNodalPointOnFace(*c, i, j);
- u_c(id, 0) = max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
+ u_c(id, 0) = std::max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
}
}
}
@@ -1338,8 +1338,8 @@ void LagrangeElasticity::EvaluateTractionStiffness(Vector &u) const {
double Un = faceElem.QNormal(q) * faceElem.VectorValue(q,u);
double w = faceElem.QWeight(q);
- s_max = max(s_max,abs(Un));
- s += w * Un * Un;
+ s_max = std::max(s_max, abs(Un));
+ s += w * Un * Un;
}
}
}
@@ -1395,8 +1395,8 @@ std::pair<double, double> LagrangeElasticity::detF(const Vector &u) const {
for (int q = 0; q < E.nQ(); ++q) {
Tensor DU = E.VectorGradient(q, u);
double J = det(One + DU);
- Jmin = min(Jmin, J);
- Jmax = max(Jmax, J);
+ Jmin = std::min(Jmin, J);
+ Jmax = std::max(Jmax, J);
}
}
Jmin = PPM->Min(Jmin);
diff --git a/src/elasticity/assemble/MixedElasticity.cpp b/src/elasticity/assemble/MixedElasticity.cpp
index 04aae5f1b7bb43b9deb48a54d2609c57f2d29723..77f2137fc96648b9e70e6bd8622c51d8d7e5fcc3 100644
--- a/src/elasticity/assemble/MixedElasticity.cpp
+++ b/src/elasticity/assemble/MixedElasticity.cpp
@@ -25,7 +25,7 @@ void MixedElasticity::Energy(const cell &c, const Vector &W, double &energy) con
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
- Scalar p = E.PressureValue(q, W);
+ double p = E.PressureValue(q, W);
VectorField u = E.VelocityField(q, W);
Tensor Du = E.VelocityFieldGradient(q, W);
Tensor F = FA_inv + Q * Du * FA_inv * QI;
@@ -80,7 +80,7 @@ void MixedElasticity::Jacobi(const cell &c, const Vector &W, Matrix &A) const {
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
- Scalar p = E.PressureValue(q, W);
+ double p = E.PressureValue(q, W);
VectorField u = E.VelocityField(q, W);
Tensor Du = E.VelocityFieldGradient(q, W);
Tensor F = FA_inv + Q * Du * FA_inv * QI;
@@ -114,7 +114,7 @@ void MixedElasticity::Jacobi(const cell &c, const Vector &W, Matrix &A) const {
}
for (int j = 0; j < E.PressureSize(); ++j) {
- Scalar p_j = E.PressureValue(q, j);
+ double p_j = E.PressureValue(q, j);
// Volumetric Piola Kirchhoff Stress
A_c(E.VelocityIndex(), E.PressureIndex(), i, j, k, 0) +=
@@ -128,9 +128,9 @@ void MixedElasticity::Jacobi(const cell &c, const Vector &W, Matrix &A) const {
}
for (int i = 0; i < E.PressureSize(); ++i) {
- Scalar p_i = E.PressureValue(q, i);
+ double p_i = E.PressureValue(q, i);
for (int j = 0; j < E.PressureSize(); ++j) {
- Scalar p_j = E.PressureValue(q, j);
+ double p_j = E.PressureValue(q, j);
A_c(E.PressureIndex(), E.PressureIndex(), i, j, 0, 0) +=
w * p_i * p_j / eProblem.GetMaterial(d_part < 2).Volumetric().Penalty();
}
@@ -158,7 +158,7 @@ void MixedElasticity::Residual(const cell &c, const Vector &W, Vector &r) const
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
- Scalar p = E.PressureValue(q, W);
+ double p = E.PressureValue(q, W);
VectorField u = E.VelocityField(q, W);
Tensor Du = E.VelocityFieldGradient(q, W);
Tensor F = FA_inv + Q * Du * FA_inv * QI;
@@ -193,7 +193,7 @@ void MixedElasticity::Residual(const cell &c, const Vector &W, Vector &r) const
}
for (int j = 0; j < E.PressureSize(); ++j) {
- Scalar p_j = E.PressureValue(q, j);
+ double p_j = E.PressureValue(q, j);
// Volumetric Energy
r_c(E.PressureIndex(), j, 0) +=
w * p_j * mat::volumetric_energy(eProblem.GetMaterial(d_part < 2), det(F));
diff --git a/src/elasticity/data_processing/CardiacLogger.h b/src/elasticity/data_processing/CardiacLogger.h
index 8f3dd8d69950a6c68dea13f20a93759d70b5732a..f71febc7beef131208c80e21d917dafd34dc61a1 100644
--- a/src/elasticity/data_processing/CardiacLogger.h
+++ b/src/elasticity/data_processing/CardiacLogger.h
@@ -6,10 +6,10 @@
class CardiacLogger {
protected:
- const string &logFileName;
+ const std::string &logFileName;
+
public:
- CardiacLogger(const string &fileName) : logFileName(fileName) {
- };
+ CardiacLogger(const std::string &fileName) : logFileName(fileName) {};
/**
* Writes all internally stored logging arguments into the file @logFileName.
diff --git a/src/elasticity/data_processing/DisplacementLogger.cpp b/src/elasticity/data_processing/DisplacementLogger.cpp
index e462474adc66dee58519bd5cc844ef1eb0db75dd..a4c6186f76bde7cf895837701190acfcfab68f2b 100644
--- a/src/elasticity/data_processing/DisplacementLogger.cpp
+++ b/src/elasticity/data_processing/DisplacementLogger.cpp
@@ -1,15 +1,15 @@
#include "DisplacementLogger.hpp"
-
-DisplacementLogger::DisplacementLogger(const string &logName, const string &dataName)
+DisplacementLogger::DisplacementLogger(const std::string &logName,
+ const std::string &dataName)
: CardiacLogger(logName), dataFileName(dataName) {
std::ifstream dataFile("../cardiacmechanics/data/" + dataFileName);
if (dataFile.is_open()) {
- string line;
+ std::string line;
while (getline(dataFile, line)) {
// Expects to points written as six consecutive doubles
double pointsInLine[3];
- string value;
+ std::string value;
std::istringstream linestream(line);
int j = 0;
while (getline(linestream, value, ' ')) {
diff --git a/src/elasticity/data_processing/DisplacementLogger.h b/src/elasticity/data_processing/DisplacementLogger.h
index 7dec347262d2a7a46838225a3e69384be54525d3..d907d34223e64d70ec3cdcb9843bcc3ba2327439 100644
--- a/src/elasticity/data_processing/DisplacementLogger.h
+++ b/src/elasticity/data_processing/DisplacementLogger.h
@@ -5,7 +5,7 @@
class DisplacementLogger : CardiacLogger {
- const string &dataFileName;
+ const std::string &dataFileName;
///
std::vector <Point> *trackedPoints = new std::vector<Point>;
@@ -20,7 +20,7 @@ class DisplacementLogger : CardiacLogger {
}
public:
- DisplacementLogger(const string &logName, const string &dataName);
+ DisplacementLogger(const std::string &logName, const std::string &dataName);
/**
* Calculates the distances between the paired points under the given displacement.
diff --git a/src/elasticity/data_processing/ThicknessLogger.cpp b/src/elasticity/data_processing/ThicknessLogger.cpp
index a7077a54ac4754b4426a81043042a03ad3b57c71..8a97340d8507941c6075d5e9deb927850f069fa3 100644
--- a/src/elasticity/data_processing/ThicknessLogger.cpp
+++ b/src/elasticity/data_processing/ThicknessLogger.cpp
@@ -1,15 +1,15 @@
#include "ThicknessLogger.hpp"
-
-ThicknessLogger::ThicknessLogger(const string &logName, const string &dataName) : CardiacLogger(
- logName), dataFileName(dataName) {
+ThicknessLogger::ThicknessLogger(const std::string &logName,
+ const std::string &dataName)
+ : CardiacLogger(logName), dataFileName(dataName) {
std::ifstream dataFile("../geo/" + dataFileName);
if (dataFile.is_open()) {
- string line;
+ std::string line;
while (getline(dataFile, line)) {
// Expects to points written as six consecutive doubles
double pointsInLine[6];
- string value;
+ std::string value;
std::istringstream linestream(line);
int j = 0;
while (getline(linestream, value, ' ')) {
diff --git a/src/elasticity/data_processing/ThicknessLogger.h b/src/elasticity/data_processing/ThicknessLogger.h
index 77851c3ab8c3586c54f213f8b323122f3af54c54..a0e6fee693431813bdc2d32c8eda98adb9a9947e 100644
--- a/src/elasticity/data_processing/ThicknessLogger.h
+++ b/src/elasticity/data_processing/ThicknessLogger.h
@@ -5,7 +5,7 @@
class ThicknessLogger : CardiacLogger {
- const string &dataFileName;
+ const std::string &dataFileName;
///
unordered_map <Point, Point> *pairedPoints = new unordered_map<Point, Point>;
@@ -20,7 +20,7 @@ class ThicknessLogger : CardiacLogger {
}
public:
- ThicknessLogger(const string &logName, const string &dataName);
+ ThicknessLogger(const std::string &logName, const std::string &dataName);
/**
* Calculates the distances between the paired points under the given displacement.
@@ -43,10 +43,11 @@ public:
};
class MultiThicknessLogger : CardiacLogger {
- vector<ThicknessLogger *> loggers;
+ std::vector<ThicknessLogger *> loggers;
public:
- MultiThicknessLogger(int loggerSize, const string &logName, const string &dataName)
+ MultiThicknessLogger(int loggerSize, const std::string &logName,
+ const std::string &dataName)
: CardiacLogger(logName), loggers(loggerSize) {
for (int l = 0; l < loggerSize; l++) {
loggers[l] = new ThicknessLogger(logName + std::to_string(l), dataName + std::to_string(l));
diff --git a/src/elasticity/data_processing/VectorLogger.cpp b/src/elasticity/data_processing/VectorLogger.cpp
index 9a2c288212b2ea369ffc4140143682019be4b414..b7715ed9cc269bd237dbe96477cfa8d1ddfc48f2 100644
--- a/src/elasticity/data_processing/VectorLogger.cpp
+++ b/src/elasticity/data_processing/VectorLogger.cpp
@@ -8,7 +8,7 @@ void VectorLogger::Log(double time) {
Point P(0.0, 0.0, 0.0);
P.t() = time;
- vector<double> valuesToSend;
+ std::vector<double> valuesToSend;
for (auto V : Vs) {
for (row r = V.rows(); r != V.rows_end(); r++) {
if (dist(r(), P) < Eps) {
@@ -25,8 +25,7 @@ void VectorLogger::Log(double time) {
E.Communicate();
if (PPM->master()) {
- auto *snapshot = new unordered_map < Point, vector<double>
- * > ;
+ auto *snapshot = new unordered_map<Point, std::vector<double> *>;
for (short q = 0; q < PPM->size(); ++q) {
while (E.Receive(q).size() < E.ReceiveSize(q)) {
@@ -35,7 +34,7 @@ void VectorLogger::Log(double time) {
size_t count = 0;
E.Receive(q) >> count;
- auto valuesRecieved = new vector<double>;
+ auto valuesRecieved = new std::vector<double>;
for (int i = 0; i < count; ++i) {
double tmp = 0.0;
E.Receive(q) >> tmp;
diff --git a/src/elasticity/data_processing/VectorLogger.h b/src/elasticity/data_processing/VectorLogger.h
index 42f5b6fb01c1c251e949fcfabe92439de78c3424..27a14be7530586973c5bd40ad7e740fbe51dfd11 100644
--- a/src/elasticity/data_processing/VectorLogger.h
+++ b/src/elasticity/data_processing/VectorLogger.h
@@ -7,16 +7,14 @@
class VectorLogger : CardiacLogger {
- list<unordered_map < Point, vector<double> *> *> *
- data = new list < unordered_map < Point, vector<double>
- *> *>;
+ list<unordered_map<Point, std::vector<double> *> *> *data =
+ new list<unordered_map<Point, std::vector<double> *> *>;
- vector <Vector> Vs;
+ std::vector<Vector> Vs;
int size;
public:
- VectorLogger(const string &filename, int s, Vector &ref) : CardiacLogger(filename), size(s),
- Vs(s, ref) {
- }
+ VectorLogger(const std::string &filename, int s, Vector &ref)
+ : CardiacLogger(filename), size(s), Vs(s, ref) {}
void Set(int i, Vector &v) {
Vs[i] = v;
diff --git a/src/elasticity/materials/Material.cpp b/src/elasticity/materials/Material.cpp
index f157135623613ef9b7fed45d3249cba056868dff..8d6f0252b18b69a5be5a6a24efb670d96093dacd 100644
--- a/src/elasticity/materials/Material.cpp
+++ b/src/elasticity/materials/Material.cpp
@@ -7,8 +7,7 @@
#include "LaplaceMaterial.hpp"
#include "MooneyRivlinMaterial.hpp"
-
-Material::Material(const string &materialName) {
+Material::Material(const std::string &materialName) {
// === Create HyperElasticMaterial ==== //
if (materialName == "Laplace") hyperelMat = std::make_unique<LaplaceMaterial>();
else if (materialName == "Linear") hyperelMat = std::make_unique<LinearMaterial>();
@@ -30,7 +29,8 @@ Material::Material(const string &materialName) {
viscoMat = GetDamping();
}
-Material::Material(const string &materialName, const vector<double> &materialParameters) {
+Material::Material(const std::string &materialName,
+ const std::vector<double> &materialParameters) {
if (materialName == "Laplace")
hyperelMat = std::make_unique<LaplaceMaterial>(materialParameters);
else if (materialName == "Linear")
@@ -80,13 +80,13 @@ std::unique_ptr<ViscoelasticDamping> Material::GetDamping(const std::string &dam
}
std::unique_ptr<VolumetricPenalty> Material::GetVolumetricPenalty() {
- string penaltyType = "None";
+ std::string penaltyType = "None";
Config::Get("QuasiCompressiblePenalty", penaltyType);
return GetVolumetricPenalty(penaltyType);
}
std::unique_ptr<ViscoelasticDamping> Material::GetDamping() {
- string dampingType = "None";
+ std::string dampingType = "None";
Config::Get("ViscoelasticDamping", dampingType);
return GetDamping(dampingType);
}
diff --git a/src/elasticity/materials/Material.hpp b/src/elasticity/materials/Material.hpp
index cdb1f6a517a7b73da6841af0fab689c505981423..8f99b8f3228ebbccc770a702b044e74697f78f69 100644
--- a/src/elasticity/materials/Material.hpp
+++ b/src/elasticity/materials/Material.hpp
@@ -141,18 +141,21 @@ class Material {
std::unique_ptr<VolumetricPenalty> GetVolumetricPenalty();
- std::unique_ptr<VolumetricPenalty> GetVolumetricPenalty(const string &penaltyType);
+ std::unique_ptr<VolumetricPenalty> GetVolumetricPenalty(
+ const std::string &penaltyType);
std::unique_ptr<ViscoelasticDamping> GetDamping();
- std::unique_ptr<ViscoelasticDamping> GetDamping(const string &dampingType);
+ std::unique_ptr<ViscoelasticDamping> GetDamping(
+ const std::string &dampingType);
public:
Material() = default;
- explicit Material(const string &materialName);
+ explicit Material(const std::string &materialName);
- Material(const string &materialName, const vector<double> &materialParameters);
+ Material(const std::string &materialName,
+ const std::vector<double> &materialParameters);
double IsoEnergy(const Tensor &F) const {
return hyperelMat->IsotropicEnergy(F);
diff --git a/src/elasticity/problem/BeamProblems.cpp b/src/elasticity/problem/BeamProblems.cpp
index 26233d4c33c96e7a570630ab5e55640e21cdc57d..44cb560139acb5891638f6f91b4cf6126673b892 100644
--- a/src/elasticity/problem/BeamProblems.cpp
+++ b/src/elasticity/problem/BeamProblems.cpp
@@ -80,7 +80,8 @@ std::string BeamProblem::Evaluate(const Vector &solution) const {
return evaluation;
}
-vector<double> BeamProblem::EvaluationResults(const Vector &solution) const {
+std::vector<double> BeamProblem::EvaluationResults(
+ const Vector &solution) const {
bool test = false;
auto displacements = atEvaluationPoints(solution);
@@ -102,11 +103,10 @@ vector<double> BeamProblem::EvaluationResults(const Vector &solution) const {
return {};
}
-string BeamProblem::Name() const {
+std::string BeamProblem::Name() const {
return "Default Beam Problem";
}
-
VectorField PressurePullBeamProblem::Deformation(double time, const Point &x) const {
return VectorField(x[0], 0, 0);
}
diff --git a/src/elasticity/problem/BeamProblems.hpp b/src/elasticity/problem/BeamProblems.hpp
index 0458df06eb54fe1e2f0540526b31b9e6dcf4ddd3..87f3048b0eb92c4fb3f06a890072bbae8ef4df8e 100644
--- a/src/elasticity/problem/BeamProblems.hpp
+++ b/src/elasticity/problem/BeamProblems.hpp
@@ -10,7 +10,7 @@
class BeamProblem : public ElasticityProblem {
protected:
int spaceDim{3};
- string is_tet{"Tet"};
+ std::string is_tet{"Tet"};
std::array<VectorField, 11> atEvaluationPoints(const Vector &solution) const;
@@ -28,9 +28,10 @@ public:
}
}
- BeamProblem(const std::string &matName, const vector<double> &matParameters,
- std::string defaultMesh = "") :
- ElasticityProblem(matName, matParameters) {
+ BeamProblem(const std::string &matName,
+ const std::vector<double> &matParameters,
+ std::string defaultMesh = "")
+ : ElasticityProblem(matName, matParameters) {
Config::Get("ProblemDimension", spaceDim);
Config::Get("ProblemGeometry", is_tet);
@@ -39,11 +40,11 @@ public:
}
}
- string Evaluate(const Vector &solution) const override;
+ std::string Evaluate(const Vector &solution) const override;
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
- string Name() const override;
+ std::string Name() const override;
};
template<int degree>
@@ -83,7 +84,9 @@ public:
return true;
}
- string Name() const override { return "DirichletBeamDegree" + std::to_string(degree); }
+ std::string Name() const override {
+ return "DirichletBeamDegree" + std::to_string(degree);
+ }
};
class PressureBeamProblem : public BeamProblem {
diff --git a/src/elasticity/problem/DynamicBoundaryProblems.cpp b/src/elasticity/problem/DynamicBoundaryProblems.cpp
index e10228d3eb5de7d3d20e56c05bb98f1712a41815..768e619eb8b389ebc45e9cf2fc726dac37a6b88c 100644
--- a/src/elasticity/problem/DynamicBoundaryProblems.cpp
+++ b/src/elasticity/problem/DynamicBoundaryProblems.cpp
@@ -1,5 +1,7 @@
#include "DynamicBoundaryProblems.hpp"
+#include <numbers>
+
VectorField DExponentialStiffnessProblem::Deformation(double time, const Point &x) const {
VectorField u(
exp(x[0]), exp(x[1]), spaceDim > 2 ? exp(x[2]) : 0.0
@@ -211,46 +213,59 @@ double DPolynomialNeumannProblem::pressure(double t, const Point &x, int bndCond
}
}
-
-VectorField DExponentialNeumannProblem::Deformation(double time, const Point &x) const {
- VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]),cos(0.5 * Pi * x[2]));
- return -2.0 / Pi * sin(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Deformation(double time,
+ const Point &x) const {
+ using std::numbers::pi;
+ VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+ cos(0.5 * pi * x[2]));
+ return -2.0 / pi * sin(0.5 * pi * time) * u;
}
-VectorField DExponentialNeumannProblem::Velocity(double time, const Point &x) const {
- VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]), cos(0.5 * Pi * x[2]));
- return -cos(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Velocity(double time,
+ const Point &x) const {
+ using std::numbers::pi;
+ VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+ cos(0.5 * pi * x[2]));
+ return -cos(0.5 * pi * time) * u;
}
-VectorField DExponentialNeumannProblem::Acceleration(double time, const Point &x) const {
- VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]), cos(0.5 * Pi * x[2]) );
- return 0.5 * Pi * sin(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Acceleration(double time,
+ const Point &x) const {
+ using std::numbers::pi;
+ VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+ cos(0.5 * pi * x[2]));
+ return 0.5 * pi * sin(0.5 * pi * time) * u;
}
-
-Tensor DExponentialNeumannProblem::DeformationGradient(double time, const Point &x) const {
- Tensor Du(- 0.5 * Pi * sin(0.5 * Pi * x[0]), 0.0, 0.0,
- 0.0, - 0.5 * Pi * sin(0.5 * Pi * x[1]), 0.0,
- 0.0, 0.0, - 0.5 * Pi * sin(0.5 * Pi * x[2])
- );
- return -2.0 / Pi * sin(0.5 * Pi * time) * Du;
+Tensor DExponentialNeumannProblem::DeformationGradient(double time,
+ const Point &x) const {
+ using std::numbers::pi;
+ Tensor Du(-0.5 * pi * sin(0.5 * pi * x[0]), 0.0, 0.0, 0.0,
+ -0.5 * pi * sin(0.5 * pi * x[1]), 0.0, 0.0, 0.0,
+ -0.5 * pi * sin(0.5 * pi * x[2]));
+ return -2.0 / pi * sin(0.5 * pi * time) * Du;
}
-VectorField DExponentialNeumannProblem::Load(double time, const Point &x, const Cell &c) const {
+VectorField DExponentialNeumannProblem::Load(double time, const Point &x,
+ const Cell &c) const {
+ using std::numbers::pi;
VectorField a = Acceleration(time, x);
- VectorField u(-0.25 * Pi * Pi * cos(0.5 * Pi * x[0]), -0.25 * Pi * Pi * cos(0.5 * Pi * x[1]),
- -0.25 * Pi * Pi * cos(0.5 * Pi * x[2]));
+ VectorField u(-0.25 * pi * pi * cos(0.5 * pi * x[0]),
+ -0.25 * pi * pi * cos(0.5 * pi * x[1]),
+ -0.25 * pi * pi * cos(0.5 * pi * x[2]));
- return a - ((-2.0 / Pi * sin(0.5 * Pi * time)) * u);
+ return a - ((-2.0 / pi * sin(0.5 * pi * time)) * u);
}
-double DExponentialNeumannProblem::pressure(double t, const Point &x, int bndCondition) const {
+double DExponentialNeumannProblem::pressure(double t, const Point &x,
+ int bndCondition) const {
+ using std::numbers::pi;
if (bndCondition == 230) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[0]) ;
+ return -sin(0.5 * pi * t) * sin(0.5 * pi * x[0]);
} else if (bndCondition == 231) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[1]);
+ return -sin(0.5 * pi * t) * sin(0.5 * pi * x[1]);
} else if (bndCondition == 232) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[2]) ;
+ return -sin(0.5 * pi * t) * sin(0.5 * pi * x[2]);
} else {
return 0;
}
diff --git a/src/elasticity/problem/DynamicBoundaryProblems.hpp b/src/elasticity/problem/DynamicBoundaryProblems.hpp
index b132fdb3cff7f74ee61f1e400c980cce6f5546fd..558d6f5d6f4628943c4b8c9e124adc410f4531b5 100644
--- a/src/elasticity/problem/DynamicBoundaryProblems.hpp
+++ b/src/elasticity/problem/DynamicBoundaryProblems.hpp
@@ -6,7 +6,8 @@
class DynamicBoundaryProblem : public ElasticityProblem {
protected:
int spaceDim{3};
- string is_tet{"Tet"};
+ std::string is_tet{"Tet"};
+
public:
explicit DynamicBoundaryProblem(std::string mName) :
ElasticityProblem("Laplace", {1.0}) {
@@ -28,7 +29,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
VectorField Load(double t, const Point &x, const Cell &c) const override;
@@ -51,7 +52,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
VectorField Load(double t, const Point &x, const Cell &c) const override;
@@ -74,7 +75,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
VectorField Load(double t, const Point &x, const Cell &c) const override;
@@ -100,7 +101,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
VectorField Load(double t, const Point &x, const Cell &c) const override;
@@ -126,7 +127,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
VectorField Load(double t, const Point &x, const Cell &c) const override;
@@ -178,7 +179,8 @@ public:
class OscillationProblem : public ElasticityProblem {
protected:
int spaceDim{3};
- string is_tet{"Tet"};
+ std::string is_tet{"Tet"};
+
public:
explicit OscillationProblem(std::string mName = "BenchmarkBeam") :
ElasticityProblem() {
@@ -195,9 +197,7 @@ public:
bool IsStatic() const override{return false;}
- string Name() const override {
- return "Oscillation Problem";
- }
+ std::string Name() const override { return "Oscillation Problem"; }
};
#endif //DYNAMICBOUNDARYPROBLEMS_HPP
diff --git a/src/elasticity/problem/ElasticityBenchmarkProblem.hpp b/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
index f72c3fe8b7bba36de83455bd6481d838449d610d..ef0b66466d7bfa5d8e3e617e48350896efd08fc1 100644
--- a/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
+++ b/src/elasticity/problem/ElasticityBenchmarkProblem.hpp
@@ -7,11 +7,12 @@
class LinearBenchmarkProblem : public ElasticityProblemOld {
bool isPlanar;
public:
- LinearBenchmarkProblem(bool planar, vector<double> linearMatParameters) : ElasticityProblemOld(
- "LinearBenchmarkProblem: u(x,y) = (-y,x)", "Linear", linearMatParameters),
- isPlanar(planar) {}
+ LinearBenchmarkProblem(bool planar, std::vector<double> linearMatParameters)
+ : ElasticityProblemOld("LinearBenchmarkProblem: u(x,y) = (-y,x)",
+ "Linear", linearMatParameters),
+ isPlanar(planar) {}
- string MeshName() const override {
+ std::string MeshName() const override {
return (meshFromConfig()? meshName : (isPlanar ? "UnitBlock2DQuad" : "UnitBlock3DQuad"));
}
@@ -29,11 +30,13 @@ public:
class QuadraticBenchmarkProblem : public ElasticityProblemOld {
bool isPlanar;
public:
- QuadraticBenchmarkProblem(bool planar, vector<double> linearMatParameters) : ElasticityProblemOld(
- "QuadraticBenchmarkProblem: u(x,y) = (x*y, 0)", "Linear", linearMatParameters),
- isPlanar(planar) {}
+ QuadraticBenchmarkProblem(bool planar,
+ std::vector<double> linearMatParameters)
+ : ElasticityProblemOld("QuadraticBenchmarkProblem: u(x,y) = (x*y, 0)",
+ "Linear", linearMatParameters),
+ isPlanar(planar) {}
- string MeshName() const override {
+ std::string MeshName() const override {
return (meshFromConfig()? meshName : (isPlanar ?"UnitSquare" : "UnitCube"));
}
diff --git a/src/elasticity/problem/ElasticityProblem.cpp b/src/elasticity/problem/ElasticityProblem.cpp
index 9be4e99f3be0888b2dcfb51f1eb90fa351875533..118df923c004ca390523139820ebbf82e04fecc9 100644
--- a/src/elasticity/problem/ElasticityProblem.cpp
+++ b/src/elasticity/problem/ElasticityProblem.cpp
@@ -6,7 +6,9 @@
#include "VentricleProblems.hpp"
#include "DGVectorFieldElement.hpp"
-std::vector<double> findDisplacements(const Vector &u, const vector<Point> &points, int index) {
+std::vector<double> findDisplacements(const Vector &u,
+ const std::vector<Point> &points,
+ int index) {
std::vector<double> displacements(points.size());
if (!contains(u.GetDisc().DiscName(),"DG")) {
@@ -35,7 +37,8 @@ std::vector<double> findDisplacements(const Vector &u, const vector<Point> &poin
return displacements;
}
-std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point> &points) {
+std::vector<VectorField> findDisplacements(const Vector &u,
+ const std::vector<Point> &points) {
std::vector<VectorField> displacements(points.size());
if (!contains(u.GetDisc().DiscName(),"DG")) {
for (row r = u.rows(); r != u.rows_end(); ++r) {
@@ -80,7 +83,6 @@ std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point>
return displacements;
}
-
ElasticityProblem::ElasticityProblem() : CardMechIProblem(), IElasticityProblem(), IPressureProblem() {
Config::Get("MechProblemVerbose", verbose);
@@ -110,7 +112,8 @@ std::unique_ptr<ElasticityProblem> GetElasticityProblem() {
return GetElasticityProblem(problemName);
}
-std::unique_ptr<ElasticityProblem> GetElasticityProblem(const string &problemName) {
+std::unique_ptr<ElasticityProblem> GetElasticityProblem(
+ const std::string &problemName) {
if (problemName == "Default" || problemName == "ReadFromConfig") {
return std::make_unique<DefaultElasticityProblem>();
}
diff --git a/src/elasticity/problem/ElasticityProblem.hpp b/src/elasticity/problem/ElasticityProblem.hpp
index 14e32b3ad53d3eb047a69982b5b8d9095e79c055..9f15414f08fb6e8b044c29332e42dcf71a54437a 100644
--- a/src/elasticity/problem/ElasticityProblem.hpp
+++ b/src/elasticity/problem/ElasticityProblem.hpp
@@ -160,13 +160,13 @@ public:
ElasticityProblem(const std::string &matName,
const std::vector<double> &matParameters);
- [[nodiscard]] const string &GetMeshName() const { return meshesName; }
+ [[nodiscard]] const std::string &GetMeshName() const { return meshesName; }
virtual const void SetMeshName(std::string meshName) {
meshesName = meshName;
}
- [[nodiscard]] const string &GetMaterialName(bool active) const {
+ [[nodiscard]] const std::string &GetMaterialName(bool active) const {
return (active ? activeMaterial.Name() : passiveMaterial.Name());
};
@@ -278,13 +278,13 @@ public:
return refProblem.DeformationGradient(0.0, x);
}
- [[nodiscard]] const string &Model() const override {
+ [[nodiscard]] const std::string &Model() const override {
return refProblem.Model();
}
[[nodiscard]] bool IsStatic() const override { return refProblem.IsStatic(); }
- [[nodiscard]] string Name() const override {
+ [[nodiscard]] std::string Name() const override {
return "Initial Prestress Problem";
}
};
diff --git a/src/elasticity/problem/EllipsoidProblems.cpp b/src/elasticity/problem/EllipsoidProblems.cpp
index f8ecf04f3835e7d00bdd3ca863356b47f3f2924d..cafcfe0fda7b7f93d7b8263bfd54a2a21badd0c2 100644
--- a/src/elasticity/problem/EllipsoidProblems.cpp
+++ b/src/elasticity/problem/EllipsoidProblems.cpp
@@ -5,9 +5,8 @@
#include "EllipsoidProblems.hpp"
#include "EllipsoidPoints.hpp"
-
-vector<string> EllipsoidProblems::EvaluationMetrics() const {
- return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
+std::vector<std::string> EllipsoidProblems::EvaluationMetrics() const {
+ return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
}
std::vector<double> EllipsoidProblems::EvaluationResults(const Vector &solution) const {
diff --git a/src/elasticity/problem/EllipsoidProblems.hpp b/src/elasticity/problem/EllipsoidProblems.hpp
index d5ed25389bb50e4e0e257dc8fc8848fc9752f366..1fc42e879c61a0050156c6c348a3a2852f492749 100644
--- a/src/elasticity/problem/EllipsoidProblems.hpp
+++ b/src/elasticity/problem/EllipsoidProblems.hpp
@@ -25,8 +25,8 @@ class EllipsoidProblems : public ElasticityProblem {
std::vector<Point> outerEvaluationPoints{};
protected:
int spaceDim{3};
- string is_tet{"Tet"};
- //string meshName{};
+ std::string is_tet{"Tet"};
+ // std::string meshName{};
std::array<double, 4> chamberPressure{};
@@ -46,14 +46,19 @@ public:
FillEvaluationPoints();
}
- EllipsoidProblems(const string &matName, const vector<double> &matPar, EllipsoidType eType = FULL, std::string const &defaultMesh="") :
- ElasticityProblem( matName, matPar), type(eType) , chamberPressure(pressureFromConfig()) {
- Config::Get("ProblemDimension", spaceDim);
- Config::Get("ProblemGeometry", is_tet);
-
- if(!meshFromConfig()){
- meshesName = defaultMesh.empty() ? meshType(eType) : defaultMesh;
- }
+ EllipsoidProblems(const std::string &matName,
+ const std::vector<double> &matPar,
+ EllipsoidType eType = FULL,
+ std::string const &defaultMesh = "")
+ : ElasticityProblem(matName, matPar),
+ type(eType),
+ chamberPressure(pressureFromConfig()) {
+ Config::Get("ProblemDimension", spaceDim);
+ Config::Get("ProblemGeometry", is_tet);
+
+ if (!meshFromConfig()) {
+ meshesName = defaultMesh.empty() ? meshType(eType) : defaultMesh;
+ }
FillEvaluationPoints();
}
@@ -61,15 +66,15 @@ public:
std::vector<double> EvaluationResults(const Vector &solution) const override;
- vector<string> EvaluationMetrics() const override;
+ std::vector<std::string> EvaluationMetrics() const override;
- std::pair<double, double> printMyocard(const Vector &u) const;
+ std::pair<double, double> printMyocard(const Vector &u) const;
- void printLineDisplacements(const Vector &u) const;
+ void printLineDisplacements(const Vector &u) const;
- double ActiveStretch(double t, const Point &x) const override;
+ double ActiveStretch(double t, const Point &x) const override;
- string Name() const override { return "Ellipsoid Problem"; }
+ std::string Name() const override { return "Ellipsoid Problem"; }
};
class LandEllipsoidProblem : public EllipsoidProblems {
@@ -77,9 +82,9 @@ public:
LandEllipsoidProblem() : EllipsoidProblems("Bonet", {2, 8, 2, 4}, ORIENTED) {
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
- string Name() const override { return "LandEllipsoid Problem"; }
+ std::string Name() const override { return "LandEllipsoid Problem"; }
};
class LandBiventricleProblem : public EllipsoidProblems {
@@ -89,7 +94,7 @@ public:
std::vector<double> EvaluationResults(const Vector &solution) const override;
- string Name() const override { return "LandBiventricle Problem"; }
+ std::string Name() const override { return "LandBiventricle Problem"; }
};
#endif //CARDMECH_ELLIPSOIDPROBLEMS_HPP
diff --git a/src/elasticity/problem/TestProblems.hpp b/src/elasticity/problem/TestProblems.hpp
index 190dba6f17cf9b55535d96f577ecda5cd3ec6ec1..a8e2d64e5141cc4c9634ef76a538a64119ee73b0 100644
--- a/src/elasticity/problem/TestProblems.hpp
+++ b/src/elasticity/problem/TestProblems.hpp
@@ -42,7 +42,9 @@ public:
return true;
}
- string Name() const override { return "LaplaceProblemDegree" + std::to_string(degree); }
+ std::string Name() const override {
+ return "LaplaceProblemDegree" + std::to_string(degree);
+ }
};
diff --git a/src/elasticity/problem/VentricleProblems.cpp b/src/elasticity/problem/VentricleProblems.cpp
index f3af98e6b6c391aa68ec0d48ce858ea13e8d1f53..5436b0886ecdb223c9ba1871a6d80eb17885c853 100644
--- a/src/elasticity/problem/VentricleProblems.cpp
+++ b/src/elasticity/problem/VentricleProblems.cpp
@@ -23,23 +23,27 @@ std::pair<double, double> VentricleProblem::printMyocard(const Vector &u) const
PPM->Sum(meanThickness) / ventricleEndocard.size()};
}
-vector<std::string> VentricleProblem::EvaluationMetrics() const{
- return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
+std::vector<std::string> VentricleProblem::EvaluationMetrics() const {
+ return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
}
-vector<double> VentricleProblem::EvaluationResults(const Vector &solution) const {
- auto myocardValues = printMyocard(solution);
-
- auto apexDisplacement = findDisplacements(solution, {ventricleApex});
- std::vector<double> apexDisplacement0 = {apexDisplacement[0][0], apexDisplacement[0][1],
- apexDisplacement[0][2]};
- std::vector<double> apexDisplacement0AndApex = {ventricleApex[0] + apexDisplacement0[0],
- ventricleApex[1] + apexDisplacement0[1],
- ventricleApex[2] + apexDisplacement0[2]};
- mout.PrintInfo("LeftVentricle", verbose,
- /*PrintInfoEntry("MyocardStrain", meanStrain, 1),*/
- PrintInfoEntry("ApexLocation", apexDisplacement0AndApex, 1),
- PrintInfoEntry("ApexDisplacement", apexDisplacement0, 1));
-
- return {myocardValues.first, myocardValues.second, apexDisplacement[0][0], apexDisplacement[0][1], apexDisplacement[0][2]};
+std::vector<double> VentricleProblem::EvaluationResults(
+ const Vector &solution) const {
+ auto myocardValues = printMyocard(solution);
+
+ auto apexDisplacement = findDisplacements(solution, {ventricleApex});
+ std::vector<double> apexDisplacement0 = {apexDisplacement[0][0],
+ apexDisplacement[0][1],
+ apexDisplacement[0][2]};
+ std::vector<double> apexDisplacement0AndApex =
+ {ventricleApex[0] + apexDisplacement0[0],
+ ventricleApex[1] + apexDisplacement0[1],
+ ventricleApex[2] + apexDisplacement0[2]};
+ mout.PrintInfo("LeftVentricle", verbose,
+ /*PrintInfoEntry("MyocardStrain", meanStrain, 1),*/
+ PrintInfoEntry("ApexLocation", apexDisplacement0AndApex, 1),
+ PrintInfoEntry("ApexDisplacement", apexDisplacement0, 1));
+
+ return {myocardValues.first, myocardValues.second, apexDisplacement[0][0],
+ apexDisplacement[0][1], apexDisplacement[0][2]};
}
diff --git a/src/elasticity/problem/VentricleProblems.hpp b/src/elasticity/problem/VentricleProblems.hpp
index ce444cad406d3189b923dba9743fa8b596e18a5e..08e2986bba7df85c1ddc8a5cfcf79849554139b4 100644
--- a/src/elasticity/problem/VentricleProblems.hpp
+++ b/src/elasticity/problem/VentricleProblems.hpp
@@ -18,12 +18,13 @@ public:
}
}
- VentricleProblem(const string &matName, const vector<double> &matPar,
+ VentricleProblem(const std::string &matName,
+ const std::vector<double> &matPar,
std::string const &defaultMesh = "")
- : VentricleProblem(matName, matPar,pressureFromConfig(0), defaultMesh) {
- }
+ : VentricleProblem(matName, matPar, pressureFromConfig(0), defaultMesh) {}
- VentricleProblem(const string &matName, const vector<double> &matPar, int pressure,
+ VentricleProblem(const std::string &matName,
+ const std::vector<double> &matPar, int pressure,
std::string const &defaultMesh = "")
: ElasticityProblem(matName, matPar), ventriclePressure(pressure) {
if (!meshFromConfig()) {
@@ -31,13 +32,13 @@ public:
}
}
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
- vector<string> EvaluationMetrics() const override;
+ std::vector<std::string> EvaluationMetrics() const override;
std::pair<double, double> printMyocard(const Vector &u) const;
- string Name() const override { return "Ventricle Problem"; }
+ std::string Name() const override { return "Ventricle Problem"; }
};
diff --git a/src/elasticity/problem/old/ElasticityProblemOld.cpp b/src/elasticity/problem/old/ElasticityProblemOld.cpp
index 10d2053ca8ffcb95f81b134fd230425f2795b2e2..7d759f69dab50913ed621e9e521ad36761c607bd 100644
--- a/src/elasticity/problem/old/ElasticityProblemOld.cpp
+++ b/src/elasticity/problem/old/ElasticityProblemOld.cpp
@@ -4,8 +4,9 @@
#include "VentricleProblems_old.hpp"
#include "../DynamicProblems.hpp"
-
-std::vector<double> findDisplacements(const Vector &u, const vector<Point> &points, int index) {
+std::vector<double> findDisplacements(const Vector &u,
+ const std::vector<Point> &points,
+ int index) {
std::vector<double> displacements(points.size());
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto[i, isNear] = isIn(r(), points);
@@ -24,7 +25,8 @@ std::vector<double> findDisplacements(const Vector &u, const vector<Point> &poin
return displacements;
}
-std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point> &points) {
+std::vector<VectorField> findDisplacements(const Vector &u,
+ const std::vector<Point> &points) {
std::vector<VectorField> displacements(points.size());
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto[i, isNear] = isIn(r(), points);
@@ -45,9 +47,9 @@ std::vector<VectorField> findDisplacements(const Vector &u, const vector<Point>
return displacements;
}
-
std::unique_ptr<ElasticityProblemOld> getElasticityProblemOld() {
- return getElasticityProblemOld(problemIndex[Config::GetWithoutCheck<string>("MechProblem")]);
+ return getElasticityProblemOld(
+ problemIndex[Config::GetWithoutCheck<std::string>("MechProblem")]);
}
std::unique_ptr<ElasticityProblemOld> getElasticityProblemOld(const ProblemType cProblem) {
@@ -56,9 +58,12 @@ std::unique_ptr<ElasticityProblemOld> getElasticityProblemOld(const ProblemType
switch (cProblem) {
case ElasticityBenchmarkLinear:
- return std::make_unique<LinearBenchmarkProblem>(false, vector<double>{1, 1});
+ return std::make_unique<LinearBenchmarkProblem>(false,
+ std::vector<double>{1,
+ 1});
case ElasticityBenchmarkQuadratic:
- return std::make_unique<QuadraticBenchmarkProblem>(false, vector<double>{1, 1});
+ return std::make_unique<
+ QuadraticBenchmarkProblem>(false, std::vector<double>{1, 1});
case FullEllipsoid:
return std::make_unique<EllipsoidProblem>(FULL);
case QuarterEllipsoid:
diff --git a/src/elasticity/problem/old/ElasticityProblemOld.hpp b/src/elasticity/problem/old/ElasticityProblemOld.hpp
index bd403e5a1a0a403dd1ce4f856877f019b32341f8..6f342fec4ecac7f570a8b2f9e97c27cfd10f433c 100644
--- a/src/elasticity/problem/old/ElasticityProblemOld.hpp
+++ b/src/elasticity/problem/old/ElasticityProblemOld.hpp
@@ -45,8 +45,9 @@ protected:
Material activeMaterial;
Material passiveMaterial;
- string aMatName = "Linear";
- string pMatName = "Linear";
+ std::string aMatName = "Linear";
+ std::string pMatName = "Linear";
+
public:
explicit ElasticityProblemOld(std::string pName) : CardiacProblem(std::move(pName), "Mech") {
Config::Get("MechDiscretization", elasticityModel);
@@ -67,17 +68,21 @@ public:
Config::Get("VolumetricSplit", useFlorySplit);
}
- ElasticityProblemOld(std::string pName, const std::string &aMat, const vector<double> &aPar,
- const std::string &pMat, const vector<double> &pPar) :
- CardiacProblem(std::move(pName), "Mech"), aMatName(aMat), pMatName(pMat),
- activeMaterial(aMat, aPar), passiveMaterial(pMat, pPar) {
+ ElasticityProblemOld(std::string pName, const std::string &aMat,
+ const std::vector<double> &aPar, const std::string &pMat,
+ const std::vector<double> &pPar)
+ : CardiacProblem(std::move(pName), "Mech"),
+ aMatName(aMat),
+ pMatName(pMat),
+ activeMaterial(aMat, aPar),
+ passiveMaterial(pMat, pPar) {
Config::Get("MechDiscretization", elasticityModel);
Config::Get("VolumetricSplit", useFlorySplit);
}
- ElasticityProblemOld(const string &pName, const string &matName, const vector<double> &matPar)
- : ElasticityProblemOld(pName, matName, matPar, matName, matPar) {
- }
+ ElasticityProblemOld(const std::string &pName, const std::string &matName,
+ const std::vector<double> &matPar)
+ : ElasticityProblemOld(pName, matName, matPar, matName, matPar) {}
/// Right hand side in domain
virtual VectorField Load(double t, const Point &x, const Cell &c) const { return zero; };
@@ -106,7 +111,7 @@ public:
virtual Tensor DeformationGradient(double time, const Point &x) const { return Zero; };
- const string &GetMaterialName(bool active) const {
+ const std::string &GetMaterialName(bool active) const {
return (active ? aMatName : pMatName);
};
diff --git a/src/elasticity/problem/old/EllipsoidProblems_Old.cpp b/src/elasticity/problem/old/EllipsoidProblems_Old.cpp
index 19e17333e35d4cf9a88640f0f4f47efb134c8cab..e0fd335877125aceb771ccb218278aa243c37569 100644
--- a/src/elasticity/problem/old/EllipsoidProblems_Old.cpp
+++ b/src/elasticity/problem/old/EllipsoidProblems_Old.cpp
@@ -1,12 +1,12 @@
#include "EllipsoidProblems_Old.hpp"
#include "EllipsoidPoints.hpp"
-
-vector<string> EllipsoidProblem::EvaluationMetrics() const {
+std::vector<std::string> EllipsoidProblem::EvaluationMetrics() const {
return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
}
-vector<double> EllipsoidProblem::EvaluationResults(const Vector &solution) const {
+std::vector<double> EllipsoidProblem::EvaluationResults(
+ const Vector &solution) const {
printLineDisplacements(solution);
auto myoCardvalues = printMyocard(solution);
diff --git a/src/elasticity/problem/old/EllipsoidProblems_Old.hpp b/src/elasticity/problem/old/EllipsoidProblems_Old.hpp
index 03d82ec78401f8e43631bd1cab40eebd67e6ee0e..5a314af3c1cc27880a4f8ec5f53a8614897d2108 100644
--- a/src/elasticity/problem/old/EllipsoidProblems_Old.hpp
+++ b/src/elasticity/problem/old/EllipsoidProblems_Old.hpp
@@ -22,22 +22,26 @@ class EllipsoidProblem : public ElasticityProblemOld {
EllipsoidType type;
protected:
int spaceDim{3};
- string meshName{};
+ std::string meshName{};
+
public:
explicit EllipsoidProblem(EllipsoidType eType = FULL) :
ElasticityProblemOld("Ellipsoid Problem"), type(eType){
}
- EllipsoidProblem(const string &matName, const vector<double> &matPar, EllipsoidType eType = FULL) :
- ElasticityProblemOld("Ellipsoid Problem", matName, matPar), type(eType) {
+ EllipsoidProblem(const std::string &matName,
+ const std::vector<double> &matPar,
+ EllipsoidType eType = FULL)
+ : ElasticityProblemOld("Ellipsoid Problem", matName, matPar),
+ type(eType) {}
+ std::string MeshName() const override {
+ return meshFromConfig() ? meshName : meshType(type);
}
- string MeshName() const override { return meshFromConfig()? meshName : meshType(type); }
-
- vector<string> EvaluationMetrics() const override;
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<std::string> EvaluationMetrics() const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
std::pair<double, double> printMyocard(const Vector &u) const;
@@ -52,7 +56,7 @@ public:
}
- string Evaluate(const Vector &solution) const override { return ""; };
+ std::string Evaluate(const Vector &solution) const override { return ""; };
};
diff --git a/src/elasticity/problem/old/VentricleProblems_old.cpp b/src/elasticity/problem/old/VentricleProblems_old.cpp
index a39b3cfbbc864343ac19280eb82ecb28b433aca1..a7b60b666c24df0c547a697598deee36de2d7051 100644
--- a/src/elasticity/problem/old/VentricleProblems_old.cpp
+++ b/src/elasticity/problem/old/VentricleProblems_old.cpp
@@ -24,10 +24,11 @@ std::pair<double, double> VentricleProblem::printMyocard(const Vector &u) const
PPM->Sum(meanThickness) / ventricleEndocard.size()};
}
-vector<std::string> VentricleProblem::EvaluationMetrics() const{
+std::vector<std::string> VentricleProblem::EvaluationMetrics() const {
return {"MyocardStrain", "MyocardThickness", "ApexX", "ApexY", "ApexZ"};
}
-vector<double> VentricleProblem::EvaluationResults(const Vector &solution) const {
+std::vector<double> VentricleProblem::EvaluationResults(
+ const Vector &solution) const {
auto myocardValues = printMyocard(solution);
auto apexDisplacement = findDisplacements(solution, {ventricleApex});
diff --git a/src/elasticity/problem/old/VentricleProblems_old.hpp b/src/elasticity/problem/old/VentricleProblems_old.hpp
index c4634035f187909449b90884151381b39d6aae67..712761fbc95d7833926401343974524a59505a04 100644
--- a/src/elasticity/problem/old/VentricleProblems_old.hpp
+++ b/src/elasticity/problem/old/VentricleProblems_old.hpp
@@ -13,20 +13,20 @@ public:
}
- VentricleProblem(const string &matName, const vector<double> &matPar) :
- ElasticityProblemOld("Ventricle Problem", matName, matPar) {
+ VentricleProblem(const std::string &matName,
+ const std::vector<double> &matPar)
+ : ElasticityProblemOld("Ventricle Problem", matName, matPar) {}
+ std::string MeshName() const override {
+ return meshFromConfig() ? meshName : "LeftVentricle";
}
- string MeshName() const override { return meshFromConfig()? meshName : "LeftVentricle"; }
-
-
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
VectorField Traction(double t, const Point &x, const VectorField &u, const VectorField &v,
const VectorField &n, int sd) const override;
- vector<string> EvaluationMetrics() const override;
+ std::vector<std::string> EvaluationMetrics() const override;
std::pair<double, double> printMyocard(const Vector &u) const;
};
@@ -51,8 +51,9 @@ public:
}
- ContractingVentricleProblem(const string &matName, const vector<double> &matPar) :
- VentricleProblem(matName, matPar) {
+ ContractingVentricleProblem(const std::string &matName,
+ const std::vector<double> &matPar)
+ : VentricleProblem(matName, matPar) {
Config::Get("TractionK", k);
Config::Get("TractionC", c);
}
diff --git a/src/elasticity/solvers/DynamicMechanicsSolver.cpp b/src/elasticity/solvers/DynamicMechanicsSolver.cpp
index 9cf43f3000bf342bbd73bfc1a2dab0ee0356b7ff..dd29934eb540ef7a067ed63ac657a22bd8a02edc 100644
--- a/src/elasticity/solvers/DynamicMechanicsSolver.cpp
+++ b/src/elasticity/solvers/DynamicMechanicsSolver.cpp
@@ -12,8 +12,8 @@ std::unique_ptr<NonLinearSolverT<IElasticity>> GetDynamicSolver() {
return GetDynamicSolver(dynamics);
}
-std::unique_ptr<NonLinearSolverT<IElasticity>>
-GetDynamicSolver(const string &solverName) {
+std::unique_ptr<NonLinearSolverT<IElasticity>> GetDynamicSolver(
+ const std::string &solverName) {
if (contains(solverName, "Euler"))
return std::make_unique<DynamicEulerSolver>(contains(solverName, "Visco"),
contains(solverName, "Rayleigh"));
diff --git a/src/elasticity/solvers/DynamicMechanicsSolver.hpp b/src/elasticity/solvers/DynamicMechanicsSolver.hpp
index 6347e5e5a5469eeccc50d1ec483df259bd2d104b..4344f73323b1d203f5d8bef86baffc9eb32f5580 100644
--- a/src/elasticity/solvers/DynamicMechanicsSolver.hpp
+++ b/src/elasticity/solvers/DynamicMechanicsSolver.hpp
@@ -8,7 +8,8 @@
std::unique_ptr<NonLinearSolverT<IElasticity>> GetDynamicSolver();
-std::unique_ptr<NonLinearSolverT<IElasticity>> GetDynamicSolver(const string &solverName);
+std::unique_ptr<NonLinearSolverT<IElasticity>> GetDynamicSolver(
+ const std::string &solverName);
class ElastodynamicTimeIntegrator {
int verbose{1};
@@ -22,10 +23,8 @@ public:
explicit ElastodynamicTimeIntegrator() : ElastodynamicTimeIntegrator(
GetDynamicSolver()) {}
- explicit ElastodynamicTimeIntegrator(
- const string &solverName) : ElastodynamicTimeIntegrator(
- GetDynamicSolver(
- solverName)) {}
+ explicit ElastodynamicTimeIntegrator(const std::string &solverName)
+ : ElastodynamicTimeIntegrator(GetDynamicSolver(solverName)) {}
explicit ElastodynamicTimeIntegrator(
std::unique_ptr<NonLinearSolverT<IElasticity>> &&nlSolver) : nonlinearSolver(
diff --git a/src/elasticity/solvers/ElasticityEuler.cpp b/src/elasticity/solvers/ElasticityEuler.cpp
index 5270d7770c87d6ea296b8a3f6c99783b749459e7..137d8f25dcb6db61dffc1525ad7e6ce99ed0d343 100644
--- a/src/elasticity/solvers/ElasticityEuler.cpp
+++ b/src/elasticity/solvers/ElasticityEuler.cpp
@@ -41,7 +41,7 @@ void DynamicEulerSolver::operator()(const IElasticity &assemble, Vector &u) {
updateOldValues(assemble.StepSize(), u);
}
-string DynamicEulerSolver::Name() const {
+std::string DynamicEulerSolver::Name() const {
return "Dynamic Euler";
}
diff --git a/src/elasticity/solvers/ElasticityNewmark.cpp b/src/elasticity/solvers/ElasticityNewmark.cpp
index 8a84a297cb31cf63d7a4d7985298cab84c6b0273..94d02c2ea51d2972f4bdcaccc97bab5474e92a09 100644
--- a/src/elasticity/solvers/ElasticityNewmark.cpp
+++ b/src/elasticity/solvers/ElasticityNewmark.cpp
@@ -55,7 +55,7 @@ void NewmarkSolver::operator()(const IElasticity &assemble, Vector &u) {
//PlotIteration(t);
}
-string NewmarkSolver::Name() const {
+std::string NewmarkSolver::Name() const {
return "Newmark";
}
@@ -124,8 +124,8 @@ void NewmarkSolver::updateOldValues(double dt, const Vector &u) {
void NewmarkSolver::PlotIteration(double time) const {
auto &plot_vel = mpp::plot("Velocity");
plot_vel.AddData("Velocity", oldValues.Velocity());
- plot_vel.PlotFile("Velocity" + to_string(time));
+ plot_vel.PlotFile("Velocity" + std::to_string(time));
auto &plot_acc = mpp::plot("Acceleration");
plot_acc.AddData("Acceleration", oldValues.Acceleration());
- plot_acc.PlotFile("Acceleration" + to_string(time));
+ plot_acc.PlotFile("Acceleration" + std::to_string(time));
}
\ No newline at end of file
diff --git a/src/electrophysiology/MainMonodomain.cpp b/src/electrophysiology/MainMonodomain.cpp
index 3213cc0e07b241ece5bd1032d57afd59980ff6cd..19f1884252e79daad871e149d5addd605ed21ca3 100644
--- a/src/electrophysiology/MainMonodomain.cpp
+++ b/src/electrophysiology/MainMonodomain.cpp
@@ -68,7 +68,7 @@ std::vector<double> MainMonodomain::Evaluate(const Vector &solution) {
mout << "\n" << "==== Monodomain Information ============== \n\n";
std::string elphyModelClass = "IBTElphyModel";
Config::Get("ElphyModelClass", elphyModelClass);
- string cmName = "FitzHughNagumo";
+ std::string cmName = "FitzHughNagumo";
if (elphyModelClass == "IBTElphyModel") {
Config::Get("IBTVentricleModel", cmName);
} else {
@@ -148,13 +148,10 @@ void MainMonodomain::ResetLevel(int newLevel) {
elevel = newLevel;
}
-double MainMonodomain::EvaluateQuantity(const Vector &solution, const string &quantity) const {
+double MainMonodomain::EvaluateQuantity(const Vector &solution,
+ const std::string &quantity) const {
if (quantity == "L2Error") { return elphyA->L2Error(solution); }
if (quantity == "L2Norm") { return elphyA->L2Norm(solution); }
if (quantity == "H1Norm") { return elphyA->H1Norm(solution); }
else { return infty; }
}
-
-
-
-
diff --git a/src/electrophysiology/assemble/IElphyAssemble.cpp b/src/electrophysiology/assemble/IElphyAssemble.cpp
index 8cbd040261922cfc89ae34b70d3fb7cfc1598d28..891c9bbbc807f7558e6dad17e4e75ea2de4893ba 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.cpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.cpp
@@ -38,7 +38,7 @@ void IElphyAssemble::IExt(const cell &c, const Vector &u, Vector &r) const {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
double iext = E.Value(q, *Iext);
r_c(j) += w * (StepSize() * 1000.0 * iext * v_j);
}
@@ -59,16 +59,16 @@ void IElphyAssemble::updateActivationTime(const Vector &V, double time) {
}
}
-void IElphyAssemble::PrintInftyNormOnCells(const Vector &V, const string &name) const {
+void IElphyAssemble::PrintInftyNormOnCells(const Vector &V, const std::string &name) const {
vout(0) << " ||" << name << "||_infty = " << InftyNormOnCells(V) << endl;
}
-void IElphyAssemble::PrintL2Norm(const Vector &V, const string &name) const {
+void IElphyAssemble::PrintL2Norm(const Vector &V, const std::string &name) const {
vout(0) << " ||" << name << "||_L2 = " << L2Norm(V) << endl;
}
-void IElphyAssemble::PrintH1Norm(const Vector &V, const string &name) const {
+void IElphyAssemble::PrintH1Norm(const Vector &V, const std::string &name) const {
vout(0) << " ||" << name << "||_H1 = " << H1Norm(V) << endl;
}
void IElphyAssemble::PrintVariable(const Vector &var,const std::string &varname) const {
@@ -199,7 +199,7 @@ void IElphyAssemble::PrintV(const Vector &V) const {
PrintH1Norm(V, "v");
vout(0) << endl;
}
-void IElphyAssemble::PlotExcitation(const string &name) const{
+void IElphyAssemble::PlotExcitation(const std::string &name) const {
if (vtkplot < 1) return;
auto &plot = mpp::plot(name);
@@ -242,11 +242,11 @@ void IElphyAssemble::PlotActivation() const {
plot.PlotFile("ActivationTime");
}
void IElphyAssemble::writeActivationTimeFile() const {
- string lLevel;
+ std::string lLevel;
Config::Get("ElphyLevel", lLevel);
- string tLevel="0";
+ std::string tLevel = "0";
Config::Get("ElphyTimeLevel", tLevel);
- string filename("data/AT_l"+lLevel+"j"+ tLevel+".txt");
+ std::string filename("data/AT_l" + lLevel + "j" + tLevel + ".txt");
std::fstream file_out;
ExchangeBuffer B;
for (row r= (*activationTime).rows(); r != (*activationTime).rows_end(); ++r) {
@@ -295,8 +295,8 @@ void IElphyAssemble::PrintActivationTime() const {
}
}
-
-void IElphyAssemble::PlotPotential(Vector &V, const string &varname, int step) const {
+void IElphyAssemble::PlotPotential(Vector &V, const std::string &varname,
+ int step) const {
if (vtkplot < 1 || step % vtksteps != 0) return;
auto &plot = mpp::plot(varname + stepAsString(step));
plot.AddData("V", V);
@@ -308,7 +308,7 @@ Vector &IElphyAssemble::ActivationTime() const { return *activationTime; }
const Vector &IElphyAssemble::ExcitationData() const { return excData; }
void IElphyAssemble::PrintMesh(const Vector &u) const {
- vector<Point> z(u.size());
+ std::vector<Point> z(u.size());
for (row r = u.rows(); r != u.rows_end(); ++r) {
z[r.Id()] = r();
}
@@ -329,9 +329,9 @@ void IElphyAssemble::PrintMesh(const Vector &u) const {
}
mout << "FACES:" << endl;
mout << "VDATA:" << endl;
- vector<double> excitation(u.size());
- vector<double> amplitude(u.size());
- vector<double> duration(u.size());
+ std::vector<double> excitation(u.size());
+ std::vector<double> amplitude(u.size());
+ std::vector<double> duration(u.size());
for (row r = u.rows(); r != u.rows_end(); ++r) {
amplitude[r.Id()] = Amplitude(excData, r, 0);
excitation[r.Id()] = Excitation(excData, r, 0);
@@ -342,9 +342,9 @@ void IElphyAssemble::PrintMesh(const Vector &u) const {
}
mout << "CDATA:" << endl;
- vector<double> fibre(u.size());
- vector<double> sheet(u.size());
- vector<double> normal(u.size());
+ std::vector<double> fibre(u.size());
+ std::vector<double> sheet(u.size());
+ std::vector<double> normal(u.size());
for (cell c = u.GetMesh().cells(); c != u.GetMesh().cells_end(); ++c) {
const auto &data = c.GetData();
VectorField fibre = Direction(data);
@@ -378,7 +378,7 @@ double IElphyAssemble::L2Norm(const Vector &V) const {
for (cell c = V.cells(); c != V.cells_end(); ++c) {
ScalarElement E(V, *c);
for (int q = 0; q < E.nQ(); ++q) {
- Scalar V_p = E.Value(q, V, 0);
+ double V_p = E.Value(q, V, 0);
double w_q = E.QWeight(q);
v_sum += w_q * V_p * V_p;
@@ -392,7 +392,8 @@ double IElphyAssemble::L2Error(const Vector &V) const {
for (cell c = V.cells(); c != V.cells_end(); ++c) {
ScalarElement E(V, *c);
for (int q = 0; q < E.nQ(); ++q) {
- Scalar V_p = E.Value(q, V, 0) - GetElphyProblem().Potential(Time(), E.QPoint(q));
+ double V_p =
+ E.Value(q, V, 0) - GetElphyProblem().Potential(Time(), E.QPoint(q));
double w_q = E.QWeight(q);
v_sum += w_q * V_p * V_p;
@@ -407,7 +408,7 @@ double IElphyAssemble::L2AvgNorm(const Vector &u) const {
ScalarElement E(u, *c);
double w = 0.0;
- Scalar U = 0.0;
+ double U = 0.0;
for (int q = 0; q < E.nQ(); ++q) {
w += E.QWeight(q);
U += E.Value(q, u);
@@ -423,7 +424,7 @@ double IElphyAssemble::H1Norm(const Vector &V) const {
for (cell c = V.cells(); c != V.cells_end(); ++c) {
ScalarElement E(V, *c);
for (int q = 0; q < E.nQ(); ++q) {
- Scalar V_p = E.Value(q, V, 0);
+ double V_p = E.Value(q, V, 0);
double w_q = E.QWeight(q);
Tensor DV = E.Derivative(q, V, 0);
h1 += w_q * (V_p * V_p + Frobenius(DV, DV));
@@ -477,9 +478,9 @@ void IElphyAssemble::PlotIteration(const Vectors &elphyValues, int step) const {
#include "Monodomain.hpp"
#include "MonodomainSplitting.hpp"
-std::unique_ptr<IElphyAssemble>
-CreateElphyAssemble(ElphyProblem &problem, MCellModel &cellModel, const string &modelName,
- int discDegree) {
+std::unique_ptr<IElphyAssemble> CreateElphyAssemble(
+ ElphyProblem &problem, MCellModel &cellModel, const std::string &modelName,
+ int discDegree) {
int hexaQuadExactupTo = 2 * discDegree;
if (problem.Name() == "HexaTestProblem") {
Config::Get("HexaQuadExactupTo", hexaQuadExactupTo);
diff --git a/src/electrophysiology/assemble/IElphyAssemble.hpp b/src/electrophysiology/assemble/IElphyAssemble.hpp
index 68425c47874c1211e3444705e46afc2b3cce66eb..675641cfc8394d6d9e6a784e860f5ddff5ef8deb 100644
--- a/src/electrophysiology/assemble/IElphyAssemble.hpp
+++ b/src/electrophysiology/assemble/IElphyAssemble.hpp
@@ -104,10 +104,10 @@ public:
THROW("Not implemented in the basis class")
};
- void PrintInftyNormOnCells(const Vector &V, const string &name ) const;
- void PrintL2Norm(const Vector &V, const string &name) const;
+ void PrintInftyNormOnCells(const Vector &V, const std::string &name ) const;
+ void PrintL2Norm(const Vector &V, const std::string &name) const;
- void PrintH1Norm(const Vector &V, const string &name) const;
+ void PrintH1Norm(const Vector &V, const std::string &name) const;
void updateActivationTime(const Vector &V, double time);
@@ -118,7 +118,7 @@ public:
void PrintMesh(const Vector &u) const;
void PlotExcitation() const;
- void PlotExcitation(const string &name) const;
+ void PlotExcitation(const std::string &name) const;
void PlotActivation() const;
diff --git a/src/electrophysiology/assemble/Monodomain.cpp b/src/electrophysiology/assemble/Monodomain.cpp
index aba40c5c7fec5936dfd9187f879d7c2d1e4c77d6..5e3db1f374d2cf8cb4df51055142370b0556079d 100644
--- a/src/electrophysiology/assemble/Monodomain.cpp
+++ b/src/electrophysiology/assemble/Monodomain.cpp
@@ -74,10 +74,10 @@ void Monodomain::SystemMatrix(const cell &c, const Vectors &VCW, Matrix &A) cons
vcw[i] = E.Value(q, VCW[i], E.IndexPart(j));
}
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField S_DV_j = S * E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
//Original:
//A_c(i, j, 0, 0) +=w * (V_i * V_j + StepSize() * ((1 / (beta * C_m)) * (DV_i * S_DV_j)+ M.TimeScale() * DVIion * V_i * V_j));
@@ -103,7 +103,7 @@ void Monodomain::addReactionSAndRHS(const cell &c, const Vectors &VCW, Matrix &A
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
double Iion = elphyModel.GetIion(vcw);
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
//Original
@@ -112,7 +112,7 @@ void Monodomain::addReactionSAndRHS(const cell &c, const Vectors &VCW, Matrix &A
(Iion - DVIion * vcw[0] - iext)) * V_j;
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
A_c(i, j, 0, 0) +=
w *
(StepSize() * elphyModel.TimeScale() * (1.0 / C_m) * (1 / cmSquareTommSquare) * DVIion *
@@ -134,9 +134,9 @@ void Monodomain::addReaction(const cell &c, const Vectors &VCW, Matrix &A) const
//vcw[i] = E.Value(q, VCW[i], E.IndexPart(j));
//}
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
A_c(i, j, 0, 0) +=
w *
(StepSize() * elphyModel.TimeScale() * (1.0 / C_m) * (1 / cmSquareTommSquare) * DVIion *
@@ -155,10 +155,10 @@ void Monodomain::MPlusK(const cell &c, const Vector &V, Matrix &A) const {
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField S_DV_j = S * E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
A_c(i, j, 0, 0) +=
w * (J * V_i * V_j + StepSize() * (1 / (beta * C_m * mircoFToF)) * (DV_i * S_DV_j));
@@ -181,7 +181,7 @@ void Monodomain::RHS(const cell &c, const Vectors &VCW, Vector &r) const {
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
double Iion = elphyModel.GetIion(vcw);
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
//Original
@@ -257,7 +257,7 @@ MonodomainLINodes::addReactionSAndRHS(const cell &c, const Vectors &VCW, Matrix
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
double Iion = E.Value(q, *IionVEC, E.IndexPart(j));
double DVIion = E.Value(q, *DVIionVEC, E.IndexPart(j));
//Original
@@ -266,7 +266,7 @@ MonodomainLINodes::addReactionSAndRHS(const cell &c, const Vectors &VCW, Matrix
(Iion - DVIion * v_old - iext)) * V_j;
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
A_c(i, j, 0, 0) +=
w *
(StepSize() * elphyModel.TimeScale() * (1.0 / C_m) * (1 / cmSquareTommSquare) * DVIion *
@@ -289,10 +289,10 @@ void MonodomainLINodes::SystemMatrix(const cell &c, const Vectors &VCW, Matrix &
for (int j = 0; j < E.size(); ++j) {
E.Values(q, VCW, vcw, E.IndexPart(j));
double DVIion = E.Value(q, *DVIionVEC, E.IndexPart(j));
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField S_DV_j = S * E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
A_c(i, j, 0, 0) +=
w * (V_i * V_j + StepSize() * ((1 / (beta * C_m)) * (DV_i * S_DV_j)
@@ -315,7 +315,7 @@ void MonodomainLINodes::RHS(const cell &c, const Vectors &VCW, Vector &r) const
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
double Iion = E.Value(q, *IionVEC, E.IndexPart(j));
double DVIion = E.Value(q, *DVIionVEC, E.IndexPart(j));
@@ -379,10 +379,10 @@ void MonodomainQuadrature::SystemMatrix(const cell &c, const Vector &V, Matrix &
}
updateCaAndW(c, vcw, (q * E.size()) + j);
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField S_DV_j = S * E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
A_c(i, j, 0, 0) +=
w * (V_i * V_j + StepSize() * ((1 / (beta * C_m)) * (DV_i * S_DV_j)
@@ -407,7 +407,7 @@ void MonodomainQuadrature::RHS(const cell &c, const Vectors &VCW, Vector &r) con
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
double Iion = elphyModel.GetIion(vcw);
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
r_c(j) +=
@@ -445,7 +445,7 @@ void MonodomainFullNewton::Residual(const cell &c, const Vector &V, Vector &r) c
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
VectorField Dv_j = E.Derivative(q, j);
double Iion = elphyModel.GetIion(vcw);
double v_old = E.Value(q, (*VACAW)[0], E.IndexPart(j));
@@ -470,10 +470,10 @@ void MonodomainFullNewton::Jacobi(const cell &c, const Vector &V, Matrix &jacobi
E.Values(q, (*VACAW), vcw, E.IndexPart(j));
vcw[0] = E.Value(q, V, E.IndexPart(j));
double DVIion = elphyModel.JacobiEntry(0, 0, vcw);
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField S_DV_j = S * E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
//A_c(i, j, 0, 0) +=w * (StepSize()* M.TimeScale() * (1.0/C_m)*(1/cmSquareTommSquare)*DVIion * V_i * V_j);
A_c(i, j, 0, 0) += w * (V_i * V_j +
@@ -515,7 +515,7 @@ void MonodomainSemiImplicit::RHS(const cell &c, const Vectors &VCW, Vector &r) c
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
r_c(j) += w * (J_old * vcw[0] -
J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
@@ -530,13 +530,13 @@ void MonodomainSemiImplicitNodes::RHS(const cell &c, const Vectors &VCW, Vector
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, VCW[0],E.IndexPart(j));
+ double v = E.Value(q, VCW[0], E.IndexPart(j));
double Iion = E.Value(q, *IionVEC, E.IndexPart(j));
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
r_c(j) += w * (v -
elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
@@ -558,8 +558,8 @@ void MonodomainSemiImplicitNodes::updateIionVecs(const Vectors &VCW) {
vout(10).EndBlock();
}
-/*void MonodomainSemiImplicitCells::RHS(const cell &c,const Vector &V,const Vectors &VCW_onCells, Vector &r) const {
- Point P(2.125, 2.75, 2.875);
+/*void MonodomainSemiImplicitCells::RHS(const cell &c,const Vector &V,const
+Vectors &VCW_onCells, Vector &r) const { Point P(2.125, 2.75, 2.875);
double J = problem.Jacobian(*c);
double J_old = problem.Jacobian_old(*c);
MultiPartScalarElementWG E(V, *c, DomainPart(*c));
@@ -572,15 +572,15 @@ void MonodomainSemiImplicitNodes::updateIionVecs(const Vectors &VCW) {
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, V,E.IndexPart(j));
+ double v = E.Value(q, V,E.IndexPart(j));
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
r_c(j) += w * (J_old * v -
- J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
- (Iion - iext)) * V_j;
+ J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1
+/ cmSquareTommSquare) * (Iion - iext)) * V_j;
}
}
}*/
@@ -594,12 +594,12 @@ void MonodomainSemiImplicitCells::RHS(const cell &c,const Vector &V, Vector &r)
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, V,E.IndexPart(j));
+ double v = E.Value(q, V,E.IndexPart(j));
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
r_c(j) += w * (J_old * v -
J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
@@ -633,8 +633,8 @@ void MonodomainSemiImplicitCellsIextPerCell::RHS(const cell &c,const Vector &V,
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, V,E.IndexPart(j));
- Scalar V_j = E.Value(q, j);
+ double v = E.Value(q, V, E.IndexPart(j));
+ double V_j = E.Value(q, j);
r_c(j) += w * (J_old * v -
J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
@@ -691,14 +691,14 @@ void MonodomainSemiImplicitCellsSVI::RHS(const cell &c,const Vector &V,const Vec
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, V,E.IndexPart(j));
+ double v = E.Value(q, V, E.IndexPart(j));
vcw[0]=v;
double Iion = elphyModel.GetIion(vcw);
double iext = 0.0;
if (!isNormIextZero) {
iext = E.Value(q, *Iext, E.IndexPart(j));
}
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
r_c(j) += w * (J_old * v -
J * elphyModel.TimeScale() * StepSize() * (1.0 / C_m) * (1 / cmSquareTommSquare) *
(Iion - iext)) * V_j;
diff --git a/src/electrophysiology/assemble/MonodomainSplitting.cpp b/src/electrophysiology/assemble/MonodomainSplitting.cpp
index 4b641905b1a50bb80f4ad60789bc7995b1a0b6e1..16ca0598090d516e2a2d1ba6a8b26637b048232d 100644
--- a/src/electrophysiology/assemble/MonodomainSplitting.cpp
+++ b/src/electrophysiology/assemble/MonodomainSplitting.cpp
@@ -56,10 +56,10 @@ void MonodomainSplitting::SystemMatrix(const cell &c, const Vector &V, Matrix &A
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField DV_j = E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
// Crank-Nicolson Scheme
A_c(i, j) += w * (J*V_i * V_j + J*theta * (StepSize() / (beta * C_m)) * DV_i * (S * DV_j));
@@ -77,10 +77,10 @@ void MonodomainSplitting::RHSMAT(const cell &c, const Vector &V, Matrix &R) cons
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar V_j = E.Value(q, j);
+ double V_j = E.Value(q, j);
VectorField DV_j = E.Derivative(q, j);
for (int i = 0; i < E.size(); ++i) {
- Scalar V_i = E.Value(q, i);
+ double V_i = E.Value(q, i);
VectorField DV_i = E.Derivative(q, i);
// Crank-Nicolson Scheme
R_c(i, j) +=
@@ -101,10 +101,10 @@ void MonodomainSplitting::RHS(const cell &c, const Vector &V, Vector &r) const {
Point xi = E.QPoint(q);
for (int j = 0; j < E.size(); ++j) {
- Scalar v = E.Value(q, V);
+ double v = E.Value(q, V);
VectorField S_Dv = S * E.Derivative(q, V);
- Scalar v_j = E.Value(q, j);
+ double v_j = E.Value(q, j);
VectorField Dv_j = E.Derivative(q, j);
r_c(j) += w * (J_old*v * v_j -J_old* (1.0 / (beta * C_m)) * (1 - theta) * StepSize() * S_Dv * Dv_j);
if (iextInPDE) {
diff --git a/src/electrophysiology/problem/DiffusionProblem.cpp b/src/electrophysiology/problem/DiffusionProblem.cpp
index 6a0bb43b58f2a6e35113d4f6ffc885fba54b9d55..1df6419a78338d87c2e4c3150ef9272a5a994db0 100644
--- a/src/electrophysiology/problem/DiffusionProblem.cpp
+++ b/src/electrophysiology/problem/DiffusionProblem.cpp
@@ -1,5 +1,6 @@
#include "DiffusionProblem.hpp"
+#include <numbers>
#include <utility>
DiffusionProblem::DiffusionProblem(std::array<double, 9> entries) :
@@ -15,8 +16,8 @@ DiffusionProblem::DiffusionProblem(std::array<double, 9> entries) :
}
}
-
-vector<double> DiffusionProblem::EvaluationResults(const Vector &solution) const {
+std::vector<double> DiffusionProblem::EvaluationResults(
+ const Vector &solution) const {
std::vector<double> eval{};
Vector error(solution);
@@ -38,7 +39,9 @@ DiffusionExampleOne::DiffusionExampleOne() : DiffusionProblem(
{1, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0}) {}
double DiffusionExampleOne::Potential(double time, const Point &x) const {
- return 0.25 * (1 - sin(time) * sin(time)) * (1 - cos(2 * Pi * x[0])) * (1 - cos(2 * Pi * x[1]));
+ using std::numbers::pi;
+ return 0.25 * (1 - sin(time) * sin(time)) * (1 - cos(2 * pi * x[0])) *
+ (1 - cos(2 * pi * x[1]));
}
double DiffusionExampleOne::IonCurrent(double time, const Point &x) const {
@@ -48,20 +51,23 @@ double DiffusionExampleOne::IonCurrent(double time, const Point &x) const {
}
double DiffusionExampleOne::timeDerivative(double time, const Point &x) const {
- return -0.5 * (sin(time) * cos(time)) * (1 - cos(2 * Pi * x[0])) * (1 - cos(2 * Pi * x[1]));
+ using std::numbers::pi;
+ return -0.5 * (sin(time) * cos(time)) * (1 - cos(2 * pi * x[0])) *
+ (1 - cos(2 * pi * x[1]));
}
double DiffusionExampleOne::divergence(int row, const Point &x) const {
+ using std::numbers::pi;
int i = row;
int j = (row + 1) % SpaceDimension;
int k = (row + 2) % SpaceDimension;
- return Pi * Pi *
- (conductivity[i][i] * cos(2 * Pi * x[i]) * (1 - cos(2 * Pi * x[j])) *
- (1 - cos(2 * Pi * x[k]))) *
- (conductivity[i][j] * sin(2 * Pi * x[j]) * (1 + sin(2 * Pi * x[i])) *
- (1 - cos(2 * Pi * x[k]))) *
- (conductivity[i][k] * sin(2 * Pi * x[k]) * (1 + sin(2 * Pi * x[i])) *
- (1 - cos(2 * Pi * x[j])));
+ return pi * pi *
+ (conductivity[i][i] * cos(2 * pi * x[i]) * (1 - cos(2 * pi * x[j])) *
+ (1 - cos(2 * pi * x[k]))) *
+ (conductivity[i][j] * sin(2 * pi * x[j]) * (1 + sin(2 * pi * x[i])) *
+ (1 - cos(2 * pi * x[k]))) *
+ (conductivity[i][k] * sin(2 * pi * x[k]) * (1 + sin(2 * pi * x[i])) *
+ (1 - cos(2 * pi * x[j])));
}
DiffusionExampleTwo::DiffusionExampleTwo() : DiffusionProblem(
@@ -71,7 +77,7 @@ double DiffusionExampleTwo::Potential(double time, const Point &x) const {
double v =
256 * (1 - sin(time) * sin(time)) * x[0] * x[0] * (x[0] - 1) * (x[0] - 1) * x[1] * x[1] *
(x[1] - 1) * (x[1] - 1) * x[2] * x[2] * (x[2] - 1) * (x[2] - 1);
- double g = 0.5 * (1 + cos(Pi * (x[1] - 0.5)));
+ double g = 0.5 * (1 + cos(std::numbers::pi * (x[1] - 0.5)));
return v * g;
}
diff --git a/src/electrophysiology/problem/DiffusionProblem.hpp b/src/electrophysiology/problem/DiffusionProblem.hpp
index 2e185210473bc10b3c5899f9483da0d4c3a3271f..a59f169f5633c714beb060f018dcf3d429932ea2 100644
--- a/src/electrophysiology/problem/DiffusionProblem.hpp
+++ b/src/electrophysiology/problem/DiffusionProblem.hpp
@@ -16,7 +16,7 @@ public:
return "Diffusion";
}
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
Tensor Conductivity(const Cell &c) const override;
};
diff --git a/src/electrophysiology/problem/ElphyProblem.cpp b/src/electrophysiology/problem/ElphyProblem.cpp
index ecc22ca74e9968239638b772765e990e61505e3f..48b5287fbb041aee2764f6f1538191bc4780d6df 100644
--- a/src/electrophysiology/problem/ElphyProblem.cpp
+++ b/src/electrophysiology/problem/ElphyProblem.cpp
@@ -23,7 +23,7 @@ GetElphyProblem() {
return GetElphyProblem(problemName);
}
-std::unique_ptr<ElphyProblem> GetElphyProblem(const string &problemName) {
+std::unique_ptr<ElphyProblem> GetElphyProblem(const std::string &problemName) {
if (problemName == "VariableMesh")
return std::make_unique<VariableMeshProblem>();
if (problemName == "DiffusionProblem1")
@@ -75,7 +75,6 @@ std::unique_ptr<ElphyProblem> GetElphyProblem(const string &problemName) {
if (problemName != "Default") Warning("Problem name unknown")
return std::make_unique<ElphyProblem>();
-
}
std::vector<double>
diff --git a/src/electrophysiology/problem/ElphyProblem.hpp b/src/electrophysiology/problem/ElphyProblem.hpp
index 0c4bc23c4b28292a16c79d59ab13654a8bfa2a04..534d0ea792bc6152c8a0caab1aa7c8d206f977c4 100644
--- a/src/electrophysiology/problem/ElphyProblem.hpp
+++ b/src/electrophysiology/problem/ElphyProblem.hpp
@@ -24,8 +24,8 @@ struct EvaluationPointData {
const int index;
};
-
-static std::vector<double> findPotential(const Vector &u, const vector<Point> &points) {
+static std::vector<double> findPotential(const Vector &u,
+ const std::vector<Point> &points) {
std::vector<double> potentials(points.size());
for (row r = u.rows(); r != u.rows_end(); ++r) {
auto [i, isNear] = isIn(r(), points);
@@ -139,7 +139,7 @@ public:
const Point &getEvaluationPointAt(int i) const { return evaluationPoints.at(i); }
- virtual void chamberSpecialization(const string &chamberName) const {}
+ virtual void chamberSpecialization(const std::string &chamberName) const {}
void
SetNumberOfEvaluationPointsBesideMesh(int number) { numberEvaluationPointsBesideMesh = number; };
@@ -166,11 +166,9 @@ public:
ElphyProblem(bool at, std::array<double, 4> sigmaParams)
: CardMechIProblem(), IElphyProblem(at, sigmaParams) {}
- [[nodiscard]] const string &GetMeshName() const { return meshesName; }
+ [[nodiscard]] const std::string &GetMeshName() const { return meshesName; }
- string Name() const override {
- return "Default Elphy Problem";
- }
+ std::string Name() const override { return "Default Elphy Problem"; }
[[nodiscard]] virtual std::string Evaluate(const Vector &solution) const { return ""; };
@@ -191,14 +189,12 @@ public:
void InitializeEvaluationPoints() override;
- string Name() const override {
- return "Monodomain";
- };
+ std::string Name() const override { return "Monodomain"; };
};
std::unique_ptr<ElphyProblem> GetElphyProblem();
-std::unique_ptr<ElphyProblem> GetElphyProblem(const string &problemName);
+std::unique_ptr<ElphyProblem> GetElphyProblem(const std::string &problemName);
#endif //ELPHYPROBLEM_H
diff --git a/src/electrophysiology/problem/OneCellProblem.hpp b/src/electrophysiology/problem/OneCellProblem.hpp
index de7a13285d905f5e3427f0d3f92a54463b5924d7..2c5e86a1ef7a19834609f54e173ca7fa0dd26161 100644
--- a/src/electrophysiology/problem/OneCellProblem.hpp
+++ b/src/electrophysiology/problem/OneCellProblem.hpp
@@ -21,7 +21,7 @@ public:
void InitializeEvaluationPoints() override;
- vector<double> EvaluationResults(const Vector &solution) const override;
+ std::vector<double> EvaluationResults(const Vector &solution) const override;
std::string Evaluate(const Vector &solution) const override;
diff --git a/src/electrophysiology/solvers/ElphySolver.cpp b/src/electrophysiology/solvers/ElphySolver.cpp
index 40eea2d4a4e2321b8cb5b78544a337540ecc34d4..6aed64ee84744324c5123f3c0040a5fba6f6e9c5 100644
--- a/src/electrophysiology/solvers/ElphySolver.cpp
+++ b/src/electrophysiology/solvers/ElphySolver.cpp
@@ -37,9 +37,10 @@ void ElphySolver::Method(IElphyAssemble &A, Vector &potential) {
}
-std::unique_ptr<ElphySolver> GetElphySolver(const string &solverName, MCellModel &cellModel,
+std::unique_ptr<ElphySolver> GetElphySolver(const std::string &solverName,
+ MCellModel &cellModel,
const Vector &excitationData) {
- string cmClass = "MElphyModel";
+ std::string cmClass = "MElphyModel";
Config::Get("ElphyModelClass", cmClass);
if (contains(cmClass, "IBT")) {
return std::make_unique<IBTSplittingSolver>();
@@ -64,7 +65,6 @@ std::unique_ptr<ElphySolver> GetElphySolver(const string &solverName, MCellModel
THROW(solverName + " No Solver for this Model implemented.")
}
-
std::unique_ptr<ElphySolver>
GetElphySolver(MCellModel &cellModel, const Vector &excitationData) {
std::string solverName{"SemiImplicit"};
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
index 178c7c52d6929e7652181c3a59237f0936a4f9b5..663bc0d711172a1860d1efab286683012bfdf527 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.cpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
@@ -228,7 +228,7 @@ void LinearImplicitSolver::PrintMaxMinGating() {
void LinearImplicitSolver::SolveConcentration(IElphyAssemble &A, double dt, const Vectors &values) {
std::vector<double> vcw(M.Size() + (M.Type() == TENSION));
//std::cout<<"M.Size() "<<M.Size()<<" ElphySize<"<<M.ElphySize()<< " size gating "<<(*gating).size()<<endl;
- std::vector<vector<double>> G(ionSteps);
+ std::vector<std::vector<double>> G(ionSteps);
for (int i = 0; i < ionSteps; ++i) {
G[i] = std::vector<double>(M.Size());
}
@@ -259,7 +259,7 @@ void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Ve
double dt = A.StepSize();
std::vector<double> vcwPerCell(M.Size() + (M.Type() == TENSION));
- std::vector<vector<double>> G(ionSteps);
+ std::vector<std::vector<double>> G(ionSteps);
for (int i = 0; i < ionSteps; ++i) {
G[i] = std::vector<double>(M.Size());
}
@@ -353,7 +353,8 @@ void LinearImplicitSolver::GatingVectorAti( Vector &values, int i)const{
values=(*gating)[i];
}
-void LinearImplicitSolver::selectIonScheme(const string &solvingIonScheme) {
+void LinearImplicitSolver::selectIonScheme(
+ const std::string &solvingIonScheme) {
if (solvingIonScheme == "ExplicitEuler") {
ionSteps = 1;
SolveConcentrationStep = ExplicitEulerIonStep;
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.hpp b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
index 02557a89f824e9eb624caaf1959fe9ba2c1f4e5b..511415a7c8adde5f0383cb936b1b6597bad543eb 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.hpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.hpp
@@ -51,7 +51,7 @@ protected:
std::shared_ptr<LagrangeDiscretization> disc0;
std::unique_ptr<Vectors> externalCurrentOnCells;
- void selectIonScheme(const string &solvingIonScheme);
+ void selectIonScheme(const std::string &solvingIonScheme);
std::function<void(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
std::vector<std::vector<double>> &g)> SolveConcentrationStep;
@@ -63,7 +63,7 @@ public:
max_estep(100000), check_conv(true), M(cellModel)//, M(cellModel.ElphyModel())
{
M.UseExpInt(true);
- string solvingIonScheme{"ExplicitEuler"};
+ std::string solvingIonScheme{"ExplicitEuler"};
Config::Get("IonScheme", solvingIonScheme);
selectIonScheme(solvingIonScheme);
diff --git a/src/electrophysiology/solvers/SplittingSolver.hpp b/src/electrophysiology/solvers/SplittingSolver.hpp
index 60d5959c819ad52b46db19c7be74e80c038a565c..ffa6d291e8733a3040b7593f4fbc9decde4c6353 100644
--- a/src/electrophysiology/solvers/SplittingSolver.hpp
+++ b/src/electrophysiology/solvers/SplittingSolver.hpp
@@ -19,8 +19,8 @@ enum SPLITTING {
STEIN = 3
};
-static SPLITTING getSplitting(const string &prefix = "Elphy") {
- string s = "Godunov";
+static SPLITTING getSplitting(const std::string &prefix = "Elphy") {
+ std::string s = "Godunov";
Config::Get(prefix + "SplittingMethod", s);
if (s == "OnlyODE") return ONLYODE;
if (s == "Godunov") return GODUNOV;
diff --git a/test/cellmodels/MTest.hpp b/test/cellmodels/MTest.hpp
index 2bd25c927a855141c08c84300ca06be60f00d82b..19783b10ac8ebec52db969259b6f22827b1da0da 100644
--- a/test/cellmodels/MTest.hpp
+++ b/test/cellmodels/MTest.hpp
@@ -1,6 +1,8 @@
#ifndef MTEST_HPP
#define MTEST_HPP
+#include <numbers>
+
#include "MCellModel.hpp"
class MTest : public MElphyModel {
@@ -33,11 +35,13 @@ public:
return t;
}
- void RHSUpdate(const vector<double> &var, vector<double> &g, double t) const override {
+ void RHSUpdate(const std::vector<double> &var, std::vector<double> &g,
+ double t) const override {
g[0] = Diff(t, var);
}
- void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override {
+ void UpdateValues(double dt, std::vector<double> &var,
+ const std::vector<double> &g) const override {
var[0] += dt * g[0];
}
@@ -54,14 +58,14 @@ public:
MExact() : MTest() {}
double Value(double t, const std::vector<double> &vw) const override {
- return sin(2 * Pi * t);
+ return sin(2 * std::numbers::pi * t);
}
double Diff(double t, const std::vector<double> &vw) const override {
- return 2 * Pi * cos(2 * Pi * t);
+ return 2 * std::numbers::pi * cos(2 * std::numbers::pi * t);
}
- void InitialValues(double t, vector<double> &vw) override {
+ void InitialValues(double t, std::vector<double> &vw) override {
vw[0] = Value(0, vw);
}
};
@@ -79,7 +83,7 @@ public:
return sin((vw[0] + t) * (vw[0] + t));
}
- void InitialValues(double t, vector<double> &vw) override {
+ void InitialValues(double t, std::vector<double> &vw) override {
vw[0] = -1.0;
}
};
diff --git a/test/cellmodels/TestElphyModels.cpp b/test/cellmodels/TestElphyModels.cpp
index 89276b673dc0a55e4b905a0e9b24b43d8e315d5c..c8d3f241a9c9391a5fcc443c89882c2760c0402a 100644
--- a/test/cellmodels/TestElphyModels.cpp
+++ b/test/cellmodels/TestElphyModels.cpp
@@ -36,8 +36,8 @@ protected:
cellSolver->Solve(TS, v);
}
-
- vector<double> GlobalError(const std::string &methodName, std::vector<double> &error) {
+ std::vector<double> GlobalError(const std::string &methodName,
+ std::vector<double> &error) {
int O = 14;
int N = 20;
@@ -66,7 +66,7 @@ protected:
double v = value[n - O][i];
double r = ref[dn * i];
- e = max(e, abs(v - r));
+ e = std::max(e, abs(v - r));
}
error[n - O] = e;
}
@@ -99,7 +99,7 @@ TYPED_TEST(MElphyTest, ConvergenceEuler) {
SelectExternalCurrentFunction("Arctan");
this->cellModel->SetExternal(0.0, 0.002,20.0);
- vector<double> error;
+ std::vector<double> error;
this->GlobalError("ExplicitEuler", error);
for (int i = 1; i < error.size(); ++i) {
if (error[i] < 1e-12) continue;
@@ -112,7 +112,7 @@ TYPED_TEST(MElphyTest, ConvergenceRK2) {
SelectExternalCurrentFunction("Arctan");
this->cellModel->SetExternal(0.0, 0.002,20.0);
- vector<double> error;
+ std::vector<double> error;
this->GlobalError("RungeKutta2", error);
for (int i = 1; i < error.size(); ++i) {
if (error[i] < 1e-12) continue;
@@ -124,7 +124,7 @@ TYPED_TEST(MElphyTest, ConvergenceRK4) {
SelectExternalCurrentFunction("Arctan");
this->cellModel->SetExternal(0.0, 0.002,20.0);
- vector<double> error;
+ std::vector<double> error;
this->GlobalError("RungeKutta4", error);
for (int i = 1; i < error.size(); ++i) {
if (error[i] < 1e-12) continue;
diff --git a/test/cellmodels/TestSolvers.cpp b/test/cellmodels/TestSolvers.cpp
index 6da3b67415ad46f732ad68fc60959d7577f75d9c..40acaf6bdf5f859fc4104ba677933acfb91818e2 100644
--- a/test/cellmodels/TestSolvers.cpp
+++ b/test/cellmodels/TestSolvers.cpp
@@ -67,7 +67,7 @@ TEST_P(MCellSolverTest, WithExactSolution) {
for (int i = 0; i < value[n - O].size(); ++i) {
double t = time[n - O][i];
double v = value[n - O][i];
- e = max(e, v - cellModel->Value(t, std::vector<double>{v}));
+ e = std::max(e, v - cellModel->Value(t, std::vector<double>{v}));
}
error[n - O] = e;
}
@@ -127,7 +127,7 @@ TEST_P(MCellSolverTest, WithoutExactSolution) {
double v = value[n - O][i];
double r = ref[dn * i];
- e = max(e, abs(v - r));
+ e = std::max(e, abs(v - r));
}
error[n - O] = e;
}
diff --git a/test/coupled/TestCoupledConsistency.cpp b/test/coupled/TestCoupledConsistency.cpp
index ccef0111663b9c83c60bb48eedf05af9d211d4ad..719423fe4d34d358df889f44f4693e8db7358988 100644
--- a/test/coupled/TestCoupledConsistency.cpp
+++ b/test/coupled/TestCoupledConsistency.cpp
@@ -27,7 +27,7 @@ protected:
std::unique_ptr<IElasticity> mechA = nullptr;
CoupledConsistencyTest() {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
initializeConfig(testConfig);
problem = GetCoupledProblem("ContractingBeam");
@@ -43,7 +43,7 @@ protected:
mechA->ResetTime(0.0, 0.1, 0.0001);
}
- void initializeConfig(std::map<string, string> &testConfig) {
+ void initializeConfig(std::map<std::string, std::string> &testConfig) {
testConfig["Model"] = "Monodomain";
testConfig["ElphyModel"] = GetParam().ElphyModel;
testConfig["ElphyLevel"] = "2";
@@ -87,7 +87,6 @@ protected:
Config::Initialize(testConfig);
}
-
void TearDown() override {
Config::Close();
}
diff --git a/test/coupled/TestDeformedDiffusion.cpp b/test/coupled/TestDeformedDiffusion.cpp
index 4438690203eff62073675778445e5ec8cd3544d7..c1fb8f3cd4d2f948e719f82e9f08b9c522bf760d 100644
--- a/test/coupled/TestDeformedDiffusion.cpp
+++ b/test/coupled/TestDeformedDiffusion.cpp
@@ -16,7 +16,7 @@ protected:
std::unique_ptr<CoupledProblem> deformedProblem = nullptr;
DeformedDiffusionTest() {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
initializeConfig(testConfig);
problem = std::make_unique<CompositeCoupledProblem>(
@@ -36,7 +36,7 @@ protected:
mechA = std::make_unique<LagrangeElasticity>(*problem, 1, false);
}
- void initializeConfig(std::map<string, string> &testConfig) {
+ void initializeConfig(std::map<std::string, std::string> &testConfig) {
testConfig["Model"] = "Monodomain";
testConfig["ElphyModel"] = "Diffusion";
testConfig["ElphyLevel"] = "2";
@@ -67,7 +67,6 @@ protected:
Config::Initialize(testConfig);
}
-
void TearDown() override {
Config::Close();
}
diff --git a/test/elasticity/TestAssembleConsistency.cpp b/test/elasticity/TestAssembleConsistency.cpp
index 5555b877dbe61587500bdf495cb1409cf2064dc8..32374a1407eb1f6a3e720340e7764a418a91d079 100644
--- a/test/elasticity/TestAssembleConsistency.cpp
+++ b/test/elasticity/TestAssembleConsistency.cpp
@@ -34,7 +34,7 @@ protected:
void SetUp() override {
- std::map<string, string> testConfig = CARDIACBEAM_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIACBEAM_TEST_CONFIG;
testConfig["MechDynamics"] = "Static";
testConfig["LAPressure"] = "0.004";
diff --git a/test/elasticity/TestCFLCondition.cpp b/test/elasticity/TestCFLCondition.cpp
index 5fafe030c67f99ac3178c8966fac934645bd7ade..e47fdb1406286362526800b2c6037c2743c762d6 100644
--- a/test/elasticity/TestCFLCondition.cpp
+++ b/test/elasticity/TestCFLCondition.cpp
@@ -44,7 +44,7 @@ protected:
}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Newmark";
testConfig["MechDiscretization"] = "Lagrange";
diff --git a/test/elasticity/TestCirculation.cpp b/test/elasticity/TestCirculation.cpp
index b28d92a1007d5318c3d74e52d17affe28e421246..82061aa107142079f206f313c42f1bfc167013f4 100644
--- a/test/elasticity/TestCirculation.cpp
+++ b/test/elasticity/TestCirculation.cpp
@@ -1,3 +1,5 @@
+#include <numbers>
+
#include <GlobalDefinitions.hpp>
#include "pressure_bnd/MCircModel.hpp"
#include "TestEnvironmentCardMech.hpp"
@@ -10,6 +12,8 @@ TEST(CirculationTest, SingleStep) {
TEST(CirculationTest, Periodic){
+ using std::numbers::pi;
+
MCircModel circModel;
TimeSeries TS(0.0, 1.0, 0.001);
@@ -25,23 +29,23 @@ TEST(CirculationTest, Periodic){
if (periodicTime >= 1.0) periodicTime -= 1.0;
if (periodicTime < 0.025) {
- pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * Pi / 0.025 * periodicTime) + 1);
+ pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * pi / 0.025 * periodicTime) + 1);
} else if (periodicTime < 0.05) {
- pressure[2] = pressure[3] = 5.0 * (cos(0.5 * Pi / 0.025 * (periodicTime - 0.025)) + 1);
+ pressure[2] = pressure[3] = 5.0 * (cos(0.5 * pi / 0.025 * (periodicTime - 0.025)) + 1);
} else if (periodicTime > 0.35) {
pressure[2] = pressure[3] = 0.0;
} else if (periodicTime > 0.3) {
- pressure[2] = pressure[3] = 5.0 * (cos(0.5 * Pi / 0.05 * (periodicTime - 0.3)) + 1);
+ pressure[2] = pressure[3] = 5.0 * (cos(0.5 * pi / 0.05 * (periodicTime - 0.3)) + 1);
} else if (periodicTime > 0.15) {
- pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * Pi / 0.15 * (periodicTime - 0.15)) + 1);
+ pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * pi / 0.15 * (periodicTime - 0.15)) + 1);
} else {
pressure[2] = pressure[3] = 0.0;
}
if (periodicTime < 0.15) {
- pressure[0] = pressure[1] = 60 * (-cos(0.5 * Pi / 0.15 * periodicTime) + 1);
+ pressure[0] = pressure[1] = 60 * (-cos(0.5 * pi / 0.15 * periodicTime) + 1);
} else if (periodicTime < 0.3) {
- pressure[0] = pressure[1] = 60.0 * (cos(0.5 * Pi / 0.15 * (periodicTime - 0.15)) + 1);
+ pressure[0] = pressure[1] = 60.0 * (cos(0.5 * pi / 0.15 * (periodicTime - 0.15)) + 1);
} else pressure[0] = pressure[1] = 0.0;
std::cout << "Volumes at t=" << t << " : " << std::endl;
diff --git a/test/elasticity/TestDirichletBeam.cpp b/test/elasticity/TestDirichletBeam.cpp
index 175dbc000cfc5a2a41dde14a469f63818db04b4f..e18216914c90fbd7570d507145ee2cd6671e20e4 100644
--- a/test/elasticity/TestDirichletBeam.cpp
+++ b/test/elasticity/TestDirichletBeam.cpp
@@ -22,7 +22,7 @@ protected:
DirichletBeamTest() = default;
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestDynamicBoundary.cpp b/test/elasticity/TestDynamicBoundary.cpp
index ee4883afa7018cc621a3e36c7f8cfffeb0f29df6..2527f761ad6476feb2776944b3711d40b4f9bc4f 100644
--- a/test/elasticity/TestDynamicBoundary.cpp
+++ b/test/elasticity/TestDynamicBoundary.cpp
@@ -38,7 +38,7 @@ protected:
}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Newmark";
testConfig["MechDiscretization"] = "Lagrange";
diff --git a/test/elasticity/TestElasticityBeam.cpp b/test/elasticity/TestElasticityBeam.cpp
index 1365a148a86ad3bb099e9d4cd7f27fe2da39cad9..620512d84d99bbf587742ab2a8cbd8674a0c06a9 100644
--- a/test/elasticity/TestElasticityBeam.cpp
+++ b/test/elasticity/TestElasticityBeam.cpp
@@ -19,7 +19,7 @@ protected:
ElasticityBeamTest() { }
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestElasticityBlock.cpp b/test/elasticity/TestElasticityBlock.cpp
index d825e191dd242abf4b296471df1a8b0b8c91adac..300116696a76e8ef8c6df62c50acbaf327fe474a 100644
--- a/test/elasticity/TestElasticityBlock.cpp
+++ b/test/elasticity/TestElasticityBlock.cpp
@@ -18,7 +18,7 @@ protected:
ElasticityBlockTest() {}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestLShape.cpp b/test/elasticity/TestLShape.cpp
index 8496adb3682e1a9652fca9a9ab48bba04f397637..2031209666830ace6466d68f8d02566a9d0eee53 100644
--- a/test/elasticity/TestLShape.cpp
+++ b/test/elasticity/TestLShape.cpp
@@ -8,7 +8,7 @@ protected:
MainElasticity *cmMain;
void SetUp() override {
- std::map<string, string> testConfig = lShapeConfig;
+ std::map<std::string, std::string> testConfig = lShapeConfig;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestLaplaceElasticity.cpp b/test/elasticity/TestLaplaceElasticity.cpp
index d2c9f65430733ca1b29e868cad7a2ecd2b964100..95bb8d655d472480f9949c4e41578ac34abce6d8 100644
--- a/test/elasticity/TestLaplaceElasticity.cpp
+++ b/test/elasticity/TestLaplaceElasticity.cpp
@@ -22,7 +22,7 @@ protected:
}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestLinearBeamProblem.cpp b/test/elasticity/TestLinearBeamProblem.cpp
index 3dc4adbc33b11592aff8bffca0edd4deaa4af514..de809a629750e95af652d3e991ec185ae5dc09f7 100644
--- a/test/elasticity/TestLinearBeamProblem.cpp
+++ b/test/elasticity/TestLinearBeamProblem.cpp
@@ -31,7 +31,7 @@ protected:
DirichletBeamTest() = default;
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestMainProgram.cpp b/test/elasticity/TestMainProgram.cpp
index bb2aa008cb83c82efad80f303a9e0656a35b75ac..695987134e61a3f6b23fdd3f9c65077c97661e2e 100644
--- a/test/elasticity/TestMainProgram.cpp
+++ b/test/elasticity/TestMainProgram.cpp
@@ -26,18 +26,20 @@ protected:
int level = mainProgram->initLevels[0];
for (auto &sampleAmount : mainProgram->initSampleAmount) {
for (int sample = 0; sample < sampleAmount; sample++) {
- string permeability = "sample_" + to_string(level) +
- "_" + to_string(sample) + "/permeability";
- string u = "sample_" + to_string(level) +
- "_" + to_string(sample) + "/u";
+ std::string permeability = "sample_" + std::to_string(level) +
+ "_" + std::to_string(sample) + "/permeability";
+ std::string u = "sample_" + std::to_string(level) +
+ "_" + std::to_string(sample) + "/u";
ASSERT_TRUE(FileExists("data/vtk/" + permeability + ".vtk"));
ASSERT_TRUE(FileExists("data/vtk/" + u + ".vtk"));
if (level != mainProgram->initLevels[0]) {
- string name_coarse = "sample_coarse_" + to_string(level) +
- "_" + to_string(sample) + "/permeability";
- string u_coarse = "sample_coarse_" + to_string(level) +
- "_" + to_string(sample) + "/u";
+ std::string name_coarse = "sample_coarse_" +
+ std::to_string(level) +
+ "_" + std::to_string(sample) + "/permeability";
+ std::string u_coarse = "sample_coarse_" + std::to_string(level)
+ +
+ "_" + std::to_string(sample) + "/u";
ASSERT_TRUE(FileExists("data/vtk/" + name_coarse + ".vtk"));
ASSERT_TRUE(FileExists("data/vtk/" + u_coarse + ".vtk"));
}
diff --git a/test/elasticity/TestPressureElasticity.hpp b/test/elasticity/TestPressureElasticity.hpp
index d4709301eef721c06ff627da24f480d539129e04..5a97e898f2b62d5aec4d0dde3daaa4b02bdb476f 100644
--- a/test/elasticity/TestPressureElasticity.hpp
+++ b/test/elasticity/TestPressureElasticity.hpp
@@ -35,7 +35,7 @@ protected:
TestPressureElasticity() = default;
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestPrestress.cpp b/test/elasticity/TestPrestress.cpp
index 2f06244d274f21f10c0aca4870ad4b50109d5fe3..e0921f19a0debd1ab70a4ce39ae9e1ef5d8876d6 100644
--- a/test/elasticity/TestPrestress.cpp
+++ b/test/elasticity/TestPrestress.cpp
@@ -20,7 +20,7 @@ protected:
PrestressTest() {}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechRuntype"] = "Default";
diff --git a/test/elasticity/TestQuadraticBeamProblem.cpp b/test/elasticity/TestQuadraticBeamProblem.cpp
index db65c17f4752e2ea6dbb0d89fc97b675858c6f2c..0057f1daf08cb264416ee966e9b3126f28947f38 100644
--- a/test/elasticity/TestQuadraticBeamProblem.cpp
+++ b/test/elasticity/TestQuadraticBeamProblem.cpp
@@ -43,7 +43,7 @@ protected:
QuadraticDirichletBeamTest() = default;
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "ActiveStrainElasticity";
testConfig["MechDynamics"] = "Static";
diff --git a/test/elasticity/TestViscoElasticity.cpp b/test/elasticity/TestViscoElasticity.cpp
index c28ef52ebd941baa29d64ed8f463776fbc3bbf51..9789068216551898b380db74fb2d3390cfc6b3d1 100644
--- a/test/elasticity/TestViscoElasticity.cpp
+++ b/test/elasticity/TestViscoElasticity.cpp
@@ -16,7 +16,7 @@ struct ViscoTestParameter {
class ViscoElasticTest : public TestWithParam<ViscoTestParameter> {
protected:
- std::map<string, string> testConfig;
+ std::map<std::string, std::string> testConfig;
std::unique_ptr<MainElasticity> cmMain;
void SetUp() override {
diff --git a/test/elasticity/materials/TestDefaultMaterial.hpp b/test/elasticity/materials/TestDefaultMaterial.hpp
index 57e27e97acea939713b98446d674da1bb624f084..1c5111d54ddcca8fffe455acf0485796ae3d7175 100644
--- a/test/elasticity/materials/TestDefaultMaterial.hpp
+++ b/test/elasticity/materials/TestDefaultMaterial.hpp
@@ -151,7 +151,8 @@ TYPED_TEST_P(MaterialTest, DerivativeConsistency) {
double isoder = mat.IsotropicDerivative(F, H) + mat.AnisotropicDerivative(F, One, H);
double numder = numMatDerivative(mat, F, One, H);
- EXPECT_NEAR(isoder, numder, 1e-6 * max(1.0, abs(max(isoder, numder))));
+ EXPECT_NEAR(isoder, numder,
+ 1e-6 * std::max(1.0, abs(std::max(isoder, numder))));
}
@@ -168,7 +169,8 @@ TYPED_TEST_P(MaterialTest, SecondDerivativeConsistency) {
mat.IsotropicSecondDerivative(F, H, G) + mat.AnisotropicSecondDerivative(F, One, H, G);
double numder = numMatSecondDerivative(mat, F, One, H, G);
- EXPECT_NEAR(isoder, numder, 1e-6 * max(1.0, abs(max(isoder, numder))));
+ EXPECT_NEAR(isoder, numder,
+ 1e-6 * std::max(1.0, abs(std::max(isoder, numder))));
}
REGISTER_TYPED_TEST_SUITE_P(MaterialTest, DefaultInstantiation, MoveInstantiation,
diff --git a/test/electrophysiology/TestDiffusion.cpp b/test/electrophysiology/TestDiffusion.cpp
index 1fc1a80184110063e993e5cfc600ca93fe5f0f9d..4495b8efb28b364ce3b2d27d9afff255220c0eed 100644
--- a/test/electrophysiology/TestDiffusion.cpp
+++ b/test/electrophysiology/TestDiffusion.cpp
@@ -25,7 +25,7 @@ protected:
DiffusionTest() {}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Model"] = "Monodomain";
testConfig["ElphyModel"] = "Diffusion";
diff --git a/test/electrophysiology/TestMonodomain.cpp b/test/electrophysiology/TestMonodomain.cpp
index 2b6f58692761572bbfaa1f892e6652dc45b7cda3..c72752b18f5c410e9d2efc73d51ff792211a522b 100644
--- a/test/electrophysiology/TestMonodomain.cpp
+++ b/test/electrophysiology/TestMonodomain.cpp
@@ -12,7 +12,7 @@ protected:
MonodomainTest() {}
void SetUp() override {
- std::map<string, string> testConfig = CARDIAC_TEST_CONFIG;
+ std::map<std::string, std::string> testConfig = CARDIAC_TEST_CONFIG;
testConfig["Distribution"] = "Stripes";
testConfig["ElphyPreconditioner"] = "SuperLU";
testConfig["Model"] = "Monodomain";
diff --git a/tools/geometries/vtu/ConvertCardiacVtu.cpp b/tools/geometries/vtu/ConvertCardiacVtu.cpp
index 97c98a297c5b126e3eb0845aee29287d808739e6..5ef2969e749cb80eb7084c8746ad482a876c719e 100644
--- a/tools/geometries/vtu/ConvertCardiacVtu.cpp
+++ b/tools/geometries/vtu/ConvertCardiacVtu.cpp
@@ -17,7 +17,8 @@ void smoothMeshData(const Mesh &mesh) {
}
if (max_exc_time >= 0.0 && max_amplitude > 0.0) {
for (int i = 0; i < c.Corners(); ++i) {
- double amplitude = min(mesh.find_vertex(c.Corner(i)).GetData()[2], max_amplitude);
+ double amplitude =
+ std::min(mesh.find_vertex(c.Corner(i)).GetData()[2], max_amplitude);
DataContainer exc_data({max_exc_time, max_duration, amplitude});
mesh.find_vertex(c.Corner(i)).SetData(exc_data);
}