Commit bcc82bef authored by thomaskummer's avatar thomaskummer
Browse files

compute von mises stress

parent 328744a9
......@@ -443,11 +443,13 @@ public:
meshPtr_Type M_fullMeshPtr;
vectorPtr_Type M_activationTimePtr;
vectorPtr_Type M_vonMisesStressPtr;
bool M_oneWayCoupling;
WallTensionEstimator<RegionMesh<LinearTetra> > M_wte;
WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
commPtr_Type M_commPtr;
......@@ -471,7 +473,9 @@ EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
M_fullMeshPtr ( ),
M_activationTimePtr ( ),
M_oneWayCoupling (true),
M_wte ( ),
M_wteTotal ( ),
M_wtePassive ( ),
M_wteActive ( ),
M_commPtr (comm),
M_data ()
{
......@@ -493,7 +497,9 @@ EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
M_fullMeshPtr ( solver.M_fullMeshPtr),
M_activationTimePtr ( solver.M_activationTimePtr),
M_oneWayCoupling ( solver.M_oneWayCoupling),
M_wte (solver.M_wte),
M_wteTotal (solver.M_wteTotal),
M_wtePassive (solver.M_wtePassive),
M_wteActive (solver.M_wteActive),
M_commPtr ( solver.M_commPtr),
M_data (solver.M_data)
......@@ -601,7 +607,9 @@ EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);
M_wte.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0);
M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "total");
M_wtePassive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "passive");
M_wteActive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "active");
}
/////////////////////
......@@ -675,11 +683,76 @@ EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
M_vonMisesStressExporterPtr.reset (new exporter_Type() );
setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
M_vonMisesStressPtr.reset (new vector_Type ( M_electroSolverPtr->potentialPtr() -> map() ) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Von Mises Stress",
"Von Mises Stress Total",
M_electroSolverPtr -> feSpacePtr(),
M_wteTotal.vonMisesStressPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"X Stress Total",
M_electroSolverPtr -> feSpacePtr(),
M_wteTotal.sigmaXPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Y Stress Total",
M_electroSolverPtr -> feSpacePtr(),
M_wteTotal.sigmaYPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Z Stress Total",
M_electroSolverPtr -> feSpacePtr(),
M_wteTotal.sigmaZPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Von Mises Stress Passive",
M_electroSolverPtr -> feSpacePtr(),
M_wtePassive.vonMisesStressPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"X Stress Passive",
M_electroSolverPtr -> feSpacePtr(),
M_wtePassive.sigmaXPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Y Stress Passive",
M_electroSolverPtr -> feSpacePtr(),
M_wtePassive.sigmaYPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Z Stress Passive",
M_electroSolverPtr -> feSpacePtr(),
M_wtePassive.sigmaZPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Von Mises Stress Active",
M_electroSolverPtr -> feSpacePtr(),
M_wteActive.vonMisesStressPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"X Stress Active",
M_electroSolverPtr -> feSpacePtr(),
M_wteActive.sigmaXPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Y Stress Active",
M_electroSolverPtr -> feSpacePtr(),
M_wteActive.sigmaYPtr(),
UInt (0) );
M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
"Z Stress Active",
M_electroSolverPtr -> feSpacePtr(),
M_wte.vonMisesStressPtr(),
M_wteActive.sigmaZPtr(),
UInt (0) );
// Mechanics
......@@ -709,6 +782,7 @@ EMSolver<Mesh, ElectroSolver>::setTimeIndex (const UInt& time)
M_electroExporterPtr -> setTimeIndex (time);
M_activationExporterPtr -> setTimeIndex (time);
M_activationTimeExporterPtr -> setTimeIndex (time);
M_vonMisesStressExporterPtr -> setTimeIndex (time);
M_mechanicsExporterPtr -> setTimeIndex (time);
}
......@@ -716,8 +790,12 @@ template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
{
M_wte.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
M_wte.analyzeTensionsRecoveryVonMisesStress();
M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
M_wtePassive.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
M_wtePassive.analyzeTensionsRecoveryVonMisesStress();
M_wteActive.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
M_wteActive.analyzeTensionsRecoveryVonMisesStress();
M_electroExporterPtr -> postProcess (time);//, restart);
//if(M_activationExporterPtr) std::cout << "\nActivation exporter available.";
......
......@@ -342,10 +342,16 @@ public:
// Assemble first piola kirchhoff tensor
firstPiola.Scale(0.0);
firstPiola += Pvol;
firstPiola += Piso;
firstPiola += Pi4;
firstPiola += Pact;
if ( invariants[5] == 0.0 || invariants[5] == 1.0 )
{
firstPiola += Pvol;
firstPiola += Piso;
firstPiola += Pi4;
}
if ( invariants[5] == 0.0 || invariants[5] == 2.0 )
{
firstPiola += Pact;
}
}
//! Get the Stiffness matrix
......
......@@ -190,7 +190,8 @@ public:
const feSpacePtr_Type& feSpace,
const feSpaceETPtr_Type& feSpaceET,
const commPtr_Type& comm,
UInt marker);
UInt marker,
std::string stressType);
//! This method computes the Cauchy stress tensor and its principal values. It uses the displacement vector that has to be set
virtual void analyzeTensions();
......@@ -306,6 +307,20 @@ public:
return *M_sigmaZ;
}
solutionVectPtr_Type sigmaXPtr() const
{
return M_sigmaX;
}
solutionVectPtr_Type sigmaYPtr() const
{
return M_sigmaY;
}
solutionVectPtr_Type sigmaZPtr() const
{
return M_sigmaZ;
}
//! Export the XX component of the stress by copying it to an external vector
/*!
* @param stressXX vector to be filled with the XX component of the stress
......@@ -527,6 +542,8 @@ protected:
//! Material class
materialPtr_Type M_material;
std::string M_stressType;
//@}
};
......@@ -587,7 +604,8 @@ WallTensionEstimator<Mesh >::setup ( const dataPtr_Type& dataMaterial,
const feSpacePtr_Type& feSpace,
const feSpaceETPtr_Type& feSpaceET,
const commPtr_Type& comm,
UInt marker)
UInt marker,
std::string stressType)
{
// Data classes & Volumes markers
M_dataMaterial = dataMaterial;
......@@ -628,7 +646,7 @@ WallTensionEstimator<Mesh >::setup ( const dataPtr_Type& dataMaterial,
M_globalEigenvalues.reset ( new solutionVect_Type (*M_FESpace->mapPtr() ) );
M_invariants.resize ( M_FESpace->fieldDim() + 1 );
M_invariants.resize ( M_FESpace->fieldDim() + 2 );
M_eigenvaluesR.resize ( M_FESpace->fieldDim() );
M_eigenvaluesI.resize ( M_FESpace->fieldDim() );
......@@ -637,6 +655,8 @@ WallTensionEstimator<Mesh >::setup ( const dataPtr_Type& dataMaterial,
M_material.reset ( new EMStructuralConstitutiveLaw<Mesh> );
//ElectrophysiologyUtility::importVectorField (M_material -> fiberVectorPtr(), fileName, fieldName, M_localMeshPtr, postDir, polynomialDegree );
M_material->setup ( M_FESpace, feSpaceET, M_FESpace->mapPtr(), M_offset, M_dataMaterial, M_displayer );
M_stressType = stressType;
}
template <typename Mesh>
......@@ -1253,6 +1273,9 @@ WallTensionEstimator<Mesh >::computeInvariantsRightCauchyGreenTensor (std::vecto
invariants[1] = I4f; // Fourth invariant C
invariants[2] = H; // Fiber activation
invariants[3] = tensorF (0, 0) * ( tensorF (1, 1) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 1) ) - tensorF (0, 1) * ( tensorF (1, 0) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 0) ) + tensorF (0, 2) * ( tensorF (1, 0) * tensorF (2, 1) - tensorF (1, 1) * tensorF (2, 0) ); // Determinant F
if (M_stressType == "total") invariants[5] = 0.0;
if (M_stressType == "passive") invariants[5] = 1.0;
if (M_stressType == "active") invariants[5] = 2.0;
//Computation of the Cofactor of F
cofactorF ( 0 , 0 ) = ( tensorF (1, 1) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 1) );
......
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