Commit 9434cd5b authored by thomaskummer's avatar thomaskummer
Browse files

automatic quad rule choice

parent 7adc9b00
......@@ -231,7 +231,7 @@ int main (int argc, char** argv)
displayer.leaderPrint ("\nSetting up EM solver ... ");
if ( dataFile ( "solid/physics/EMPassiveMaterialType", "4pt") == "PHO" )
if ( dataFile ( "solid/physics/EMPassiveMaterialType", "4pt") == "PHO" && dataFile ( "solid/space_discretization/order", "P1") == "P2" )
{
EMAssembler::quadRule.setQuadRule("15pt");
}
......
......@@ -51,6 +51,8 @@
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>
#include <lifev/structure/solver/WallTensionEstimator.hpp>
//#include <lifev/em/solver/activation/activeStressModels/ActiveStressActivation.hpp>
#include <lifev/em/solver/activation/ActivationModelsList.hpp>
......@@ -88,7 +90,7 @@ public:
typedef boost::shared_ptr<Epetra_Comm> commPtr_Type;
typedef VectorEpetra vector_Type;
typedef StructuralConstitutiveLawData structureData_Type;
typedef boost::shared_ptr<structureData_Type> structureDataPtr_Type;
......@@ -175,6 +177,11 @@ public:
EMUtility::setupExporter<Mesh> (*M_activationTimeExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
void setupVonMisesStressExporter ( std::string problemFolder = "./", std::string outputFileName = "VonMisesStress" )
{
EMUtility::setupExporter<Mesh> (*M_vonMisesStressExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
}
void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
{
if (M_mechanicsExporterPtr)
......@@ -187,7 +194,8 @@ public:
std::string electroFileName = "ElectroSolution",
std::string activationFileName = "ActivationSolution",
std::string activationTimeFileName = "ActivationTimeSolution",
std::string mechanicsFileName = "MechanicalSolution");
std::string mechanicsFileName = "MechanicalSolution",
std::string vonMisesStressFileName = "VonMisesStress");
void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
{
......@@ -269,7 +277,17 @@ public:
const std::string& postDir = "./",
const std::string& polynomialDegree = "P1" )
{
ElectrophysiologyUtility::importVectorField (getMechanicsFibers(), fileName, fieldName, M_localMeshPtr, postDir, polynomialDegree );
// if ( polynomialDegree == "P1" )
// {
// FESpace<RegionMesh<LinearTetra> , MapEpetra > p2FESpace ( fullMeshPtr(), "P2", 1, fullMeshPtr() -> comm() );
// vector_Type p2FibersRep (*getMechanicsFibers(), Repeated);
// vector_Type p1FibersRep = M_EMStructuralOperatorPtr -> dispFESpacePtr() -> feToFEInterpolate(p2FESpace, p2FibersRep);
// getMechanicsFibers().reset ( new vector_Type (p1FibersRep, Unique) );
// }
}
void setupMechanicalSheetVector( const std::string& fileName,
......@@ -419,13 +437,17 @@ public:
exporterPtr_Type M_activationExporterPtr;
exporterPtr_Type M_mechanicsExporterPtr;
exporterPtr_Type M_activationTimeExporterPtr;
exporterPtr_Type M_vonMisesStressExporterPtr;
meshPtr_Type M_localMeshPtr;
meshPtr_Type M_fullMeshPtr;
vectorPtr_Type M_activationTimePtr;
vectorPtr_Type M_vonMisesStressPtr;
bool M_oneWayCoupling;
WallTensionEstimator<RegionMesh<LinearTetra> > M_wte;
commPtr_Type M_commPtr;
......@@ -449,6 +471,7 @@ EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
M_fullMeshPtr ( ),
M_activationTimePtr ( ),
M_oneWayCoupling (true),
M_wte ( ),
M_commPtr (comm),
M_data ()
{
......@@ -470,6 +493,7 @@ 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_commPtr ( solver.M_commPtr),
M_data (solver.M_data)
......@@ -577,7 +601,8 @@ 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_wte.analyzeTensionsRecoveryVonMisesStress();
}
/////////////////////
......@@ -607,7 +632,8 @@ EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
std::string electroFileName,
std::string activationFileName,
std::string activationTimeFileName,
std::string mechanicsFileName)
std::string mechanicsFileName,
std::string vonMisesStressFileName)
{
if (M_commPtr -> MyPID() == 0)
{
......@@ -646,6 +672,17 @@ EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
M_activationTimePtr,
UInt (0) );
// Von Mises stress
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",
M_electroSolverPtr -> feSpacePtr(),
M_wte.vonMisesStressPtr(),
UInt (0) );
// Mechanics
M_mechanicsExporterPtr.reset (new exporter_Type() );
setupMechanicsExporter (problemFolder, mechanicsFileName);
......@@ -685,6 +722,8 @@ EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
//if(M_activationModelPtr -> fiberActivationPtr()) std::cout << "\nFiber activation exporter available.";
M_activationExporterPtr -> postProcess (time);//, restart);
M_activationTimeExporterPtr -> postProcess (time);
M_vonMisesStressExporterPtr -> postProcess (time);
//if(M_mechanicsExporterPtr) std::cout << "\nMechanics exporter available.";
M_mechanicsExporterPtr -> postProcess (time);//, restart);
}
......@@ -749,6 +788,145 @@ EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
} // namespace LifeV
//void constructGlobalStressVector()
//{
// //Chrono
// LifeChrono chrono;
// chrono.start();
//
// // Reset stress components
// boost::shared_ptr<Epetra_SerialDenseMatrix> sigma ( new Epetra_SerialDenseMatrix ( M_FESpace->fieldDim(), M_FESpace->fieldDim() ) );
// boost::shared_ptr<VectorEpetra> sigmaX ( new VectorEpetra (*M_FESpace->mapPtr() ) );
// boost::shared_ptr<VectorEpetra> sigmaY ( new VectorEpetra (*M_FESpace->mapPtr() ) );
// boost::shared_ptr<VectorEpetra> sigmaZ ( new VectorEpetra (*M_FESpace->mapPtr() ) );
//// *sigmaX *= 0.;
//// *sigmaY *= 0.;
//// *sigmaZ *= 0.;
//
// //Constructing the patch area vector for reconstruction purposes
// VectorEpetra patchArea (*M_displacement, Unique, Add);
// patchArea *= 0.0;
//
// constructPatchAreaVector ( patchArea );
//
// //Before assembling the reconstruction process is done
// solutionVect_Type patchAreaR (patchArea, Repeated);
//
// //Compute the area of reference element
// Real refElemArea (0);
// for (UInt iq (0); iq < M_FESpace->qr().nbQuadPt(); ++iq)
// {
// refElemArea += M_FESpace->qr().weight (iq);
// }
//
// //Setting the quadrature Points = DOFs of the element and weight = 1
// Real wQuad (refElemArea / M_FESpace->refFE().nbDof() );
// std::vector<Real> weights (M_FESpace->fe().nbFEDof(), wQuad);
// std::vector<GeoVector> coords = M_FESpace->refFE().refCoor();
//
// QuadratureRule fakeQuadratureRule;
// fakeQuadratureRule.setDimensionShape ( shapeDimension (M_FESpace->refFE().shape() ), M_FESpace->refFE().shape() );
// fakeQuadratureRule.setPoints (coords, weights);
//
// //Creating a copy of the FESpace
// feSpace_Type fakeFESpace ( M_FESpace->mesh(), M_FESpace->refFE(), M_FESpace->qr(), M_FESpace->bdQr(), 3, M_FESpace->map().commPtr() );
//
// //Set the new quadrature rule
// fakeFESpace.setQuadRule (fakeQuadratureRule);
//
// //Preliminary variables
// UInt totalDof = fakeFESpace.dof().numTotalDof();
// VectorElemental dk_loc (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//
// //Vectors for the deformation tensor
// (*M_deformationF).Scale (0.0);
// std::vector<matrix_Type> vectorDeformationF (fakeFESpace.fe().nbFEDof(), *M_deformationF);
//
// //Copying the displacement field into a vector with repeated map for parallel computations
// solutionVect_Type dRep (*M_displacement, Repeated);
//
// //Creating the local stress tensors
// VectorElemental elVecSigmaX (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
// VectorElemental elVecSigmaY (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
// VectorElemental elVecSigmaZ (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//
// //Loop on each volume
// for ( UInt i (0); i < fakeFESpace.mesh()->numVolumes(); ++i )
// {
// //fakeFESpace.fe().update ( fakeFESpace.mesh()->volumeList ( i ), UPDATE_DPHI | UPDATE_WDET );
// fakeFESpace.fe().updateFirstDerivQuadPt ( fakeFESpace.mesh()->volumeList ( i ) );
//
// elVecSigmaX.zero();
// elVecSigmaY.zero();
// elVecSigmaZ.zero();
//
// M_marker = fakeFESpace.mesh()->volumeList ( i ).markerID();
//
// UInt eleID = fakeFESpace.fe().currentLocalId();
//
// //Extracting the local displacement
// for ( UInt iNode (0); iNode < ( UInt ) fakeFESpace.fe().nbFEDof(); ++iNode )
// {
// UInt iloc = fakeFESpace.fe().patternFirst ( iNode );
//
// for ( UInt iComp = 0; iComp < fakeFESpace.fieldDim(); ++iComp )
// {
// UInt ig = fakeFESpace.dof().localToGlobalMap ( eleID, iloc ) + iComp * fakeFESpace.dim() + M_offset;
// dk_loc[iloc + iComp * fakeFESpace.fe().nbFEDof() ] = dRep[ig];
// }
// }
//
// //Compute the element tensor F
// AssemblyElementalStructure::computeLocalDeformationGradient ( dk_loc, vectorDeformationF, fakeFESpace.fe() );
//
// //Compute the local vector of the principal stresses
// for ( UInt nDOF (0); nDOF < ( UInt ) fakeFESpace.fe().nbFEDof(); ++nDOF )
// {
// UInt iloc = fakeFESpace.fe().patternFirst ( nDOF );
//
// sigma->Scale (0.0);
// M_firstPiola->Scale (0.0);
// M_cofactorF->Scale (0.0);
//
// //Compute the rightCauchyC tensor
// AssemblyElementalStructure::computeInvariantsRightCauchyGreenTensor (M_invariants, vectorDeformationF[nDOF], *M_cofactorF);
//
// //Compute the first Piola-Kirchhoff tensor
// M_material->computeLocalFirstPiolaKirchhoffTensor (*M_firstPiola, vectorDeformationF[nDOF], *M_cofactorF, M_invariants, M_marker);
//
// //Compute the Cauchy tensor
// AssemblyElementalStructure::computeCauchyStressTensor (*sigma, *M_firstPiola, M_invariants[3], vectorDeformationF[nDOF]);
//
// //Assembling the local vectors for local tensions Component X, Y, Z
// for ( UInt coor (0); coor < fakeFESpace.fieldDim(); ++coor )
// {
// (elVecSigmaX) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 0);
// (elVecSigmaY) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 1);
// (elVecSigmaZ) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 2);
// }
// }
//
// reconstructElementaryVector ( elVecSigmaX, patchAreaR, fakeFESpace );
// reconstructElementaryVector ( elVecSigmaY, patchAreaR, fakeFESpace );
// reconstructElementaryVector ( elVecSigmaZ, patchAreaR, fakeFESpace );
//
// //Assembling the three elemental vector in the three global
// for ( UInt ic = 0; ic < fakeFESpace.fieldDim(); ++ic )
// {
// assembleVector (*sigmaX, elVecSigmaX, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset + ic * totalDof );
// assembleVector (*sigmaY, elVecSigmaY, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset + ic * totalDof );
// assembleVector (*sigmaZ, elVecSigmaZ, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset + ic * totalDof );
// }
// }
//
// sigmaX->globalAssemble();
// sigmaY->globalAssemble();
// sigmaZ->globalAssemble();
//
// //Chrono
// chrono.stop();
// M_displayer->leaderPrint (" S- Cauchy stresses recovered in: ", chrono.globalDiff ( *M_displayer->comm() ), " s\n" );
//}
......
......@@ -322,7 +322,12 @@ public:
const Epetra_SerialDenseMatrix& tensorF,
const Epetra_SerialDenseMatrix& cofactorF,
const std::vector<Real>& invariants,
const UInt material) {}
const UInt material)
{
firstPiola = tensorF; //.Scale(2.);
}
//! Get the Stiffness matrix
matrixPtr_Type const stiffMatrix() const
{
......
......@@ -396,6 +396,11 @@ public:
stressVonMises = *M_sigmaVonMises;
}
solutionVectPtr_Type vonMisesStressPtr()
{
return M_sigmaVonMises;
}
//! Get the global vector for the eigenvalues
solutionVect_Type principalStresses() const
{
......
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