diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 4d4d1b6c..a15d2f88 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -6,6 +6,7 @@ add_subdirectory("curvedflow_with_preinlet")
add_subdirectory("flowaroundsphere")
add_subdirectory("microcontraction")
add_subdirectory("oneCellShear")
+add_subdirectory("oneCellShear_fedosov")
add_subdirectory("parachuting")
add_subdirectory("parallelplanes")
add_subdirectory("pipeflow")
diff --git a/examples/oneCellShear_fedosov/CMakeLists.txt b/examples/oneCellShear_fedosov/CMakeLists.txt
new file mode 100644
index 00000000..34fda455
--- /dev/null
+++ b/examples/oneCellShear_fedosov/CMakeLists.txt
@@ -0,0 +1,14 @@
+# executable will have the same name as its directory
+get_filename_component(EXEC_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+# write the resulting executable in the _current_ directory
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+
+# add executable from source files in current directory assuming
+# to add more source files, manually register them using `add_executable`
+get_filename_component(MAIN ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+add_executable(${EXEC_NAME} "${CMAKE_CURRENT_SOURCE_DIR}/${MAIN}.cpp")
+
+# link the executable to `hemocell` and `HDF5` dependencies
+target_link_libraries(${EXEC_NAME} ${PROJECT_NAME})
+target_link_libraries(${EXEC_NAME} ${HDF5_C_HL_LIBRARIES} ${HDF5_LIBRARIES})
diff --git a/examples/oneCellShear_fedosov/RBC.pos b/examples/oneCellShear_fedosov/RBC.pos
new file mode 100644
index 00000000..ea9fa124
--- /dev/null
+++ b/examples/oneCellShear_fedosov/RBC.pos
@@ -0,0 +1,2 @@
+1
+9.5 9.5 4.5 90 0 0
diff --git a/examples/oneCellShear_fedosov/RBC.xml b/examples/oneCellShear_fedosov/RBC.xml
new file mode 100644
index 00000000..275cb667
--- /dev/null
+++ b/examples/oneCellShear_fedosov/RBC.xml
@@ -0,0 +1,20 @@
+
+
+
+ Parameters for the Fedosov RBC constitutive model.
+ RBC
+ 0.0
+ 1.0
+ 0
+ 640
+ 3.91e-6
+ 81
+ 0.4545
+ 2
+ 6.3e-6
+ 0.29
+ 0.01
+ 1000.0
+ 67.0
+
+
\ No newline at end of file
diff --git a/examples/oneCellShear_fedosov/clean.sh b/examples/oneCellShear_fedosov/clean.sh
new file mode 100644
index 00000000..19f7a1e8
--- /dev/null
+++ b/examples/oneCellShear_fedosov/clean.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+echo "Warning! This will remove all files produced during the CMake build, except the binary!"
+read -p "Are you sure? [y/n]" -n 1 -r
+echo # (optional) move to a new line
+if [[ $REPLY =~ ^[Yy]$ ]]
+then
+ echo '** Clearing palabos and hemocell compiled libraries...'
+ find ../../build/hemocell/. -type f ! -name 'CMakeLists.txt' -delete
+ echo '** Clearing local build...'
+ rm -r ./build
+fi
+
diff --git a/examples/oneCellShear_fedosov/compile.sh b/examples/oneCellShear_fedosov/compile.sh
new file mode 100644
index 00000000..2d09f6c1
--- /dev/null
+++ b/examples/oneCellShear_fedosov/compile.sh
@@ -0,0 +1,26 @@
+#!/bin/bash
+trap "exit" INT
+
+# This script invokes the compilation of the example present in the current
+# directory.
+
+echo "=========== Building =========="
+date
+
+example=${PWD}
+
+if [ ! -d "../../build" ]; then
+ echo "* Running CMake..."
+ mkdir ../../build
+ cd ../../build || exit 1
+ cmake ..
+ cd "$example" || exit 1
+fi
+
+echo "* Compiling..."
+cd ../../build || exit 1
+cmake --build . --target "${example##*/}"
+cd "$example" || exit 1
+
+date
+echo "=========== Done ==========="
diff --git a/examples/oneCellShear_fedosov/config.xml b/examples/oneCellShear_fedosov/config.xml
new file mode 100644
index 00000000..e2f02201
--- /dev/null
+++ b/examples/oneCellShear_fedosov/config.xml
@@ -0,0 +1,29 @@
+
+
+
+
+ 0
+
+
+
+
+ 3.91e-6
+
+
+
+ 111.0
+ 1025
+ 1.1e-6
+ 0.5e-6
+ 0.5e-7
+ 25
+ 4.100531391e-21
+
+
+
+ 100000
+ 2000
+ 500000
+
+
+
diff --git a/examples/oneCellShear_fedosov/oneCellShear_fedosov.cpp b/examples/oneCellShear_fedosov/oneCellShear_fedosov.cpp
new file mode 100644
index 00000000..7e99477c
--- /dev/null
+++ b/examples/oneCellShear_fedosov/oneCellShear_fedosov.cpp
@@ -0,0 +1,177 @@
+/*
+This file is part of the HemoCell library
+
+HemoCell is developed and maintained by the Computational Science Lab
+in the University of Amsterdam. Any questions or remarks regarding this library
+can be sent to: info@hemocell.eu
+
+When using the HemoCell library in scientific work please cite the
+corresponding paper: https://doi.org/10.3389/fphys.2017.00563
+
+The HemoCell library is free software: you can redistribute it and/or
+modify it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+The library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see .
+*/
+#include "hemocell.h"
+#include "rbcFedosovModel.h"
+#include "helper/hemocellInit.hh"
+#include "helper/cellInfo.h"
+
+#include "palabos3D.h"
+#include "palabos3D.hh"
+
+using namespace hemo;
+
+int main(int argc, char* argv[])
+{
+ if(argc < 2)
+ {
+ cout << "Usage: " << argv[0] << " " << endl;
+ return -1;
+ }
+
+ HemoCell hemocell(argv[1], argc, argv);
+ Config * cfg = hemocell.cfg;
+
+
+
+
+// ----------------- Read in config file & calc. LBM parameters ---------------------------
+ pcout << "(OneCellShear) (Parameters) calculating shear flow parameters" << endl;
+ plint nz = 10.0*(1e-6/(*cfg)["domain"]["dx"].read());
+ plint nx = 2*nz;
+ plint ny = 2*nz;
+ param::lbm_shear_parameters((*cfg),ny);
+ param::printParameters();
+
+ // ------------------------ Init lattice --------------------------------
+
+ pcout << "(CellStretch) Initializing lattice: " << nx <<"x" << ny <<"x" << nz << " [lu]" << std::endl;
+
+ plint extendedEnvelopeWidth = 2; // Because we might use ibmKernel with with 2.
+
+ hemocell.lattice = new MultiBlockLattice3D(
+ defaultMultiBlockPolicy3D().getMultiBlockManagement(nx, ny, nz, extendedEnvelopeWidth),
+ defaultMultiBlockPolicy3D().getBlockCommunicator(),
+ defaultMultiBlockPolicy3D().getCombinedStatistics(),
+ defaultMultiBlockPolicy3D().getMultiCellAccess(),
+ new GuoExternalForceBGKdynamics(1.0/param::tau));
+
+ pcout << "(OneCellShear) Re corresponds to u_max = " << (param::re * param::nu_p)/(hemocell.lattice->getBoundingBox().getNy()*param::dx) << " [m/s]" << endl;
+ // -------------------------- Define boundary conditions ---------------------
+
+ OnLatticeBoundaryCondition3D* boundaryCondition
+ = createLocalBoundaryCondition3D();
+
+ hemocell.lattice->toggleInternalStatistics(false);
+
+ iniLatticeSquareCouette_Y(*hemocell.lattice, nx, ny, nz, *boundaryCondition, param::shearrate_lbm);
+
+ hemocell.lattice->initialize();
+ hemocell.outputInSiUnits = true;
+
+ // ----------------------- Init cell models --------------------------
+
+ hemocell.initializeCellfield();
+ hemocell.addCellType("RBC", RBC_FROM_SPHERE);
+ vector outputs = {OUTPUT_POSITION,OUTPUT_TRIANGLES,OUTPUT_FORCE,OUTPUT_FORCE_VOLUME,OUTPUT_FORCE_BENDING,OUTPUT_FORCE_LINK,OUTPUT_FORCE_AREA, OUTPUT_FORCE_VISC};
+ hemocell.setOutputs("RBC", outputs);
+
+ outputs = {OUTPUT_VELOCITY};
+ hemocell.setFluidOutputs(outputs);
+
+ // -------------------------- Setting periodicity after initializing lattice and cell fields ---------------------
+ hemocell.setSystemPeriodicity(0, true);
+ hemocell.setSystemPeriodicity(2, true);
+
+// ---------------------- Initialise particle positions if it is not a checkpointed run ---------------
+
+ //loading the cellfield
+ if (not cfg->checkpointed) {
+ hemocell.loadParticles();
+ hemocell.writeOutput();
+ } else {
+ pcout << "(OneCellShear) CHECKPOINT found!" << endl;
+ hemocell.loadCheckPoint();
+ }
+
+
+ if (hemocell.iter == 0) {
+ pcout << "(OneCellShear) fresh start: warming up cell-free fluid domain for " << (*cfg)["parameters"]["warmup"].read() << " iterations..." << endl;
+ for (plint itrt = 0; itrt < (*cfg)["parameters"]["warmup"].read(); ++itrt) {
+ hemocell.lattice->collideAndStream();
+ }
+ }
+
+ pcout << "(OneCellShea) Shear rate: " << (*cfg)["domain"]["shearrate"].read() << " s^-1." << endl;
+
+ unsigned int tmax = (*cfg)["sim"]["tmax"].read();
+ unsigned int tmeas = (*cfg)["sim"]["tmeas"].read();
+ unsigned int tcheckpoint = (*cfg)["sim"]["tcheckpoint"].read();
+
+ // Get undeformed cell values
+ CellInformationFunctionals::calculateCellVolume(&hemocell);
+ CellInformationFunctionals::calculateCellArea(&hemocell);
+ T volume_eq = (CellInformationFunctionals::info_per_cell[0].volume)/pow(1e-6/param::dx,3);
+ T surface_eq = (CellInformationFunctionals::info_per_cell[0].area)/pow(1e-6/param::dx,2);
+
+ T D0 = 2.0 * (*cfg)["ibm"]["radius"].read() * 1e6;
+
+ // Creating output log file
+ plb_ofstream fOut;
+ if(cfg->checkpointed)
+ fOut.open("stretch.log", std::ofstream::app);
+ else
+ fOut.open("stretch.log");
+
+
+ while (hemocell.iter < tmax ) {
+
+ hemocell.iterate();
+
+ if (hemocell.iter % tmeas == 0) {
+ hemocell.writeOutput();
+
+ // Fill up the static info structure with desired data
+ CellInformationFunctionals::calculateCellVolume(&hemocell);
+ CellInformationFunctionals::calculateCellArea(&hemocell);
+ CellInformationFunctionals::calculateCellPosition(&hemocell);
+ CellInformationFunctionals::calculateCellStretch(&hemocell);
+ CellInformationFunctionals::calculateCellBoundingBox(&hemocell);
+
+ T volume = (CellInformationFunctionals::info_per_cell[0].volume)/pow(1e-6/param::dx,3);
+ T surface = (CellInformationFunctionals::info_per_cell[0].area)/pow(1e-6/param::dx,2);
+ hemo::Array position = CellInformationFunctionals::info_per_cell[0].position/(1e-6/param::dx);
+ hemo::Array bbox = CellInformationFunctionals::info_per_cell[0].bbox/(1e-6/param::dx);
+ T largest_diam = (CellInformationFunctionals::info_per_cell[0].stretch)/(1e-6/param::dx);
+ T rel_D2 = (largest_diam/D0)*(largest_diam/D0);
+ T def_idx = (rel_D2 - 1.0) / (rel_D2 + 1.0) * 100.0;
+
+ pcout << "\t Cell center at: {" <& lattice,
lattice.initialize();
}
+template class Descriptor>
+void iniLatticeSquareCouette_Y(plb::MultiBlockLattice3D& lattice,
+ plint nx, plint ny, plint nz,
+ plb::OnLatticeBoundaryCondition3D& boundaryCondition, T shearRate)
+{
+ plb::Box3D top = plb::Box3D(0, nx-1, ny-1, ny-1, 0, nz-1);
+ plb::Box3D bottom = plb::Box3D(0, nx-1, 0, 0, 0, nz-1);
+
+ boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top );
+ boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, bottom );
+
+ T vHalf = (ny-1)*shearRate*0.5;
+ setBoundaryVelocity(lattice, top, plb::Array(vHalf,0.0,0.0));
+ setBoundaryVelocity(lattice, bottom, plb::Array(-vHalf,0.0,0.0));
+
+ setExternalVector( lattice, lattice.getBoundingBox(),
+ Descriptor::ExternalField::forceBeginsAt, plb::Array::d>(0.0,0.0,0.0));
+
+ lattice.initialize();
+}
+
#endif // HEMOCELLINIT_HH
diff --git a/mechanics/cellMechanics.h b/mechanics/cellMechanics.h
index cbebbc78..a3c9d0c7 100644
--- a/mechanics/cellMechanics.h
+++ b/mechanics/cellMechanics.h
@@ -76,6 +76,36 @@ class CellMechanics {
T calculate_etaM(Config & cfg ){
return cfg["MaterialModel"]["eta_m"].read() * param::dx / param::dt / param::df;
};
+
+ T calculate_kd(Config & cfg, plb::MeshMetrics & meshmetric){
+ T kd = cfg["MaterialModel"]["kd"].read();
+ return kd * param::dx / param::df; // converting to LBM units
+ };
+
+ T calculate_ka(Config & cfg, plb::MeshMetrics & meshmetric){
+ T ka = cfg["MaterialModel"]["ka"].read();
+ return ka * param::dx / param::df; // converting to LBM units
+ };
+
+ T calculate_kv(Config & cfg, plb::MeshMetrics & meshmetric){
+ T kv = cfg["MaterialModel"]["kv"].read();
+ return kv * param::dx * param::dx / param::df; // converting to LBM units
+ };
+
+ T calculate_kb(Config & cfg, plb::MeshMetrics & meshmetric){
+ T kb = cfg["MaterialModel"]["kb"].read();
+ return kb * param::kBT_lbm; // converting to LBM units
+ };
+
+ T calculate_x0(Config & cfg, plb::MeshMetrics & meshmetric){
+ T x0 = cfg["MaterialModel"]["x0"].read();
+ return x0;
+ };
+
+ T calculate_m(Config & cfg, plb::MeshMetrics & meshmetric){
+ T m = cfg["MaterialModel"]["m"].read();
+ return m;
+ };
};
}
#endif
diff --git a/mechanics/commonCellConstants.cpp b/mechanics/commonCellConstants.cpp
index 27be5978..68146bee 100644
--- a/mechanics/commonCellConstants.cpp
+++ b/mechanics/commonCellConstants.cpp
@@ -40,10 +40,15 @@ CommonCellConstants::CommonCellConstants(HemoCellField & cellField_,
vector vertex_n_vertexes_,
vector,6>> vertex_outer_edges_per_vertex_,
vector,6>> vertex_outer_edges_per_vertex_sign_,
- T volume_eq_, T area_mean_eq_,
+ T volume_eq_, T surface_eq_, T area_mean_eq_,
T edge_mean_eq_, T angle_mean_eq_,
vector> inner_edge_list_,
- vector inner_edge_length_eq_list_) :
+ vector inner_edge_length_eq_list_,
+ vector edge_length_max_list_,
+ vector edge_kp_list_,
+ vector edge_kLinkWLC_list_,
+ vector edge_sin_angle_eq_list_,
+ vector edge_cos_angle_eq_list_) :
cellField(cellField_),
triangle_list(triangle_list_),
edge_list(edge_list_),
@@ -60,12 +65,17 @@ CommonCellConstants::CommonCellConstants(HemoCellField & cellField_,
vertex_outer_edges_per_vertex(vertex_outer_edges_per_vertex_),
vertex_outer_edges_per_vertex_sign(vertex_outer_edges_per_vertex_sign_),
volume_eq(volume_eq_),
+ surface_eq(surface_eq_),
area_mean_eq(area_mean_eq_),
edge_mean_eq(edge_mean_eq_),
angle_mean_eq(angle_mean_eq_),
inner_edge_list(inner_edge_list_),
- inner_edge_length_eq_list(inner_edge_length_eq_list_)
- {};
+ inner_edge_length_eq_list(inner_edge_length_eq_list_),
+ edge_length_max_list(edge_length_max_list_),
+ edge_kp_list(edge_kp_list_),
+ edge_kLinkWLC_list(edge_kLinkWLC_list_),
+ edge_sin_angle_eq_list(edge_sin_angle_eq_list_),
+ edge_cos_angle_eq_list(edge_cos_angle_eq_list_) {};
CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCellField & cellField_, Config & modelCfg_) {
HemoCellField & cellField = cellField_;
@@ -77,7 +87,6 @@ CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCell
//TODO, every edge is in here twice, clean that?
vector> edge_list_;
-
for (const hemo::Array & triangle : triangle_list_) {
//FIX by allowing only incremental edges
@@ -98,6 +107,41 @@ CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCell
edge_length_eq_list_.push_back(cellField.meshElement->computeEdgeLength(edge[0],edge[1]));
}
+ //Calculate max edges lengths
+ T x0;
+ vector edge_length_max_list_;
+ try {
+ x0 = modelCfg_["MaterialModel"]["x0"].read();
+ for (pluint iEdge = 0 ; iEdge < edge_list_.size(); iEdge++) {
+ edge_length_max_list_.push_back(edge_length_eq_list_[iEdge] / x0);
+ }
+ } catch (std::invalid_argument & exeption) {}
+
+ //Calculate kp and k_link_WLC for each edge
+ T m;
+ T mu0;
+ vector edge_kp_list_;
+ vector edge_kLinkWLC_list_;
+ try {
+ m = modelCfg_["MaterialModel"]["m"].read();
+ mu0 = modelCfg_["MaterialModel"]["mu0"].read();
+ for (pluint i = 0 ; i < edge_list_.size(); i++) {
+ T A = std::sqrt(3) * (m + 1) / 4.0 / std::pow(edge_length_eq_list_[i] * param::dx, m + 1);
+ T B = std::sqrt(3) * param::kBT_p * (x0 / 2.0 / std::pow(1 - x0, 3) - 1.0 / 4.0 / std::pow(1 - x0, 2) + 0.25) / 4.0 / (edge_length_max_list_[i] * param::dx) / x0;
+ T C = param::kBT_p * (1 / (4.0 * std::pow(1 - x0, 2)) - 0.25 + x0) * std::pow(edge_length_eq_list_[i] * param::dx, m);
+ T persistenceLengthCoarse = (B + A * C) / mu0;
+
+ T kp = C / persistenceLengthCoarse; // in pyhical units
+ T kLinkWLC = param::kBT_p / persistenceLengthCoarse; // in pyhical units
+
+ kp *= (1/param::df) * (1/param::dx) * (1/param::dx); // Converting to LBM units
+ kLinkWLC *= (1/param::df); // Converting to LBM units
+
+ edge_kp_list_.push_back(kp);
+ edge_kLinkWLC_list_.push_back(kLinkWLC);
+ }
+ } catch (std::invalid_argument & exeption) {}
+
//Calculate eq edges angles
vector edge_angle_eq_list_;
hemo::Array x2 = {0.0,0.0,0.0};
@@ -139,6 +183,15 @@ CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCell
edge_angle_eq_list_.push_back(angle);
}
+ //Calculate eq edges angles sin and cos
+ std::vector edge_sin_angle_eq_list_;
+ std::vector edge_cos_angle_eq_list_;
+ for (size_t i = 0; i < edge_angle_eq_list_.size(); ++i) {
+ T angle = edge_angle_eq_list_[i];
+ edge_sin_angle_eq_list_.push_back(std::sin(angle));
+ edge_cos_angle_eq_list_.push_back(std::cos(angle));
+ }
+
// Define opposite points TODO: make it automatic / or add the 33 links version ;)
vector> inner_edge_list_;
try {
@@ -168,6 +221,7 @@ CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCell
T volume_eq_ = MeshMetrics(*cellField.meshElement).getVolume();
+ T surface_eq_ = MeshMetrics(*cellField.meshElement).getSurface();
//store important points for bending calculation
@@ -400,11 +454,17 @@ CommonCellConstants CommonCellConstants::CommonCellConstantsConstructor(HemoCell
vertex_outer_edges_per_vertex,
vertex_outer_edges_per_vertex_sign,
volume_eq_,
+ surface_eq_,
mean_area_eq_,
mean_edge_eq_,
mean_angle_eq_,
inner_edge_list_,
- inner_edge_length_eq_list_);
+ inner_edge_length_eq_list_,
+ edge_length_max_list_,
+ edge_kp_list_,
+ edge_kLinkWLC_list_,
+ edge_sin_angle_eq_list_,
+ edge_cos_angle_eq_list_);
return CCC;
};
diff --git a/mechanics/commonCellConstants.h b/mechanics/commonCellConstants.h
index 99cb0ff1..90717423 100644
--- a/mechanics/commonCellConstants.h
+++ b/mechanics/commonCellConstants.h
@@ -53,10 +53,15 @@ class CommonCellConstants {
std::vector vertex_n_vertexes_,
std::vector,6>> vertex_outer_edges_per_vertex_,
std::vector,6>> vertex_outer_edges_per_vertex_sign_,
- T volume_eq_, T area_mean_eq_,
+ T volume_eq_, T surface_eq_, T area_mean_eq_,
T edge_mean_eq_, T angle_mean_eq_,
std::vector> inner_edge_list_,
- std::vector inner_edge_length_eq_list_);
+ std::vector inner_edge_length_eq_list_,
+ std::vector edge_length_max_list_,
+ std::vector edge_kp_list_,
+ std::vector edge_kLinkWLC_list_,
+ std::vector edge_sin_angle_eq_list_,
+ std::vector edge_cos_angle_eq_list_);
public:
static CommonCellConstants CommonCellConstantsConstructor(HemoCellField &, hemo::Config & modelCfg_);
@@ -78,12 +83,18 @@ class CommonCellConstants {
const std::vector,6>> vertex_outer_edges_per_vertex_sign;
const T volume_eq;
+ const T surface_eq;
const T area_mean_eq;
const T edge_mean_eq;
const T angle_mean_eq;
const std::vector> inner_edge_list;
const std::vector inner_edge_length_eq_list;
+ const std::vector edge_length_max_list;
+ const std::vector edge_kp_list;
+ const std::vector edge_kLinkWLC_list;
+ const std::vector edge_sin_angle_eq_list;
+ const std::vector edge_cos_angle_eq_list;
};
}
#endif
diff --git a/mechanics/rbcFedosovModel.cpp b/mechanics/rbcFedosovModel.cpp
new file mode 100644
index 00000000..6da445f3
--- /dev/null
+++ b/mechanics/rbcFedosovModel.cpp
@@ -0,0 +1,263 @@
+/*
+This file is part of the HemoCell library
+
+HemoCell is developed and maintained by the Computational Science Lab
+in the University of Amsterdam. Any questions or remarks regarding this library
+can be sent to: info@hemocell.eu
+
+When using the HemoCell library in scientific work please cite the
+corresponding paper: https://doi.org/10.3389/fphys.2017.00563
+
+The HemoCell library is free software: you can redistribute it and/or
+modify it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+The library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see .
+*/
+#include "rbcFedosovModel.h"
+#include "logfile.h"
+
+namespace hemo {
+
+RbcFedosovModel::RbcFedosovModel(Config & modelCfg_, HemoCellField & cellField_) : CellMechanics(cellField_, modelCfg_),
+ cellField(cellField_),
+ eta_m( RbcFedosovModel::calculate_etaM(modelCfg_) ),
+ kd( RbcFedosovModel::calculate_kd(modelCfg_,*cellField_.meshmetric) ),
+ ka( RbcFedosovModel::calculate_ka(modelCfg_,*cellField_.meshmetric) ),
+ kv( RbcFedosovModel::calculate_kv(modelCfg_,*cellField_.meshmetric) ),
+ kb( RbcFedosovModel::calculate_kb(modelCfg_,*cellField_.meshmetric) ),
+ x0( RbcFedosovModel::calculate_x0(modelCfg_,*cellField_.meshmetric) ),
+ m( RbcFedosovModel::calculate_m(modelCfg_,*cellField_.meshmetric)) {};
+
+void RbcFedosovModel::ParticleMechanics(map> & particles_per_cell, const map & lpc, size_t ctype) {
+
+ for (const auto & pair : lpc) { //For all cells with at least one lsp in the local domain.
+
+ const int & cid = pair.first;
+
+ vector & cell = particles_per_cell[cid];
+
+ if (cell.size() == 0) continue;
+ if (cell[0]->sv.celltype != ctype) continue; //only execute on correct particle
+
+ //Calculate Cell Values that need all particles
+ T volume = 0.0;
+ vector triangle_areas;
+ triangle_areas.reserve(cellConstants.triangle_list.size());
+ vector> triangle_normals;
+ triangle_normals.reserve(cellConstants.triangle_list.size());
+
+ // Volume calculation
+ for (const hemo::Array & triangle : cellConstants.triangle_list) {
+ const hemo::Array & v0 = cell[triangle[0]]->sv.position;
+ const hemo::Array & v1 = cell[triangle[1]]->sv.position;
+ const hemo::Array & v2 = cell[triangle[2]]->sv.position;
+
+ //Volume
+ const T v210 = v2[0]*v1[1]*v0[2];
+ const T v120 = v1[0]*v2[1]*v0[2];
+ const T v201 = v2[0]*v0[1]*v1[2];
+ const T v021 = v0[0]*v2[1]*v1[2];
+ const T v102 = v1[0]*v0[1]*v2[2];
+ const T v012 = v0[0]*v1[1]*v2[2];
+ volume += (-v210+v120+v201-v021-v102+v012);
+ }
+ volume *= (1.0/6.0);
+
+ // Calculating the current surface area
+ T surface = 0.0;
+ for (const hemo::Array & triangle : cellConstants.triangle_list) {
+ const hemo::Array & v0 = cell[triangle[0]]->sv.position;
+ const hemo::Array & v1 = cell[triangle[1]]->sv.position;
+ const hemo::Array & v2 = cell[triangle[2]]->sv.position;
+
+ //Area
+ T area_tmp = 0.0;
+ area_tmp = computeTriangleArea(v0, v1, v2);
+ surface += area_tmp;
+ }
+
+ // Per-triangle calculations
+ int triangle_n = 0;
+ for (const hemo::Array & triangle : cellConstants.triangle_list) {
+ const hemo::Array & v0 = cell[triangle[0]]->sv.position;
+ const hemo::Array & v1 = cell[triangle[1]]->sv.position;
+ const hemo::Array & v2 = cell[triangle[2]]->sv.position;
+
+ //Area
+ T area;
+ hemo::Array t_normal;
+ computeTriangleAreaAndUnitNormal(v0, v1, v2, area, t_normal);
+
+ // ======================================= AREA =======================================
+ T alpha_a = -ka * (surface - cellConstants.surface_eq) / (4 * cellConstants.surface_eq * area);
+ T alpha_d = -kd * (area - cellConstants.triangle_area_eq_list[triangle_n]) / (4 * cellConstants.triangle_area_eq_list[triangle_n] * area);
+
+ hemo::Array v10 = v1 - v0;
+ hemo::Array v02 = v0 - v2;
+ hemo::Array v21 = v2 - v1;
+
+ hemo::Array si_k = crossProduct(v10, -v02);
+
+ hemo::Array force0localGlobal = (alpha_d + alpha_a) * crossProduct(si_k, v21);
+ hemo::Array force1localGlobal = (alpha_d + alpha_a) * crossProduct(si_k, v02);
+ hemo::Array force2localGlobal = (alpha_d + alpha_a) * crossProduct(si_k, v10);
+
+ *cell[triangle[0]]->force_area += force0localGlobal;
+ *cell[triangle[1]]->force_area += force1localGlobal;
+ *cell[triangle[2]]->force_area += force2localGlobal;
+
+ //Store values necessary later
+ triangle_areas.push_back(area);
+ triangle_normals.push_back(t_normal);
+
+ // ======================================= VOLUME =======================================
+ // Calculate centroid of the triangle
+ hemo::Array tc_k = (v0 + v1 + v2) / 3.0;
+
+ T beta_v = -kv * (volume - cellConstants.volume_eq) / cellConstants.volume_eq;
+ hemo::Array force0v = beta_v * (si_k / 3.0 + crossProduct(tc_k, v21));
+ hemo::Array force1v = beta_v * (si_k / 3.0 + crossProduct(tc_k, v02));
+ hemo::Array force2v = beta_v * (si_k / 3.0 + crossProduct(tc_k, v10));
+
+ *cell[triangle[0]]->force_volume += force0v;
+ *cell[triangle[1]]->force_volume += force1v;
+ *cell[triangle[2]]->force_volume += force2v;
+
+#ifdef INTERIOR_VISCOSITY
+ // Add the normal direction here, always pointing outward
+ // const hemo::Array local_normal_dir = (triangle_normals[triangle_n])*(triangle_areas[triangle_n]/cellConstants.area_mean_eq);
+ const hemo::Array local_normal_dir = t_normal * (area / cellConstants.area_mean_eq);
+ cell[triangle[0]]->normalDirection += local_normal_dir;
+ cell[triangle[1]]->normalDirection += local_normal_dir;
+ cell[triangle[2]]->normalDirection += local_normal_dir;
+#endif
+
+ triangle_n++;
+ }
+
+ // Per-edge calculations ===================== LINK ==================================================
+
+ int edge_n=0;
+ for (const hemo::Array & edge : cellConstants.edge_list) {
+
+ hemo::Array & p0 = cell[edge[0]]->sv.position;
+ hemo::Array & p1 = cell[edge[1]]->sv.position;
+
+ hemo::Array v32 = p1 - p0;
+ T v32_mag = norm(v32);
+ hemo::Array v32_norm = v32 / v32_mag;
+
+ // Fedosov model
+ T edge_frac = (v32_mag - cellConstants.edge_length_eq_list[edge_n]) / cellConstants.edge_length_eq_list[edge_n];
+ T x_frac = x0 * (edge_frac + 1);
+
+ // Cutting off x_frac to avoid numerical problems
+ if (x_frac > 0.999) {
+ x_frac = 0.99;
+ }
+
+ T edge_force_WLC_mag = cellConstants.edge_kLinkWLC_list[edge_n] * (1.0 / (4.0 * (1.0 - x_frac) * (1.0 - x_frac)) - 1.0 / 4.0 + x_frac);
+ T edge_force_POW_mag = cellConstants.edge_kp_list[edge_n] / std::pow(x_frac * cellConstants.edge_length_max_list[edge_n], m);
+ T edge_force_mag = edge_force_WLC_mag - edge_force_POW_mag;
+ hemo::Array edge_force = v32_norm * edge_force_mag;
+ *cell[edge[0]]->force_link += edge_force;
+ *cell[edge[1]]->force_link -= edge_force;
+
+ // ======================================= Viscosity =======================================
+
+ if (eta_m != 0.0) {
+ // Membrane viscosity of bilipid layer
+ // F = eta * (dv/l) * l.
+ const hemo::Array rel_vel = cell[edge[1]]->sv.v - cell[edge[0]]->sv.v;
+ const hemo::Array rel_vel_projection = dot(rel_vel, v32_norm) * v32_norm;
+ hemo::Array Fvisc_memb = eta_m * rel_vel_projection;
+
+ // Limit membrane viscosity
+ const T Fvisc_memb_mag = norm(Fvisc_memb);
+ if (Fvisc_memb_mag > FORCE_LIMIT / 4.0) {
+ Fvisc_memb *= (FORCE_LIMIT / 4.0) / Fvisc_memb_mag;
+ }
+
+ *cell[edge[0]]->force_visc += Fvisc_memb;
+ *cell[edge[1]]->force_visc -= Fvisc_memb;
+ }
+
+ // ======================================= Bending =======================================
+ hemo::Array v1 = cell[cellConstants.edge_bending_triangles_outer_points[edge_n][0]]->sv.position;
+ hemo::Array v2 = p0;
+ hemo::Array v3 = p1;
+ hemo::Array v4 = cell[cellConstants.edge_bending_triangles_outer_points[edge_n][1]]->sv.position;
+
+ // hemo::Array v32 = v3 - v2; // Already calculated above
+ hemo::Array v13 = v1 - v3;
+ hemo::Array v34 = v3 - v4;
+ hemo::Array v42 = v4 - v2;
+ hemo::Array v21 = v2 - v1;
+
+ // to be used for theta sign determination
+ hemo::Array v14 = v1 - v4;
+
+ hemo::Array si = crossProduct(v21, -v13);
+ hemo::Array zeta = crossProduct(v34, -v42);
+
+ T si_mag = norm(si);
+ T zeta_mag = norm(zeta);
+
+ hemo::Array si_norm = si / si_mag;
+ hemo::Array zeta_norm = zeta / zeta_mag;
+
+ // USERMESO implementation for stability
+ T cosine_theta = dot(si_norm, zeta_norm);
+ if (cosine_theta > 1.0) cosine_theta = 1.0;
+ if (cosine_theta < -1.0) cosine_theta = -1.0;
+ T sine_theta = std::sqrt(1 - cosine_theta * cosine_theta);
+ if (sine_theta < T(1.0e-6)) sine_theta = 1.0e-6;
+ T mx = dot((si_norm - zeta_norm), v14);
+ if (mx < 0.0) sine_theta = -sine_theta;
+ T theta = std::atan2(sine_theta, cosine_theta);
+
+ // Fedosov method
+ T beta = kb * (sine_theta * cellConstants.edge_cos_angle_eq_list[edge_n] - cosine_theta * cellConstants.edge_sin_angle_eq_list[edge_n]) /
+ sine_theta;
+ T dTheta = theta - cellConstants.edge_angle_eq_list[edge_n];
+
+ T b11 = -beta * cosine_theta / (si_mag * si_mag);
+ T b12 = beta / (si_mag * zeta_mag);
+ T b22 = -beta * cosine_theta / (zeta_mag * zeta_mag);
+
+ hemo::Array force1 = b11 * crossProduct(si, v32) + b12 * crossProduct(zeta, v32);
+ hemo::Array force2 = b11 * crossProduct(si, v13) + b12 * (crossProduct(si, v34) + crossProduct(zeta, v13)) + b22 * crossProduct(zeta, v34);
+ hemo::Array force3 = b11 * crossProduct(si, v21) + b12 * (crossProduct(si, v42) + crossProduct(zeta, v21)) + b22 * crossProduct(zeta, v42);
+ hemo::Array force4 = b12 * crossProduct(si, -v32) + b22 * crossProduct(zeta, -v32);
+
+ // Apply forces to the vertices
+ *cell[cellConstants.edge_bending_triangles_outer_points[edge_n][0]]->force_bending += force1;
+ *cell[edge[0]]->force_bending += force2;
+ *cell[edge[1]]->force_bending += force3;
+ *cell[cellConstants.edge_bending_triangles_outer_points[edge_n][1]]->force_bending += force4;
+
+ edge_n++;
+ }
+ }
+};
+
+void RbcFedosovModel::statistics() {
+ hlog << "(Cell-mechanics model) Fedosov model parameters for " << cellField.name << " cellfield" << std::endl;
+ hlog << "\t eta_m: " << eta_m << std::endl;
+ hlog << "\t kd: " << kd << std::endl;
+ hlog << "\t ka: " << ka << std::endl;
+ hlog << "\t kv: " << kv << std::endl;
+ hlog << "\t kb: " << kb << std::endl;
+ hlog << "\t mean_edge:" << cellConstants.edge_mean_eq << std::endl;
+ hlog << "\t N faces: " << cellConstants.triangle_list.size() << std::endl;
+};
+
+}
diff --git a/mechanics/rbcFedosovModel.h b/mechanics/rbcFedosovModel.h
new file mode 100644
index 00000000..cf3dee24
--- /dev/null
+++ b/mechanics/rbcFedosovModel.h
@@ -0,0 +1,57 @@
+/*
+This file is part of the HemoCell library
+
+HemoCell is developed and maintained by the Computational Science Lab
+in the University of Amsterdam. Any questions or remarks regarding this library
+can be sent to: info@hemocell.eu
+
+When using the HemoCell library in scientific work please cite the
+corresponding paper: https://doi.org/10.3389/fphys.2017.00563
+
+The HemoCell library is free software: you can redistribute it and/or
+modify it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+The library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see .
+*/
+#ifndef HEMOCELL_RBCFEDOSOVMODEL_H
+#define HEMOCELL_RBCFEDOSOVMODEL_H
+
+#include "config.h"
+#include "cellMechanics.h"
+#include "hemoCellField.h"
+#include "constant_defaults.h"
+
+namespace hemo {
+
+class RbcFedosovModel : public CellMechanics {
+
+ public:
+ //Variables
+ HemoCellField & cellField;
+
+ const T eta_m;
+ const T kd;
+ const T ka;
+ const T kv;
+ const T kb;
+ const T x0;
+ const T m;
+
+ public:
+ RbcFedosovModel(Config & modelCfg_, HemoCellField & cellField_);
+
+ private:
+ void ParticleMechanics(map> & particles_per_cell, const map &lpc, size_t ctype) ;
+ void statistics();
+
+};
+}
+#endif