Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/deltaMechanicalForces.cpp
b7305f7b5a8bmaster

Last Change
Raw File

deltaMechanicalForces.cpp
View Options

#include "deltaMechanicalForces.hpp"
using namespace std;
void DeltaMechanicalForces::initParameters(){
we_min= 0.01;
we_r = 10.0;
hill_co = 4.0;
h_pressure = 2;
HC_pressure = 10;
k_sig=0.3;
k_pressure=0;
k_tension=1;
max_pert = 0.01;
small_cell=false;
gamma_d = 0.1;
h_dc = 10;
dprod=3;
del_sig = 10;
hfdp = 3;
#ifdef REAL_LATTICE
h_we = 15;//220;//200
small_cell_thresh= 0.8;//0.3;
tau_sig= 0;//20;
small_bound_thresh = 1.0;
#else
h_we = 600;//220;//200
small_cell_thresh= 0.25;
tau_sig= 10;
small_bound_thresh = 0.2;
#endif
tau_mech=0.05;
k_pressure*=tau_mech;
k_tension*=tau_mech;
}
void DeltaMechanicalForces::operator()(const state_type &y,state_type &dydt,const double t){
fill(dydt.begin(), dydt.end(), 0); //necessary
// cout<<"updating delta"<<endl;
updateDelta(y,dydt);
// cout<<"updating signal"<<endl;
updateSignal(y);
// cout<<"updating pressure"<<endl;
updatePressure(y);
uint nbcells = lattice->getNbCells();
// cout<<"computing forces"<<endl;
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty() && ! ce.isFrozen()){
dTensionAndPressure(&ce,y,dydt);
}
}
double n2= norm2(dydt);
// cout<<t<<" "<<n2<<endl;
addPerturbations(dydt);
checkIntegrity(y,dydt);
//#endif
}
void DeltaMechanicalForces::updateDelta(const state_type &y, state_type &dydt){
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isFrozen()){
uint idx = ce.stateVectorIndex();
dydt[idx] = tau_sig*(getDeltaProduction(m_signal[c])-gamma_d *y[idx]);
}
}
}
void DeltaMechanicalForces::updateSignal(const state_type &y){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
m_signal[c] = 0;
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){ // need to also update frozen cells
double dcis = y[ce.stateVectorIndex()];
double dciseff = deltacisfunc(dcis);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(b->oppositeCell()){ // tissue boundary contributes nothing
double l = length(b,y);
double p = 0;
p=perimeter(b->oppositeCell());
int idx = b->oppositeCell()->stateVectorIndex();
m_signal[c] += y[idx]*l/p;
}
}
m_signal[c]*=dciseff;
}
// cout<<"signal "<<c<<" "<<m_signal[c]<<endl;;
}
}
void DeltaMechanicalForces::initState(state_type &s, double *delta){
uint idx = MechanicalForces::initState(s);
s.resize(state_dim);
int n= lattice->getNbCells();
uint m = cellStateDim();
std::cout<<"filling in..."<<m<<" "<<n<<std::endl;
for (int c=0;c<n;c++){
// std::cout<<c<<std::endl;
uint idx = lattice->getCell(c).stateVectorIndex();
for (uint i=0;i<m;i++){
s[idx++] = delta[c*m+i];
}
}
}
void DeltaMechanicalForces::initState(state_type &s){
MechanicalForces::initState(s,false);
}
void DeltaMechanicalForces::setLattice(Lattice *lat){
lattice=CellInteraction::lattice = MechanicalForces::lattice = lat;
uint n= lat->getNbCells();
m_pressure.resize(n);
m_signal.resize(n);
m_area.resize(n);
m_perimeter.resize(n);
uint m = setStateVectorIndex();
state_dim= m;
// observer.resize(m);
}
uint DeltaMechanicalForces::setStateVectorIndex(){
uint n = MechanicalForces::setStateVectorIndex();
return CellInteraction::setStateVectorIndex(n);
}
float DeltaMechanicalForces::getWeight(Boundary *b){
// return MechanicalForces::getWeight(b);
#ifdef REAL_LATTICE
if(isHC(b->cell())){
return isHC(b->oppositeCell())?100.0:weightfunc(m_signal[b->oppositeCell()->id()]) ;
}
else{
return isSC(b->oppositeCell())?50.0:weightfunc(m_signal[b->cell()->id()]);
}
#else
double s1 = m_signal[b->oppositeCell()->id()];
double s2 = m_signal[b->cell()->id()];
double comb = 0.04*pow(s1-s2,4)+1*pow(s1+s2,2);
return weightfunc(comb);
#endif
}
void DeltaMechanicalForces::fillDeltaVector(double *v,state_type& s){
uint nbcells = lattice->getNbCells();
uint m = cellStateDim();
uint cnt=0;
for (uint c=0;c<nbcells;c++){
uint idx = lattice->getCell(c).stateVectorIndex();
for (uint i=0;i<m;i++){
v[cnt++] = s[idx+i];
}
}
}
// void DeltaMechanicalForces::checkT1Transition(Boundary *b){
// // std::cout<<"check"<<std::endl;
// if(isHC(b->cell()) && isHC(b->oppositeCell())){
// double len = b->scaledLength(Vec2(getScaleX(),getScaleY()));
// if(len < 0.2){
// // cout<<"correcting between "<<b->cell()->id()<<" "<< b->oppositeCell()->id()<<endl;
// b->doT1Transition();
// }
// }
// }
// template void checkT1TransitionWrapper<DeltaMechanicalForces>(Boundary *b, void *arg);
// void myWrapper(Boundary *b, void *args){
// checkT1TransitionWrapper<DeltaMechanicalForces>(b,args);
// }
// void DeltaMechanicalForces::findT1Transitions(){
// // lattice->forAllBoundaries(MechanicalForces::checkT1TransitionWrapper,(void *)this);
// lattice->forAllBoundaries(myWrapper,(void *)this);
// }
// void DeltaMechanicalForces::checkT1TransitionWrapper(Boundary *b, void *arg){
// MechanicalForces *forces = (MechanicalForces *)arg;
// if(!b->isGhost()){
// forces->checkT1Transition(b);
// }
// }
// void MechanicalForces::checkT1Transition(Boundary *b){
// if(isHC(b->cell()) && isHC(b->oppositeCell())){
// double len = b->scaledLength(Vec2(getScaleX(),getScaleY()));
// if(len < 0.2){
// // cout<<"correcting between "<<b->cell()->id()<<" "<< b->oppositeCell()->id()<<endl;
// b->doT1Transition();
// }
// }
// }
// void DeltaMechanicalForces::findT1Transitions(){
// lattice->forAllBoundaries(DeltaMechanicalForces::checkT1TransitionWrapper,(void *)this);
// }
c4science · Help