Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/lattice.cpp
b7305f7b5a8bmaster

Last Change
Raw File

lattice.cpp
View Options

#include <iostream>
#include "lattice.hpp"
#include <boost/math/constants/constants.hpp>
using namespace std;
bool Boundary::s_periodic_bc=false;
Lattice::Lattice(LatticeType type){
uint ncells=0;
switch(type){
case TWO_CELLS: ncells=2;break;
case FIVE_CELLS: ncells=5;break;
case SEVEN_CELLS: ncells=7;break;
default:;
}
m_cells.resize(ncells);
m_bonds.resize(ncells*6);
#ifdef WITH_EXTRACELLULAR_MATRIX
m_extraCells.resize(ncells*6);//max number
m_extracell_cnt=0;
#endif
for(uint i=0;i<ncells;i++){
Cell& c = m_cells[i];
c.setId(i);
double ori=-PI/6;
for(uint b=0;b<6;b++){
Boundary& bo = m_bonds[i*6+b];
bo.setOrientation(ori);
c.insert(bo);
ori+= PI/3;
}
}
switch(type){
case SEVEN_CELLS:
makeAdjacentHex(0,4,5,1);
makeAdjacentHex(1,3,5,0);
makeAdjacentHex(1,4,6,1);
makeAdjacentHex(5,5,6,2);
makeAdjacentHex(2,3,6,0);
case FIVE_CELLS:
makeAdjacentHex(1,5,2,2);
makeAdjacentHex(1,1,3,4);
makeAdjacentHex(1,0,4,3);
makeAdjacentHex(0,0,3,3);
makeAdjacentHex(3,5,4,2);
makeAdjacentHex(2,1,4,4);
case TWO_CELLS:
makeAdjacentHex(0,5,1,2);
default:;
}
}
Lattice::Lattice(double *cells, double* bonds, uint& ncells,
double* nbonds, int& totalbonds,int& maxbonds,
double *verts, uint& nverts, double *scale):
m_cells(ncells), m_bonds(totalbonds),m_verts(nverts)
#ifdef WITH_EXTRACELLULAR_MATRIX
,m_extraCells(totalbonds)//max number
#endif
{
#ifdef WITH_EXTRACELLULAR_MATRIX
m_extracell_cnt=0;
#endif
// filling up vertices list
for(uint i=0;i<nverts;i++){
m_verts[i].setPos(verts[2*i],verts[2*i+1]);
m_verts[i].setId(i);
}
// setting scaling factors
if(scale){
m_scale=Vec2(scale[0],scale[1]);
}else{
m_scale=Vec2(1.0,1.0);
}
for (uint i=0;i<ncells;i++){
Cell& c=m_cells[i];
c.setId(i);
Vec2 center(0,0);
int nv =0;
// cout<<"building cell "<<i<< " with " <<nbonds[i]<<endl;
for (uint j=0;j<nbonds[i];j++){
int bidx = (int)cells[i*maxbonds+j];
if(bidx <0) continue;
// cout<<"building boundary "<<bidx<<endl;
Boundary& b=m_bonds[bidx];//check index -1 is done in matlab
b.setId(bidx);
int vidx = (int)bonds[bidx*4];
Vertex *from = &(m_verts[vidx]);
// from->setPos(verts[vidx*2],verts[vidx*2+1]);
// from->setId(vidx);
vidx = (int)bonds[bidx*4+1]; //-1 is done in matlab
Vertex *to = &(m_verts[vidx]);
// to->setPos(verts[vidx*2],verts[vidx*2+1]);
// to->setId(vidx);
b.setVertices(from,to);
c.insert(b);
// center += from->pos();nv++;
int ocidx = (int) bonds[bidx*4+3];
// cout<<"building boundary "<<bidx<< " with " <<ocidx<<endl;
if(ocidx<0) {
#ifdef WITH_EXTRACELLULAR_MATRIX
ExtraCell *ec = addExtracellularMatrix();
b.setExtracellularMatrix(ec);
#endif
continue;
}
Cell& oc = m_cells[ocidx];
// cout<<"assignments"<<endl;
b.setOppositeCell(&oc);
if(oc.initialized()){
// cout<<"mutual bond "<<i<<" "<<idx<<endl;
b.setOppositeBoundary(&oc);
#ifdef WITH_EXTRACELLULAR_MATRIX
ExtraCell *ec = addExtracellularMatrix();
b.setExtracellularMatrix(ec);
b.opposite()->setExtracellularMatrix(ec);
#endif
}
}
// center /= nv;
// c.setCenter(center);
}
}
// can't use this for changing the lattice...
uint Lattice::forAllCellsAndBoundaries(cellVisitorFunction cfun, void *cellargs, boundaryVisitorFunction bfun, void *boundargs ){
uint cnt=0;
for (uint c=0;c<getNbCells();c++){
Cell& ce = getCell(c);
// if(!ce.isEmpty()){ // ignore dead cells
if(cfun!=NULL)cfun(&ce,cellargs);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(bfun!=NULL) bfun(b,boundargs);
cnt++;
}
// }
}
return cnt;
}
void Lattice::getBoundaryMatrix(double *bonds){
uint idx=0;
cout<<"# of boundaries "<<m_bonds.size()<<endl;
for(uint i=0;i<m_bonds.size();i++){
Boundary& b=m_bonds[i];
if(!b.isGhost()){
bonds[idx++] = b.from()->id();
bonds[idx++] = b.to()->id();
bonds[idx++] = b.cell()->id();
bonds[idx++] = b.oppositeCell()?b.oppositeCell()->id():-1;
}
else{//cout<<"unspecified boundary "<<i<<endl;
bonds[idx++] = -1;
bonds[idx++] = -1;
bonds[idx++] = -1;
bonds[idx++] = -1;
}
}
}
void Lattice::getCellMatrix(double *cells, uint maxbonds){
uint idx=0;
for(uint i=0;i<m_cells.size();i++){
Cell& c=m_cells[i];
uint nb=c.getNbBoundaries();
if(nb>maxbonds)
{cout<<"cell has too many boundaries"<<endl;}
else{
BoundaryIterator it(&c);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
cells[idx++] = b->id();
}
for(uint j=nb;j<maxbonds;j++){
cells[idx++] = -2;
}
}
}
}
void Lattice::getCellCenterMatrix(double *cent){
uint idx=0;
for(uint i=0;i<m_cells.size();i++){
Cell& c=m_cells[i];
if(c.isEmpty()){
cent[idx++]= cent[idx++]=0.0;
}
else{
cent[idx++]=c.center()[0];
cent[idx++]=c.center()[1];
}
}
}
void Lattice::makeAdjacentHex(uint cell1, uint bound1, uint cell2, uint bound2){
uint b1 = cell1*6+bound1;
uint b2 = cell2*6+bound2;
m_bonds[b1].setOppositeCell(&(m_cells[cell2]));
m_bonds[b2].setOppositeCell(&(m_cells[cell1]));
m_bonds[b1].setOppositeBoundary(&(m_bonds[b2]));
m_bonds[b2].setOppositeBoundary(&(m_bonds[b1]));
#ifdef WITH_EXTRACELLULAR_MATRIX
ExtraCell *ec = addExtracellularMatrix();
m_bonds[b1].setExtracellularMatrix(ec);
m_bonds[b2].setExtracellularMatrix(ec);
#endif
}
void Boundary::relativePosition(Vec2& from,Vec2& to){
const Vec2& f=m_from->pos();
const Vec2& t=m_to->pos();
from = f;to=t;
relativePos(from,to);
// for (uint i=0;i<2;i++){
// if(abs(t[i]-f[i]-2*PI)<abs(t[i]-f[i])){to[i]-= 2*PI;}
// if(abs(t[i]-f[i]+2*PI)<abs(t[i]-f[i])){to[i]+= 2*PI;}
// }
}
void Boundary::relativePosition(double fromx, double fromy, double tox, double toy,Vec2& from,Vec2& to){
Vec2 f(fromx,fromy);
Vec2 t(tox,toy);
from = f;to=t;
relativePos(from,to);
// for (uint i=0;i<2;i++){
// if(abs(t[i]-f[i]-2*PI)<abs(t[i]-f[i])){to[i]-= 2*PI;}
// if(abs(t[i]-f[i]+2*PI)<abs(t[i]-f[i])){to[i]+= 2*PI;}
// }
}
void Boundary::relativePos(Vec2& from,Vec2& to){
if(s_periodic_bc){
for (uint i=0;i<2;i++){
if(abs(to[i]-from[i]-2*PI)<abs(to[i]-from[i])){to[i]-= 2*PI;}
if(abs(to[i]-from[i]+2*PI)<abs(to[i]-from[i])){to[i]+= 2*PI;}
}
}
}
Vec2 Boundary::positionRelativeToCell(bool update){
if(update){
cell()->updateCenter();
}
Vec2 f,t;
Vec2 cent = cell()->center();
relativePosition(f,t);
Vec2 mid = f+t;
mid *= 0.5;
relativePos(cent,mid);
return mid-cent;
}
Vec2 Boundary::headPositionRelativeToCell(bool update){
Vec2 center = cell()->center();
Vec2 t=to()->pos();
relativePosition(t,center);
return t-center;
}
Vec2 Boundary::tailPositionRelativeToCell(bool update){
Vec2 center = cell()->center();
Vec2 f=from()->pos();
relativePosition(f,center);
return f-center;
}
void Boundary::updateOrientation(){
if(m_from!=NULL && m_to!=NULL){
Vec2 f,t;
relativePosition(f,t);
setOrientation((double)(t-f).direction());
// cout<<f<<" "<<t<<":"<<m_orient*180/PI<<endl;
}
}
void Boundary::doT1Transition(){
Boundary *prec,*nex,*op,*oprec,*onex;
Cell *c,*cop,*cprec,*cnex;
// cout<<"removing boundary between cells "<<cell()->id()<<" and "<<oppositeCell()->id()<<endl;
prec = prev();
nex = next();
// cout<<"op"<<endl;
op = opposite();
if(op==NULL){ cout<<"no opposite boundary..."<<endl;return;}
oprec =op->prev();
// cout<<"next"<<endl;
onex = op->next();
c= cell();
// cout<<"neig1"<<endl;
cop = neighbor();
// cout<<"neigh2"<<endl;
cprec = prec->neighbor();
// cout<<"neigh3"<<endl;
cnex =nex->neighbor();
// cout<<"removing"<<endl;
c->remove(this);
cop->remove(op);
// cout<<"inserting"<<endl;
insertBefore(oprec->opposite());
// cout<<"after"<<endl;
op->insertAfter(onex->opposite());
// cout<<"opposing"<<endl;
// op->setCell(cnex);
op->setOppositeCell(cnex);
// setCell(cprec);
setOppositeCell(cprec);
// cout<<"updating orientation"<<endl;
// // update orientation of edge
// double twopi = 2*PI;
// double ori = orientation()+PI;
// ori -= ori>PI?twopi:0;
// ori += ori<-PI?twopi:0; //probably not needed
// setOrientation(ori);
// ori = op->orientation()+PI;
// ori -= ori>PI?twopi:0;
// ori += ori<-PI?twopi:0; // probably not neede
// op->setOrientation(ori);
// cout<<"updating position"<<endl;
// update Vertex position
// const Vec2& f = from()->pos();
// const Vec2& t = to()->pos();
Vec2 f;
Vec2 t;
relativePosition(f,t);
// Vec2 mid =(Vec2) (f+t)*0.5;
Vec2 mid = t+f;
mid *= 0.5;
Vec2 dir= t-f;
Vec2 perp = Vec2(dir[1],-dir[0]);
perp.Normalize();
perp*= 0.01;
*from() = mid+perp;
*to() = mid-perp;
// update orientation of edge: done just when needed
// updateOrientation();
// op->updateOrientation();
// prec->updateOrientation();
// prev->updateOrientation();
// cout<<"done"<<endl;
}
// void Lattice::initHexLatticeVertices(uint nrow, uint ncol, bool vertical=true){
// double sq3 = sqrt(3);
// uint nvert = 2*nrow*(ncol+1)+2*ncol;
// m_verts.resize(nvert);
// uint idx = 0;
// for (uint i=0;i<=ncol;i++){
// for(uint j=0;j<nrow;j++){
// m_verts[idx++] = Vec2(i*sq3-0.5,sq3*(0.5*(j%2)-0.5));
// m_verts[idx++] = Vec2(i*sq3-1,sq3*(0.5*(j%2)));
// }
// // if(i%2){
// // m_verts[
// // }
// }
// }
// void Lattice::setPositions(uint cell){
// m_cells[cell].setPosition(Vec2(0,0));
// m_cells[cell].
// set<uint> ok;
// while(ok.size()<getNbCells()){
// }
// }
bool Boundary::setOppositeBoundary(Cell *oc){
BoundaryIterator it(oc);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
if(b->oppositeCell()==m_cell){
m_opposite = b;
b->setOppositeBoundary(this);
return true;
}
}
return false;
}
void Boundary::remove(){
prev()->setNext(next());
next()->setPrevious(prev());
next()->setFrom(prev()->to());
if(next()->opposite()==NULL){
cout<<"boundary "<<next()->id()<<" has no opposite"<<endl;
}
else{
next()->opposite()->setTo(prev()->to());
}
m_prev=m_next=NULL;
m_cell=m_neighbor=NULL;
toGhost();
}
void Boundary::updateBoundaryOrientation(Boundary *b,void *arg){
b->updateOrientation();
}
void Cell::updateMaxBoundaries(Cell *c,void *arg){
uint *n = (uint *) arg;
uint m = c->getNbBoundaries();
if(*n<m){*n=m;}
}
void Cell::remove(Boundary *b){
if(b->cell()==this){
if(m_empty=!(--m_nboundaries)){
m_first=NULL;
}
else{
if(m_first==b){
m_first=b->next();
}
// done in Boundary::remove
// b->prev()->setTo(b->next->from())
// b->next()->setFrom(b->prev()->to());
b->remove();
}
}
}
void Cell::doT2Transition(){
vector<Boundary *> mem;
while(getNbBoundaries()>3){
double min_ori=2;
Boundary *min_b=NULL;
BoundaryIterator it(this);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
// b->updateOrientation();
// double ori = abs(cos(b->orientation()));
Vec2 rpos = b->opposite()->positionRelativeToCell(true);
double ori = abs(rpos.sinAngle(Vec2(1,0)));//TODO find something smarter
if(ori<min_ori){
min_ori=ori;
min_b = b;
}
// cout<< "checking "<<b->oppositeCell()->id()<<" "<<ori<<endl;
}
if(min_b){
// cout<<"trans "<<min_b->oppositeCell()->id()<<endl;
min_b->doT1Transition();
mem.push_back(min_b);
}
}
cout<<"removing cell "<<id()<<endl;
BoundaryIterator it(this);
Vertex *v = it.first()->from();
// Vec2 mp(0,0);
// BoundaryIterator it(this);
Vec2 f = it.first()->from()->pos();
Cell *tc=it.first()->next()->oppositeCell();
bool merge=false;
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// if(b->oppositeCell()->getNbBoundaries()==3){
// cout<<"bad lattice"<<endl;
// return;
// }
// }
updateCenter();
*v=center();
// v->setPos(mp[0],mp[1]);
Boundary *b = firstBoundary();
Boundary *nb;
while(!isEmpty()){
cout<<"ok "<<b->id()<<endl;
b->opposite()->next()->setFrom(v);
b->opposite()->prev()->setTo(v);
nb=b->next(); // because also empties next slot
b->oppositeCell()->remove(b->opposite());
remove(b);//do something with ghost boundaries??
b=nb;
}
cout<<"done"<<endl;
}
void Cell::updateCenter(){
m_center = Vec2(0,0);
Vec2 to;
BoundaryIterator it(this);
Vec2 from = it.first()->from()->pos();
for(Boundary *b=it.first();!it.isLast();b=it.next()){
to = b->to()->pos();
b->relativePos(from,to);
m_center += to;
from = to;
}
m_center /= getNbBoundaries();
}
void Cell::updatePerimeter(){
m_perimeter = 0.0;
BoundaryIterator it(this);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
m_perimeter += b->length();
}
}
void Cell::printPosition(){
BoundaryIterator it(this);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
cout<<b->from()->pos()<<endl;
}
}
c4science · Help