Commit 9eaaa6ab authored by thomaskummer's avatar thomaskummer
Browse files

heart geometry with patches

parent 3b3047a4
......@@ -279,7 +279,7 @@ int main (int argc, char** argv)
//============================================//
// Setup exporter for EMSolver
// Setup exporters for EMSolver
//============================================//
displayer.leaderPrint ("\nSetting up exporters ... ");
......@@ -289,20 +289,6 @@ int main (int argc, char** argv)
displayer.leaderPrint ("\ndone!");
//============================================//
// Setup vector & exporter for activation time
//============================================//
vectorPtr_Type activationTimeVector ( new vector_Type ( solver.electroSolverPtr()->potentialPtr() -> map() ) );
*activationTimeVector = -1.0;
ExporterHDF5< RegionMesh <LinearTetra> > activationTimeExporter;
activationTimeExporter.setMeshProcId (solver.localMeshPtr(), solver.comm()->MyPID() );
activationTimeExporter.addVariable (ExporterData<mesh_Type>::ScalarField, "Activation Time", solver.electroSolverPtr()->feSpacePtr(), activationTimeVector, UInt (0) );
activationTimeExporter.setPrefix ("ActivationTime");
activationTimeExporter.setPostDir (problemFolder);
//============================================//
// Electric stimulus function
//============================================//
......@@ -497,7 +483,6 @@ int main (int argc, char** argv)
}
// Set time exporter time index
activationTimeExporter.setTimeIndex(restartInputStr + 1);
solver.setTimeIndex(restartInputStr + 1);
// Load restart solutions from output files
......@@ -511,7 +496,7 @@ int main (int argc, char** argv)
}
ElectrophysiologyUtility::importScalarField (solver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
ElectrophysiologyUtility::importScalarField (activationTimeVector, "ActivationTime" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "ActivationTimeSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
circulationSolver.restartFromFile ( restartDir + "solution.dat" , nIter );
......@@ -536,7 +521,6 @@ int main (int argc, char** argv)
};
solver.saveSolution (-1.0);
activationTimeExporter.postProcess(-1.0);
LifeChrono chronoPreload;
chronoPreload.start();
......@@ -556,8 +540,6 @@ int main (int argc, char** argv)
// Solve mechanics
solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
solver.solveMechanics();
//solver.saveSolution ( i );
//activationTimeExporter.postProcess( i );
}
if ( 0 == comm->MyPID() )
......@@ -598,7 +580,6 @@ int main (int argc, char** argv)
if ( ! restart )
{
solver.saveSolution (t);
activationTimeExporter.postProcess(t);
circulationSolver.exportSolution( circulationOutputFile );
}
......@@ -616,13 +597,7 @@ int main (int argc, char** argv)
//============================================//
// Solve electrophysiology and activation
//============================================//
// if ( fmod(t, 800.) < 4 && fmod(t, 800.) > 2)
// {
// ElectrophysiologyUtility::setValueOnBoundary ( * (solver.electroSolverPtr()->potentialPtr() ), solver.fullMeshPtr(), 1.0, lvendo );
// }
solver.electroSolverPtr() -> registerActivationTime (*activationTimeVector, t, 0.9);
solver.solveElectrophysiology (stim, t);
solver.solveActivation (dt_activation);
......@@ -801,11 +776,7 @@ int main (int argc, char** argv)
// Export FE-solution
//============================================//
bool save ( std::abs(std::remainder(t, 10.)) < 0.01 );
if ( save )
{
solver.saveSolution(t, restart);
activationTimeExporter.postProcess(t);
}
if ( save ) solver.saveSolution(t, restart);
}
}
......@@ -814,9 +785,7 @@ int main (int argc, char** argv)
//============================================//
// Close all exporters
//============================================//
solver.closeExporters();
activationTimeExporter.closeFile();
#ifdef HAVE_MPI
......
......@@ -169,6 +169,11 @@ public:
{
EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
void setupActivationTimeExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationTimeSolution" )
{
EMUtility::setupExporter<Mesh> (*M_activationTimeExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
{
......@@ -181,6 +186,7 @@ public:
void setupExporters (std::string problemFolder = "./",
std::string electroFileName = "ElectroSolution",
std::string activationFileName = "ActivationSolution",
std::string activationTimeFileName = "ActivationTimeSolution",
std::string mechanicsFileName = "MechanicalSolution");
void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
......@@ -324,6 +330,11 @@ public:
{
return M_activationModelPtr;
}
vectorPtr_Type activationTimePtr()
{
return M_activationTimePtr;
}
void saveSolution (Real time, const bool& restart = 0);
......@@ -595,6 +606,7 @@ void
EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
std::string electroFileName,
std::string activationFileName,
std::string activationTimeFileName,
std::string mechanicsFileName)
{
if (M_commPtr -> MyPID() == 0)
......@@ -621,28 +633,21 @@ EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
UInt (0) );
M_mechanicsExporterPtr.reset (new exporter_Type() );
// // Activation time
// M_activationTimeExporterPtr.reset (new exporter_Type() );
// setupActivationExporter (problemFolder, activationFileName );
//
// vectorPtr_Type M_activationTime ( new vector_Type ( M_electroSolverPtr()->potentialPtr() -> map() ) );
// *activationTimeVector = -1.0;
//
// M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
// "Activation Time",
// M_electroSolverPtr -> feSpacePtr(),
// M_activationModelPtr -> fiberActivationPtr(),
// UInt (0) );
// M_mechanicsExporterPtr.reset (new exporter_Type() );
//
//
// activationTimeExporter.addVariable (ExporterData<mesh_Type>::ScalarField, "Activation Time", solver.electroSolverPtr()->feSpacePtr(), activationTimeVector, UInt (0) );
//
// activationTimeExporter.setMeshProcId (solver.localMeshPtr(), solver.comm()->MyPID() );
// activationTimeExporter.setPrefix ("ActivationTime");
// activationTimeExporter.setPostDir (problemFolder);
// Activation time
M_activationTimeExporterPtr.reset (new exporter_Type() );
setupActivationTimeExporter (problemFolder, activationTimeFileName );
M_activationTimePtr.reset (new vector_Type ( M_electroSolverPtr->potentialPtr() -> map() ) );
*M_activationTimePtr = -1.0;
M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Activation Time",
M_electroSolverPtr -> feSpacePtr(),
M_activationTimePtr,
UInt (0) );
// Mechanics
M_mechanicsExporterPtr.reset (new exporter_Type() );
setupMechanicsExporter (problemFolder, mechanicsFileName);
M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"displacement",
......@@ -667,6 +672,7 @@ EMSolver<Mesh, ElectroSolver>::setTimeIndex (const UInt& time)
{
M_electroExporterPtr -> setTimeIndex (time);
M_activationExporterPtr -> setTimeIndex (time);
M_activationTimeExporterPtr -> setTimeIndex (time);
M_mechanicsExporterPtr -> setTimeIndex (time);
}
......@@ -678,6 +684,7 @@ EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
//if(M_activationExporterPtr) std::cout << "\nActivation exporter available.";
//if(M_activationModelPtr -> fiberActivationPtr()) std::cout << "\nFiber activation exporter available.";
M_activationExporterPtr -> postProcess (time);//, restart);
M_activationTimeExporterPtr -> postProcess (time);
//if(M_mechanicsExporterPtr) std::cout << "\nMechanics exporter available.";
M_mechanicsExporterPtr -> postProcess (time);//, restart);
}
......@@ -688,6 +695,7 @@ EMSolver<Mesh, ElectroSolver>::closeExporters()
{
M_electroExporterPtr -> closeFile();
M_activationExporterPtr -> closeFile();
M_activationTimeExporterPtr -> closeFile();
M_mechanicsExporterPtr -> closeFile();
}
......@@ -728,6 +736,7 @@ EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus,
M_electroSolverPtr -> solveOneStepGatingVariablesFE();
// M_electroSolverPtr -> solveOneStepGatingVariablesRL();
M_electroSolverPtr -> solveOneICIStep();
M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
}
template<typename Mesh , typename ElectroSolver>
......
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