Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/collierSignaling.cpp
b7305f7b5a8bmaster

Last Change
Raw File

collierSignaling.cpp
View Options

#ifndef __COLLIER_SIGNALING_CPP__
#define __COLLIER_SIGNALING_CPP__
#include "collierSignaling.hpp"
#include <assert.h>
using namespace std;
//#define ADIABATIC
template<class Coupling>
CollierSignaling<Coupling>::CollierSignaling(){
v = 1.0;
coupling.init();
// set_hill_params(1.0,0.2,3.0);
}
// notch: y[cell()->stateVectorIndex()]
// delta: y[cell()->stateVectorIndex()+1]
template<class Coupling>
void CollierSignaling<Coupling>::operator()(const state_type &y, state_type &dydt,double t){
updateSignal(y);
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty() && !ce.isFrozen()){
double ndelta = 0.0;
double neighcnt=0.0;
double fac = 1.0;
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(b->oppositeCell() && !b->oppositeCell()->isFrozen()){
fac = b->length()/sqrt(b->oppositeCell()->perimeter());
#ifndef ADIABATIC
ndelta += MAX(0.0,y[b->oppositeCell()->stateVectorIndex()+1])*fac;
#else
ndelta += m_delta[b->oppositeCell()->id()]*fac;
#endif
// neighcnt++;
neighcnt += fac;
}
}
if(neighcnt>0) ndelta /= neighcnt;
uint idx = ce.stateVectorIndex();
dydt[idx] = coupling.f_func(ndelta) - y[idx];
assert(!isnan(ndelta));
assert(!(isnan(coupling.f_func(ndelta)) && cout<<ndelta));
assert(!isnan(dydt[idx]));
#ifndef ADIABATIC
dydt[idx+1] = v*(coupling.g_func(y[idx]) - y[idx+1]);
assert(!isnan(dydt[idx+1]));
#endif
}
}
}
template<class Coupling>
void CollierSignaling<Coupling>::updateSignal(const state_type& y){
#ifdef ADIABATIC
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty() && !ce.isFrozen()){
uint idx = ce.stateVectorIndex();
m_delta[c] = MAX(0,coupling.g_func(y[idx]));
}
}
#endif
}
template<class Coupling>
void CollierSignaling<Coupling>::fromState(const state_type& state){
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
uint idx = ce.stateVectorIndex();
m_notch[c] = state[idx];
#ifndef ADIABATIC
m_delta[c] = state[idx+1];
#endif
}
}
template<class Coupling>
uint CollierSignaling<Coupling>::cellStateDim()const{
#ifdef ADIABATIC
return 1;
#else
return 2;
#endif
}
template<class Coupling>
float CollierSignaling<Coupling>::getWeight(Boundary *b,const state_type &y){
double diff = fabs(y[b->cell()->stateVectorIndex()] -y[b->oppositeCell()->stateVectorIndex()]);
//double diff = MAX(y[b->cell()->stateVectorIndex()],y[b->oppositeCell()->stateVectorIndex()]);
double h = 2;// was 2
double aa = pow(h,3);
double bb = pow(diff,3);
return 0.1 + aa/(aa+bb);
// if(isHC(b->cell())){
// return isHC(b->oppositeCell())?50.0:weightfunc(m_signal[b->oppositeCell()->id()]) ;
// }
// else{
// return isSC(b->oppositeCell())?30.0:weightfunc(m_signal[b->cell()->id()]);
// }
//return 1.0;
}
template<class Coupling>
void CollierSignaling<Coupling>::initState(state_type& state, double *signal){
uint nbcells = lattice->getNbCells();
// m_notch.resize(nbcells);// already initialized
// m_delta.resize(nbcells);
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
m_notch[c] = state[ce.stateVectorIndex()]=signal[2*c+1];
m_delta[c] =signal[2*c];
#ifndef ADIABATIC
state[ce.stateVectorIndex()+1]=signal[2*c];
#endif
// cout<<"init "<<c<< " " <<m_notch[c]<<" "<<m_delta[c]<<" "<<ce.stateVectorIndex()<<endl;
}
}
template<class Coupling>
void CollierSignaling<Coupling>::updateState(state_type& state){
uint nbcells = lattice->getNbCells();
// m_notch.resize(nbcells);// already initialized
// m_delta.resize(nbcells);
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(ce.isEmpty()){
m_notch[c] = m_delta[c] = state[ce.stateVectorIndex()]=-1.0;
#ifndef ADIABATIC
state[ce.stateVectorIndex()+1]=-1.0;
#endif
}
}
}
// first delta and then notch
template<class Coupling>
void CollierSignaling<Coupling>::fillVector(double *v, int i)const{
uint cnt=0;
uint nbcells = lattice->getNbCells();
for(uint i=0;i<nbcells;i++){
v[cnt++] = m_delta[i];
v[cnt++] = m_notch[i];
}
}
template<class Coupling>
double CollierSignaling<Coupling>::delaminationScore(Boundary *bo){
double no1= m_notch[bo->cell()->id()];
double no2 = m_notch[bo->oppositeCell()->id()];
double diff = fabs(no1-no2);
return 1.0/diff;
}
#endif
#ifdef TEST_COLLIER
int main(int argc, char *argv[]){
HillCauchyCoupling sig;
sig.init();
double x= 0.6;
cout<<x<<":"<<sig.f_func(x)<<endl;
return 0;
}
#endif
c4science · Help