Commit 115a7056 authored by thomaskummer's avatar thomaskummer
Browse files

intrepid lifev update

parent c4e8250c
......@@ -342,10 +342,51 @@ public:
return m_exporter;
}
Real
externalPower ( const VectorEpetra& dispCurrent,
const VectorEpetra& dispPrevious,
const boost::shared_ptr<ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > > dispETFESpace,
Real pressure,
Real dt,
const unsigned int bdFlag) const
{
VectorEpetra traction ( dispCurrent.map() );
VectorEpetra velocity ( (dispCurrent - dispPrevious) / dt );
MatrixSmall<3,3> Id;
Id(0,0) = 1.; Id(0,1) = 0., Id(0,2) = 0.;
Id(1,0) = 0.; Id(1,1) = 1., Id(1,2) = 0.;
Id(2,0) = 0.; Id(2,1) = 0., Id(2,2) = 1.;
{
using namespace ExpressionAssembly;
auto I = value(Id);
auto Grad_u = grad( dispETFESpace, dispCurrent, 0);
auto F = Grad_u + I;
auto FmT = minusT(F);
auto J = det(F);
auto p = value(pressure);
QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
integrate ( boundary ( dispETFESpace->mesh(), bdFlag),
myBDQR,
dispETFESpace,
p * J * dot( FmT * Nface, phi_i)
//p * J * dot( FmT * Nface, phi_i)
//value(-1.0) * J * dot (vE1, FmT * Nface) * phi_i) >> intergral
) >> traction;
traction.globalAssemble();
}
return traction.dot(velocity);
}
protected:
EmSolver& M_emSolver;
Circulation& M_circulationSolver;
......@@ -455,50 +496,6 @@ protected:
}
Real
externalPower ( const VectorEpetra& dispCurrent,
const VectorEpetra& dispPrevious,
const boost::shared_ptr<ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > > dispETFESpace,
Real pressure,
Real dt,
const unsigned int bdFlag) const
{
VectorEpetra traction ( dispCurrent.map() );
VectorEpetra velocity ( (dispCurrent - dispPrevious) / dt );
MatrixSmall<3,3> Id;
Id(0,0) = 1.; Id(0,1) = 0., Id(0,2) = 0.;
Id(1,0) = 0.; Id(1,1) = 1., Id(1,2) = 0.;
Id(2,0) = 0.; Id(2,1) = 0., Id(2,2) = 1.;
{
using namespace ExpressionAssembly;
auto I = value(Id);
auto Grad_u = grad( dispETFESpace, dispCurrent, 0);
auto F = Grad_u + I;
auto FmT = minusT(F);
auto J = det(F);
auto p = value(pressure);
QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
integrate ( boundary ( dispETFESpace->mesh(), bdFlag),
myBDQR,
dispETFESpace,
p * J * dot( FmT * Nface, phi_i)
//p * J * dot( FmT * Nface, phi_i)
//value(-1.0) * J * dot (vE1, FmT * Nface) * phi_i) >> intergral
) >> traction;
traction.globalAssemble();
}
return traction.dot(velocity);
}
};
}
......
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