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

New EM solver

parent 00c8b290
......@@ -44,8 +44,16 @@
#ifndef _EMSOLVER_H_
#define _EMSOLVER_H_
#include <lifev/core/mesh/MeshLoadingUtility.hpp>
#include <lifev/em/solver/electrophysiology/EMMonodomainSolver.hpp>
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>
#include <lifev/em/solver/mechanics/EMStructuralConstitutiveLaw.hpp>
#include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>
#include <lifev/em/solver/activation/activeStressModels/ActiveStressActivation.hpp>
//
//#include <lifev/em/solver/EMStructuralOperator.hpp>
//#include <lifev/em/solver/EMGeneralizedActiveHolzapfelOgdenMaterial.hpp>
......@@ -57,7 +65,7 @@
////#include <lifev/core/interpolation/RBFscalar.hpp>
//#include <lifev/core/interpolation/RBFvectorial.hpp>
//
//#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
//
//
//#include <lifev/em/solver/EMEvaluate.hpp>
......@@ -67,47 +75,499 @@ namespace LifeV
//! EMSolver - Class featuring the solution of the electromechanical problem with monodomain equation
template<typename Mesh , typename IonicModel, typename ActivationModel>
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
class EMSolver
{
public:
typedef Mesh mesh_Type;
typedef boost::shared_ptr<mesh_Type> meshPtr_Type;
typedef Epetra_Comm comm_Type;
typedef boost::shared_ptr<Epetra_Comm> commPtr_Type;
typedef VectorEpetra vector_Type;
typedef StructuralConstitutiveLawData structureData_Type;
typedef boost::shared_ptr<structureData_Type> structureDataPtr_Type;
typedef EMStructuralOperator< mesh_Type > structuralOperator_Type;
typedef boost::shared_ptr< structuralOperator_Type > structuralOperatorPtr_Type;
typedef BCHandler bc_Type;
typedef boost::shared_ptr< bc_Type > bcPtr_Type;
typedef StructuralOperator< mesh_Type > physicalSolver_Type;
typedef BCInterface3D< bc_Type, physicalSolver_Type > bcInterface_Type;
EMSolver();
typedef boost::shared_ptr< bcInterface_Type > bcInterfacePtr_Type;
typedef boost::shared_ptr<ActivationModel> activationModelPtr_Type;
typedef ElectroSolver electroSolver_Type;
typedef boost::shared_ptr<electroSolver_Type> electroSolverPtr_Type;
typedef ElectroIonicModel ionicModel_Type;
typedef boost::shared_ptr<ionicModel_Type> ionicModelPtr_Type;
typedef ExporterHDF5< Mesh > exporter_Type;
typedef boost::shared_ptr<ExporterHDF5< Mesh > > exporterPtr_Type;
typedef FESpace< RegionMesh<LinearTetra>, MapEpetra > solidFESpace_Type;
typedef boost::shared_ptr<solidFESpace_Type> solidFESpacePtr_Type;
typedef ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > solidETFESpace_Type;
typedef boost::shared_ptr<solidETFESpace_Type> solidETFESpacePtr_Type;
typedef boost::function < Real (const Real& t,
const Real & x,
const Real & y,
const Real& z,
const ID& /*i*/ ) > function_Type;
EMSolver();
EMSolver(const EMSolver& solver);
EMSolver(EMSolver&& solver);
public:
EMMonodomainSolver<Mesh, IonicModel> M_monodomain;
ActivationModel M_activationModel;
inline void loadMesh(std::string meshName, std::string meshPath)
{
std::cout << "EMS - Loading mesh\n";
MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
}
inline void setupElectroExporter( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
{
M_electroSolverPtr -> setupExporter(*M_electroExporterPtr, outputFileName, problemFolder);
}
inline void setupActivationExporter( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
{
EMUtility::setupExporter<Mesh>(*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
inline void setupMechanicsExporter( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
{
if(M_mechanicsExporterPtr)
EMUtility::setupExporter<Mesh>(*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
void setupExporters(std::string problemFolder = "./",
std::string electroFileName = "ElectroSolution",
std::string activationFileName = "ActivationSolution",
std::string mechanicsFileName = "MechanicalSolution");
void setupElectroSolver( GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize);
void setupElectroSolver( GetPot& dataFile, Teuchos::ParameterList& list, std::string ionicModelName);
void setupMechanicalSolver( GetPot& dataFile);
void setupMechanicalBC(std::string data_file_name,
std::string section,
solidFESpacePtr_Type dFESpace);
inline void setupActivation(const MapEpetra& map)
{
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - setting up activation solver\n";
}
M_activationModelPtr.reset( new ActivationModel(map) );
}
void setup(GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize, commPtr_Type commPtr);
void setup(GetPot& dataFile, Teuchos::ParameterList& list, std::string ionicModelName, commPtr_Type commPtr);
inline void buildMechanicalSystem()
{
M_EMStructuralOperatorPtr -> buildSystem (1.0);
}
inline void buildElectroSystem()
{
M_electroSolverPtr -> setupMatrices();
}
inline void buildSystem()
{
buildMechanicalSystem(); buildElectroSystem();
}
inline void initializeElectroVariables()
{
M_electroSolverPtr -> setInitialConditions();
}
inline void initialize()
{
initializeElectroVariables();
}
inline bcInterfacePtr_Type bcInterfacePtr()
{
return M_bcInterfacePtr;
}
inline void setupMechanicalFiberVector( Real fx, Real fy, Real fz )
{
M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector( fx, fy, fz);
}
inline void setupElectroFiberVector( VectorSmall<3>& fibers)
{
M_electroSolverPtr -> setupFibers (fibers);
}
inline void setupFiberVector( Real fx, Real fy, Real fz )
{
VectorSmall<3> f; f[0]=fx; f[1]=fy; f[2]=fz;
setupElectroFiberVector(f);
M_EMStructuralOperatorPtr -> EMMaterial() -> setFiberVectorPtr ( M_electroSolverPtr -> fiberPtr() );
}
inline electroSolverPtr_Type electroSolverPtr()
{
return M_electroSolverPtr;
}
inline structuralOperator_Type structuralOperatorPtr()
{
return M_EMStructuralOperatorPtr;
}
inline activationModelPtr_Type activationModelPtr()
{
return M_activationModelPtr;
}
void saveSolution(Real time);
void closeExporters();
void oneWayCoupling();
void twoWayCoupling();
inline void setAppliedCurrent(function_Type& stimulus, Real time = 0.0)
{
M_electroSolverPtr -> setAppliedCurrentFromFunction(stimulus, time);
}
inline void solveMechanics()
{
M_EMStructuralOperatorPtr -> iterate( M_bcInterfacePtr -> handler() );
}
inline void solveElectrophysiology(function_Type& stimulus, Real time = 0.0);
inline void solveActivation(UInt index, Real dt);
// inline bcInterface_Type bcInterface()
// {
// return *M_bcInterfacePtr;
// }
protected:
electroSolverPtr_Type M_electroSolverPtr;
activationModelPtr_Type M_activationModelPtr;
bcInterfacePtr_Type M_bcInterfacePtr;
structuralOperatorPtr_Type M_EMStructuralOperatorPtr;
exporterPtr_Type M_electroExporterPtr;
exporterPtr_Type M_activationExporterPtr;
exporterPtr_Type M_mechanicsExporterPtr;
meshPtr_Type M_localMeshPtr;
meshPtr_Type M_fullMeshPtr;
vectorPtr_Type M_activationTime;
bool M_oneWayCoupling;
commPtr_Type M_commPtr;
};
/////////////////////
// CONSTRUCTORS
template<typename Mesh , typename IonicModel, typename ActivationModel>
EMSolver<Mesh, IonicModel, ActivationModel>::EMSolver()
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
EMSolver<Mesh, ElectroSolver, ActivationModel>::EMSolver() :
M_electroSolverPtr ( ),
M_activationModelPtr ( ),
M_bcInterfacePtr ( ),
M_EMStructuralOperatorPtr(),
M_electroExporterPtr ( ),
M_activationExporterPtr ( ),
M_mechanicsExporterPtr ( ),
M_localMeshPtr ( ),
M_fullMeshPtr ( ),
M_activationTime ( ),
M_oneWayCoupling (true),
M_commPtr ( )
{
}
/////////////////////
// COPY CONSTRUCTORS
template<typename Mesh , typename IonicModel, typename ActivationModel>
EMSolver<Mesh, IonicModel, ActivationModel>::EMSolver(const EMSolver& solver) :
M_monodomain(solver.M_monodomain),
M_activationModel(solver.M_activationModel)
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
EMSolver<Mesh, ElectroSolver, ActivationModel>::EMSolver(const EMSolver& solver) :
M_electroSolverPtr(solver.M_electroSolverPtr),
M_activationModelPtr(solver.M_activationModelPtr),
M_bcInterfacePtr ( solver.M_bcInterfacePtr),
M_EMStructuralOperatorPtr(solver.M_EMStructuralOperatorPtr),
M_electroExporterPtr ( solver.M_electroExporterPtr),
M_activationExporterPtr ( solver.M_activationExporterPtr),
M_mechanicsExporterPtr ( solver.M_mechanicsExporterPtr),
M_localMeshPtr ( solver.M_localMeshPtr),
M_fullMeshPtr ( solver.M_fullMeshPtr),
M_activationTime ( solver.M_activationTime),
M_oneWayCoupling ( solver.M_oneWayCoupling),
M_commPtr ( solver.M_commPtr)
{
}
/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setup( GetPot& dataFile,
Teuchos::ParameterList& list,
short int ionicModelSize,
commPtr_Type commPtr)
{
M_commPtr = commPtr;
setupElectroSolver( dataFile, list, ionicModelSize );
setupMechanicalSolver( dataFile );
setupActivation( M_electroSolverPtr -> potentialPtr() ->map() );
}
/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setup( GetPot& dataFile,
Teuchos::ParameterList& list,
std::string ionicModelName,
commPtr_Type commPtr)
{
M_commPtr = commPtr;
setupElectroSolver( dataFile, list, ionicModelName );
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - electro solver setup done! ";
}
setupMechanicalSolver( dataFile );
setupActivation( M_electroSolverPtr -> potentialPtr() ->map() );
}
/////////////////////
// MOVE CONSTRUCTORS
template<typename Mesh , typename IonicModel, typename ActivationModel>
EMSolver<Mesh, IonicModel, ActivationModel>::EMSolver(EMSolver&& solver) :
M_monodomain(solver.M_monodomain)
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupElectroSolver( GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize)
{
}
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupElectroSolver( GetPot& dataFile, Teuchos::ParameterList& list, std::string ionicModelName)
{
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - creating ionic model ";
}
ionicModelPtr_Type ionicModelPtr;
ionicModelPtr.reset(ionicModel_Type::IonicModelFactory::instance().createObject ( ionicModelName ));
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - setting up electrophysiology solver ";
}
M_electroSolverPtr.reset( new ElectroSolver() );
M_electroSolverPtr -> setIonicModelPtr(ionicModelPtr);
if(M_localMeshPtr)
{
M_electroSolverPtr -> setLocalMeshPtr(M_localMeshPtr);
if(M_fullMeshPtr)
{
M_electroSolverPtr -> setFullMeshPtr(M_fullMeshPtr);
}
M_electroSolverPtr -> setup(dataFile, ionicModelPtr->Size() );
}
M_electroSolverPtr -> setParameters ( list );
if(M_commPtr -> MyPID()==0)
{
std::cout << "... `Done\n";
}
}
/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupMechanicalSolver( GetPot& dataFile)
{
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - setting up mechanical solver\n";
}
boost::shared_ptr<StructuralConstitutiveLawData> dataStructure (new StructuralConstitutiveLawData( ) );
dataStructure->setup (dataFile);
std::string dOrder = dataFile ( "solid/space_discretization/order", "P1");
solidFESpacePtr_Type dFESpace ( new solidFESpace_Type (M_localMeshPtr, dOrder, 3, M_commPtr) );
solidETFESpacePtr_Type dETFESpace ( new solidETFESpace_Type ( M_localMeshPtr,
& (dFESpace->refFE() ),
& (dFESpace->fe().geoMap() ),
M_commPtr) );
std::string data_file_name = dataFile.get(0,"NO_DATA_FILENAME_FOUND");
setupMechanicalBC(data_file_name, "solid", dFESpace);
M_EMStructuralOperatorPtr.reset(new structuralOperator_Type() );
M_EMStructuralOperatorPtr->setup( dataStructure,
dFESpace,
dETFESpace,
M_bcInterfacePtr->handler(),
M_commPtr);
M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
}
/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupMechanicalBC(std::string data_file_name,
std::string section,
solidFESpacePtr_Type dFESpace)
{
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - setting up bc interface\n";
}
M_bcInterfacePtr.reset(new bcInterface_Type() );
M_bcInterfacePtr->createHandler();
M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
M_bcInterfacePtr->handler()->bcUpdate( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
}
/////////////////////
//Setup exporters
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupExporters(std::string problemFolder,
std::string electroFileName,
std::string activationFileName,
std::string mechanicsFileName)
{
if(M_commPtr -> MyPID()==0)
{
std::cout << "EMS - setting up exporters\n";
}
M_electroExporterPtr.reset(new exporter_Type() );
setupElectroExporter(problemFolder, electroFileName);
M_activationExporterPtr.reset(new exporter_Type() );
setupActivationExporter(problemFolder, activationFileName );
M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Activation",
M_electroSolverPtr -> feSpacePtr(),
M_activationModelPtr -> activationPtr(),
UInt (0) );
M_mechanicsExporterPtr.reset(new exporter_Type() );
setupMechanicsExporter(problemFolder, mechanicsFileName);
M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"displacement",
M_EMStructuralOperatorPtr -> dispFESpacePtr(),
M_EMStructuralOperatorPtr -> displacementPtr(),
UInt (0) );
M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"fibers",
M_EMStructuralOperatorPtr -> dispFESpacePtr(),
M_EMStructuralOperatorPtr -> EMMaterial() -> fiberVectorPtr(),
UInt (0) );
}
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::saveSolution(Real time)
{
M_electroExporterPtr -> postProcess(time);
M_activationExporterPtr -> postProcess(time);
M_mechanicsExporterPtr -> postProcess(time);
}
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::closeExporters()
{
M_electroExporterPtr -> closeFile();
M_activationExporterPtr -> closeFile();
M_mechanicsExporterPtr -> closeFile();
}
///////////////////////////////
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::oneWayCoupling()
{
M_electroSolverPtr -> setMechanicsModifiesConductivity(false);
M_electroSolverPtr -> displacementPtr().reset();
M_EMStructuralOperatorPtr -> EMMaterial() -> setActivationPtr( M_activationModelPtr -> activationPtr() );
}
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::twoWayCoupling()
{
M_electroSolverPtr -> setMechanicsModifiesConductivity(true);
M_electroSolverPtr -> setDisplacementPtr( M_EMStructuralOperatorPtr -> displacementPtr() );
M_EMStructuralOperatorPtr -> EMMaterial() -> setActivationPtr( M_activationModelPtr -> activationPtr() );
}
////////////////////////////
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::solveElectrophysiology(function_Type& stimulus, Real time )
{
setAppliedCurrent ( stimulus, time );
M_electroSolverPtr -> solveOneStepGatingVariablesFE();
M_electroSolverPtr -> solveOneICIStep();
}
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::solveActivation(UInt index, Real dt)
{
M_activationModelPtr -> solveModel( *( (M_electroSolverPtr -> globalSolution())[index] ), dt);
}
} // namespace LifeV
......
Supports Markdown
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