Commit cc664c34 authored by thomaskummer's avatar thomaskummer
Browse files

compute von mises stress

parent b4bcb6d5
......@@ -576,6 +576,13 @@ int main (int argc, char** argv)
LifeChrono chronoPreload;
chronoPreload.start();
for (int j (0); j < 25; ++j) {
solver.solveElectrophysiology (stim, j*0.1);
solver.solveActivation (0.1);
solver.saveSolution (j);
}
for (int i (1); i <= preloadSteps; i++)
{
if ( 0 == comm->MyPID() )
......
......@@ -621,7 +621,7 @@ EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);
M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "total");
M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
// M_wtePassive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "passive");
// M_wteActive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "active");
}
......@@ -816,6 +816,8 @@ EMSolver<Mesh, ElectroSolver>::setTimeIndex (const UInt& time)
M_activationExporterPtr -> setTimeIndex (time);
M_activationTimeExporterPtr -> setTimeIndex (time);
M_vonMisesStressExporterPtr -> setTimeIndex (time);
M_vonMisesStressExporterPtrP -> setTimeIndex (time);
M_vonMisesStressExporterPtrA -> setTimeIndex (time);
M_mechanicsExporterPtr -> setTimeIndex (time);
}
......
......@@ -339,6 +339,9 @@ public:
// Active stress part
auto Pact = tensorF;
Pact.Scale( invariants[1] * 0.5 * std::pow(invariants[2], 2.0) * 1300000 );
//if ( invariants[2] > 0. ) {
// std::cout << invariants[2] << std::endl;
//Pact.Print(std::cout); }
// Assemble first piola kirchhoff tensor
firstPiola.Scale(0.0);
......
......@@ -191,7 +191,7 @@ public:
const feSpaceETPtr_Type& feSpaceET,
const commPtr_Type& comm,
UInt marker,
std::string stressType);
const boost::shared_ptr< EMStructuralConstitutiveLaw<Mesh> >& EMMaterial);
//! This method computes the Cauchy stress tensor and its principal values. It uses the displacement vector that has to be set
virtual void analyzeTensions();
......@@ -610,7 +610,7 @@ WallTensionEstimator<Mesh >::setup ( const dataPtr_Type& dataMaterial,
const feSpaceETPtr_Type& feSpaceET,
const commPtr_Type& comm,
UInt marker,
std::string stressType)
const boost::shared_ptr< EMStructuralConstitutiveLaw<Mesh> >& EMMaterial )
{
// Data classes & Volumes markers
M_dataMaterial = dataMaterial;
......@@ -657,11 +657,11 @@ WallTensionEstimator<Mesh >::setup ( const dataPtr_Type& dataMaterial,
// Materials
//M_material.reset ( material_Type::StructureMaterialFactory::instance().createObject ( M_dataMaterial->solidType() ) );
M_material.reset ( new EMStructuralConstitutiveLaw<Mesh> );
M_material = EMMaterial;
//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_material->setup ( M_FESpace, feSpaceET, M_FESpace->mapPtr(), M_offset, M_dataMaterial, M_displayer );
M_stressType = stressType;
M_stressType = "total";
}
template <typename Mesh>
......@@ -1115,8 +1115,9 @@ WallTensionEstimator<Mesh >::constructGlobalStressVector ()
//Copying the displacement field into a vector with repeated map for parallel computations
solutionVect_Type dRep (*M_displacement, Repeated);
solutionVect_Type fRep (*(M_material -> fiberVectorPtr()), Repeated);
solutionVect_Type ARep (*(M_material -> fiberActivationPtr()), Repeated);
solutionVect_Type fARep (*(M_material -> fiberActivationPtr()), Repeated);
//if (fARep.maxValue() > 0. ) std::cout << fARep.maxValue() << "maxV" << std::endl;
//Creating the local stress tensors
VectorElemental elVecSigmaX (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
VectorElemental elVecSigmaY (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
......@@ -1147,9 +1148,13 @@ WallTensionEstimator<Mesh >::constructGlobalStressVector ()
dk_loc[iloc + iComp * fakeFESpace.fe().nbFEDof() ] = dRep[ig + iComp * fakeFESpace.dim()];
fk_loc[iloc + iComp * fakeFESpace.fe().nbFEDof() ] = fRep[ig + iComp * fakeFESpace.dim()];
}
fAk_loc[iloc] = fRep[ig];
//if ( fARep[ig] > 0. ) std::cout << fARep[ig] << "iiiiiii" << std::endl;
fAk_loc[iloc] = fARep[ig];
}
//if (dRep.maxValue() > 0. ) std::cout << dRep.maxValue() << "maxd" << std::endl;
//if (fARep.maxValue() > 0. ) std::cout << fARep.maxValue() << "maxV" << std::endl;
//if (fRep.maxValue() > 0. ) std::cout << fRep.maxValue() << "maxf" << std::endl;
//Compute the element tensor F
AssemblyElementalStructure::computeLocalDeformationGradient ( dk_loc, vectorDeformationF, fakeFESpace.fe() );
......
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