Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F102271948
runCollierSignaling.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Feb 18, 23:38
Size
5 KB
Mime Type
text/x-c
Expires
Thu, Feb 20, 23:38 (2 d)
Engine
blob
Format
Raw Data
Handle
24321000
Attached To
R9411 tisue modeling
runCollierSignaling.cpp
View Options
#define WITH_OBSERVER_MICHA
#include <signal.h>
#include <iostream>
#include <fstream>
#include "matrix.h"
#include "mex.h"
#include <boost/array.hpp>
#include <boost/range/algorithm.hpp>
#include <boost/numeric/odeint.hpp>
#include "collierDoubleSignaling.hpp"
// includes the cpp files containing template class methods
#include "collierSignaling.cpp"
#include "signalAndMechanicalForces.cpp"
using namespace std;
using namespace boost::numeric::odeint;
void INThandler(int sig)
{
printf("Ctrl-C pressed\n");
exit(0);
}
//SignalAndMechanicalForces <CollierSignaling <HillCoupling> > forces;
SignalAndMechanicalForces <CollierSignaling <HillCauchyCoupling> > forces;
//SignalAndMechanicalForces <CollierDoubleSignaling> forces;
//bool stop(const state_type &x){return forces.findSmallCell(x)!=-1;}
bool stop(const state_type &x){return 0;}
//typedef runge_kutta_dopri5< double , double , double , double , vector_space_algebra > stepper_type;
typedef runge_kutta_cash_karp54< state_type , double > error_stepper_type;
//typedef runge_kutta_fehlberg78 <state_type> error_stepper_type;
//typedef rosenbrock4_controller< state_type , double > error_stepper_type;
//typedef runge_kutta4< state_type> stepper_type;
void mexFunction(int nlhs,mxArray* plhs[],int nrhs, const mxArray *prhs[]){
//bool with_scale = (nrhs>6);
const mxArray *cells = prhs[0];
const mxArray *bonds = prhs[1];
const mxArray *nbonds = prhs[2];
const mxArray *par = prhs[3];
const mxArray *verts =prhs[4];
const mxArray *state = prhs[5];
const mxArray *scale = prhs[6];
const mxArray *frozen_cells = prhs[7];
const mxArray *bc = prhs[8];
double *params = mxGetPr(par);
uint ncells = mxGetN(cells);
int maxbonds = mxGetM(cells);
int totbonds = mxGetN(bonds);
int nbparams = mxGetN(par);
uint nverts = mxGetN(verts);
double *scale_pt = mxGetPr(scale);
uint steps =0;
bool with_delamination= false;
signal(SIGINT, INThandler);
// cout<<"c running..."<<params[0]<<" "<<ncells<<endl;
Lattice lat(mxGetPr(cells), mxGetPr(bonds), ncells, mxGetPr(nbonds),totbonds,maxbonds, mxGetPr(verts), nverts,scale_pt);
// cout<<"conversion done"<<endl;
lat.freeze(mxGetPr(frozen_cells));
// cout<<"freezing done"<<endl;
if(*mxGetPr(bc) > 0.5){
lat.setPeriodicBC();
}
srand(time(NULL));
// MechanicalForces forces;
forces.setLattice(&lat);
//forces.setLattice(&lat,mxGetPr(state));
// cout<<"lattice set"<<endl;
double duration = params[0];
// state_type& y = obs.getState();
state_type y;
//forces.initState(y);
cout<<"initialization..."<<endl;
forces.initState(y,mxGetPr(state));
//forces.noMechanics();
cout<<"initialization done "<<y.size()<<endl;
KeepLast obs;
obs.init(y.size());
// steps = integrate(forces,y,0.0,duration,0.001,obs);
auto stepper = make_controlled<error_stepper_type>(0.001,0.001);
// steps = integrate_adaptive(stepper,forces,y,0.0,duration,0.01,obs);
uint nmstep=10;
double curtime=0.0;
for (uint msteps=0;msteps<nmstep;msteps++){
double mstep_dur = duration/nmstep*(msteps+1);
while(mstep_dur-curtime>0.001){
// cout<<"stop condition"<<endl;
stop(y);
cout<<"integrating..."<<endl;
// if(msteps)
auto iter= boost::find_if(make_adaptive_range(stepper,forces,y,curtime,mstep_dur,0.01,obs),stop);
// forces.fillSignalVector(mxGetPr(plhs[4]));
cout<<"integration done"<<endl;
int smallcell = forces.findSmallCell(y);
forces.fromState(y);
if(smallcell >= 0 && with_delamination){
cout<<"delaminating..." <<smallcell<<endl;
forces.delaminate(smallcell);
// lat.delaminate(smallcell);
cout<<"T1s..."<<endl;
// forces.findT1Transitions();// not very efficient...
cout<<"reinit..."<<endl;
forces.updateState(y);
cout<<"cont"<<endl;
}
// cout<<"initializing..."<<endl;
curtime = obs.getTime();
cout<<"time: "<<curtime<<endl;
}
}
// steps = integrate_adaptive(make_controlled<error_stepper_type>(0.001,0.001),forces,y,0.0,duration,0.01,obs);
// stepper_type stepper;
//steps = integrate_const(stepper,forces,y,0.0,duration,0.0001,obs);
// cout<<"done in "<<iter<< " steps"<<endl;
// y = obs.getState();
// cout<<"preparing output... "<<totbonds<<" "<<lat.getNbBoundaries()<<endl;
uint newmaxbonds= lat.maxNbBoundariesInCell();
forces.updateSignal(y);
// fill up return values
plhs[0] = mxCreateDoubleMatrix(2*nverts,1,mxREAL);// vertices
double *ver0= mxGetPr(plhs[0]);
for(uint i=0;i<2*nverts;i++){ver0[i]= y[i];}
plhs[1] = mxCreateDoubleMatrix(4,totbonds,mxREAL);// boundaries
lat.getBoundaryMatrix(mxGetPr(plhs[1]));
plhs[2] = mxCreateDoubleMatrix(newmaxbonds,lat.getNbCells(),mxREAL);// cells
lat.getCellMatrix(mxGetPr(plhs[2]),newmaxbonds);
plhs[3] = mxCreateDoubleMatrix(2,lat.getNbCells(),mxREAL);//cell centers
lat.getCellCenterMatrix(mxGetPr(plhs[3]));
uint nso = forces.getSignalingOutputNumber(); // signals
for(uint sigo=0;sigo<nso;sigo++){
pair<int,int> si = forces.getSignalingOutputSize(sigo);
plhs[sigo+4] = mxCreateDoubleMatrix(si.first,si.second,mxREAL);
forces.fillSignalingVector(mxGetPr(plhs[4+sigo]),sigo);
}
}
Event Timeline
Log In to Comment