Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/polarSignaling.cpp
b7305f7b5a8bmaster

Last Change
Raw File

polarSignaling.cpp
View Options

#include "polarSignaling.hpp"
using namespace std;
PolarSignaling::PolarSignaling(){
tau_pol = 0.01;
}
void PolarSignaling::operator()(const state_type &y, state_type &dydt,double t){
// only allow exchange with neighboring boundaries
// position-dependent pull factor
// cout<<"starting us"<<endl;
updateSignal(y);
// cout<<"uss done"<<endl;
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty() && !ce.isFrozen()){
double sum=0.0;
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
// cout<<"idx 0 "<<b->id()<<endl;
uint idx = b->stateVectorIndex();
uint cidx = b->cell()->id();
if(b->oppositeCell()){
uint oidx = b->oppositeCell()->id();
// cout<<"idx "<<idx<<" "<<oidx<<endl;
// cout<<"idx2 "<<idx<<" "<<oidx<<endl;
Vec2 v = b->positionRelativeToCell();
try{v.Normalize();}
catch(...){v=Vec2(0,0);}
// cout<<"position "<<v<<endl;
// cout<<"vi "<<v<<endl;
// if (isnan(y[idx])) {cout<<"cell "<<c<<" vec "<<v<<endl;exit(0);}
// v *= fabs(y[idx])/v.Norm();
v *= fabs(y[idx]);///v.Norm(); //why not 1/v.Norm()???
double dv = tau_pol*sin(v.Dot(m_cellPolarity[oidx]+m_cellPolarity[cidx]));
// cout<<"dv "<<dv<<" / "<<v.Dot(m_cellPolarity[oidx]) << ":"<< m_cellPolarity[oidx] << " ++ " << v<<endl;
dydt[idx] = dv;
sum += dv;
// cout<<"dv out "<<sum<<endl;
}
}
sum /= ce.getNbBoundaries();
// cout<<"sum "<<ce.id()<<": "<<sum<<endl;
// remove average dydt
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint idx = b->stateVectorIndex();
dydt[idx] -= sum;
}
}
}
}
/*
void PolarSignaling::operator()(const state_type &y, state_type &dydt,double t){
// only allow exchange with neighboring boundaries
// position-dependent pull factor
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isFrozen()){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint idx = b->stateVectorIndex();
uint oidx = b->opposite()->stateVectorIndex();
uint pidx = b->prev()->stateVectorIndex();
uint nidx = b->next()->stateVectorIndex();
uint opidx = b->prev()->opposite()->stateVectorIndex();
uint onidx = b->next()->opposite()->stateVectorIndex();
Vec2 v = b->positionRelativeToCell();
Vec2 pv = b->prev()->positionRelativeToCell();
Vec2 nv = b->next()->positionRelativeToCell();
double pfac = fabs(fabs(v.sinAngle(Vec2(1,0)))-fabs(pv.sinAngle(Vec2(1,0))));
double nfac = fabs(fabs(v.sinAngle(Vec2(1,0)))-fabs(nv.sinAngle(Vec2(1,0))));
dydt[idx] = tau_pol*y[idx]*
(pfac*y[pidx]*(y[oidx]-y[opidx])+nfac*y[nidx]*(y[oidx]-y[onidx]));
}
}
}
}
*/
float PolarSignaling::getWeight(Boundary *b){
return 1.0;
}
void PolarSignaling::initBoundary(Boundary *bo, state_type *state){
// compute position and share across boundaries
uint idx = bo->stateVectorIndex();
uint n= bo->cell()->getNbBoundaries();
double pol = (RND(0.001)+0.9995)/n;
(*state)[idx]=pol;
m_polarity[bo->id()]=pol;
cout<<pol<<endl;
}
void PolarSignaling::updateSignal(const state_type& state){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
double sum =0.0;
m_cellPolarity[ce.id()]=Vec2(0,0);
if(!ce.isEmpty() && ! ce.isFrozen()){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
double vi = state[b->stateVectorIndex()];
m_polarity[b->id()] =vi;
Vec2 v = b->positionRelativeToCell(true);
// cout<<"us "<< vi<< " : "<<v<<endl;
try{
v.Normalize();
}catch(...){v=Vec2(0,0);}
m_cellPolarity[ce.id()] +=v*vi;
sum += vi;
}
if(sum>EPSILON) m_cellPolarity[ce.id()] /=sum;
}
}
// cout<<"us done"<<endl;
}
void PolarSignaling::initState(state_type& state, double *signal){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
m_polarity[b->id()] = state[b->stateVectorIndex()]=signal[b->id()];
}
}
}
}
void PolarSignaling::initState(state_type& state){
CellInteraction::initState(state);
// uint nbcells = lattice->getNbCells();
// for (uint c=0;c<nbcells;c++){
// Cell& ce = lattice->getCell(c);
// if(!ce.isEmpty()){
// BoundaryIterator it(&ce);
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// state[b->stateVectorIndex()]=m_polarity[b->id()];
// }
// }
// }
}
void PolarSignaling::fillVector(double *v, int i)const{
uint cnt=0;
switch(i){
case 0:
for(uint i=0;i<m_polarity.size();i++){
v[i]=m_polarity[i];
}
break;
case 1:
for(uint i=0;i<m_cellPolarity.size();i++){
v[cnt++]=m_cellPolarity[i][0];
v[cnt++]=m_cellPolarity[i][1];
}
break;
default:break;
}
}
c4science · Help