Commit 7e0c979f authored by thomaskummer's avatar thomaskummer
Browse files

euler lifev update

parent 28a109c3
......@@ -41,7 +41,7 @@ monodomain_xml_file = ParamListE.xml
[./physics]
ActivationModel = ActiveStrainRossi14 #ActiveStressRossi14
CalciumIndex = 3 # potential for activstress, 2 from minimal for active strain.
ActiveForceCoefficient = -500 #-700.0 #-7.0 #-8.0
ActiveForceCoefficient = -500 #-700.0 #-7.0 #-8.0
InverseViscosity = 0.000002
Tmax = 10000
ActiveStress_Beta = 2.279
......@@ -52,6 +52,13 @@ ActiveForceCoefficient = -500 #-700.0 #-7.0 #-8.0
endtime = 100000
timestep = 0.05
[../pathology]
strength = 0.1
infarctX = 0
infarctY = 0
infarctZ = 0
radius = 0.0
[../]
......
......@@ -143,7 +143,7 @@ EMData::setupActivationParameters(GetPot& dataFile, const std::string& section)
double timestep = dataFile ( ( section + "/time_discretization/timestep" ).data(), 0.1 );
M_activationParametersList.set ("timestep", timestep);
}
......
......@@ -628,7 +628,7 @@ EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
M_bcInterfacePtr->handler(),
M_commPtr);
M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);
M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data, dataFile);
M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
// M_wtePassive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "passive");
......
......@@ -517,7 +517,7 @@ public:
void setParameters(EMData& data);
void setParameters(EMData& data, GetPot& dataFile);
void showMaterialParameters()
......@@ -568,6 +568,11 @@ protected:
vectorPtr_Type M_FirstPiolaKStress;
VectorSmall<3> M_PathologyCenter;
Real M_PathologyRadius;
Real M_PathologyStrength;
class OrthonormalizeVector
{
......@@ -960,7 +965,7 @@ protected:
auto X = vectors[2];
auto gf = g;
if ( X [2] < -3 ) gf *= 0.1;
if ( (X - M_Pathology).norm <= M_PathologyRadius ) gf *= M_PathologyStrength;
auto gn = 4 * gf;
auto gs = 1 / ( (gf + 1) * (gn + 1) ) - 1;
......@@ -1346,7 +1351,7 @@ protected:
return_Type operator() (const Real& g, const VectorSmall<3>& X)
{
auto gf = g;
if ( X[2] < -3 ) gf *= 0.1;
if ( (X - M_Pathology).norm <= M_PathologyRadius ) gf *= M_PathologyStrength;
return gf;
}
......@@ -1639,13 +1644,22 @@ void EMStructuralConstitutiveLaw<MeshType>::computeStiffness ( const vector_Type
}
template <typename MeshType>
void EMStructuralConstitutiveLaw<MeshType>::setParameters(EMData& data)
void EMStructuralConstitutiveLaw<MeshType>::setParameters(EMData& data, GetPot& dataFile)
{
if (M_activeStressMaterialPtr)
{
M_activeStressMaterialPtr-> setParameters(data);
}
M_PathologyCenter[0] = dataFile("activation/pathology/infarctX", 0.0);
M_PathologyCenter[1] = dataFile("activation/pathology/infarctY", 0.0);
M_PathologyCenter[2] = dataFile("activation/pathology/infarctZ", 0.0);
M_PathologyRadius = dataFile("activation/pathology/radius", 0.0);
M_PathologyStrength = dataFile("activation/pathology/strength", 1.0);
if (M_passiveMaterialPtr)
{
M_passiveMaterialPtr-> setParameters(data);
......
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