Commit 65e6d079 authored by thomaskummer's avatar thomaskummer
Browse files

intrepid lifev update

parent 34fb8dcd
......@@ -425,6 +425,8 @@ public:
void computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
void computeDeformedFiberDirection (VectorEpetra& f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
vectorPtr_Type getElectroFibers()
{
return M_electroSolverPtr -> fiberPtr();
......@@ -1043,6 +1045,48 @@ EMSolver<Mesh, ElectroSolver>::computeI4f (VectorEpetra& i4f, VectorEpetra& f0_,
}
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::computeDeformedFiberDirection (VectorEpetra& f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr)
{
VectorEpetra dUdx (disp);
VectorEpetra dUdy (disp);
VectorEpetra dUdz (disp);
dUdx = GradientRecovery::ZZGradient (feSpacePtr, disp, 0);
dUdy = GradientRecovery::ZZGradient (feSpacePtr, disp, 1);
dUdz = GradientRecovery::ZZGradient (feSpacePtr, disp, 2);
int n = f.epetraVector().MyLength();
int i (0); int j (0); int k (0);
MatrixSmall<3,3> F; VectorSmall<3> f0;
for (int p (0); p < n; p++)
{
i = dUdx.blockMap().GID (p);
j = dUdx.blockMap().GID (p + n);
k = dUdx.blockMap().GID (p + 2 * n);
F(0,0) = 1.0 + dUdx[i];
F(0,1) = dUdy[i];
F(0,2) = dUdz[i];
F(1,0) = dUdx[j];
F(1,1) = 1.0 + dUdy[j];
F(1,2) = dUdz[j];
F(2,0) = dUdx[k];
F(2,1) = dUdy[k];
F(2,2) = 1.0 + dUdz[k];
f0(0) = f0_[i];
f0(1) = f0_[j];
f0(2) = f0_[k];
f0.normalize();
f[i] = F * f0;
}
}
} // namespace LifeV
......
......@@ -269,7 +269,8 @@ public:
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"Fibers",
M_emSolver.structuralOperatorPtr()->dispFESpacePtr(),
M_emSolver.structuralOperatorPtr()->EMMaterial()->fiberVectorPtr(),
M_emSolver.structuralOperatorPtr()->fPtr(),
// M_emSolver.structuralOperatorPtr()->EMMaterial()->fiberVectorPtr(),
UInt (0) );
m_exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
......@@ -314,9 +315,15 @@ public:
void postProcess(const Real& time)
{
// Compute Von Mises stress
M_emSolver.tensionEstimator().setDisplacement ( M_emSolver.structuralOperatorPtr()->displacement() );
M_emSolver.tensionEstimator().analyzeTensionsRecoveryVonMisesStress();
// Compute deformed fiber direction
M_emSolver.computeI4f (M_emSolver.structuralOperatorPtr()->f(), *M_emSolver.structuralOperatorPtr()->EMMaterial()->fiberVectorPtr(), *M_emSolver.structuralOperatorPtr()->displacementPtr(), M_emSolver.structuralOperatorPtr()->dispFESpacePtr());
// Write on hdf5 output file
m_exporter->postProcess(time);
}
......
......@@ -178,7 +178,6 @@ public:
M_I4fPtr = ptr;
}
void setI4f(vector_Type& i4f)
{
*M_I4fPtr = i4f;
......@@ -188,6 +187,16 @@ public:
{
return M_I4fPtr;
}
vectorPtr_Type fPtr()
{
return M_fPtr;
}
VectorEpetra& f()
{
return *M_fPtr;
}
protected:
......@@ -204,6 +213,8 @@ protected:
vectorPtr_Type M_I4fPtr;
vectorPtr_Type M_fPtr;
};
//====================================
......@@ -234,8 +245,15 @@ EMStructuralOperator<Mesh>::setup (boost::shared_ptr<data_Type> data,
this->super::setup (data, dFESpace, dETFESpace, BCh, comm);
M_EMMaterial = boost::dynamic_pointer_cast<material_Type> (this -> material() );
M_boundaryVectorPtr.reset(new vector_Type ( this->M_disp->map(), Repeated ) );
// Setup I4 vector
M_I4fPtr.reset(new vector_Type ( M_EMMaterial->scalarETFESpacePtr()->map() ) );
*M_I4fPtr += 1.0;
// Setup deformed fiber vector
M_fPtr.reset(new vector_Type ( this->M_disp->map() ));
*M_fPtr *= 0.;
}
......
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