Commit 46309652 authored by Thomas Kummer's avatar Thomas Kummer
Browse files

End visit

parent 10cc253d
No preview for this file type
......@@ -5394,6 +5394,9 @@
<Group
location = "group:example_XBridge"
name = "example_XBridge">
<FileRef
location = "group:data">
</FileRef>
<FileRef
location = "group:CMakeLists.txt">
</FileRef>
......@@ -5832,6 +5835,12 @@
<Group
location = "group:activeStressModels"
name = "activeStressModels">
<FileRef
location = "group:ActiveStressRossiModel14XB.cpp">
</FileRef>
<FileRef
location = "group:ActiveStressRossiModel14XB.hpp">
</FileRef>
<FileRef
location = "group:ActiveStressActivation.hpp">
</FileRef>
......
No preview for this file type
......@@ -221,6 +221,17 @@ public:
ElectroETAMonodomainSolver (std::string meshName, std::string meshPath,
GetPot& dataFile, ionicModelPtr_Type model);
//! Constructor
/*!
* @param meshName file name of the mesh
* @param meshPath path to the mesh
* @param datafile GetPot file to setup the preconditioner
* @param model shared pointer to the chosen ionic model
* @param list xml Parameter list data
*/
ElectroETAMonodomainSolver (std::string meshName, std::string meshPath,
GetPot& dataFile, ionicModelPtr_Type model, list_Type list);
//! Constructor
/*!
* @param string file name of the mesh
......@@ -1450,6 +1461,18 @@ ElectroETAMonodomainSolver<Mesh>::ElectroETAMonodomainSolver (
setup (meshName, meshPath, dataFile, M_ionicModelPtr->Size() );
}
//!constructor
template<typename Mesh>
ElectroETAMonodomainSolver<Mesh>::ElectroETAMonodomainSolver (
std::string meshName, std::string meshPath, GetPot& dataFile,
ionicModelPtr_Type model, list_Type list)
{
M_verbose = false;
setParameters(list);
init (model);
setup (meshName, meshPath, dataFile, M_ionicModelPtr->Size() );
}
//!constructor with communicator
template<typename Mesh>
ElectroETAMonodomainSolver<Mesh>::ElectroETAMonodomainSolver (
......@@ -1719,7 +1742,7 @@ void ElectroETAMonodomainSolver<Mesh>::setupLumpedMassMatrix()
}
{
using namespace ExpressionAssembly;
// Todo: Check here whether p1 or p2, change quadRule to maybe 10pt.
integrate (elements (M_localMeshPtr), quadRuleTetra4ptNodal,
M_ETFESpacePtr, M_ETFESpacePtr, phi_i * phi_j)
>> M_massMatrixPtr;
......@@ -1804,11 +1827,11 @@ template<typename Mesh>
void ElectroETAMonodomainSolver<Mesh>::setupGlobalSolution (
short int ionicSize)
{
M_globalSolution.push_back (M_potentialPtr);
M_globalSolution.resize(ionicSize);
M_globalSolution[0] = M_potentialPtr;
for (int k = 1; k < ionicSize; ++k)
{
M_globalSolution.push_back (
* (new vectorPtr_Type (new VectorEpetra (M_ETFESpacePtr->map() ) ) ) );
M_globalSolution[k].reset (new VectorEpetra (M_ETFESpacePtr->map() ) ) ;
}
}
......@@ -1816,11 +1839,11 @@ template<typename Mesh>
void ElectroETAMonodomainSolver<Mesh>::setupGlobalRhs (
short int ionicSize)
{
M_globalRhs.push_back (M_rhsPtrUnique);
M_globalRhs.resize(ionicSize);
M_globalRhs[0] = M_rhsPtrUnique;
for (int k = 1; k < ionicSize; ++k)
{
M_globalRhs.push_back (
* (new vectorPtr_Type (new VectorEpetra (M_ETFESpacePtr->map() ) ) ) );
M_globalRhs[k].reset ( new VectorEpetra (M_ETFESpacePtr->map() ) );
}
}
......
......@@ -115,7 +115,7 @@ IonicTenTusscher06::IonicTenTusscher06() :
if (flag == Endo)
{
Gks = 0.073;
Gto = 0.073;
}
else
{
......@@ -151,37 +151,37 @@ IonicTenTusscher06::IonicTenTusscher06() :
inversevssF2 = 1. / (2.*Vss * F);
//V
M_restingConditions.at (0) = -86.2;
M_restingConditions.at (0) = -85.4392; //-86.2;
//m
M_restingConditions.at (1) = 0.0;
M_restingConditions.at (1) = 0.00139553; //0.0;
//h
M_restingConditions.at (2) = 0.75;
M_restingConditions.at (2) = 0.750059; //0.75;
//j
M_restingConditions.at (3) = 0.75;
M_restingConditions.at (3) = 0.750006; //0.75;
//d
M_restingConditions.at (4) = 0.;
M_restingConditions.at (4) = 1.48661e-6; //0.;
//f
M_restingConditions.at (5) = 1.;
//f2
M_restingConditions.at (6) = 1;
M_restingConditions.at (6) = 0.999999; //1;
//fCass
M_restingConditions.at (7) = 1.0;
//r
M_restingConditions.at (8) = 0.;
M_restingConditions.at (8) = 1.10759e-10; //0.;
//s
M_restingConditions.at (9) = 1.0;
//Xr1
M_restingConditions.at (10) = 0.0;
M_restingConditions.at (10) = 8.59227e-8; //0.0;
//Xr2
M_restingConditions.at (11) = 1.0;
M_restingConditions.at (11) = 0.985661; //1.0;
//Xs
M_restingConditions.at (12) = 0.0;
M_restingConditions.at (12) = 7.47962e-7; //0.0;
//Nai
M_restingConditions.at (13) = 7.67;
//Ki
M_restingConditions.at (14) = 138.3;
//Cai
M_restingConditions.at (15) = 0.00007;
M_restingConditions.at (15) = 6.99985e-5; //0.00007;
//Cass
M_restingConditions.at (16) = 0.00007;
//Casr
......
......@@ -2,14 +2,14 @@
<Parameter name="surfaceVolumeRatio" type="double" value="1400.0"/><!-- cm ^ -1-->
<Parameter name="timeStep" type="double" value="0.1" /><!-- ms ^ -1-->
<Parameter name="endTime" type="double" value="55." />
<Parameter name="saveStep" type="double" value="1.0" />
<Parameter name="saveStep" type="double" value="0.1" />
<Parameter name="longitudinalDiffusion" type="double" value="1.3342" /><!-- k Ohm ^-1 cm ^ -1-->
<Parameter name="transversalDiffusion" type="double" value ="0.17606" /><!--k Ohm ^-1 cm ^ -1-->
<Parameter name="elementsOrder" type="string" value="P1" />
<Parameter name="ionic_model" type="string" value="TenTusscher06" />
<Parameter name="elementsOrder" type="string" value="P2" />
<Parameter name="ionic_model" type="string" value="MinimalModel" />
<Parameter name="solutionMethod" type="string" value="ICI" />
<Parameter name="subiter" type="int" value="1" />
<Parameter name="LumpedMass" type="bool" value="true"/>
<Parameter name="LumpedMass" type="bool" value="false"/>
<Parameter name="mesh_name" type="string" value="benchmark_05mm.mesh" />
<Parameter name="mesh_path" type="string" value="./" />
<Parameter name="Reuse Preconditioner" type="bool" value="true"/>
......
......@@ -167,7 +167,7 @@ Real PacingProtocolMM ( const Real& t, const Real& x, const Real& y, const Real&
Real pacingSite_Y = 0.0;
Real pacingSite_Z = 0.0;
Real stimulusRadius = 0.15;
Real stimulusValue = 10;
Real stimulusValue = 100;
Real returnValue;
......
......@@ -275,7 +275,7 @@ Int main ( Int argc, char** argv )
std::cout << "\nBuilding Monodomain Solver ... \n";
}
monodomainSolverPtr_Type solver ( new monodomainSolver_Type ( meshName, meshPath, dataFile , model ) );
monodomainSolverPtr_Type solver ( new monodomainSolver_Type ( meshName, meshPath, dataFile , model, monodomainList) );
if ( Comm->MyPID() == 0 )
{
std::cout << "\t... solver created!";
......@@ -316,8 +316,9 @@ Int main ( Int argc, char** argv )
// be around 1e-1 - 1e0.
// ---------------------------------------------------------------
solver -> setParameters ( monodomainList );
//solver -> setParameters ( monodomainList );
//solver -> setup( dataFile, model -> Size() );
// ---------------------------------------------------------------
// Cardiac tissue is typically model as transversely isotropic.
// We need to define the preferred direction (or fiber direction).
......
......@@ -321,7 +321,7 @@ inline void setupFibers ( VectorEpetra& fiberVector, Real fx, Real fy, Real fz)
*/
inline void setValueOnBoundary ( VectorEpetra& vec, boost::shared_ptr< RegionMesh<LinearTetra> > fullMesh, Real value, UInt flag)
{
// fetofespace...
for ( Int j (0); j < vec.epetraVector().MyLength() ; ++j )
{
if ( fullMesh -> point ( vec.blockMap().GID (j) ).markerID() == flag )
......
No preview for this file type
......@@ -34,6 +34,7 @@ SET(HEADERS ${HEADERS}
${solver_mechanics_materials_passive_HEADERS}
${solver_mechanics_materials_activeStress_HEADERS}
${solver_activation_HEADERS}
${solver_activation_xBridgeModels_HEADERS}
${solver_activation_activeStressModels_HEADERS}
${solver_activation_activeStrainModels_HEADERS})
......@@ -46,6 +47,7 @@ SET(SOURCES ${SOURCES}
${solver_mechanics_materials_passive_SOURCES}
${solver_mechanics_materials_activeStress_SOURCES}
${solver_activation_SOURCES}
${solver_activation_xBridgeModels_SOURCES}
${solver_activation_activeStressModels_SOURCES}
${solver_activation_activeStrainModels_SOURCES}
)
......
......@@ -16,18 +16,18 @@ monodomain_xml_file = ParamList.xml
[./physics]
IonicModel = MinimalModel
fiberDiffusion = 10.0
fiberDiffusion = 40.0
sheetDiffusion = 3.0
normalDiffusion = 3.0
[../time_discretization]
endtime = 40
endtime = 10
timestep = 0.05
[../discretization]
LumpedMass = false
LumpedMass = false // use when too coarse mesh
elementsOrder = P1
elementsOrder = P2
[../flags]
lvendo = 36
......@@ -41,27 +41,27 @@ monodomain_xml_file = ParamList.xml
[./physics]
#ActivationModels available and relative required parameters
# ActiveStressRossi14:
# ActiveStressRossi14: # only correct with AlievPa., wrong parameters for others.
# ActiveStress_Beta
# ActiveStressMu
# Tmax ( MaxActiveTension )
# ActiveStressNashPanfilov04
# kTa
# epsilon0
# ActiveStrainRossi14
# ActiveStrainRossi14 # only use with minimalmodel and neoHoekan and Holzapfel and exponential
# calciumIndex
# InverseViscosity
# ActiveForceCoefficient
# ChemicalThreshold
ActivationModel = ActiveStressRossi14
CalciumIndex = 2
CalciumIndex = 0 # potential for activstress, 2 from minimal for active strain.
ActiveForceCoefficient = -4.0
# Tmax = 250 # Default: 50
EMActiveStrainOrthotropicParameter = -666.0
#Active Strain types Available
# Anisotropic
# Orthotropic
# Orthotropic # set kappa = 3 or 4 positiv to be orthotropic
# TransverselyIsotropic
EMActiveStrainType = TransverselyIsotropic
......@@ -98,11 +98,15 @@ monodomain_xml_file = ParamList.xml
#PMR = Passive Mooney Rivlin
#SimpleActiveStress = Simple Fiber Active Stress
EMPassiveMaterialType = PNH
EMPassiveMaterialType = PHO
EMActiveStressMaterialType = SimpleActiveStress #SFAS
Tmax = 30000 #500000
Tmax = 30000 #500000 # 500000 value is thought for Holzapfel
#NeoHokean
mu = 4960
# Holzapfel
a = 3330
af = 185350
as = 25640
......@@ -111,7 +115,9 @@ monodomain_xml_file = ParamList.xml
bf = 15.972
bs = 10.446
bfs = 11.602
BulkModulus = 35000
BulkModulus = 3500000
# Fang type materials
C = 20000.0
bff = 8.0
bss = 2.0
......@@ -122,11 +128,12 @@ monodomain_xml_file = ParamList.xml
[../boundary_conditions]
# list = 'Top Bottom Front Back Left Right' # Cube
list = 'Base'
LvPreloadPressure = 100
numPreloadSteps = 100
LvFlag = 36
[./Base]
type = Essential
flag = 34
......@@ -134,55 +141,13 @@ monodomain_xml_file = ParamList.xml
component = 3
function = '[0.0, 0.0, 0.0]'
[../Top]
type = Natural
flag = 400
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
[../Bottom]
type = Essential
flag = 200
mode = Component
component = 2
function = '0.0'
[../Front]
type = Natural
flag = 500
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
[../Back]
type = Essential
flag = 300
mode = Component
component = 1
function = '0.0'
[../Left]
type = Natural
flag = 100
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
[../Right]
type = Essential
flag = 900
mode = Component
component = 0
function = '0.0'
[../]
[../time_discretization]
initialtime = 0.
endtime = 1000.
timestep = 1
endtime = 100.
timestep = 0.5
theta = 0.35
zeta = 0.75
BDF_order = 2
......@@ -193,12 +158,12 @@ mesh_type = .mesh
mesh_dir = ./
mesh_file = ventFine.mesh
mesh_scaling = 0.1
order = P2
order = P1
anisotropic = true
fiber_name = FiberDirectionP2
fiber_name = FiberDirectionP1
fiber_fieldname = fibers
fiber_dir = ./
sheet_name = SheetsDirectionP2
sheet_name = SheetsDirectionP1
sheet_fieldname = sheets
sheet_dir = ./
......
......@@ -43,12 +43,13 @@ using namespace LifeV;
Real Iapp (const Real& t, const Real& X, const Real& Y, const Real& Z, const ID& /*i*/)
{
return 0; // ( X < 0.25 && Y < 0.25 && Z < 0.25 && t < 2 ? 10 : 0 );
return ( Y > 1.5 && Y < 3 /*std::abs(X) < 1 && std::abs(Z-3) < 1 && Y < 0*/ /*&& Y < 0.25 && Z < 0.25 */ && t < 2 ? 30 : 0 );
// setAppliedCurrent in electrophys. module.
}
Real potentialMultiplyerFcn (const Real& t, const Real& X, const Real& Y, const Real& Z, const ID& /*i*/)
{
return 1; // ( Y < 2.5 && Y > 0.5 ? 1.0 : 0.0 );
return 0; // ( Y < 2.5 && Y > 0.5 ? 1.0 : 0.0 );
}
......@@ -318,9 +319,9 @@ int main (int argc, char** argv)
// Get b.c. flags
ID LvFlag = dataFile ( "solid/boundary_conditions/LvFlag", 0);
// Get preload pressure // Todo: loop for ramp until prel. p. is reached
Real pPreloadLvBC = 0.000001; // ( "solid/boundary_conditions/LvPreloadPressure", 0.0);
Real preloadSteps = 5;
// Get preload pressure and steps
Real pPreloadLvBC = dataFile ( "solid/boundary_conditions/LvPreloadPressure", 0.0);
int preloadSteps = dataFile ( "solid/boundary_conditions/numPreloadSteps", 0);
// Declare pressure b.c. obtained from circulation
Real pLvBC;
......@@ -334,12 +335,7 @@ int main (int argc, char** argv)
solver.bcInterfacePtr() -> handler() -> addBC("LvPressure", LvFlag, Natural, Full, *pLvBCVectorPtr, 3);
solver.bcInterfacePtr() -> handler() -> bcUpdate( *solver.structuralOperatorPtr() -> dispFESpacePtr() -> mesh(), solver.structuralOperatorPtr() -> dispFESpacePtr() -> feBd(), solver.structuralOperatorPtr() -> dispFESpacePtr() -> dof() );
// ComputeBC<solidETFESpace_Type>(localSolidMesh, solid.displacement(), boundaryVectorPtr, dETFESpace, flag);
// Real LVPreloadPressure = dataFile ( "solid/boundary_conditions/LV_preload_pressure", 0.0);
// bool deformedPressure = dataFile ( "solid/boundary_conditions/deformed_pressure", 1 );
//********************************************//
// Preload
//********************************************//
......@@ -351,7 +347,7 @@ int main (int argc, char** argv)
std::cout << "\n*********************\n";
// Update pressure b.c.
*endoLvVectorPtr = - pPreloadLvBC * ( i / preloadSteps );
*endoLvVectorPtr = - pPreloadLvBC * i / preloadSteps;
pLvBCVectorPtr.reset ( ( new bcVector_Type (*endoLvVectorPtr, solver.structuralOperatorPtr() -> dispFESpacePtr() -> dof().numTotalDof(), 1) ) );
solver.bcInterfacePtr() -> handler() -> modifyBC(LvFlag, *pLvBCVectorPtr);
......@@ -360,7 +356,7 @@ int main (int argc, char** argv)
solver.solveMechanics();
}
//********************************************//
// Time loop
//********************************************//
......@@ -372,7 +368,7 @@ int main (int argc, char** argv)
UInt maxiter = static_cast<UInt>( endtime / dt_activation ) ;
Real t = 0;
solver.saveSolution (0.0);
for (int k (1); k <= maxiter; k++)
{
solver.electroSolverPtr() -> registerActivationTime (*activationTimeVector, t, 0.9);
......@@ -389,7 +385,7 @@ int main (int argc, char** argv)
if ( k % saveIter == 0 )
{
// Update pressure b.c.
pLvBC = Circulation::computePressure(0.1);
pLvBC = Circulation::computePressure(1/pPreloadLvBC);
*endoLvVectorPtr = - pLvBC;
pLvBCVectorPtr.reset ( ( new bcVector_Type (*endoLvVectorPtr, solver.structuralOperatorPtr() -> dispFESpacePtr() -> dof().numTotalDof(), 1) ) );
solver.bcInterfacePtr() -> handler() -> modifyBC(LvFlag, *pLvBCVectorPtr);
......
......@@ -72,6 +72,7 @@
#include <lifev/em/solver/activation/xBridgeModels/XBridge4SM.hpp>
using namespace LifeV;
#define SolutionTestNorm -2476.4745158656560307
......@@ -134,13 +135,20 @@ Int main ( Int argc, char** argv )
//********************************************//
Real Iapp (0.0);
//********************************************//
// Setup X-Bridge Model. //
//********************************************//
XBridge4SM xb;
std::vector<Real> statesXb ( xb.xbVarVec() );
//********************************************//
// Simulation starts on t=0 and ends on t=TF. //
// The timestep is given by dt //
//********************************************//
Real TF (500.0);
Real dt (0.1);
Real TF (15500.0);
Real dt (0.001);
//********************************************//
......@@ -158,20 +166,29 @@ Int main ( Int argc, char** argv )
std::ofstream output ("output.txt");
output << 0 << "\t";
for ( int j (0); j < ionicModel.Size() - 1; j++)
for ( int j (0); j < ionicModel.Size(); j++)
{
output << states[j] << "\t";
}
output << states[ ionicModel.Size() - 1 ] << "\n";
for ( int jXb (0); jXb < xb.size(); ++jXb )
{
output << statesXb[jXb] << "\t";
}
output << "\n";
//********************************************//
// Time loop starts. //
//********************************************//
std::cout << "Time loop starts...\n";
int savestep ( 1.0 / dt );
int savestep ( 10.0 / dt );
int iter (0);
int stimulusStep (1000 / dt);
int stimulusStepBreak (1 / dt);
std::vector<Real> v (states);
......@@ -182,11 +199,12 @@ Int main ( Int argc, char** argv )
// Compute the applied current. This is a //
// simple switch. //
//********************************************//
if ( t > 50 && t < 52 )
// if ( t > 0 && t < 1 && iter % stimulusStep == 0)
if ( iter % stimulusStep == 0)
{
Iapp = 20.;
Iapp = 52.;
}
else
if( ( iter - stimulusStepBreak ) % stimulusStep == 0)
{
Iapp = 0;
}
......@@ -201,13 +219,60 @@ Int main ( Int argc, char** argv )
rhs[0] = ionicModel.computeLocalPotentialRhs (states);
ionicModel.addAppliedCurrent (rhs);
ionicModel.computeGatingVariablesWithRushLarsen ( states, dt);
states[0] = states[0] + dt * rhs[0]; //This should be divided by the membrane capacitance ionicModel.membraneCapacitance();
//ionicModel.computeGatingVariablesWithRushLarsen ( states, dt);
Real V = states[0];
Real m = states[1];
Real h = states[2];
Real j = states[3];
Real d = states[4];
Real f = states[5];
Real f2 = states[6];
Real fcass = states[7];
Real r = states[8];
Real s = states[9];
Real xr1 = states[10];
Real xr2 = states[11];
Real xs = states[12];
Real Nai = states[13];
Real Ki = states[14];
Real Cai = states[15];
Real CaSS = states[16];
Real CaSR = states[17];
Real RR = states[18];
states[18] = ionicModel.solveRR(CaSR, CaSS, RR, dt);
states[17] = ionicModel.solveCaSR(Cai, CaSR, CaSS, RR, dt);
states[16] = ionicModel.solveCaSS(Cai, CaSR, CaSS, RR, V, d, f,