Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/runMechanicalForces.cpp
b7305f7b5a8bmaster

Last Change
Raw File

runMechanicalForces.cpp
View Options

#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 "mechanicalForces.hpp"
//#include "deltaMechanicalForces.hpp"
#include "collierSignaling.cpp"
#include "signalAndMechanicalForces.cpp"
using namespace std;
using namespace boost::numeric::odeint;
MechanicalForces forces;
//SignalAndMechanicalForces <CollierSignaling <HillCauchyCoupling> > forces;
bool stop(const state_type &x){return false;//forces.findSmallCell(x)!=-1;
}
//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 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;
//cout<<"c running..."<<params[0]<<" "<<ncells<<endl;
Lattice lat(mxGetPr(cells), mxGetPr(bonds), ncells, mxGetPr(nbonds),totbonds,maxbonds, mxGetPr(verts), nverts,scale_pt);
lat.freeze(mxGetPr(frozen_cells));
if(*mxGetPr(bc) > 0.5){
lat.setPeriodicBC();
}
// cout<<"conversion done"<<endl;
srand(time(NULL));
// MechanicalForces forces;
forces.setLattice(&lat,mxGetPr(state));
double duration = params[0];
// KeepLast& obs =forces.observeLast();
// state_type& y = obs.getState();
state_type y;
forces.initState(y);
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;
for (uint msteps=0;msteps<nmstep;msteps++){
double mstep_dur = duration/nmstep;
// while(mstep_dur>0){
stop(y);
// cout<<"integrating..."<<endl;
// auto iter= boost::find_if(make_adaptive_range(stepper,forces,y,0.0,duration/nmstep,0.01,obs),stop);
// cout<<"before "<<y[2]<<" "<<y[3]<<endl;
integrate_adaptive(stepper, forces, y,0.0,duration/nmstep,0.01,obs);
int smallcell = forces.findSmallCell(y);
// cout<<"after "<<duration<<" " <<y[2]<<" "<<y[3]<<endl;
forces.fromState(y);
if(smallcell >= 0 && with_delamination){
// cout<<"delaminating..."<<endl;
lat.delaminate(smallcell);
// cout<<"T1transitions..."<<endl;
forces.findT1Transitions();// not very efficient...
forces.initState(y); // why not working in if block?
}
// cout<<"time: "<<obs.getTime()<<endl;
// y = obs.getState();
// mstep_dur -= obs.getTime();
// }
}
// 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;
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];}
// cout<<"boundaries.."<<endl;
plhs[1] = mxCreateDoubleMatrix(4,totbonds,mxREAL);// boundaries
lat.getBoundaryMatrix(mxGetPr(plhs[1]));
// cout<<"cells.."<<endl;
uint newmaxbonds= lat.maxNbBoundariesInCell();
plhs[2] = mxCreateDoubleMatrix(newmaxbonds,lat.getNbCells(),mxREAL);// cells
lat.getCellMatrix(mxGetPr(plhs[2]),newmaxbonds);
}
c4science · Help