Commit f8a15981 authored by Simone Rossi's avatar Simone Rossi
Browse files

Global update in the EM branch

parent 6cf20dc9
......@@ -1323,11 +1323,12 @@ public:
//@}
private:
public:
//! Set default parameters
void setParameters();
public:
//! initialization in constructor
void init();
......@@ -1349,6 +1350,20 @@ private:
*/
void init (ionicModelPtr_Type model);
void showParameters()
{
std::cout << "\nM_surfaceVolumeRatio =" << M_surfaceVolumeRatio;
std::cout << "\nM_diffusionTensor[0] =" << M_diffusionTensor[0];
std::cout << "\nM_diffusionTensor[1] =" << M_diffusionTensor[1];
std::cout << "\nM_diffusionTensor[2] =" << M_diffusionTensor[2];
std::cout << "\nM_initialTime =" << M_initialTime;
std::cout << "\nM_endTime =" << M_endTime;
std::cout << "\nM_timeStep =" << M_timeStep;
std::cout << "\nM_elementsOrder =" << M_elementsOrder;
std::cout << "\nM_lumpedMassMatrix =" << M_lumpedMassMatrix;
}
protected:
//surface to volume ratio
Real M_surfaceVolumeRatio;
......@@ -1765,8 +1780,7 @@ void ElectroETAMonodomainSolver<Mesh, IonicModel>::setupGlobalMatrix()
}
template<typename Mesh, typename IonicModel>
void ElectroETAMonodomainSolver<Mesh, IonicModel>::setupLinearSolver (
GetPot dataFile)
void ElectroETAMonodomainSolver<Mesh, IonicModel>::setupLinearSolver ( GetPot dataFile )
{
prec_Type* precRawPtr;
basePrecPtr_Type precPtr;
......
<ParameterList> <!-- LinearSolver parameters -->
<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="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="timeStep" type="double" value="0.01" /><!-- ms ^ -1-->
<Parameter name="endTime" type="double" value="10." />
<Parameter name="saveStep" type="double" value="0.01" />
<Parameter name="longitudinalDiffusion" type="double" value="3.3342" /> <!-- 1.3342k Ohm ^-1 cm ^ -1-->
<Parameter name="transversalDiffusion" type="double" value ="3.17606" /><!-- 0.17606 k Ohm ^-1 cm ^ -1-->
<Parameter name="elementsOrder" type="string" value="P1" />
<Parameter name="ionic_model" type="string" value="TenTusscher06" />
<Parameter name="ionic_model" type="string" value="AlievPanfilov" />
<Parameter name="solutionMethod" type="string" value="ICI" />
<Parameter name="subiter" type="int" value="1" />
<Parameter name="LumpedMass" type="bool" value="true"/>
<Parameter name="mesh_name" type="string" value="benchmark_05mm.mesh" />
<Parameter name="mesh_name" type="string" value="cube4.mesh" />
<Parameter name="mesh_path" type="string" value="./" />
<Parameter name="Reuse Preconditioner" type="bool" value="true"/>
<Parameter name="Max Iterations For Reuse" type="int" value="80"/>
......
......@@ -34,7 +34,8 @@ SET(HEADERS ${HEADERS}
${solver_mechanics_materials_passive_HEADERS}
${solver_mechanics_materials_activeStress_HEADERS}
${solver_activation_HEADERS}
${solver_activation_activeStressModels_HEADERS})
${solver_activation_activeStressModels_HEADERS}
${solver_activation_activeStrainModels_HEADERS})
SET(SOURCES ${SOURCES}
${solver_SOURCES}
......@@ -46,6 +47,7 @@ SET(SOURCES ${SOURCES}
${solver_mechanics_materials_activeStress_SOURCES}
${solver_activation_SOURCES}
${solver_activation_activeStressModels_SOURCES}
${solver_activation_activeStrainModels_SOURCES}
)
ADD_SUBDIRECTORY(util)
......
<ParameterList>
<!-- LinearSolver parameters -->
<Parameter name="timeStep" type="double" value="0.02" />
<Parameter name="endTime" type="double" value="1.0" />
<Parameter name="timeStep" type="double" value="0.01" />
<Parameter name="endTime" type="double" value="1." />
<Parameter name="saveStep" type="double" value="1.0" />
<Parameter name="meth" type="double" value="1.0" />
<Parameter name="emdt" type="double" value="0.1" />
<Parameter name="longitudinalDiffusion" type="double" value="3.10" />
<Parameter name="transversalDiffusion" type="double" value="3.010" />
<Parameter name="longitudinalDiffusion" type="double" value="3.0" />
<Parameter name="transversalDiffusion" type="double" value="3.00" />
<Parameter name="elementsOrder" type="string" value="P1" />
<Parameter name="ionic_model" type="string" value="AlievPanfilov" />
<Parameter name="solid_mesh_name" type="string" value="cube4.mesh" />
......
......@@ -23,6 +23,56 @@ multimesh = false
start = 0
save = 1
[electrophysiology]
[./physics]
IonicModel = AlievPanfilov
fiberDiffusion = 3.0
sheetDiffusion = 3.0
normalDiffusion = 3.0
[../time_discretization]
endtime = 50.
timestep = 0.01
[../discretization]
LumpedMass = false
[../]
[activation]
[./physics]
#ActivationModels available and relative required parameters
# ActiveStressRossi14:
# ActiveStress_Beta
# ActiveStressMu
# Tmax ( MaxActiveTension )
# ActiveStressNashPanfilov04
# kTa
# epsilon0
# ActiveStrainRossi14
# calciumIndex
# InverseViscosity
# ActiveForceCoefficient
# ChemicalThreshold
ActivationModel = ActiveStressRossi14
CalciumIndex = 0
EMActiveStrainOrthotropicParameter = -666.0
#Active Strain types Available
# Anisotropic
#Orthotropic
# TransverselyIsotropic
EMActiveStrainType = TransverselyIsotropic
[../time_discretization]
endtime = 50.
timestep = 0.01
[../]
[solid]
......@@ -38,122 +88,92 @@ gamma = 1.0
gammaf = 0.0
solidType = EMMaterial
lawType = nonlinear
#Passive material available types for EM
#PNH = Passive NeoHookean
#PNHL = Passive NeoHookean Linearized
#PIE = Passive Isoptropic Exponential
#PTIE = Passive Transversely Isoptropic Exponential
EMPassiveMaterialType = PTIE
#Passive material available types for EM
#PNH = Passive NeoHookean
#PNHL = Passive NeoHookean Linearized
#PIE = Passive Isoptropic Exponential
#PTIE = Passive Transversely Isoptropic Exponential
#PIEWS = Passive Isoptropic Exponential With Shear
#PHO = Passive Holzapfel Ogden
#POF = Passive Orthotropic Fung
#PTIF = Passive Transversely Isotropic Fung
#PV = Passive Volumetric
#ANH = Active (Strain) Neo Hookean
#SimpleActiveStress = Simple Fiber Active Stress
EMPassiveMaterialType = POF
EMActiveStressMaterialType = SimpleActiveStress #SFAS
Tmax = 497000
mu = 4960
a = 3330
af = 185350
as = 25640
afs = 4170
b = 9.242
bf = 15.972
bs = 10.446
bfs = 11.602
BulkModulus = 500000
C = 20000.0
bff = 8.0
bss = 2.0
bnn = 2.0
bfs = 4.0
bfn = 4.0
bsn = 2.0
[../boundary_conditions]
list = 'Top Bottom Front Back Left Right'
# list = 'Top InnerRing Epicardium'
[./Top]
type = Natural
flag = 400
mode = Full
component = 3
function = '[0.0, 0.0, -0.0]'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '[0.0, 0.0, 0.0]'
[../Bottom]
type = Essential
flag = 200
mode = Component
component = 2
function = '0.0'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '0.0'
[../Front]
[../Front]
type = Natural
flag = 200
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '[0.0, 0.0, 0.0]'
[../Back]
type = Essential
flag = 300
mode = Component
component = 1
function = '0.0'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '0.0'
[../Left]
type = Natural
flag = 100
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '[0.0, 0.0, 0.0]'
[../Right]
type = Essential
flag = 900
mode = Component
component = 0
function = '0.0'
# functionSD = RobinWall
[./RobinAlpha]
function = ' 5.00 * 1000' # D
[../RobinBeta]
function = ' 5.00 * 1000' # D
[../]
function = '0.0'
[../]
# [../]
[../time_discretization]
initialtime = 0.
endtime = 0.4
endtime = 50
timestep = 0.1
theta = 0.35
zeta = 0.75
......@@ -161,8 +181,8 @@ BDF_order = 2
[../space_discretization]
mesh_type = .mesh
mesh_dir = /usr/scratch/srossi/meshes/ #./ #/usr/scratch/srossi/meshes/
mesh_file = CoarseSymmCube.mesh #SymmCube2.mesh #SymmCube1.mesh #SuperCoarseSymmCube.mesh #cubeHolzapfel.mesh #CoarseSymmCube.mesh #cubeHolzapfel.mesh #cyl.mesh #StructuredCube4_test_structuralsolver.mesh #cyl.mesh
mesh_dir = /usr/scratch/srossi/meshes/
mesh_file = CoarseSymmCube.mesh
order = P1
......
#include <lifev/core/LifeV.hpp>
#include <lifev/electrophysiology/solver/ElectroETAMonodomainSolver.hpp>
#include <lifev/electrophysiology/solver/IonicModels/IonicMinimalModel.hpp>
#include <lifev/electrophysiology/solver/IonicModels/ElectroIonicModel.hpp>
#include <lifev/electrophysiology/solver/IonicModels/IonicAlievPanfilov.hpp>
//#include <lifev/electrophysiology/solver/IonicModels/IonicAlievPanfilov.hpp>
#include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>
#include <lifev/structure/solver/StructuralConstitutiveLaw.hpp>
#include <lifev/structure/solver/StructuralOperator.hpp>
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>
#include <lifev/em/solver/mechanics/EMStructuralConstitutiveLaw.hpp>
#include <lifev/em/solver/EMETAFunctors.hpp>
#include <lifev/em/solver/EMSolver.hpp>
#include <lifev/em/util/EMUtility.hpp>
#include <lifev/em/solver/electrophysiology/EMMonodomainSolver.hpp>
#include <lifev/core/filter/ExporterEnsight.hpp>
#include <lifev/core/filter/ExporterVTK.hpp>
#ifdef HAVE_HDF5
#include <lifev/core/filter/ExporterHDF5.hpp>
#endif
#include <lifev/core/filter/ExporterEmpty.hpp>
#include <lifev/eta/fem/ETFESpace.hpp>
#include <lifev/eta/expression/Integrate.hpp>
//#include <lifev/eta/fem/ETFESpace.hpp>
//#include <lifev/eta/expression/Integrate.hpp>
//
////#include <lifev/em/solver/mechanics/materials/EMMaterial.hpp>
//#include <lifev/em/solver/mechanics/materials/EMMaterial.hpp>
//#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
#include <lifev/em/solver/activation/activeStressModels/ActiveStressRossiModel14.hpp>
#include <lifev/em/solver/activation/activeStressModels/ActiveStressNashPanfilovModel04.hpp>
#include <sys/stat.h>
//#include <sys/stat.h>
using namespace LifeV;
Real Iapp (const Real& t, const Real& X, const Real& Y, const Real& Z, const ID& /*i*/)
{
Real r = 0.1;
Real t0 = 2;
if (X < r && Y < r && Z < r && t < t0)
// Real r = 0.1;
// Real t0 = 2;
// if (X < r && Y < r && Z < r && t < t0)
// {
// return 5.0;
// }
// else
// {
// return 0.0;
// }
Real pacingSite_X = 0.0;
Real pacingSite_Y = 0.0;
Real pacingSite_Z = 0.0;
Real stimulusRadius = 0.15;
Real stimulusValue = 10;
Real returnValue;
if ( std::abs ( X - pacingSite_X ) <= stimulusRadius
&&
std::abs ( Y - pacingSite_Z ) <= stimulusRadius
&&
std::abs ( Z - pacingSite_Y ) <= stimulusRadius
&&
t <= 0.5)
{
return 5.0;
returnValue = stimulusValue;
}
else
{
return 0.0;
returnValue = 0.;
}
return returnValue;
}
int main (int argc, char** argv)
{
......@@ -77,7 +89,7 @@ int main (int argc, char** argv)
typedef BCHandler bc_Type;
typedef boost::shared_ptr< bc_Type > bcPtr_Type;
typedef StructuralOperator< RegionMesh<LinearTetra> > physicalSolver_Type;
typedef EMStructuralOperator< RegionMesh<LinearTetra> > physicalSolver_Type;
typedef BCInterface3D< bc_Type, physicalSolver_Type > bcInterface_Type;
typedef boost::shared_ptr< bcInterface_Type > bcInterfacePtr_Type;
......@@ -98,13 +110,15 @@ int main (int argc, char** argv)
//===========================================================
//===========================================================
boost::shared_ptr<Epetra_Comm> comm ( new Epetra_MpiComm (MPI_COMM_WORLD) );
if ( comm->MyPID() == 0 )
{
cout << "% using MPI" << endl;
}
EMSolver<mesh_Type, monodomain_Type, ActiveStressRossiModel14> solver;
EMSolver<mesh_Type, monodomain_Type> solver(comm);
//********************************************//
// Import parameters from an xml list. Use //
......@@ -169,6 +183,8 @@ int main (int argc, char** argv)
const std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
GetPot dataFile (data_file_name);
// std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$";
// std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$";
// std::cout << "read dataFile name: " << data_file_name << ", GetPot: " << dataFile.get(0,"data000000");
......@@ -260,14 +276,17 @@ int main (int argc, char** argv)
// solidBC->handler()->bcUpdate( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
//solver.setupMechanicalBC(data_file_name, "solid", dFESpace);
ionicModelPtr_Type ionicModel (new IonicAlievPanfilov() );
std::string ionicModelType = parameterList.get ("ionic_model", "");
if ( comm->MyPID() == 0 )
if( 0 == comm->MyPID() )
{
std::cout << "Setting up EM solver ... ";
}
solver.setup (dataFile, parameterList);
if( 0 == comm->MyPID() )
{
std::cout << "\nsetup EM solver" << std::endl;
std::cout << " done!" << std::endl;
}
solver.setup (dataFile, parameterList, ionicModelType, comm);
// if ( comm->MyPID() == 0 )
// {
// std::cout << "\nsetup structural operator" << std::endl;
......@@ -283,7 +302,19 @@ int main (int argc, char** argv)
// solid.setup(dataStructure,dFESpace, dETFESpace, solver.bcInterfacePtr()->handler(), comm);
// BCh -> showMe();
solver.setupFiberVector (1, 0, 0);
if( 0 == comm->MyPID() )
{
std::cout << "Setting up anisotropy vectors ...";
}
solver.setupFiberVector (1., 0., 0.);
solver.setupSheetVector (0., 1., 0.);
if( 0 == comm->MyPID() )
{
std::cout << " done!" << std::endl;
}
// solid.EMMaterial() -> setupFiberVector(1, 0, 0);
......@@ -421,13 +452,53 @@ int main (int argc, char** argv)
//
if( 0 == comm->MyPID() )
{
std::cout << "Setting solver properties ... ";
}
solver.oneWayCoupling();
if( 0 == comm->MyPID() )
{
std::cout << " done!" << std::endl;
}
//Here we initialize the electrophysiology
//with given intial conditions
if( 0 == comm->MyPID() )
{
std::cout << "Initialize electrophysiology ... ";
}
solver.initialize();
if( 0 == comm->MyPID() )
{
std::cout << " done!" << std::endl;
}
//Here we are building the matrices
//mass matrix for mechanic and the others for electrophysiology
if( 0 == comm->MyPID() )
{
std::cout << "Buildin matrices ... ";
}
solver.buildSystem();
if( 0 == comm->MyPID() )
{
std::cout << " done!" << std::endl;
}
// solid.buildSystem (1.0);
// vectorPtr_Type rhs (new vector_Type (solid.displacement(), Unique) );
// vectorPtr_Type disp (new vector_Type (solid.displacement(), Unique) );
// vectorPtr_Type initialDisplacement (new vector_Type (solid.displacement(), Unique) );
// solid.initialize ( initialDisplacement );
solver.initialize();
//
//
......@@ -438,14 +509,17 @@ int main (int argc, char** argv)
// std::cout << "\nsetup solid exporter" << std::endl;
// }
//
// boost::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporter;
// exporter.reset ( new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, "structure" ) );
//
// exporter->setPostDir ( problemFolder );
// exporter->setMeshProcId ( localSolidMesh, comm->MyPID() );
//
//// vectorPtr_Type solidDisp ( new vector_Type (solid.displacement(), exporter->mapType() ) );
// exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "displacement", dFESpace, solid.displacementPtr(), UInt (0) );
// boost::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporter;
// exporter.reset ( new ExporterVTK<RegionMesh<LinearTetra> > ( dataFile, "structure" ) );
// //
// exporter->setPostDir ( problemFolder );
// exporter->setMeshProcId (solver.M_localMeshPtr, comm->MyPID() );
//
// //
// vectorPtr_Type solidDisp ( new vector_Type (solver.M_EMStructuralOperatorPtr->displacement(), exporter->mapType() ) );
// exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "displacement",
// solver.M_EMStructuralOperatorPtr->dispFESpacePtr(),
// solidDisp, UInt (0) );
//
// exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "Fibers", dFESpace, solid.EMMaterial() -> fiberVectorPtr(), UInt (0) );
//
......@@ -765,82 +839,125 @@ int main (int argc, char** argv)
function_Type stim = &Iapp;
Real dt = 0.02;
Real dt_activation = solver.data().activationParameter<Real>("timestep");
Real endtime = solver.data().activationParameter<Real>("endtime");
Real dt_mechanics = solver.data().solidParameter<Real>("timestep");
std::cout << "\nEndTime: " << endtime;
UInt maxiter = static_cast<UInt>( endtime / dt_activation ) ;
UInt saveIter = static_cast<UInt>( dt_mechanics / dt_activation );
Real t = 0;
// RossiModel14 activationModel(monodomain.potentialPtr()->map());
// NashPanfilovModel04 activationModel(monodomain.potentialPtr()->map());
// std::cout << "Activation Model set!\n";
//
// vectorPtr_Type activationPtr(&activationModel.activation());
//
// solid.EMMaterial()->setActivationPtr(activationPtr);