Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/mechanicalForces.cpp
22354a679068master

Last Change
Raw File

mechanicalForces.cpp
View Options

#include "mechanicalForces.hpp"
#include <typeinfo>
// void MechanicalForces::dArea(Cell *c, state_type &y,state_type &dydt){
// double ar = c->getArea();
// double diff = ar-Cell::restArea();
// BoundaryIterator it(c);
// Vec2 from1,to1;
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
// dydt[fidx] += to1[2]*lattice->scaleY()*diff;
// dydt[tidx] -= from1[2]*lattice->scaleY()*diff;
// dydt[fidx+1] -= to1[1]*lattice->scaleX()*diff;
// dydt[tidx+1] += from1[1]*lattice->scaleX()*diff;
// }
// }
//#define PI (boost::math::constants::pi<double>())
using namespace std;
void MechanicalForces::initParameters(){
we_min= 0.1;
we_r = 5.0;//5.0
h_we = 0.5;
hill_co = 4.0;
h_pressure = 2;
HC_pressure = 3;
k_sig=1;
k_pressure=0.001;
k_tension=0.01;
max_pert = 0.1;
small_cell_thresh= 0.2;
small_bound_thresh= 0.2;
small_cell=false;
}
void MechanicalForces::operator()(const state_type &y,state_type &dydt,const double t){
// updateGeometry(y); //done in find_if
// updateSignal(y);
fromState(y);
updatePressure(y);
fill(dydt.begin(), dydt.end(), 0);
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty() && ! ce.isFrozen()){
dTensionAndPressure(&ce,y,dydt);
}
}
// if(norm2(dydt)<0.0001){
double n2= norm2(dydt);
// cout<<t<<" "<<n2<<endl;
// if(n2<0.0001){
addPerturbations(dydt);
//}
// checkIntegrity(y,dydt);
checkConvexity(y,dydt);
//observer.setTime(t); // not very clean but simple
}
double MechanicalForces::norm2(const state_type &s){
uint n= s.size();
double res=0.0;
for(uint i=0;i<n;i++){res+=s[i]*s[i];}
return res;
}
// not so good, each vertex is perturbed twice.
void MechanicalForces::addPerturbations(state_type &s){
// for(uint i=0;i<s.size();i++){
// s[i]+=RND(max_pert)-0.5*max_pert;
// }
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
BoundaryIterator it(&ce);
uint n= ce.getNbBoundaries();
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint tidx = b->to()->vertexStateIndex();
s[tidx] +=RND(max_pert)-0.5*max_pert;
s[tidx+1] +=RND(max_pert)-0.5*max_pert;
}
}
}
}
void MechanicalForces::checkConvexity(const state_type &y,state_type &dydt){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
BoundaryIterator it(&ce);
uint n= ce.getNbBoundaries();
Boundary *b=it.first();
uint fidx = b->from()->vertexStateIndex();
uint tidx = b->to()->vertexStateIndex();
Vec2 vf,vt;
b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],vf,vt);
// b->relativePos(vf,vt);
Vec2 v1 = vt-vf;
Vec2 v2;
while(!it.isLast()){
Boundary *n=it.next();
uint fidx = n->from()->vertexStateIndex();
uint tidx = n->to()->vertexStateIndex();
n->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],vf,vt);
// n->relativePos(vf,vt);
v2 = vt-vf;
if(v2.sinAngle(v1) < EPSILON){ // if two edges are aligned
cout<<"bad "<<c<<" "<<n->oppositeCell()->id()<<endl;
cout<<v1<<" "<<v2<<endl;
Vec2 v3(dydt[fidx],dydt[fidx+1]);
if (v1.sinAngle(v3)<0){ // if vertex moves "inside" cell, project on edges
v1.Normalize();
double c = v1.Dot(v3);
dydt[fidx] = c*v1[0];
dydt[fidx+1] = c*v1[1];
}
}
v1 = v2;
}
}
}
}
void MechanicalForces::checkIntegrity(const state_type &y,state_type &dydt){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
BoundaryIterator it(&ce);
uint n= ce.getNbBoundaries();
Vec2 verts[20];
Vec2 from,to;
uint fidx = it.first()->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// from=Vec2(y[fidx],y[fidx+1]);
verts[0]=Vec2(y[fidx],y[fidx+1]);
Vec2 cent =Vec2(0.0,0.0);// verts[0];
uint i=1;
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint tidx = b->to()->vertexStateIndex();
// to = Vec2(y[tidx],y[tidx+1]);
verts[i] = Vec2(y[tidx],y[tidx+1]);
// b->relativePos(from,to);
b->relativePos(verts[i-1],verts[i]);
cent += verts[i];
i++;
}//verts[n]= verts[0];
cent /=((float) n);
for(i=0;i<n;i++){
verts[i] -= cent;
}
i =0;
for(Boundary *b=it.first();!it.isLast();b=it.next()){
// if(!(verts[i].clockwise(verts[(i+1)%n]))){
int j=(i+1)%n;
Vec2 bd = verts[j]-verts[i];
uint tidx = b->to()->vertexStateIndex();
Vec2 dj = Vec2(dydt[tidx],dydt[tidx+1]);
// if(bd.Dot(dj)<0.0 && bd.Norm2()<0.001){
// try{
// Vec2 perp(bd[1],-bd[0]);
// perp.Normalize();
// double p= perp.Dot(dj);
// dydt[tidx]=p*perp[0];
// dydt[tidx+1]=p*perp[1];
// }catch(char const *e){std::cout<<e<<"while checking integrity "<<std::endl;}
// }
// uint fidx = b->from()->vertexStateIndex();
// Vec2 di = Vec2(dydt[fidx],dydt[fidx+1]);
try{
verts[i].Normalize();
if(verts[i].sinAngle(verts[j],true)>0.0){
// std::cout<<c<<"->"<<i<<":"<<verts[i]<<"-"<<verts[(i+1)%n]<<std::endl;
uint tidx = b->to()->vertexStateIndex();
if(verts[j].sinAngle(Vec2(dydt[tidx],dydt[tidx+1]))>0){//if dydt does not point away.
double p = verts[j][0]*dydt[tidx]+ verts[j][1]*dydt[tidx+1]; // dot product
dydt[tidx]=p*verts[j][0]; //projection on the radial axis
dydt[tidx+1]=p*verts[j][1];
// uint fidx = b->from()->vertexStateIndex();
// p = verts[i][0]*dydt[fidx]+ verts[i][1]*dydt[fidx+1]; // dot product
// dydt[fidx]=p*verts[i][0]; //actually project on the perp
// dydt[fidx+1]=p*verts[i][1];
}
}
}catch(char const *e){std::cout<<e<<"while checking integrity "<<std::endl;}
i++;
}
}
}
}
float MechanicalForces::length(Boundary *b,const state_type &y){
Vec2 from,to,v,s;
uint fidx = b->from()->vertexStateIndex();
uint tidx = b->to()->vertexStateIndex();
b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from, to);
v= to-from;
v *= Vec2(getScaleX(y),getScaleY(y));
return v.Norm();
}
float MechanicalForces::getWeight(Boundary *b){
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()]);
}
}
// returns the index of the smallest cell if below delamination threshold, -1 otherwise
// also updates center position
int MechanicalForces::updateGeometry(const state_type &y){
int small =-1;
uint nbcells = lattice->getNbCells();
double thresh = small_cell_thresh;
// cout<<"updating geometry ..."<<endl;
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
// if(c==80)cout<<"computing area "<<c;
m_perimeter[c] =0;
m_area[c] =0;
BoundaryIterator it(&ce);
Vec2 from,to,cent(0,0);
uint fidx = it.first()->from()->vertexStateIndex();
uint n= ce.getNbBoundaries();
// uint tidx = b->to()->vertexStateIndex();
from=Vec2(y[fidx],y[fidx+1]);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint tidx = b->to()->vertexStateIndex();
to = Vec2(y[tidx],y[tidx+1]);
b->relativePos(from,to);
m_area[c] += from[1]*to[0] - from[0]*to[1];
m_perimeter[c] += length(b,y); //one should remove this one probably
from=to;
cent += to;
// if(c==80) cout<<".";
}
ce.setPerimeter(m_perimeter[c]);//stored twice for convenience
cent/=n;
ce.setCenter(cent);
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from, to);
// }
m_area[c] *= getScaleX(y)*getScaleY(y);
//cout<<"area "<<c<<" "<<m_area[c]<<" < "<<thresh<<endl;
if(m_area[c]<thresh){
small_cell=true;
small=c;
thresh = m_area[c];
// std::cout<<"small cell "<<c<<" "<< m_area[c]<<std::endl;
}
}
}
// cout<<"geometry done"<<endl;
return small;
}
void MechanicalForces::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()){
if(isSC(&ce)){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(isHC(b->oppositeCell())){
double l = length(b,y);
double p = perimeter(b->oppositeCell());
m_signal[c] += k_sig*l/p;
}
}
}
}
}
}
void MechanicalForces::updatePressure(const state_type &y){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()){
if(isSC(&ce)){
double s = pow(m_signal[c],hill_co);
double h = pow(h_pressure,hill_co);
m_pressure[c] = s/(s+h);
}
if(isHC(&ce)){
m_pressure[c] = HC_pressure;
}
m_pressure[c] /= m_area[c];
}
}
}
void MechanicalForces::dTensionAndPressure(Cell *c, const state_type &y,state_type &dydt){
BoundaryIterator it(c);
Vec2 from1,to1;
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(!b->oppositeCell()->isFrozen() && !b->next()->oppositeCell()->isFrozen()){
uint fidx = b->from()->vertexStateIndex();
uint tidx = b->to()->vertexStateIndex();
double we = getWeight(b,y);
b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
// b->relativePosition(from1, to1);
Vec2 v= to1-from1;
double delta_press = pressure(b->cell())- pressure(b->oppositeCell());
Vec2 vperp(v[1],-v[0]);
dydt[tidx] -= k_tension*we*v[0] + k_pressure*delta_press* vperp[0];
dydt[tidx+1] -= k_tension*we*v[1]+ k_pressure*delta_press* vperp[1];
}
}
//TODO update scaling parameters
}
void MechanicalForces::setLattice(Lattice *lat,double *celltypes){
lattice=lat;
uint n= lat->getNbCells();
m_pressure.resize(n);
m_signal.resize(n);
m_area.resize(n);
m_perimeter.resize(n);
if(celltypes){
m_celltype.resize(n);
for (uint i=0;i<n;i++){
m_celltype[i] = celltypes[i]>15?'h':'s';
}
}
uint m =setStateVectorIndex();
// observer.resize(m);
// state.resize(2*m+2);
}
uint MechanicalForces::setStateVectorIndex(){
uint n =lattice->getNbVertices();
for(uint i=0;i<n;i++){
lattice->getVertex(i).setVertexStateIndex(2*i);
}
return 2*n+2; //to include scaling factors
}
int MechanicalForces::initState(state_type &s,bool resize){
uint n =lattice->getNbVertices();
if(resize){
s.resize(2*n+2);
}
for(uint i=0;i<n;i++){
// std::cout<<i<<"/"<<n<<std::endl;
const Vec2& v=lattice->getVertex(i).pos();
s[2*i] =v[0];
s[2*i+1] =v[1];
}
s[2*n]=lattice->getScale()[0];
s[2*n+1]=lattice->getScale()[1];
//observer(s,0);//maybe not the best place to put it
state=&s;
return 2*n+2;
}
bool MechanicalForces::fromState(const state_type &y){
uint n =lattice->getNbVertices();
if(y.size()<2*n+2){return false;}
for(uint i=0;i<n;i++){
Vertex &v= lattice->getVertex(i);
v.setPos(y[v.vertexStateIndex()],y[v.vertexStateIndex()+1]);
}
lattice->setScale(y[2*n],y[2*n+1]);
lattice->updateAllOrientations();
updateGeometry(y);
updateSignal(y);
return true;
}
// static method
// void MechanicalForces::checkT1TransitionWrapper(Boundary *b, void *arg){
// MechanicalForces *forces = (MechanicalForces *)arg;
// if(!b->isGhost()){
// forces->checkT1Transition(b);
// }
// }
//should ulitmately be replaced by a random process
void MechanicalForces::checkT1Transition(Boundary *b){
// std::cout<<"check"<<std::endl;
//boundaries between two HC are unstable
if(isHC(b->cell()) && isHC(b->oppositeCell())){
double len = b->scaledLength(Vec2(getScaleX(),getScaleY()));
if(len < small_bound_thresh){
if(b->cell()->getNbBoundaries()>3){
cout<<"1. correcting between "<<b->cell()->id()<<" "<< b->oppositeCell()->id()<<endl;
cout<<" sig "<<m_signal[b->cell()->id()]<<" and "<<m_signal[b->oppositeCell()->id()]<<endl;
b->doT1Transition();
// cout<<"done"<<endl;
}
}
}
else{
if(isSC(b->cell()) && isSC(b->oppositeCell())){
if(isHC(b->next()->oppositeCell()) && isSC(b->prev()->oppositeCell())){
Vec2 rpos1 = b->prev()->opposite()->tailPositionRelativeToCell();
Vec2 rpos2 = b->next()->opposite()->headPositionRelativeToCell();
float ori1 = fabs(rpos1.sinAngle(Vec2(1,0)));
float ori2 = fabs(rpos2.sinAngle(Vec2(1,0)));
if(ori1 > 0.8 && ori2>0.8){
double len = b->scaledLength(Vec2(getScaleX(),getScaleY()));
if(len < 0.2){
cout<<"2. correcting between "<<b->cell()->id()<<" "<< b->oppositeCell()->id()<<endl;
if(b->cell()->getNbBoundaries()>3){
b->doT1Transition();
// cout<<"done"<<endl;
}
}
}
}
}
}
}
void MechanicalForces::findT1Transitions(){
// uint nbound = lattice->getNbBoundaries();
// for(uint i=0;i<nbound;i++){
// Boundary& b = lattice->getBoundary(i);
// if(!b.isGhost()){
// checkT1Transition(&b);
// cout<<" ok with "<<b.id()<<endl;
// }
// }
vector<Boundary *> cand;
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isEmpty()&& ! ce.isFrozen()){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
cand.push_back(b);
}
while(cand.size()){
checkT1Transition(cand.back());
cand.pop_back();
}
}
}
}
void MechanicalForces::fillSignalVector(double *v){
uint nbcells = lattice->getNbCells();
for (uint c=0;c<nbcells;c++){
v[c] = m_signal[c];
// cout<<"signal "<<c<<" "<<m_signal[c]<<endl;;
}
}
// not tested
// double MechanicalForces::dTension(Cell *c, state_type &y,state_type &dydt){
// BoundaryIterator it(c);
// Vec2 from1,to1;
// double sum=0;
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
// Vec2 diff = to1-from1;
// double nor = diff.weightedNorm(Vec(scaleX(),scaleY()));
// if(nor<EPSILON){
// dydt[fidx] -= diff[0]*scaleX()/nor;
// dydt[tidx] += diff[0]*scaleX()/nor;
// dydt[fidx+1] -= diff[1]*scaleY()/nor;
// dydt[tidx+1] += diff[1]*scaleY()/nor;
// }
// sum+=nor;
// }
// return sum;
// }
// bool MechanicalForces::stopIntegration(const state_type &y)
// {return getInstance()->delaminationTime();}
c4science · Help