Commit 37edab4e authored by thomaskummer's avatar thomaskummer
Browse files

cervin lifev update

parent 723519f6
......@@ -269,7 +269,7 @@ monodomain_xml_file = ParamListE.xml
[../newton]
maxiter = 50
reltol = 1.e-7
abstol = 0.8e-6
abstol = 0.9e-6
NonLinearLineSearch = 0
[../solver]
......
......@@ -19,34 +19,71 @@
namespace LifeV
{
template <class EmSolver>
class HeartSolver {
public:
HeartSolver(EmSolver& emSolver, Circulation& circulationSolver) :
M_emSolver (emSolver),
M_circulationSolver (circulationSolver)
{};
class HeartData
{
public:
HeartSolver(const HeartSolver& heartSolver) :
M_emSolver (heartSolver.M_emSolver),
M_circulationSolver (heartSolver.M_circulationSolver)
{};
HeartData(const GetPot& datafile)
{}
virtual ~HeartData() {};
//UserQueues(const UserQueues&) = delete;
protected:
// const Real dt_activation = solver.data().electroParameter<Real>("timestep");
// const Real dt_loadstep = dataFile ( "solid/time_discretization/dt_loadstep", 1.0 );
// const Real activationLimit_loadstep = dataFile ( "solid/time_discretization/activation_limit_loadstep", 0.0 );
// const Real dt_mechanics = solver.data().solidParameter<Real>("timestep");
// const Real dt_save = dataFile ( "exporter/save", 10. );
// const Real endtime = solver.data().electroParameter<Real>("endtime");
// const UInt mechanicsLoadstepIter = static_cast<UInt>( dt_loadstep / dt_activation );
// const UInt mechanicsCouplingIter = static_cast<UInt>( dt_mechanics / dt_activation );
// const UInt maxiter = static_cast<UInt>( endtime / dt_activation ) ;
//
// const Real pPerturbationFe = dataFile ( "solid/coupling/pPerturbationFe", 1e-2 );
// const Real pPerturbationCirc = dataFile ( "solid/coupling/pPerturbationCirc", 1e-3 );
// const Real couplingError = dataFile ( "solid/coupling/couplingError", 1e-6 );
// const UInt couplingJFeSubIter = dataFile ( "solid/coupling/couplingJFeSubIter", 1 );
// const UInt couplingJFeSubStart = dataFile ( "solid/coupling/couplingJFeSubStart", 1 );
// const UInt couplingJFeIter = dataFile ( "solid/coupling/couplingJFeIter", 1 );
//
// const Real dpMax = dataFile ( "solid/coupling/dpMax", 0.1 );
//
// std::vector<std::vector<std::string> > bcNames { { "lv" , "p" } , { "rv" , "p" } };
// std::vector<double> bcValues { p ( "lv" ) , p ( "rv") };
// std::vector<double> bcValuesPre ( bcValues );
//
// VectorSmall<4> ABdplv, ABdprv, ABcoef;
// ABcoef (0) = 55/24; ABcoef (1) = -59/24; ABcoef (2) = 37/24; ABcoef (3) = -3/8;
//
// VectorSmall<2> VCirc, VCircNew, VCircPert, VFe, VFeNew, VFePert, R, dp;
// MatrixSmall<2,2> JFe, JCirc, JR;
//
// UInt iter (0);
// Real t (0);
};
virtual ~HeartSolver() {};
template <class EmSolver>
class HeartSolver {
public:
HeartSolver(EmSolver& emSolver, Circulation& circulationSolver) :
M_emSolver (emSolver),
M_circulationSolver (circulationSolver)
{}
virtual ~HeartSolver() {}
void preloadHeart(const VectorSmall<2>& endocardiaBC);
// void setupHeartData(const GetPot& datafile)
// {
// M_heartData = HeartData(datafile);
// }
......@@ -59,15 +96,14 @@ protected:
EmSolver M_emSolver;
Circulation M_circulationSolver;
//HeartData M_heartData;
VectorSmall<2> M_pressure;
VectorSmall<2> M_volume;
Real patchForce (const Real& t, const Real& Tmax, const Real& tmax, const Real& tduration) const
{
......
......@@ -247,7 +247,15 @@ public:
if ( vertex1Idx < dofh.size() ) u[1] = M_u(vertex1Idx);
if ( vertex2Idx < dofh.size() ) u[2] = M_u(vertex2Idx);
element -> initRestart ( u , M_time );
std::vector<double> uPrev0 {M_uPrev0(elementIdx) , 0.0, 0.0 };
if ( vertex1Idx < dofh.size() ) uPrev0[1] = M_uPrev0(vertex1Idx);
if ( vertex2Idx < dofh.size() ) uPrev0[2] = M_uPrev0(vertex2Idx);
std::vector<double> uPrev1 {M_uPrev1(elementIdx) , 0.0, 0.0 };
if ( vertex1Idx < dofh.size() ) uPrev1[1] = M_uPrev1(vertex1Idx);
if ( vertex2Idx < dofh.size() ) uPrev1[2] = M_uPrev1(vertex2Idx);
element -> initRestart ( u , uPrev0, uPrev1, M_time );
}
}
......
......@@ -32,7 +32,7 @@ public:
const std::string node(const unsigned int idx) const {return M_nodes[idx];} // todo: return pointer to Vertex object
const std::vector<std::string> nodes() const {return M_nodes;}
const double init() const {return M_init;}
virtual void initRestart(const std::vector<double>& u, const double& time){}
virtual void initRestart(const std::vector<double>& u, const std::vector<double>& uPrev0, const std::vector<double>& uPrev1, const double& time){}
virtual const double rhs(const std::vector<double>& u, const double& time) = 0;
virtual const double lhs(const unsigned int& idx, const std::vector<double>& u, const double& time) = 0; // u containing (Q p1 p2 R)
......@@ -312,12 +312,12 @@ public:
}
}
virtual void initRestart(const std::vector<double>& u, const double& time)
virtual void initRestart(const std::vector<double>& u, const std::vector<double>& uPrev0, const std::vector<double>& uPrev1, const double& time)
{
M_time = time;
M_D = M_param[0] * u[0] / ( u[1] - u[2] );
M_Dprev0 = M_D;
M_Dprev1 = M_D;
M_Dprev0 = M_param[0] * uPrev0[0] / ( uPrev0[1] - uPrev0[2] );
M_Dprev1 = M_param[0] * uPrev1[0] / ( uPrev1[1] - uPrev1[2] );
}
......@@ -363,12 +363,12 @@ public:
};
}
virtual void initRestart(const std::vector<double>& u, const double& time)
virtual void initRestart(const std::vector<double>& u, const std::vector<double>& uPrev0, const std::vector<double>& uPrev1, const double& time)
{
M_time = time;
M_D = M_param[0] * u[0] * u[1] / ( u[1] - u[2] );
M_Dprev0 = M_D;
M_Dprev1 = M_D;
M_Dprev0 = M_param[0] * uPrev0[0] / ( uPrev0[1] - uPrev0[2] );
M_Dprev1 = M_param[0] * uPrev1[0] / ( uPrev1[1] - uPrev1[2] );
}
};
......
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