Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82272040
runDeltaMechanicalForces.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, Sep 10, 12:24
Size
4 KB
Mime Type
text/x-c
Expires
Thu, Sep 12, 12:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20674577
Attached To
R9411 tisue modeling
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
];
}
}
Event Timeline
Log In to Comment