...
 
Commits (2)
......@@ -122,6 +122,8 @@ BaseReaction::createReaction(std::vector<double> &paraValue,
else if(idValue=="WallMechanics::SpringEpidermalCell" ||
idValue=="VertexFromEpidermalCellWallSpring")
return new WallMechanics::SpringEpidermalCell(paraValue,indValue);
else if(idValue=="WallMechanics::ViscoElastic")
return new WallMechanics::ViscoElastic(paraValue,indValue);
else if(idValue=="VertexFromWallSpringMTConcentrationHill")
return new VertexFromWallSpringMTConcentrationHill(paraValue,indValue);
else if(idValue=="VertexFromDoubleWallSpringMTConcentrationHill")
......
......@@ -554,6 +554,139 @@ namespace WallMechanics {
}
}
}
ViscoElastic::
ViscoElastic(std::vector<double> &paraValue,
std::vector< std::vector<size_t> >
&indValue )
{
// Do some checks on the parameters and variable indeces
if( paraValue.size()!=2 ) {
std::cerr << "WallMechanics::ViscoElastic::"
<< "ViscoElastic() "
<< "Uses two parameters viscosity and k_spring." << std::endl;
exit(EXIT_FAILURE);
}
if( (indValue.size() < 1 ||
indValue.size() >2 ||
indValue[0].size() != 1 ||
(indValue.size() == 2 && indValue[1].size() != 1 )) ) {
std::cerr << "WallMechanics::ViscoElastic::"
<< "ViscoElastic() "
<< "Wall length index given in first level. "
<< "If two levels are given, wall index for saving strain is in second level. " << std::endl;
exit(EXIT_FAILURE);
}
// Set the variable values
setId("WallMechanics::ViscoElastic");
setParameter(paraValue);
setVariableIndex(indValue);
// Set the parameter identities
std::vector<std::string> tmp( numParameter() );
tmp[0] = "eta";
tmp[1] = "k_spring";
setParameterId( tmp );
}
void ViscoElastic::
derivs(Tissue &T,
DataMatrix &cellData,
DataMatrix &wallData,
DataMatrix &vertexData,
DataMatrix &cellDerivs,
DataMatrix &wallDerivs,
DataMatrix &vertexDerivs ) {
// Set parameters
double h = 0.1; // ad hoc value (works with constant step size for solver set to this value?)
double eta = parameter(0);
double k_spring = parameter(1);
// Do the update for each wall
size_t numWalls = T.numWall();
size_t wallLengthIndex = variableIndex(0,0);
// relaxation time
double B = -k_spring/eta;
for( size_t i=0 ; i<numWalls ; ++i ) {
size_t v1 = T.wall(i).vertex1()->index();
size_t v2 = T.wall(i).vertex2()->index();
size_t dimension = vertexData[v1].size();
assert( vertexData[v2].size()==dimension );
// Calculate shared factors
double distance=0.0;
for( size_t d=0 ; d<dimension ; d++ )
distance += (vertexData[v1][d]-vertexData[v2][d])*
(vertexData[v1][d]-vertexData[v2][d]);
distance = std::sqrt(distance);
double wallLength=wallData[i][wallLengthIndex];
// Calculate integral over old times
double integral = 0.0;
size_t numStep = historyTime_.size();
size_t finalStep = numStep-1;
size_t t = historyTime_[finalStep];
for (size_t step=0; step<numStep; ++step) {
integral += k_spring*historyData_[i][step]*std::exp(B*(t-historyTime_[step])); //should have dt
}
// Add to the integral (assuming du/dt constant)
integral += k_spring*(std::exp(B*(h)) + 1.0)*historyData_[i][finalStep] ; //should have dt
// Elastic + Viscous contributions
double coeff = parameter(1)*((1.0/wallLength)-(1.0/distance))
+ integral;
//Save force in wall variable if appropriate
if( (numVariableIndexLevel()==2 && numVariableIndex(1)>0) ) {
wallData[i][variableIndex(1,0)] = coeff;
}
//Update both vertices for each dimension
for(size_t d=0 ; d<dimension ; d++ ) {
double div = (vertexData[v1][d]-vertexData[v2][d])*coeff;
vertexDerivs[v1][d] -= div;
vertexDerivs[v2][d] += div;
}
}
}
void ViscoElastic::update(Tissue &T,
DataMatrix &cellData,
DataMatrix &wallData,
DataMatrix &vertexData,
double h) {
size_t lengthIndex = variableIndex(0,0);
historyTime_.push_back(h);
size_t numWall = T.numWall();
for (size_t wallIndex=0; wallIndex<numWall; ++wallIndex) {
historyData_[wallIndex].push_back(wallData[lengthIndex][wallIndex]);
// To Do! Should be...
//historyData_[wallIndex].push_back(wallDerivs[lengthIndex][wallIndex]);
}
}
void ViscoElastic::initiate(Tissue &T,
DataMatrix &cellData,
DataMatrix &walldata,
DataMatrix &vertexData,
DataMatrix &cellderivs,
DataMatrix &wallderivs,
DataMatrix &vertexDerivs )
{
historyData_.reserve(10000);
historyTime_.reserve(10000);
}
} // end namespace WallMechanics
VertexFromWallSpringMTnew::
......
......@@ -314,6 +314,82 @@ namespace WallMechanics {
DataMatrix &wallDerivs,
DataMatrix &vertexDerivs);
};
///
/// @brief A viscoelastic update of edges
///
/// The update (in all dimensions) are given by
///
/// @f[ \frac{dx_i}{dt} = @f]
///
/// In a model file the reaction is defined as:
///
/// @verbatim
/// WallMechanics::ViscoElastic 2 1 1
/// eta # viscosity
/// k # spring constant
/// L_ij-index # index of the resting lenth
/// @endverbatim
/// or
/// @verbatim
/// WallMechanics::ViscoElastic 2 2 1 1
/// eta # viscosity
/// k # spring constant
/// L_ij-index # index of the resting lenth
/// strain_index # optional wall index to store strain
/// @endverbatim
///
class ViscoElastic : public BaseReaction {
private:
DataMatrix historyData_;
std::vector<double> historyTime_;
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,...)
///
ViscoElastic(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 );
void initiate(Tissue &T,
DataMatrix &cellData,
DataMatrix &walldata,
DataMatrix &vertexData,
DataMatrix &cellderivs,
DataMatrix &wallderivs,
DataMatrix &vertexDerivs );
void update(Tissue &T,
DataMatrix &cellData,
DataMatrix &wallData,
DataMatrix &vertexData,
double h);
};
} // end namespace WallMechanics
///
......
......@@ -15,6 +15,157 @@
#include <iostream>
#include <sstream>
namespace TRBS {
double AreaFromEdges(std::vector<double> &edgeLength) {
assert(edgeLength.size()==3);
return std::sqrt( ( edgeLength[0]+edgeLength[1]+edgeLength[2])*
(-edgeLength[0]+edgeLength[1]+edgeLength[2])*
( edgeLength[0]-edgeLength[1]+edgeLength[2])*
( edgeLength[0]+edgeLength[1]-edgeLength[2]) )*0.25;
}
void CosFromEdges(std::vector<double> &restingLength, std::vector<double> &cosAngle) {
assert(restingLength.size()==3 && cosAngle.size()==3);
cosAngle[0] = (restingLength[0]*restingLength[0]+restingLength[2]*restingLength[2]
-restingLength[1]*restingLength[1])
/ (restingLength[0]*restingLength[2]*2);
cosAngle[1] = (restingLength[0]*restingLength[0]+restingLength[1]*restingLength[1]
-restingLength[2]*restingLength[2])
/ (restingLength[0]*restingLength[1]*2);
cosAngle[2] = (restingLength[1]*restingLength[1]+restingLength[2]*restingLength[2]
-restingLength[0]*restingLength[0])
/ (restingLength[1]*restingLength[2]*2);
}
void CotanFromCos(std::vector<double> &cosAngle, std::vector<double> &cotanAngle) {
assert(cosAngle.size()==cotanAngle.size());
for (size_t i=0;i<cotanAngle.size(); ++i) {
cotanAngle[i] = cosAngle[i]/std::sqrt(1-cosAngle[i]*cosAngle[i]);
}
}
void Stiffness(double lambda, double mio, double areaFactor, std::vector<double> &cotan,
std::vector<double> &tensileStiffness, std::vector<double> &angularStiffness) {
assert( cotan.size()==tensileStiffnes.size() && cotan.size()==angularStiffness.size());
double fac = lambda+2.0*mio;
size_t numVertex=cotan.size();
for (size_t i=0; i<numVertex; ++i) {
size_t ip1 = (i+1)%numVertex;
size_t im1 = (i+2)%numVertex; //assumes a triangle
tensileStiffness[i] = (2*cotan[im1]*cotan[im1]*fac+2*mio)*areaFactor;
angularStiffness[i] = (2*cotan[im1]*cotan[ip1]*fac-2*mio)*areaFactor;
}
}
void BiQuadraticStrain(std::vector<double> &length, std::vector<double> &restingLength,
std::vector<double> &Delta) {
size_t numVertex=length.size();
assert(restingLength.size()==numVertex && Delta.size()==numVertex);
for (size_t i=0; i<numVertex; ++i) {
Delta[i] = length[i]*length[i] - restingLength[i]*restingLength[i];
}
}
void ShapeVector(double Xa, double Xb, double Xc, std::vector< std::vector<double> > &shapeVector) {
shapeVector[0][0] = shapeVector[0][2] = shapeVector[2][2] = 0.0;
shapeVector[0][1] = 1/Xc;
shapeVector[1][0] = -1.0/Xb;
shapeVector[1][1] = (Xa-Xb)/(Xb*Xc);
shapeVector[1][2] = 1.0;
shapeVector[2][0] = 1.0/Xb;
shapeVector[2][1] = -Xa/(Xb*Xc);
}
void AddForceIsotropic(std::vector< std::vector<double> > &pos,
std::vector<double> &tS, std::vector<double> &aS,
std::vector<double> &Delta,std::vector< std::vector<double> > &Force) {
size_t dimension = pos.size();
assert(tS.size()==dimension && aS.size()==dimension && Delta.size()==dimension && Force.size()==dimension);
//for (size_t i=0; i<dimension; ++i) {
//for (size_t j=0; j<dimension; ++j) {
// Force[i][j]= (tS[0]*Delta[0]+aS[1]*Delta[1]+aS[0]*Delta[2])*(pos[1][0]-pos[i][j])
// + (tS[2]*Delta[2]+aS[2]*Delta[1]+aS[0]*Delta[0])*(pos[2][0]-pos[i][j]);
//}
//}
Force[0][0]= (tS[0]*Delta[0]+aS[1]*Delta[1]+aS[0]*Delta[2])*(pos[1][0]-pos[0][0])
+ (tS[2]*Delta[2]+aS[2]*Delta[1]+aS[0]*Delta[0])*(pos[2][0]-pos[0][0]);
Force[0][1]= (tS[0]*Delta[0]+aS[1]*Delta[1]+aS[0]*Delta[2])*(pos[1][1]-pos[0][1])
+ (tS[2]*Delta[2]+aS[2]*Delta[1]+aS[0]*Delta[0])*(pos[2][1]-pos[0][1]);
Force[0][2]= (tS[0]*Delta[0]+aS[1]*Delta[1]+aS[0]*Delta[2])*(pos[1][2]-pos[0][2])
+ (tS[2]*Delta[2]+aS[2]*Delta[1]+aS[0]*Delta[0])*(pos[2][2]-pos[0][2]);
Force[1][0]= (tS[0]*Delta[0]+aS[0]*Delta[2]+aS[1]*Delta[1])*(pos[0][0]-pos[1][0])
+ (tS[1]*Delta[1]+aS[2]*Delta[2]+aS[1]*Delta[0])*(pos[2][0]-pos[1][0]);
Force[1][1]= (tS[0]*Delta[0]+aS[0]*Delta[2]+aS[1]*Delta[1])*(pos[0][1]-pos[1][1])
+ (tS[1]*Delta[1]+aS[2]*Delta[2]+aS[1]*Delta[0])*(pos[2][1]-pos[1][1]);
Force[1][2]= (tS[0]*Delta[0]+aS[0]*Delta[2]+aS[1]*Delta[1])*(pos[0][2]-pos[1][2])
+ (tS[1]*Delta[1]+aS[2]*Delta[2]+aS[1]*Delta[0])*(pos[2][2]-pos[1][2]);
Force[2][0]= (tS[2]*Delta[2]+aS[0]*Delta[0]+aS[2]*Delta[1])*(pos[0][0]-pos[2][0])
+ (tS[1]*Delta[1]+aS[1]*Delta[0]+aS[2]*Delta[2])*(pos[1][0]-pos[2][0]);
Force[2][1]= (tS[2]*Delta[2]+aS[0]*Delta[0]+aS[2]*Delta[1])*(pos[0][1]-pos[2][1])
+ (tS[1]*Delta[1]+aS[1]*Delta[0]+aS[2]*Delta[2])*(pos[1][1]-pos[2][1]);
Force[2][2]= (tS[2]*Delta[2]+aS[0]*Delta[0]+aS[2]*Delta[1])*(pos[0][2]-pos[2][2])
+ (tS[1]*Delta[1]+aS[1]*Delta[0]+aS[2]*Delta[2])*(pos[1][2]-pos[2][2]);
}
void RotationMatrix(std::vector< std::vector<double> > &pos,std::vector< std::vector<double> > &rotation) {
size_t dimension=pos.size();
assert( rot.size()==dimension);
std::vector<double> xCurrent(dimension), bCurrent(dimension), zCurrent(dimension), yCurrent(dimension);
double temp=std::sqrt((pos[2][0]-pos[1][0])*(pos[2][0]-pos[1][0])+
(pos[2][1]-pos[1][1])*(pos[2][1]-pos[1][1])+
(pos[2][2]-pos[1][2])*(pos[2][2]-pos[1][2]) );
double tempB=std::sqrt((pos[0][0]-pos[1][0])*(pos[0][0]-pos[1][0])+
(pos[0][1]-pos[1][1])*(pos[0][1]-pos[1][1])+
(pos[0][2]-pos[1][2])*(pos[0][2]-pos[1][2]) );
for (size_t i=0; i<dimension; ++i) {
xCurrent[i] = (pos[2][i]-pos[1][i])/temp;
bCurrent[i] = (pos[0][i]-pos[1][i])/tempB;
}
zCurrent[0] = xCurrent[1]*bCurrent[2]-xCurrent[2]*bCurrent[1];
zCurrent[1] = xCurrent[2]*bCurrent[0]-xCurrent[0]*bCurrent[2];
zCurrent[2] = xCurrent[0]*bCurrent[1]-xCurrent[1]*bCurrent[0];
//normalise
temp = 1.0/std::sqrt(zCurrent[0]*zCurrent[0]+zCurrent[1]*zCurrent[1]+zCurrent[2]*zCurrent[2]);
for (size_t i=0; i<dimension; ++i) {
zCurrent[i] *= temp;
}
yCurrent[0] = zCurrent[1]*xCurrent[2]-zCurrent[2]*xCurrent[1];
yCurrent[1] = zCurrent[2]*xCurrent[0]-zCurrent[0]*xCurrent[2];
yCurrent[2] = zCurrent[0]*xCurrent[1]-zCurrent[1]*xCurrent[0];
for (size_t i=0; i<dimension; ++i) {
rotation[i][0] = xCurrent[i];
rotation[i][1] = yCurrent[i];
rotation[i][2] = zCurrent[i];
}
}
void Rotate(std::vector<double> &inVector, std::vector< std::vector<double> > &rotation,
std::vector<double> &outVector) {
size_t dimension = rotation.size();
for (size_t i=0; i<dimension; ++i) {
outVector[i]=0.0;
for (size_t j=0; j<dimension; ++j) {
outVector[i] += rotation[j][i]*inVector[j];
}
}
}
} // end namespace TRBS
VertexFromTRBS::
VertexFromTRBS(std::vector<double> &paraValue,
std::vector< std::vector<size_t> >
......
......@@ -12,6 +12,76 @@
#include"baseReaction.h"
#include<cmath>
///
/// @brief Mechanical models using the Triangular Biquadratic Spring (TRBS) for triangular plates
///
/// @details These reactions update the vertex positions from the mechanical feedback on 2D triangular
/// elements.
/// The theory of the mechanical model comes from H. Delingette,
/// Triangular springs for modelling non-linear membranes, IEEE Trans
/// Vis Comput Graph 14, 329-41 (2008)
///
/// @see https://gitlab.com/slcu/teamhj/publications/bozorg_etal_2014 (and bozorg_etal_2016) for examples
///
namespace TRBS {
// Some functions for general triangle manipulations
///
/// @brief Calculates the triangular area from edges using Heron's formula
///
double AreaFromEdges(std::vector<double> &edgeLength);
///
/// @brief Calculates the cosine of the triangle angles from the edge lengths
///
void CosFromEdges(std::vector<double> &restingLength, std::vector<double> &cosAngle);
///
/// @brief Calculates the Cotanges of the angles of the triangle from the cosines
///
void CotanFromCos(std::vector<double> &cosAngle, std::vector<double> &cotanAngle);
///
/// @brief Calculates the tensile and angular stiffness from cotan vector
///
void Stiffness(double lambda, double mio, double areaFactor, std::vector<double> &cotan,
std::vector<double> &tensilStiffness, std::vector<double> &angularStiffness);
///
/// @brief Calculates the BiQuadratic strain
///
void BiQuadraticStrain(std::vector<double> &length, std::vector<double> &restingLength,
std::vector<double> &Delta);
///
/// @brief Shape vector matrix in local coordinate system from shape factors (P or Q = X)
///
/// @details This matrix is the inverse of coordinate matrix. Only first two elements are used in calculations
/// i.e. shapeVector[3][2] although implemented as 3x3.
///
void ShapeVector(double Xa, double Xb, double Xc, std::vector< std::vector<double> > &shapeVector);
///
/// @brief Calculates Forces from isotropic contribution of the material (youngT,poissonT)
///
/// The calculation uses positions, tensile and angular stiffnesses and Delta and adds to the Force matrix.
///
void AddForceIsotropic(std::vector< std::vector<double> > &position,
std::vector<double> &tensileStiffness, std::vector<double> &angularStiffness,
std::vector<double> &Delta,
std::vector< std::vector<double> > &Force);
///
/// @brief Calculates the rotation matrix from positions
///
void RotationMatrix(std::vector< std::vector<double> > &position, std::vector< std::vector<double> > &rotation);
///
/// @brief Rotates inVector into outVector using the rotation matrix
///
void Rotate(std::vector<double> &inVector, std::vector< std::vector<double> > &rotation,
std::vector<double> &outVector);
namespace CenterTriangulation {
} // end namespace CenterTriangulation
} // end namespace TRBS
///
/// @brief Triangular spring model for plates (2D walls) assuming
/// triangular faces (cells).
......