Page MenuHomec4science
Files F92037553

runCollierSignaling.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 16, 20:13

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

c4science · Help