Skip to content

Commit

Permalink
Modified definition of contact material point in several sliding and …
Browse files Browse the repository at this point in the history
…tied interfaced to facilitate plotting of contact gap and pressure gap in tied-biphasic and tied-multiphasic interfaces
  • Loading branch information
gateshian committed Sep 2, 2024
1 parent d725351 commit 7d99494
Show file tree
Hide file tree
Showing 19 changed files with 417 additions and 537 deletions.
2 changes: 1 addition & 1 deletion FEBioMech/FEBioMechPlot.h
Original file line number Diff line number Diff line change
Expand Up @@ -686,7 +686,7 @@ class FEPlotShellDirector : public FEPlotDomainData
class FEPlotElementElasticity : public FEPlotDomainData
{
public:
FEPlotElementElasticity(FEModel* pfem) : FEPlotDomainData(pfem, PLT_TENS4FS, FMT_ITEM){}
FEPlotElementElasticity(FEModel* pfem) : FEPlotDomainData(pfem, PLT_TENS4FS, FMT_ITEM) { SetUnits(UNIT_PRESSURE); }
bool Save(FEDomain& dom, FEDataStream& a);
};

Expand Down
2 changes: 2 additions & 0 deletions FEBioMix/FEBioMix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,7 @@ void FEBioMix::InitModule()
REGISTER_FECORE_CLASS(FEPlotActualFluidPressure , "fluid pressure" );
REGISTER_FECORE_CLASS(FEPlotFluidFlux , "fluid flux" );
REGISTER_FECORE_CLASS(FEPlotNodalFluidFlux , "nodal fluid flux");
REGISTER_FECORE_CLASS(FEPlotContactGapMP , "contact gap" );
REGISTER_FECORE_CLASS(FEPlotPressureGap , "pressure gap" );
REGISTER_FECORE_CLASS(FEPlotFluidForce , "fluid force" );
REGISTER_FECORE_CLASS(FEPlotFluidForce2 , "fluid force2" );
Expand Down Expand Up @@ -331,6 +332,7 @@ void FEBioMix::InitModule()
REGISTER_FECORE_CLASS(FEPlotEffectiveSoluteConcentration , "effective solute concentration");
REGISTER_FECORE_CLASS(FEPlotEffectiveShellSoluteConcentration, "effective shell solute concentration");
REGISTER_FECORE_CLASS(FEPlotActualSoluteConcentration , "solute concentration");
REGISTER_FECORE_CLASS(FEPlotConcentrationGap , "concentration gap" );
REGISTER_FECORE_CLASS(FEPlotPartitionCoefficient , "partition coefficient");
REGISTER_FECORE_CLASS(FEPlotSoluteFlux , "solute flux" );
REGISTER_FECORE_CLASS(FEPlotSoluteVolumetricFlux , "solute volumetric flux" );
Expand Down
102 changes: 100 additions & 2 deletions FEBioMix/FEBioMixPlot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,15 @@ SOFTWARE.*/
#include "FEMultiphasicSolidDomain.h"
#include "FEMultiphasicShellDomain.h"
#include <FEBioMech/FEElasticSolidDomain.h>
#include <FEBioMech/FESlidingInterface.h>
#include "FEBiphasic.h"
#include "FEBiphasicSolute.h"
#include "FETriphasic.h"
#include "FEMultiphasic.h"
#include "FEBiphasicContactSurface.h"
#include "FESlidingInterfaceMP.h"
#include "FETiedBiphasicInterface.h"
#include "FETiedMultiphasicInterface.h"
#include "FEBioMech/FEDonnanEquilibrium.h"
#include "FEBioMech/FEElasticMixture.h"
#include <FECore/FEModel.h>
Expand Down Expand Up @@ -135,13 +139,34 @@ bool FEPlotMixtureFluidFlowRate::Save(FESurface &surf, FEDataStream &a)

//-----------------------------------------------------------------------------
// Plot contact gap
bool FEPlotContactGapMP::Save(FESurface& surf, FEDataStream& a)
{
FEContactSurface* pcs = dynamic_cast<FEContactSurface*>(&surf);
if (pcs == 0) return false;

writeAverageElementValue<double>(surf, a, [](const FEMaterialPoint& mp) {
const FEContactMaterialPoint* pt = dynamic_cast<const FEContactMaterialPoint*>(&mp);
double d = (pt ? pt->m_gap : 0);
if (d == 0) {
const FETiedBiphasicContactPoint* pd = dynamic_cast<const FETiedBiphasicContactPoint*>(&mp);
vec3d vd = (pd ? pd->m_dg : vec3d(0,0,0));
d = (pd ? vd.unit() : 0);
}
return d;
});

return true;
}

//-----------------------------------------------------------------------------
// Plot pressure gap
bool FEPlotPressureGap::Save(FESurface& surf, FEDataStream& a)
{
FEBiphasicContactSurface* pcs = dynamic_cast<FEBiphasicContactSurface*>(&surf);
if (pcs == 0) return false;

writeNodalProjectedElementValues<double>(surf, a, [](const FEMaterialPoint& mp) {
const FEBiphasicContactPoint* pt = mp.ExtractData<FEBiphasicContactPoint>();
writeAverageElementValue<double>(surf, a, [](const FEMaterialPoint& mp) {
const FEBiphasicContactPoint* pt = dynamic_cast<const FEBiphasicContactPoint*>(&mp);
return (pt ? pt->m_pg : 0);
});

Expand Down Expand Up @@ -246,6 +271,79 @@ bool FEPlotFluidLoadSupport::Save(FESurface &surf, FEDataStream &a)
return true;
}

//-----------------------------------------------------------------------------
// Plot concentration gap
FEPlotConcentrationGap::FEPlotConcentrationGap(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_ARRAY, FMT_ITEM)
{
DOFS& dofs = pfem->GetDOFS();
int nsol = dofs.GetVariableSize("concentration");
SetArraySize(nsol);

// collect the names
int ndata = pfem->GlobalDataItems();
vector<string> s;
for (int i = 0; i<ndata; ++i)
{
FESoluteData* ps = dynamic_cast<FESoluteData*>(pfem->GetGlobalData(i));
if (ps)
{
s.push_back(ps->GetName());
m_sol.push_back(ps->GetID());
}
}
assert(nsol == (int)s.size());
SetArrayNames(s);
SetUnits(UNIT_CONCENTRATION);
}

bool FEPlotConcentrationGap::Save(FESurface& surf, FEDataStream& a)
{
FEContactSurface* pcs = dynamic_cast<FEContactSurface*>(&surf);
if (pcs == 0) return false;

for (int i=0; i<surf.Elements(); ++i) {
FESurfaceElement& el = surf.Element(i);

FEElement* se = (el.m_elem[0]).pe;
FEMaterial* mat = GetFEModel()->GetMaterial(se->GetMatID());
FESoluteInterface* pm = dynamic_cast<FESoluteInterface*>(mat);
if ((pm == 0) || (pm->Solutes() == 0)) return false;

// figure out the local solute IDs. This depends on the material
int nsols = (int)m_sol.size();
vector<int> lid(nsols, -1);
int nsc = 0;
for (int i = 0; i<(int)m_sol.size(); ++i)
{
lid[i] = pm->FindLocalSoluteID(m_sol[i]);
if (lid[i] != -1) nsc++;
}
if (nsc == 0) return false;

for (int k=0; k<nsols; ++k)
{
int nsid = lid[k];
if (nsid == -1) a << 0.f;
else
{
// calculate average concentration gp
double ew = 0;
for (int j = 0; j<el.GaussPoints(); ++j)
{
FEMaterialPoint& mp = *el.GetMaterialPoint(j);
const FEMultiphasicContactPoint* pt = dynamic_cast<const FEMultiphasicContactPoint*>(&mp);
const FETiedMultiphasicContactPoint* tt = dynamic_cast<const FETiedMultiphasicContactPoint*>(&mp);
if (pt) ew += pt->m_cg[nsid];
else if (tt) ew += tt->m_cg[nsid];
}
ew /= el.GaussPoints();
a << ew;
}
}
}
return true;
}

//=============================================================================
// D O M A I N D A T A
//=============================================================================
Expand Down
25 changes: 24 additions & 1 deletion FEBioMix/FEBioMixPlot.h
Original file line number Diff line number Diff line change
Expand Up @@ -363,14 +363,23 @@ class FEPlotFluidForce2 : public FEPlotSurfaceData
bool Save(FESurface& surf, FEDataStream& a);
};

//-----------------------------------------------------------------------------
//! Contact gap for multiphasic interfaces
//!
class FEPlotContactGapMP : public FEPlotSurfaceData
{
public:
FEPlotContactGapMP(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_FLOAT, FMT_ITEM) { SetUnits(UNIT_LENGTH); }
bool Save(FESurface& surf, FEDataStream& a);
};

//-----------------------------------------------------------------------------
//! Fluid pressure gap
//!
class FEPlotPressureGap : public FEPlotSurfaceData
{
public:
FEPlotPressureGap(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_FLOAT, FMT_MULT) { SetUnits(UNIT_PRESSURE); }
FEPlotPressureGap(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_FLOAT, FMT_ITEM) { SetUnits(UNIT_PRESSURE); }
bool Save(FESurface& surf, FEDataStream& a);
};

Expand All @@ -383,3 +392,17 @@ class FEPlotFluidLoadSupport : public FEPlotSurfaceData
FEPlotFluidLoadSupport(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_FLOAT, FMT_REGION){}
bool Save(FESurface& surf, FEDataStream& a);
};

//-----------------------------------------------------------------------------
//! concentration gap
//!
class FEPlotConcentrationGap : public FEPlotSurfaceData
{
public:
FEPlotConcentrationGap(FEModel* pfem);
bool Save(FESurface& surf, FEDataStream& a);

protected:
vector<int> m_sol;
};

12 changes: 12 additions & 0 deletions FEBioMix/FEBiphasicContactSurface.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,18 @@ SOFTWARE.*/
void FEBiphasicContactPoint::Serialize(DumpStream& ar)
{
FEContactMaterialPoint::Serialize(ar);
ar & m_dg;
ar & m_Lmd;
ar & m_Lmt;
ar & m_epsn;
ar & m_epsp;
ar & m_p1;
ar & m_nu;
ar & m_s1;
ar & m_tr;
ar & m_rs;
ar & m_rsp;
ar & m_bstick;
ar & m_Lmp & m_pg & m_mueff & m_fls;
}

Expand Down
28 changes: 26 additions & 2 deletions FEBioMix/FEBiphasicContactSurface.h
Original file line number Diff line number Diff line change
Expand Up @@ -34,14 +34,38 @@ SOFTWARE.*/
class FEBIOMIX_API FEBiphasicContactPoint : public FEContactMaterialPoint
{
public:
double m_Lmp; //!< lagrange multipliers for fluid pressures
double m_pg; //!< pressure "gap" for biphasic contact
vec3d m_dg; //!< vector gap
double m_Lmd; //!< Lagrange multipliers for normal traction
vec3d m_Lmt; //!< Lagrange multipliers for vector traction
double m_epsn; //!< penalty factor
double m_epsp; //!< pressure penalty factor
double m_p1; //!< fluid pressure
vec3d m_nu; //!< normal at integration points
vec3d m_s1; //!< tangent along slip direction
vec3d m_tr; //!< contact traction
vec2d m_rs; //!< natural coordinates of projection
vec2d m_rsp; //!< m_rs at the previous time step
bool m_bstick; //!< stick flag
double m_Lmp; //!< lagrange multipliers for fluid pressures
double m_pg; //!< pressure "gap" for biphasic contact
double m_mueff; //!< effective friction coefficient
double m_fls; //!< local fluid load support

void Init() override
{
FEContactMaterialPoint::Init();
m_Lmd = 0.0;
m_Lmt = m_tr = vec3d(0,0,0);
m_Lmp = 0.0;
m_epsn = 1.0;
m_epsp = 1.0;
m_pg = 0.0;
m_p1 = 0.0;
m_mueff = 0.0;
m_fls = 0.0;
m_nu = m_s1 = m_dg = vec3d(0,0,0);
m_rs = m_rsp = vec2d(0,0);
m_bstick = false;
m_Lmp = 0.0;
m_pg = 0.0;
m_mueff = 0.0;
Expand Down
Loading

0 comments on commit 7d99494

Please sign in to comment.