Commit 372508ba authored by Ross Carter's avatar Ross Carter

Added new reaction for generating a growth force in the direction of

another cell, not implemented yet
parent 1df7751c
......@@ -256,6 +256,9 @@ BaseReaction::createReaction(std::vector<double> &paraValue,
idValue=="CenterTriangulation::GrowthForce::Radial" ||
idValue == "MoveVertexRadiallycenterTriangulation")
return new GrowthForce::CenterTriangulation::Radial(paraValue, indValue);
else if (idValue=="GrowthForce::centerTriangulation::ForceToCell" ||
idValue=="centerTriangulation::GrowthForce::ForceToCell")
return new GrowthForce::CenterTriangulation::ForceToCell(paraValue, indValue);
else if(idValue=="GrowthForce::EpidermalRadial" ||
idValue == "MoveEpidermalVertexRadially")
return new GrowthForce::EpidermalRadial(paraValue, indValue);
......
......@@ -142,39 +142,126 @@ namespace GrowthForce {
// Move vertices / apply force to vertices
for( size_t i=0 ; i<numVertices ; ++i ) {
double fac=parameter(0);
if( parameter(1)==0.0 ) {
double r=0.0;
for( size_t d=0 ; d<dimension ; ++d )
r += vertexData[i][d]*vertexData[i][d];
if( r>0.0 )
r = std::sqrt(r);
if( r>0.0 )
fac /= r;
else
fac=0.0;
}
for( size_t d=0 ; d<dimension ; ++d )
vertexDerivs[i][d] += fac*vertexData[i][d];
double fac=parameter(0);
if( parameter(1)==0.0 ) {
double r=0.0;
for( size_t d=0 ; d<dimension ; ++d )
r += vertexData[i][d]*vertexData[i][d];
if( r>0.0 )
r = std::sqrt(r);
if( r>0.0 )
fac /= r;
else
fac=0.0;
}
for( size_t d=0 ; d<dimension ; ++d )
vertexDerivs[i][d] += fac*vertexData[i][d];
}
// Move / apply force to vertices defined in cell centers
for( size_t i=0 ; i<numCells ; ++i ) {
double fac=parameter(0);
if( parameter(1)==0.0 ) {
double r=0.0;
for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
r += cellData[i][d]*cellData[i][d];
if( r>0.0 )
r = std::sqrt(r);
if( r>0.0 )
fac /= r;
else
fac=0.0;
}
for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
cellDerivs[i][d] += fac*cellData[i][d];
double fac=parameter(0);
if( parameter(1)==0.0 ) {
double r=0.0;
for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
r += cellData[i][d]*cellData[i][d];
if( r>0.0 )
r = std::sqrt(r);
if( r>0.0 )
fac /= r;
else
fac=0.0;
}
for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
cellDerivs[i][d] += fac*cellData[i][d];
}
}
ForceToCell::
ForceToCell(std::vector<double> &paraValue,
std::vector< std::vector<size_t> >
&indValue ) {
// // Do some checks on the parameters and variable indeces
// //
// if( paraValue.size()!=2 || ( paraValue[1]!=0 && paraValue[1]!=1) ) {
// std::cerr << "GrowthForce::CenterTriangulation::ForceToCell::"
// << "ForceToCell() " << std::endl
// << "Uses two parameters k_growth and r_pow (0,1)" << std::endl;
// exit(EXIT_FAILURE);
// }
// if( indValue.size() != 1 || indValue[0].size() != 1 ) {
// std::cerr << "GrowthForce::CenterTriangulation::ForceToCell::"
// << "ForceToCell() " << std::endl
// << "Start of additional Cell variable indices (center(x,y,z) "
// << "L_1,...,L_n, n=num vertex) is given in first level. "
// << "See Documentation for namespace CenterTriangulation."
// << std::endl;
// exit(EXIT_FAILURE);
// }
// // Set the variable values
// //
// setId("GrowthForce::CenterTriangulation::ForceToCell");
// setParameter(paraValue);
// setVariableIndex(indValue);
// // Set the parameter identities
// //
// std::vector<std::string> tmp( numParameter() );
// tmp.resize( numParameter() );
// tmp[0] = "k_growth";
// tmp[0] = "r_pow";
// setParameterId( tmp );
}
void ForceToCell::
derivs(Tissue &T,
DataMatrix &cellData,
DataMatrix &wallData,
DataMatrix &vertexData,
DataMatrix &cellDerivs,
DataMatrix &wallDerivs,
DataMatrix &vertexDerivs ) {
// size_t numVertices = T.numVertex();
// size_t numCells = T.numCell();
// size_t dimension=vertexData[0].size();
// // Move vertices / apply force to vertices
// for( size_t i=0 ; i<numVertices ; ++i ) {
// double fac=parameter(0);
// if( parameter(1)==0.0 ) {
// double r=0.0;
// for( size_t d=0 ; d<dimension ; ++d )
// r += vertexData[i][d]*vertexData[i][d];
// if( r>0.0 )
// r = std::sqrt(r);
// if( r>0.0 )
// fac /= r;
// else
// fac=0.0;
// }
// for( size_t d=0 ; d<dimension ; ++d )
// vertexDerivs[i][d] += fac*vertexData[i][d];
// }
// // Move / apply force to vertices defined in cell centers
// for( size_t i=0 ; i<numCells ; ++i ) {
// double fac=parameter(0);
// if( parameter(1)==0.0 ) {
// double r=0.0;
// for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
// r += cellData[i][d]*cellData[i][d];
// if( r>0.0 )
// r = std::sqrt(r);
// if( r>0.0 )
// fac /= r;
// else
// fac=0.0;
// }
// for( size_t d=variableIndex(0,0) ; d<variableIndex(0,0)+dimension ; ++d )
// cellDerivs[i][d] += fac*cellData[i][d];
// }
}
} // end namespace CenterTriangulation
EpidermalRadial::
......
......@@ -142,6 +142,62 @@ namespace GrowthForce {
DataMatrix &wallDerivs,
DataMatrix &vertexDerivs );
};
///
/// @brief Growth via a force acting towards the centroid of a cell whose
/// index is stored in cell data
///
/// UPDATE BELOW
///
/// The tissue grows from vertex movement radially outwards, and also
/// includes moving the vertex defining the 'center' of the cells in the
/// center triangulated mesh. The update is given by
/// @f[ \frac{dr}{dt} = p_{0} @f] (if @f$ p_1=0 @f$) or
/// @f[ \frac{dr}{dt} = p_{0} r @f] (if @f$ p_{1}=1 @f$)
///
/// @f$ p_{0} @f$ is the rate (@f$ k_{growth} @f$),
/// @f$ p_{1} @f$ {0,1} is a flag determining which function to be used (@f$ r_{pow} @f$).
/// In a model file the reaction is defined as
/// @verbatim
/// GrowthForce::CenterTriangulation::Radial 2 1 1
/// p_0 p_1
/// InternalVarStartIndex
/// @endverbatim
/// @see Force::Radial (same but without moving the central vertices)
/// @note Used to be named MoveVertexRadiallycenterTriangulation (still allowed).
///
class ForceToCell : public BaseReaction {
public:
///
/// @brief Main constructor
///
/// This is the main constructor which sets the parameters and variable
/// indices that defines the reaction.
///
/// @param paraValue vector with parameters
///
/// @param indValue vector of vectors with variable indices
///
/// @see BaseReaction::createReaction(std::vector<double> &paraValue,...)
///
ForceToCell(std::vector<double> &paraValue,
std::vector< std::vector<size_t> >
&indValue );
///
/// @brief Derivative function for this reaction class
///
/// @see BaseReaction::derivs(Tissue &T,...)
///
void derivs(Tissue &T,
DataMatrix &cellData,
DataMatrix &wallData,
DataMatrix &vertexData,
DataMatrix &cellDerivs,
DataMatrix &wallDerivs,
DataMatrix &vertexDerivs );
};
} // end namespace CenterTriangulation
///
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment