Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/runDeltaMechanicalForces.cpp
b7305f7b5a8bmaster

Last Change
Raw File

runDeltaMechanicalForces.cpp
View Options

#define WITH_OBSERVER_MICHA
#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"
using namespace std;
using namespace boost::numeric::odeint;
DeltaMechanicalForces forces;
bool stop(const state_type &x){return 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;
// 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,mxGetPr(state));
// cout<<"initialization done"<<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){
// cout<<"delaminating..."<<endl;
lat.delaminate(smallcell); //problem here
// cout<<"T1s..."<<endl;
forces.findT1Transitions();// not very efficient...
}
// cout<<"initializing..."<<endl;
forces.initState(y); // why not working in if block?
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();
plhs[0] = mxCreateDoubleMatrix(2*nverts,1,mxREAL);// vertices
plhs[1] = mxCreateDoubleMatrix(4,totbonds,mxREAL);// boundaries
plhs[2] = mxCreateDoubleMatrix(newmaxbonds,lat.getNbCells(),mxREAL);// cells
plhs[3] = mxCreateDoubleMatrix(forces.cellStateDim(),lat.getNbCells(),mxREAL);// delta
plhs[4] = mxCreateDoubleMatrix(1,lat.getNbCells(),mxREAL);// signal
lat.getBoundaryMatrix(mxGetPr(plhs[1]));
lat.getCellMatrix(mxGetPr(plhs[2]),newmaxbonds);
forces.fillDeltaVector(mxGetPr(plhs[3]),y);
// forces.updateSignal(y);
forces.fillSignalVector(mxGetPr(plhs[4]));
double *ver0= mxGetPr(plhs[0]);
// double *ver1= mxGetPr(plhs[1]);
for(uint i=0;i<2*nverts;i++){
ver0[i]= y[i];
}
}
c4science · Help