Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77744266
mechanicalForces.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
Fri, Aug 16, 04:30
Size
14 KB
Mime Type
text/x-c
Expires
Sun, Aug 18, 04:30 (2 d)
Engine
blob
Format
Raw Data
Handle
19913438
Attached To
R9411 tisue modeling
mechanicalForces.cpp
View Options
#include "mechanicalForces.hpp"
#include <typeinfo>
// void MechanicalForces::dArea(Cell *c, state_type &y,state_type &dydt){
// double ar = c->getArea();
// double diff = ar-Cell::restArea();
// BoundaryIterator it(c);
// Vec2 from1,to1;
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
// dydt[fidx] += to1[2]*lattice->scaleY()*diff;
// dydt[tidx] -= from1[2]*lattice->scaleY()*diff;
// dydt[fidx+1] -= to1[1]*lattice->scaleX()*diff;
// dydt[tidx+1] += from1[1]*lattice->scaleX()*diff;
// }
// }
//#define PI (boost::math::constants::pi<double>())
using
namespace
std
;
void
MechanicalForces
::
initParameters
(){
we_min
=
0.1
;
we_r
=
5.0
;
//5.0
h_we
=
0.5
;
hill_co
=
4.0
;
h_pressure
=
2
;
HC_pressure
=
3
;
k_sig
=
1
;
k_pressure
=
0.001
;
k_tension
=
0.01
;
max_pert
=
0.1
;
small_cell_thresh
=
0.2
;
small_bound_thresh
=
0.2
;
small_cell
=
false
;
}
void
MechanicalForces
::
operator
()(
const
state_type
&
y
,
state_type
&
dydt
,
const
double
t
){
// updateGeometry(y); //done in find_if
// updateSignal(y);
fromState
(
y
);
updatePressure
(
y
);
fill
(
dydt
.
begin
(),
dydt
.
end
(),
0
);
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()
&&
!
ce
.
isFrozen
()){
dTensionAndPressure
(
&
ce
,
y
,
dydt
);
}
}
// if(norm2(dydt)<0.0001){
double
n2
=
norm2
(
dydt
);
// cout<<t<<" "<<n2<<endl;
// if(n2<0.0001){
addPerturbations
(
dydt
);
//}
// checkIntegrity(y,dydt);
checkConvexity
(
y
,
dydt
);
//observer.setTime(t); // not very clean but simple
}
double
MechanicalForces
::
norm2
(
const
state_type
&
s
){
uint
n
=
s
.
size
();
double
res
=
0.0
;
for
(
uint
i
=
0
;
i
<
n
;
i
++
){
res
+=
s
[
i
]
*
s
[
i
];}
return
res
;
}
// not so good, each vertex is perturbed twice.
void
MechanicalForces
::
addPerturbations
(
state_type
&
s
){
// for(uint i=0;i<s.size();i++){
// s[i]+=RND(max_pert)-0.5*max_pert;
// }
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
BoundaryIterator
it
(
&
ce
);
uint
n
=
ce
.
getNbBoundaries
();
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
s
[
tidx
]
+=
RND
(
max_pert
)
-
0.5
*
max_pert
;
s
[
tidx
+
1
]
+=
RND
(
max_pert
)
-
0.5
*
max_pert
;
}
}
}
}
void
MechanicalForces
::
checkConvexity
(
const
state_type
&
y
,
state_type
&
dydt
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
BoundaryIterator
it
(
&
ce
);
uint
n
=
ce
.
getNbBoundaries
();
Boundary
*
b
=
it
.
first
();
uint
fidx
=
b
->
from
()
->
vertexStateIndex
();
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
Vec2
vf
,
vt
;
b
->
relativePosition
(
y
[
fidx
],
y
[
fidx
+
1
],
y
[
tidx
],
y
[
tidx
+
1
],
vf
,
vt
);
// b->relativePos(vf,vt);
Vec2
v1
=
vt
-
vf
;
Vec2
v2
;
while
(
!
it
.
isLast
()){
Boundary
*
n
=
it
.
next
();
uint
fidx
=
n
->
from
()
->
vertexStateIndex
();
uint
tidx
=
n
->
to
()
->
vertexStateIndex
();
n
->
relativePosition
(
y
[
fidx
],
y
[
fidx
+
1
],
y
[
tidx
],
y
[
tidx
+
1
],
vf
,
vt
);
// n->relativePos(vf,vt);
v2
=
vt
-
vf
;
if
(
v2
.
sinAngle
(
v1
)
<
EPSILON
){
// if two edges are aligned
cout
<<
"bad "
<<
c
<<
" "
<<
n
->
oppositeCell
()
->
id
()
<<
endl
;
cout
<<
v1
<<
" "
<<
v2
<<
endl
;
Vec2
v3
(
dydt
[
fidx
],
dydt
[
fidx
+
1
]);
if
(
v1
.
sinAngle
(
v3
)
<
0
){
// if vertex moves "inside" cell, project on edges
v1
.
Normalize
();
double
c
=
v1
.
Dot
(
v3
);
dydt
[
fidx
]
=
c
*
v1
[
0
];
dydt
[
fidx
+
1
]
=
c
*
v1
[
1
];
}
}
v1
=
v2
;
}
}
}
}
void
MechanicalForces
::
checkIntegrity
(
const
state_type
&
y
,
state_type
&
dydt
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
BoundaryIterator
it
(
&
ce
);
uint
n
=
ce
.
getNbBoundaries
();
Vec2
verts
[
20
];
Vec2
from
,
to
;
uint
fidx
=
it
.
first
()
->
from
()
->
vertexStateIndex
();
// uint tidx = b->to()->vertexStateIndex();
// from=Vec2(y[fidx],y[fidx+1]);
verts
[
0
]
=
Vec2
(
y
[
fidx
],
y
[
fidx
+
1
]);
Vec2
cent
=
Vec2
(
0.0
,
0.0
);
// verts[0];
uint
i
=
1
;
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
// to = Vec2(y[tidx],y[tidx+1]);
verts
[
i
]
=
Vec2
(
y
[
tidx
],
y
[
tidx
+
1
]);
// b->relativePos(from,to);
b
->
relativePos
(
verts
[
i
-
1
],
verts
[
i
]);
cent
+=
verts
[
i
];
i
++
;
}
//verts[n]= verts[0];
cent
/=
((
float
)
n
);
for
(
i
=
0
;
i
<
n
;
i
++
){
verts
[
i
]
-=
cent
;
}
i
=
0
;
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
// if(!(verts[i].clockwise(verts[(i+1)%n]))){
int
j
=
(
i
+
1
)
%
n
;
Vec2
bd
=
verts
[
j
]
-
verts
[
i
];
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
Vec2
dj
=
Vec2
(
dydt
[
tidx
],
dydt
[
tidx
+
1
]);
// if(bd.Dot(dj)<0.0 && bd.Norm2()<0.001){
// try{
// Vec2 perp(bd[1],-bd[0]);
// perp.Normalize();
// double p= perp.Dot(dj);
// dydt[tidx]=p*perp[0];
// dydt[tidx+1]=p*perp[1];
// }catch(char const *e){std::cout<<e<<"while checking integrity "<<std::endl;}
// }
// uint fidx = b->from()->vertexStateIndex();
// Vec2 di = Vec2(dydt[fidx],dydt[fidx+1]);
try
{
verts
[
i
].
Normalize
();
if
(
verts
[
i
].
sinAngle
(
verts
[
j
],
true
)
>
0.0
){
// std::cout<<c<<"->"<<i<<":"<<verts[i]<<"-"<<verts[(i+1)%n]<<std::endl;
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
if
(
verts
[
j
].
sinAngle
(
Vec2
(
dydt
[
tidx
],
dydt
[
tidx
+
1
]))
>
0
){
//if dydt does not point away.
double
p
=
verts
[
j
][
0
]
*
dydt
[
tidx
]
+
verts
[
j
][
1
]
*
dydt
[
tidx
+
1
];
// dot product
dydt
[
tidx
]
=
p
*
verts
[
j
][
0
];
//projection on the radial axis
dydt
[
tidx
+
1
]
=
p
*
verts
[
j
][
1
];
// uint fidx = b->from()->vertexStateIndex();
// p = verts[i][0]*dydt[fidx]+ verts[i][1]*dydt[fidx+1]; // dot product
// dydt[fidx]=p*verts[i][0]; //actually project on the perp
// dydt[fidx+1]=p*verts[i][1];
}
}
}
catch
(
char
const
*
e
){
std
::
cout
<<
e
<<
"while checking integrity "
<<
std
::
endl
;}
i
++
;
}
}
}
}
float
MechanicalForces
::
length
(
Boundary
*
b
,
const
state_type
&
y
){
Vec2
from
,
to
,
v
,
s
;
uint
fidx
=
b
->
from
()
->
vertexStateIndex
();
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
b
->
relativePosition
(
y
[
fidx
],
y
[
fidx
+
1
],
y
[
tidx
],
y
[
tidx
+
1
],
from
,
to
);
v
=
to
-
from
;
v
*=
Vec2
(
getScaleX
(
y
),
getScaleY
(
y
));
return
v
.
Norm
();
}
float
MechanicalForces
::
getWeight
(
Boundary
*
b
){
if
(
isHC
(
b
->
cell
())){
return
isHC
(
b
->
oppositeCell
())
?
50.0
:
weightfunc
(
m_signal
[
b
->
oppositeCell
()
->
id
()])
;
}
else
{
return
isSC
(
b
->
oppositeCell
())
?
30.0
:
weightfunc
(
m_signal
[
b
->
cell
()
->
id
()]);
}
}
// returns the index of the smallest cell if below delamination threshold, -1 otherwise
// also updates center position
int
MechanicalForces
::
updateGeometry
(
const
state_type
&
y
){
int
small
=-
1
;
uint
nbcells
=
lattice
->
getNbCells
();
double
thresh
=
small_cell_thresh
;
// cout<<"updating geometry ..."<<endl;
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
// if(c==80)cout<<"computing area "<<c;
m_perimeter
[
c
]
=
0
;
m_area
[
c
]
=
0
;
BoundaryIterator
it
(
&
ce
);
Vec2
from
,
to
,
cent
(
0
,
0
);
uint
fidx
=
it
.
first
()
->
from
()
->
vertexStateIndex
();
uint
n
=
ce
.
getNbBoundaries
();
// uint tidx = b->to()->vertexStateIndex();
from
=
Vec2
(
y
[
fidx
],
y
[
fidx
+
1
]);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
to
=
Vec2
(
y
[
tidx
],
y
[
tidx
+
1
]);
b
->
relativePos
(
from
,
to
);
m_area
[
c
]
+=
from
[
1
]
*
to
[
0
]
-
from
[
0
]
*
to
[
1
];
m_perimeter
[
c
]
+=
length
(
b
,
y
);
//one should remove this one probably
from
=
to
;
cent
+=
to
;
// if(c==80) cout<<".";
}
ce
.
setPerimeter
(
m_perimeter
[
c
]);
//stored twice for convenience
cent
/=
n
;
ce
.
setCenter
(
cent
);
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from, to);
// }
m_area
[
c
]
*=
getScaleX
(
y
)
*
getScaleY
(
y
);
//cout<<"area "<<c<<" "<<m_area[c]<<" < "<<thresh<<endl;
if
(
m_area
[
c
]
<
thresh
){
small_cell
=
true
;
small
=
c
;
thresh
=
m_area
[
c
];
// std::cout<<"small cell "<<c<<" "<< m_area[c]<<std::endl;
}
}
}
// cout<<"geometry done"<<endl;
return
small
;
}
void
MechanicalForces
::
updateSignal
(
const
state_type
&
y
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
m_signal
[
c
]
=
0
;
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
if
(
isSC
(
&
ce
)){
BoundaryIterator
it
(
&
ce
);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
if
(
isHC
(
b
->
oppositeCell
())){
double
l
=
length
(
b
,
y
);
double
p
=
perimeter
(
b
->
oppositeCell
());
m_signal
[
c
]
+=
k_sig
*
l
/
p
;
}
}
}
}
}
}
void
MechanicalForces
::
updatePressure
(
const
state_type
&
y
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()){
if
(
isSC
(
&
ce
)){
double
s
=
pow
(
m_signal
[
c
],
hill_co
);
double
h
=
pow
(
h_pressure
,
hill_co
);
m_pressure
[
c
]
=
s
/
(
s
+
h
);
}
if
(
isHC
(
&
ce
)){
m_pressure
[
c
]
=
HC_pressure
;
}
m_pressure
[
c
]
/=
m_area
[
c
];
}
}
}
void
MechanicalForces
::
dTensionAndPressure
(
Cell
*
c
,
const
state_type
&
y
,
state_type
&
dydt
){
BoundaryIterator
it
(
c
);
Vec2
from1
,
to1
;
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
if
(
!
b
->
oppositeCell
()
->
isFrozen
()
&&
!
b
->
next
()
->
oppositeCell
()
->
isFrozen
()){
uint
fidx
=
b
->
from
()
->
vertexStateIndex
();
uint
tidx
=
b
->
to
()
->
vertexStateIndex
();
double
we
=
getWeight
(
b
,
y
);
b
->
relativePosition
(
y
[
fidx
],
y
[
fidx
+
1
],
y
[
tidx
],
y
[
tidx
+
1
],
from1
,
to1
);
// b->relativePosition(from1, to1);
Vec2
v
=
to1
-
from1
;
double
delta_press
=
pressure
(
b
->
cell
())
-
pressure
(
b
->
oppositeCell
());
Vec2
vperp
(
v
[
1
],
-
v
[
0
]);
dydt
[
tidx
]
-=
k_tension
*
we
*
v
[
0
]
+
k_pressure
*
delta_press
*
vperp
[
0
];
dydt
[
tidx
+
1
]
-=
k_tension
*
we
*
v
[
1
]
+
k_pressure
*
delta_press
*
vperp
[
1
];
}
}
//TODO update scaling parameters
}
void
MechanicalForces
::
setLattice
(
Lattice
*
lat
,
double
*
celltypes
){
lattice
=
lat
;
uint
n
=
lat
->
getNbCells
();
m_pressure
.
resize
(
n
);
m_signal
.
resize
(
n
);
m_area
.
resize
(
n
);
m_perimeter
.
resize
(
n
);
if
(
celltypes
){
m_celltype
.
resize
(
n
);
for
(
uint
i
=
0
;
i
<
n
;
i
++
){
m_celltype
[
i
]
=
celltypes
[
i
]
>
15
?
'h'
:
's'
;
}
}
uint
m
=
setStateVectorIndex
();
// observer.resize(m);
// state.resize(2*m+2);
}
uint
MechanicalForces
::
setStateVectorIndex
(){
uint
n
=
lattice
->
getNbVertices
();
for
(
uint
i
=
0
;
i
<
n
;
i
++
){
lattice
->
getVertex
(
i
).
setVertexStateIndex
(
2
*
i
);
}
return
2
*
n
+
2
;
//to include scaling factors
}
int
MechanicalForces
::
initState
(
state_type
&
s
,
bool
resize
){
uint
n
=
lattice
->
getNbVertices
();
if
(
resize
){
s
.
resize
(
2
*
n
+
2
);
}
for
(
uint
i
=
0
;
i
<
n
;
i
++
){
// std::cout<<i<<"/"<<n<<std::endl;
const
Vec2
&
v
=
lattice
->
getVertex
(
i
).
pos
();
s
[
2
*
i
]
=
v
[
0
];
s
[
2
*
i
+
1
]
=
v
[
1
];
}
s
[
2
*
n
]
=
lattice
->
getScale
()[
0
];
s
[
2
*
n
+
1
]
=
lattice
->
getScale
()[
1
];
//observer(s,0);//maybe not the best place to put it
state
=&
s
;
return
2
*
n
+
2
;
}
bool
MechanicalForces
::
fromState
(
const
state_type
&
y
){
uint
n
=
lattice
->
getNbVertices
();
if
(
y
.
size
()
<
2
*
n
+
2
){
return
false
;}
for
(
uint
i
=
0
;
i
<
n
;
i
++
){
Vertex
&
v
=
lattice
->
getVertex
(
i
);
v
.
setPos
(
y
[
v
.
vertexStateIndex
()],
y
[
v
.
vertexStateIndex
()
+
1
]);
}
lattice
->
setScale
(
y
[
2
*
n
],
y
[
2
*
n
+
1
]);
lattice
->
updateAllOrientations
();
updateGeometry
(
y
);
updateSignal
(
y
);
return
true
;
}
// static method
// void MechanicalForces::checkT1TransitionWrapper(Boundary *b, void *arg){
// MechanicalForces *forces = (MechanicalForces *)arg;
// if(!b->isGhost()){
// forces->checkT1Transition(b);
// }
// }
//should ulitmately be replaced by a random process
void
MechanicalForces
::
checkT1Transition
(
Boundary
*
b
){
// std::cout<<"check"<<std::endl;
//boundaries between two HC are unstable
if
(
isHC
(
b
->
cell
())
&&
isHC
(
b
->
oppositeCell
())){
double
len
=
b
->
scaledLength
(
Vec2
(
getScaleX
(),
getScaleY
()));
if
(
len
<
small_bound_thresh
){
if
(
b
->
cell
()
->
getNbBoundaries
()
>
3
){
cout
<<
"1. correcting between "
<<
b
->
cell
()
->
id
()
<<
" "
<<
b
->
oppositeCell
()
->
id
()
<<
endl
;
cout
<<
" sig "
<<
m_signal
[
b
->
cell
()
->
id
()]
<<
" and "
<<
m_signal
[
b
->
oppositeCell
()
->
id
()]
<<
endl
;
b
->
doT1Transition
();
// cout<<"done"<<endl;
}
}
}
else
{
if
(
isSC
(
b
->
cell
())
&&
isSC
(
b
->
oppositeCell
())){
if
(
isHC
(
b
->
next
()
->
oppositeCell
())
&&
isSC
(
b
->
prev
()
->
oppositeCell
())){
Vec2
rpos1
=
b
->
prev
()
->
opposite
()
->
tailPositionRelativeToCell
();
Vec2
rpos2
=
b
->
next
()
->
opposite
()
->
headPositionRelativeToCell
();
float
ori1
=
fabs
(
rpos1
.
sinAngle
(
Vec2
(
1
,
0
)));
float
ori2
=
fabs
(
rpos2
.
sinAngle
(
Vec2
(
1
,
0
)));
if
(
ori1
>
0.8
&&
ori2
>
0.8
){
double
len
=
b
->
scaledLength
(
Vec2
(
getScaleX
(),
getScaleY
()));
if
(
len
<
0.2
){
cout
<<
"2. correcting between "
<<
b
->
cell
()
->
id
()
<<
" "
<<
b
->
oppositeCell
()
->
id
()
<<
endl
;
if
(
b
->
cell
()
->
getNbBoundaries
()
>
3
){
b
->
doT1Transition
();
// cout<<"done"<<endl;
}
}
}
}
}
}
}
void
MechanicalForces
::
findT1Transitions
(){
// uint nbound = lattice->getNbBoundaries();
// for(uint i=0;i<nbound;i++){
// Boundary& b = lattice->getBoundary(i);
// if(!b.isGhost()){
// checkT1Transition(&b);
// cout<<" ok with "<<b.id()<<endl;
// }
// }
vector
<
Boundary
*>
cand
;
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
Cell
&
ce
=
lattice
->
getCell
(
c
);
if
(
!
ce
.
isEmpty
()
&&
!
ce
.
isFrozen
()){
BoundaryIterator
it
(
&
ce
);
for
(
Boundary
*
b
=
it
.
first
();
!
it
.
isLast
();
b
=
it
.
next
()){
cand
.
push_back
(
b
);
}
while
(
cand
.
size
()){
checkT1Transition
(
cand
.
back
());
cand
.
pop_back
();
}
}
}
}
void
MechanicalForces
::
fillSignalVector
(
double
*
v
){
uint
nbcells
=
lattice
->
getNbCells
();
for
(
uint
c
=
0
;
c
<
nbcells
;
c
++
){
v
[
c
]
=
m_signal
[
c
];
// cout<<"signal "<<c<<" "<<m_signal[c]<<endl;;
}
}
// not tested
// double MechanicalForces::dTension(Cell *c, state_type &y,state_type &dydt){
// BoundaryIterator it(c);
// Vec2 from1,to1;
// double sum=0;
// for(Boundary *b=it.first();!it.isLast();b=it.next()){
// uint fidx = b->from()->vertexStateIndex();
// uint tidx = b->to()->vertexStateIndex();
// b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
// Vec2 diff = to1-from1;
// double nor = diff.weightedNorm(Vec(scaleX(),scaleY()));
// if(nor<EPSILON){
// dydt[fidx] -= diff[0]*scaleX()/nor;
// dydt[tidx] += diff[0]*scaleX()/nor;
// dydt[fidx+1] -= diff[1]*scaleY()/nor;
// dydt[tidx+1] += diff[1]*scaleY()/nor;
// }
// sum+=nor;
// }
// return sum;
// }
// bool MechanicalForces::stopIntegration(const state_type &y)
// {return getInstance()->delaminationTime();}
Event Timeline
Log In to Comment