Commit 3f06907e authored by thomaskummer's avatar thomaskummer
Browse files

euler lifev update

parent cd9d17d4
......@@ -126,6 +126,11 @@ public:
solver.bcInterfacePtr()->handler()->modifyBC(m_patchFlag, *m_patchDispBCPtr);
}
vectorPtr_Type patchDisplacement()
{
return m_patchDispPtr;
}
protected:
......@@ -194,7 +199,7 @@ public:
~EssentialPatchBCHandler(){}
addPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver)
void addPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver)
{
for ( UInt i (0) ; i < m_patchNumber ; ++i )
{
......@@ -206,7 +211,7 @@ public:
}
}
applyPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver)
void applyPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver)
{
for (auto& patch : m_patchBCPtrVec)
{
......@@ -214,21 +219,41 @@ public:
}
}
modifyPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver, const Real& time)
void modifyPatchBC(EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver, const Real& time)
{
for (auto& patch : m_patchBCPtrVec)
{
patch->modifyPatchBC(solver, time);
}
updatePatchDisplacementSum();
}
vectorPtr_Type patchDisplacement()
{
return m_patchDisplacementSum;
}
private:
void updatePatchDisplacementSum()
{
m_patchDisplacementSum *= 0.0;
for (auto& patch : m_patchBCPtrVec)
{
m_patchDisplacementSum += patch->patchDisplacement();
}
}
const std::string m_patchListName;
const GetPot& m_dataFile;
const int m_patchNumber;
vectorPtr_Type m_patchDisplacementSum;
std::vector<EssentialPatchBC*> m_patchBCPtrVec;
};
......
......@@ -58,25 +58,25 @@ protected:
UInt jGID = vectorField->blockMap().GID (j + nCompLocalDof);
UInt kGID = vectorField->blockMap().GID (j + 2 * nCompLocalDof);
// Vector3D coord;
//
// coord(0) = dFeSpace->mesh()->point(iGID).x() + (*m_dispPtr)[iGID];
// coord(1) = dFeSpace->mesh()->point(iGID).y() + (*m_dispPtr)[jGID];
// coord(2) = dFeSpace->mesh()->point(iGID).z() + (*m_dispPtr)[kGID];
//
// // Radial and axial distance to center line
// Vector3D currentPatchCenter = m_Center + activationFunction(time) * direction;
// auto radialDistance = ( (coord - m_Center).cross(coord - currentPatchCenter) ).norm() / (m_Center - currentPatchCenter).norm();
// auto axialDistance = (coord - currentPatchCenter).dot(direction) * direction;
//
// // If coordiantes inside or outside of a certain radius
// auto displacement = (m_EdgeDispFactor * disp - disp) * dispDistributionWeight(coord) + disp;
Vector3D coord;
coord(0) = dFeSpace->mesh()->point(iGID).x() + (*m_dispPtr)[iGID];
coord(1) = dFeSpace->mesh()->point(iGID).y() + (*m_dispPtr)[jGID];
coord(2) = dFeSpace->mesh()->point(iGID).z() + (*m_dispPtr)[kGID];
// Radial and axial distance to center line
Vector3D currentPatchCenter = m_Center + activationFunction(time) * direction;
auto radialDistance = ( (coord - m_Center).cross(coord - currentPatchCenter) ).norm() / (m_Center - currentPatchCenter).norm();
auto axialDistance = (coord - currentPatchCenter).dot(direction) * direction;
// If coordiantes inside or outside of a certain radius
auto displacement = (m_EdgeDispFactor * disp - disp) * dispDistributionWeight(coord) + disp;
// If patch inside or outside the structure
// Scale the direction vector
auto displacementVec = direction; //* displacement;
auto displacementVec = direction * displacement;
(*vectorField)[iGID] = displacementVec[0];
(*vectorField)[jGID] = displacementVec[1];
(*vectorField)[kGID] = displacementVec[2];
......@@ -124,7 +124,7 @@ protected:
Real dispWeight (0);
if ( nodeInsideEllipsoid(ellipsoidCoord) )
{
dispWeight = = ellipsoidFuncEval(ellipsoidCoord);
dispWeight = ellipsoidFuncEval(ellipsoidCoord);
}
return dispWeight;
......
......@@ -192,6 +192,8 @@ int main (int argc, char** argv)
EssentialPatchBCHandler patchHandler ("listEssentialPatchBC", dataFile);
patchHandler.addPatchBC(solver);
heartSolver.setPatchDisplacement(patchHandler.patchDisplacement());
//============================================
// Setup solver (including fe-spaces & b.c.)
......
......@@ -205,35 +205,41 @@ public:
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Von Mises Stress",
"Von Mises stress",
M_emSolver.electroSolverPtr()->feSpacePtr(),
M_emSolver.tensionEstimator().vonMisesStressPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Principal Stress",
"Principal stress",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.tensionEstimator().principalStressesPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"X Stress Total",
"X stress",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.tensionEstimator().sigmaXPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Y Stress Total",
"Y stress",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.tensionEstimator().sigmaYPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Z Stress Total",
"Z stress",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.tensionEstimator().sigmaZPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Patch displacement",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.structuralOperatorPtr()->displacementPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Fibers",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
......@@ -254,11 +260,11 @@ public:
M_emSolver.activationModelPtr()->fiberActivationPtr(),
UInt (0) );
m_exporter -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Activation Time",
M_emSolver.electroSolverPtr() -> feSpacePtr(),
M_emSolver.activationTimePtr(),
UInt (0) );
m_exporter -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Activation time",
M_emSolver.electroSolverPtr() -> feSpacePtr(),
M_emSolver.activationTimePtr(),
UInt (0) );
for (int i = 0; i < M_emSolver.electroSolverPtr()->ionicModelPtr()->Size(); ++i)
{
......@@ -353,7 +359,12 @@ public:
return traction.dot(velocity);
}
void setPatchDisplacement(boost::shared_ptr<VectorEpetra>& patchDisplacement)
{
m_patchDisplacement = patchDisplacement;
}
protected:
......@@ -369,6 +380,7 @@ protected:
VectorSmall<4> m_ABdplv, m_ABdprv;
boost::shared_ptr<VectorEpetra> m_patchDisplacement;
......
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