Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F106925801
polarSignaling.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
Wed, Apr 2, 14:49
Size
4 KB
Mime Type
text/x-c
Expires
Fri, Apr 4, 14:49 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
25307431
Attached To
R9411 tisue modeling
polarSignaling.cpp
View Options
#include "polarSignaling.hpp"
using
namespace
std
;
PolarSignaling
::
PolarSignaling
(){
tau_pol
=
0.01
;
}
void
PolarSignaling
::
operator
()(
const
state_type
&
y
,
state_type
&
dydt
,
double
t
){
// only allow exchange with neighboring boundaries
// position-dependent pull factor
// cout<<"starting us"<<endl;
updateSignal
(
y
);
// cout<<"uss done"<<endl;
int
n
=
lattice
->
getNbCells
();
for
(
int
c
=
0
;
c
<
n
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()
&&
!
ce
.
isFrozen
()){
double
sum
=
0.0
;
BoundaryIterator
it
(
&
ce
);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
// cout<<"idx 0 "<<b->id()<<endl;
uint
idx
=
b
->
stateVectorIndex
();
uint
cidx
=
b
->
cell
()
->
id
();
if
(
b
->
oppositeCell
()){
uint
oidx
=
b
->
oppositeCell
()
->
id
();
// cout<<"idx "<<idx<<" "<<oidx<<endl;
// cout<<"idx2 "<<idx<<" "<<oidx<<endl;
Vec2
v
=
b
->
positionRelativeToCell
();
try
{
v
.
Normalize
();}
catch
(...){
v
=
Vec2
(
0
,
0
);}
// cout<<"position "<<v<<endl;
// cout<<"vi "<<v<<endl;
// if (isnan(y[idx])) {cout<<"cell "<<c<<" vec "<<v<<endl;exit(0);}
// v *= fabs(y[idx])/v.Norm();
v
*=
fabs
(
y
[
idx
]);
///v.Norm(); //why not 1/v.Norm()???
double
dv
=
tau_pol
*
sin
(
v
.
Dot
(
m_cellPolarity
[
oidx
]
+
m_cellPolarity
[
cidx
]));
// cout<<"dv "<<dv<<" / "<<v.Dot(m_cellPolarity[oidx]) << ":"<< m_cellPolarity[oidx] << " ++ " << v<<endl;
dydt
[
idx
]
=
dv
;
sum
+=
dv
;
// cout<<"dv out "<<sum<<endl;
}
}
sum
/=
ce
.
getNbBoundaries
();
// cout<<"sum "<<ce.id()<<": "<<sum<<endl;
// remove average dydt
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
uint
idx
=
b
->
stateVectorIndex
();
dydt
[
idx
]
-=
sum
;
}
}
}
}
/*
void PolarSignaling::operator()(const state_type &y, state_type &dydt,double t){
// only allow exchange with neighboring boundaries
// position-dependent pull factor
int n= lattice->getNbCells();
for (int c=0;c<n;c++){
Cell& ce = lattice->getCell(c);
if(!ce.isFrozen()){
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint idx = b->stateVectorIndex();
uint oidx = b->opposite()->stateVectorIndex();
uint pidx = b->prev()->stateVectorIndex();
uint nidx = b->next()->stateVectorIndex();
uint opidx = b->prev()->opposite()->stateVectorIndex();
uint onidx = b->next()->opposite()->stateVectorIndex();
Vec2 v = b->positionRelativeToCell();
Vec2 pv = b->prev()->positionRelativeToCell();
Vec2 nv = b->next()->positionRelativeToCell();
double pfac = fabs(fabs(v.sinAngle(Vec2(1,0)))-fabs(pv.sinAngle(Vec2(1,0))));
double nfac = fabs(fabs(v.sinAngle(Vec2(1,0)))-fabs(nv.sinAngle(Vec2(1,0))));
dydt[idx] = tau_pol*y[idx]*
(pfac*y[pidx]*(y[oidx]-y[opidx])+nfac*y[nidx]*(y[oidx]-y[onidx]));
}
}
}
}
*/
float
PolarSignaling
::
getWeight
(
Boundary
*
b
){
return
1.0
;
}
void
PolarSignaling
::
initBoundary
(
Boundary
*
bo
,
state_type
*
state
){
// compute position and share across boundaries
uint
idx
=
bo
->
stateVectorIndex
();
uint
n
=
bo
->
cell
()
->
getNbBoundaries
();
double
pol
=
(
RND
(
0.001
)
+
0.9995
)
/
n
;
(
*
state
)[
idx
]
=
pol
;
m_polarity
[
bo
->
id
()]
=
pol
;
cout
<<
pol
<<
endl
;
}
void
PolarSignaling
::
updateSignal
(
const
state_type
&
state
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
double
sum
=
0.0
;
m_cellPolarity
[
ce
.
id
()]
=
Vec2
(
0
,
0
);
if
(
!
ce
.
isEmpty
()
&&
!
ce
.
isFrozen
()){
BoundaryIterator
it
(
&
ce
);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
double
vi
=
state
[
b
->
stateVectorIndex
()];
m_polarity
[
b
->
id
()]
=
vi
;
Vec2
v
=
b
->
positionRelativeToCell
(
true
);
// cout<<"us "<< vi<< " : "<<v<<endl;
try
{
v
.
Normalize
();
}
catch
(...){
v
=
Vec2
(
0
,
0
);}
m_cellPolarity
[
ce
.
id
()]
+=
v
*
vi
;
sum
+=
vi
;
}
if
(
sum
>
EPSILON
)
m_cellPolarity
[
ce
.
id
()]
/=
sum
;
}
}
// cout<<"us done"<<endl;
}
void
PolarSignaling
::
initState
(
state_type
&
state
,
double
*
signal
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
BoundaryIterator
it
(
&
ce
);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
m_polarity
[
b
->
id
()]
=
state
[
b
->
stateVectorIndex
()]
=
signal
[
b
->
id
()];
}
}
}
}
void
PolarSignaling
::
initState
(
state_type
&
state
){
CellInteraction
::
initState
(
state
);
// uint nbcells = lattice->getNbCells();
// for (uint c=0;c<nbcells;c++){
// Cell& ce = lattice->getCell(c);
// if(!ce.isEmpty()){
// BoundaryIterator it(&ce);
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// state[b->stateVectorIndex()]=m_polarity[b->id()];
// }
// }
// }
}
void
PolarSignaling
::
fillVector
(
double
*
v
,
int
i
)
const
{
uint
cnt
=
0
;
switch
(
i
){
case
0
:
for
(
uint
i
=
0
;
i
<
m_polarity
.
size
();
i
++
){
v
[
i
]
=
m_polarity
[
i
];
}
break
;
case
1
:
for
(
uint
i
=
0
;
i
<
m_cellPolarity
.
size
();
i
++
){
v
[
cnt
++
]
=
m_cellPolarity
[
i
][
0
];
v
[
cnt
++
]
=
m_cellPolarity
[
i
][
1
];
}
break
;
default
:
break
;
}
}
Event Timeline
Log In to Comment