From 7d99494a9710857fb41901d20016c70ac0319ae2 Mon Sep 17 00:00:00 2001 From: Gerard Ateshian Date: Sun, 1 Sep 2024 21:00:34 -0400 Subject: [PATCH] Modified definition of contact material point in several sliding and tied interfaced to facilitate plotting of contact gap and pressure gap in tied-biphasic and tied-multiphasic interfaces --- FEBioMech/FEBioMechPlot.h | 2 +- FEBioMix/FEBioMix.cpp | 2 + FEBioMix/FEBioMixPlot.cpp | 102 ++++++++++++++- FEBioMix/FEBioMixPlot.h | 25 +++- FEBioMix/FEBiphasicContactSurface.cpp | 12 ++ FEBioMix/FEBiphasicContactSurface.h | 28 ++++- FEBioMix/FESlidingInterface2.cpp | 62 +++------ FEBioMix/FESlidingInterface2.h | 18 --- FEBioMix/FESlidingInterface3.cpp | 126 +++++++++---------- FEBioMix/FESlidingInterfaceBiphasic.cpp | 103 ++++----------- FEBioMix/FESlidingInterfaceBiphasic.h | 26 ---- FEBioMix/FESlidingInterfaceBiphasicMixed.cpp | 84 ++++--------- FEBioMix/FESlidingInterfaceBiphasicMixed.h | 24 ---- FEBioMix/FESlidingInterfaceMP.cpp | 108 +++++----------- FEBioMix/FESlidingInterfaceMP.h | 43 +++---- FEBioMix/FETiedBiphasicInterface.cpp | 47 +++---- FEBioMix/FETiedBiphasicInterface.h | 44 ++++--- FEBioMix/FETiedMultiphasicInterface.cpp | 56 +++------ FEBioMix/FETiedMultiphasicInterface.h | 42 +++---- 19 files changed, 417 insertions(+), 537 deletions(-) diff --git a/FEBioMech/FEBioMechPlot.h b/FEBioMech/FEBioMechPlot.h index 18a529e1b..9bbf372b7 100644 --- a/FEBioMech/FEBioMechPlot.h +++ b/FEBioMech/FEBioMechPlot.h @@ -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); }; diff --git a/FEBioMix/FEBioMix.cpp b/FEBioMix/FEBioMix.cpp index c55eb22ab..915eb86e7 100644 --- a/FEBioMix/FEBioMix.cpp +++ b/FEBioMix/FEBioMix.cpp @@ -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" ); @@ -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" ); diff --git a/FEBioMix/FEBioMixPlot.cpp b/FEBioMix/FEBioMixPlot.cpp index 276a09363..dccb91ae6 100644 --- a/FEBioMix/FEBioMixPlot.cpp +++ b/FEBioMix/FEBioMixPlot.cpp @@ -37,11 +37,15 @@ SOFTWARE.*/ #include "FEMultiphasicSolidDomain.h" #include "FEMultiphasicShellDomain.h" #include +#include #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 @@ -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(&surf); + if (pcs == 0) return false; + + writeAverageElementValue(surf, a, [](const FEMaterialPoint& mp) { + const FEContactMaterialPoint* pt = dynamic_cast(&mp); + double d = (pt ? pt->m_gap : 0); + if (d == 0) { + const FETiedBiphasicContactPoint* pd = dynamic_cast(&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(&surf); if (pcs == 0) return false; - writeNodalProjectedElementValues(surf, a, [](const FEMaterialPoint& mp) { - const FEBiphasicContactPoint* pt = mp.ExtractData(); + writeAverageElementValue(surf, a, [](const FEMaterialPoint& mp) { + const FEBiphasicContactPoint* pt = dynamic_cast(&mp); return (pt ? pt->m_pg : 0); }); @@ -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 s; + for (int i = 0; i(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(&surf); + if (pcs == 0) return false; + + for (int i=0; iGetMaterial(se->GetMatID()); + FESoluteInterface* pm = dynamic_cast(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 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(&mp); + const FETiedMultiphasicContactPoint* tt = dynamic_cast(&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 //============================================================================= diff --git a/FEBioMix/FEBioMixPlot.h b/FEBioMix/FEBioMixPlot.h index d1d9e7e27..b7cce2e78 100644 --- a/FEBioMix/FEBioMixPlot.h +++ b/FEBioMix/FEBioMixPlot.h @@ -363,6 +363,15 @@ 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 @@ -370,7 +379,7 @@ class FEPlotFluidForce2 : public FEPlotSurfaceData 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); }; @@ -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 m_sol; +}; + diff --git a/FEBioMix/FEBiphasicContactSurface.cpp b/FEBioMix/FEBiphasicContactSurface.cpp index 0eccacfb1..833aae59d 100644 --- a/FEBioMix/FEBiphasicContactSurface.cpp +++ b/FEBioMix/FEBiphasicContactSurface.cpp @@ -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; } diff --git a/FEBioMix/FEBiphasicContactSurface.h b/FEBioMix/FEBiphasicContactSurface.h index fa3a0eeff..adbb89c27 100644 --- a/FEBioMix/FEBiphasicContactSurface.h +++ b/FEBioMix/FEBiphasicContactSurface.h @@ -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; diff --git a/FEBioMix/FESlidingInterface2.cpp b/FEBioMix/FESlidingInterface2.cpp index d4fa32f89..60e81d412 100644 --- a/FEBioMix/FESlidingInterface2.cpp +++ b/FEBioMix/FESlidingInterface2.cpp @@ -59,30 +59,6 @@ BEGIN_FECORE_CLASS(FESlidingInterface2, FEContactInterface) ADD_PARAMETER(m_bdupr , "dual_proj" ); END_FECORE_CLASS(); -//----------------------------------------------------------------------------- -FESlidingSurface2::Data::Data() -{ - m_Lmd = 0.0; - m_Lmp = 0.0; - m_epsn = 1.0; - m_epsp = 1.0; - m_pg = 0.0; - m_p1 = 0.0; - m_nu = vec3d(0,0,0); - m_rs = vec2d(0,0); -} - -void FESlidingSurface2::Data::Serialize(DumpStream& ar) -{ - FEBiphasicContactPoint::Serialize(ar); - ar & m_Lmd; - ar & m_epsn; - ar & m_epsp; - ar & m_p1; - ar & m_nu; - ar & m_rs; -} - //----------------------------------------------------------------------------- // FESlidingSurface2 //----------------------------------------------------------------------------- @@ -132,7 +108,7 @@ bool FESlidingSurface2::Init() //! create material point data FEMaterialPoint* FESlidingSurface2::CreateMaterialPoint() { - return new FESlidingSurface2::Data; + return new FEBiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -250,7 +226,7 @@ vec3d FESlidingSurface2::GetContactForceFromElementStress() // evaluate the contact force for that element for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors vec3d g[2]; @@ -286,7 +262,7 @@ double FESlidingSurface2::GetContactArea() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors @@ -326,7 +302,7 @@ vec3d FESlidingSurface2::GetFluidForce() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors vec3d g[2]; @@ -377,7 +353,7 @@ vec3d FESlidingSurface2::GetFluidForceFromElementPressure() // evaluate the fluid force for that element for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors vec3d g[2]; @@ -423,7 +399,7 @@ void FESlidingSurface2::GetContactPressure(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_Ln; } pg /= ni; @@ -437,7 +413,7 @@ void FESlidingSurface2::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt -= data.m_nu*data.m_Ln; } pt /= ni; @@ -545,7 +521,7 @@ void FESlidingInterface2::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = pt.m_pme; if (pe != 0) { @@ -630,7 +606,7 @@ void FESlidingInterface2::CalcAutoPenalty(FESlidingSurface2& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -652,7 +628,7 @@ void FESlidingInterface2::CalcAutoPressurePenalty(FESlidingSurface2& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -800,7 +776,7 @@ void FESlidingInterface2::ProjectSurface(FESlidingSurface2& ss, FESlidingSurface for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point r = ss.Local2Global(el, j); @@ -1087,7 +1063,7 @@ void FESlidingInterface2::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) for (j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -1308,7 +1284,7 @@ void FESlidingInterface2::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& for (j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -1696,7 +1672,7 @@ void FESlidingInterface2::UpdateContactPressures() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(i)); gap = pt.m_gap; eps = m_epsn*pt.m_epsn; @@ -1707,7 +1683,7 @@ void FESlidingInterface2::UpdateContactPressures() int mint = pme->GaussPoints(); double ti[FEElement::MAX_NODES]; for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& mp = static_cast(*el.GetMaterialPoint(j)); gap = mp.m_gap; eps = m_epsn*mp.m_epsn; ti[j] = MBRACKET(mp.m_Lmd + m_epsn*mp.m_epsn*mp.m_gap); @@ -1746,7 +1722,7 @@ bool FESlidingInterface2::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); normL0 += ds.m_Lmd*ds.m_Lmd; } } @@ -1755,7 +1731,7 @@ bool FESlidingInterface2::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); normL0 += dm.m_Lmd*dm.m_Lmd; } } @@ -1772,7 +1748,7 @@ bool FESlidingInterface2::Augment(int naug, const FETimeInfo& tp) vec3d tn[FEElement::MAX_INTPOINTS]; if (m_bsmaug) m_ss.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface if (m_bsmaug) { // replace this multiplier with a smoother version @@ -1808,7 +1784,7 @@ bool FESlidingInterface2::Augment(int naug, const FETimeInfo& tp) vec3d tn[FEElement::MAX_INTPOINTS]; if (m_bsmaug) m_ms.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface if (m_bsmaug) { // replace this multiplier with a smoother version diff --git a/FEBioMix/FESlidingInterface2.h b/FEBioMix/FESlidingInterface2.h index 7b8441f06..43251dc04 100644 --- a/FEBioMix/FESlidingInterface2.h +++ b/FEBioMix/FESlidingInterface2.h @@ -33,24 +33,6 @@ SOFTWARE.*/ //----------------------------------------------------------------------------- class FEBIOMIX_API FESlidingSurface2 : public FEBiphasicContactSurface { -public: - //! Integration point data - class Data : public FEBiphasicContactPoint - { - public: - Data(); - - void Serialize(DumpStream& ar) override; - - public: - double m_Lmd; //!< lagrange multipliers for displacement - double m_epsn; //!< penalty factor - double m_epsp; //!< pressure penatly factor - double m_p1; //!< fluid pressure - vec3d m_nu; //!< normal at integration points - vec2d m_rs; //!< natrual coordinates of projection - }; - public: //! constructor FESlidingSurface2(FEModel* pfem); diff --git a/FEBioMix/FESlidingInterface3.cpp b/FEBioMix/FESlidingInterface3.cpp index f4b29c17f..4cca61c30 100644 --- a/FEBioMix/FESlidingInterface3.cpp +++ b/FEBioMix/FESlidingInterface3.cpp @@ -28,6 +28,7 @@ SOFTWARE.*/ #include "stdafx.h" #include "FESlidingInterface3.h" +#include "FESlidingInterfaceMP.h" #include "FEBiphasic.h" #include "FEBiphasicSolute.h" #include "FECore/FEModel.h" @@ -64,35 +65,6 @@ BEGIN_FECORE_CLASS(FESlidingInterface3, FEContactInterface) ADD_PARAMETER(m_ambc , "ambient_concentration"); END_FECORE_CLASS(); -//----------------------------------------------------------------------------- -FESlidingSurface3::Data::Data() -{ - m_nu = vec3d(0,0,0); - m_rs = vec2d(0,0); - m_Lmd = 0.0; - m_Lmc = 0.0; - m_Lmp = 0.0; - m_epsn = 1.0; - m_epsc = 1.0; - m_epsp = 1.0; - m_pg = 0.0; - m_cg = 0.0; -} - -//----------------------------------------------------------------------------- -void FESlidingSurface3::Data::Serialize(DumpStream& ar) -{ - FEBiphasicContactPoint::Serialize(ar); - ar & m_Lmd; - ar & m_Lmc; - ar & m_epsn; - ar & m_epsp; - ar & m_epsc; - ar & m_cg; - ar & m_nu; - ar & m_rs; -} - //----------------------------------------------------------------------------- // FESlidingSurface3 //----------------------------------------------------------------------------- @@ -113,7 +85,7 @@ FESlidingSurface3::~FESlidingSurface3() //! create material point data FEMaterialPoint* FESlidingSurface3::CreateMaterialPoint() { - return new FESlidingSurface3::Data; + return new FEMultiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -190,6 +162,28 @@ bool FESlidingSurface3::Init() } } } + + // allocate data structures + int NE = Elements(); + for (int i=0; i(*el.GetMaterialPoint(j)); + data.m_Lmc.resize(nsol); + data.m_epsc.resize(nsol); + data.m_cg.resize(nsol); + data.m_c1.resize(nsol); + data.m_epsc.assign(nsol,1); + data.m_c1.assign(nsol,0.0); + data.m_cg.assign(nsol,0.0); + data.m_Lmc.assign(nsol,0.0); + } + } + } return true; } @@ -292,7 +286,7 @@ double FESlidingSurface3::GetContactArea() for (int i=0; i(*el.GetMaterialPoint(i)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); double s = (data.m_Ln > 0) ? 1 : 0; // get the base vectors @@ -381,7 +375,7 @@ void FESlidingSurface3::GetContactPressure(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_Ln; } pg /= ni; @@ -395,7 +389,7 @@ void FESlidingSurface3::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt -= data.m_nu*data.m_Ln; } pt /= ni; @@ -514,7 +508,7 @@ void FESlidingInterface3::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = data.m_pme; if (pe != 0) @@ -605,7 +599,7 @@ void FESlidingInterface3::CalcAutoPenalty(FESlidingSurface3& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -627,7 +621,7 @@ void FESlidingInterface3::CalcAutoPressurePenalty(FESlidingSurface3& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -728,8 +722,8 @@ void FESlidingInterface3::CalcAutoConcentrationPenalty(FESlidingSurface3& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); - pt.m_epsc = eps; + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); + pt.m_epsc[0] = eps; } } } @@ -889,7 +883,7 @@ void FESlidingInterface3::ProjectSurface(FESlidingSurface3& ss, FESlidingSurface for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point r = ss.Local2Global(el, j); @@ -962,7 +956,7 @@ void FESlidingInterface3::ProjectSurface(FESlidingSurface3& ss, FESlidingSurface double cm[FEElement::MAX_NODES]; for (int k=0; kNodes(); ++k) cm[k] = mesh.Node(pme->m_node[k]).get(m_dofC + sid); double c2 = pme->eval(cm, rs[0], rs[1]); - pt.m_cg = c1 - c2; + pt.m_cg[0] = c1 - c2; } } else @@ -975,8 +969,8 @@ void FESlidingInterface3::ProjectSurface(FESlidingSurface3& ss, FESlidingSurface pt.m_pg = 0; } if (ssolu) { - pt.m_Lmc = 0; - pt.m_cg = 0; + pt.m_Lmc[0] = 0; + pt.m_cg[0] = 0; } } } @@ -990,8 +984,8 @@ void FESlidingInterface3::ProjectSurface(FESlidingSurface3& ss, FESlidingSurface pt.m_pg = 0; } if (ssolu) { - pt.m_Lmc = 0; - pt.m_cg = 0; + pt.m_Lmc[0] = 0; + pt.m_cg[0] = 0; } } } @@ -1207,7 +1201,7 @@ void FESlidingInterface3::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) // note that we are integrating over the current surface for (int j = 0; j(*se.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; if (pme) @@ -1341,9 +1335,9 @@ void FESlidingInterface3::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) int ndof = nseln + nmeln; // calculate the flow rate - double epsc = m_epsc*pt.m_epsc*psf; + double epsc = m_epsc*pt.m_epsc[0]*psf; - double jn = pt.m_Lmc + epsc*pt.m_cg; + double jn = pt.m_Lmc[0] + epsc*pt.m_cg[0]; // fill the LM LM.resize(ndof); @@ -1468,7 +1462,7 @@ void FESlidingInterface3::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& // loop over all integration points for (j=0; j(*se.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -1793,7 +1787,7 @@ void FESlidingInterface3::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& double dt = fem.GetTime().timeIncrement; double epsp = (tn > 0) ? m_epsp*pt.m_epsp*psf : 0.; - double epsc = (tn > 0) ? m_epsc*pt.m_epsc*psf : 0.; + double epsc = (tn > 0) ? m_epsc*pt.m_epsc[0]*psf : 0.; // --- S O L I D - P R E S S U R E / S O L U T E C O N T A C T --- @@ -1824,7 +1818,7 @@ void FESlidingInterface3::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& } double wn = pt.m_Lmp + epsp*pt.m_pg; - double jn = pt.m_Lmc + epsc*pt.m_cg; + double jn = pt.m_Lmc[0] + epsc*pt.m_cg[0]; // b. A-term //------------------------------------- @@ -1986,7 +1980,7 @@ void FESlidingInterface3::UpdateContactPressures() double gap, eps; for (int i=0; i(*el.GetMaterialPoint(i)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(i)); gap = pt.m_gap; eps = m_epsn*pt.m_epsn; @@ -1997,7 +1991,7 @@ void FESlidingInterface3::UpdateContactPressures() int mint = pme->GaussPoints(); double ti[FEElement::MAX_NODES]; for (int j=0; j(*pme->GetMaterialPoint(j)); + FEMultiphasicContactPoint& md = static_cast(*pme->GetMaterialPoint(j)); gap = md.m_gap; eps = m_epsn*md.m_epsn; @@ -2039,7 +2033,7 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); normL0 += ds.m_Lmd*ds.m_Lmd; } } @@ -2048,7 +2042,7 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); normL0 += dm.m_Lmd*dm.m_Lmd; } } @@ -2064,7 +2058,7 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) vec3d tn[FEElement::MAX_INTPOINTS]; if (m_bsmaug) m_ss.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface if (m_bsmaug) { // replace this multiplier with a smoother version @@ -2094,12 +2088,12 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) if (m_ss.m_bsolu) { Lc = 0; if (Ln > 0) { - epsc = m_epsc*data.m_epsc; - Lc = data.m_Lmc + epsc*data.m_cg; - maxcg = max(maxcg,fabs(data.m_cg)); - normDC += data.m_cg*data.m_cg; + epsc = m_epsc*data.m_epsc[0]; + Lc = data.m_Lmc[0] + epsc*data.m_cg[0]; + maxcg = max(maxcg,fabs(data.m_cg[0])); + normDC += data.m_cg[0]*data.m_cg[0]; } - data.m_Lmc = Lc; + data.m_Lmc[0] = Lc; } if (Ln > 0) maxgap = max(maxgap,fabs(data.m_gap)); @@ -2111,7 +2105,7 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) vec3d tn[FEElement::MAX_INTPOINTS]; if (m_bsmaug) m_ms.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface if (m_bsmaug) { // replace this multiplier with a smoother version @@ -2141,12 +2135,12 @@ bool FESlidingInterface3::Augment(int naug, const FETimeInfo& tp) if (m_ms.m_bsolu) { Lc = 0; if (Ln > 0) { - epsc = m_epsc*data.m_epsc; - Lc = data.m_Lmc + epsc*data.m_cg; - maxcg = max(maxcg,fabs(data.m_cg)); - normDC += data.m_cg*data.m_cg; + epsc = m_epsc*data.m_epsc[0]; + Lc = data.m_Lmc[0] + epsc*data.m_cg[0]; + maxcg = max(maxcg,fabs(data.m_cg[0])); + normDC += data.m_cg[0]*data.m_cg[0]; } - data.m_Lmc = Lc; + data.m_Lmc[0] = Lc; } if (Ln > 0) maxgap = max(maxgap,fabs(data.m_gap)); diff --git a/FEBioMix/FESlidingInterfaceBiphasic.cpp b/FEBioMix/FESlidingInterfaceBiphasic.cpp index 3b91c1036..6441db46d 100644 --- a/FEBioMix/FESlidingInterfaceBiphasic.cpp +++ b/FEBioMix/FESlidingInterfaceBiphasic.cpp @@ -65,59 +65,6 @@ BEGIN_FECORE_CLASS(FESlidingInterfaceBiphasic, FEContactInterface) ADD_PARAMETER(m_bshellbm , "shell_bottom_secondary"); END_FECORE_CLASS(); -//----------------------------------------------------------------------------- -FESlidingSurfaceBiphasic::Data::Data() -{ - 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; -} - -void FESlidingSurfaceBiphasic::Data::Init() -{ - FEBiphasicContactPoint::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; -} - -//----------------------------------------------------------------------------- -void FESlidingSurfaceBiphasic::Data::Serialize(DumpStream& ar) -{ - FEBiphasicContactPoint::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; -} - //----------------------------------------------------------------------------- // FESlidingSurfaceBiphasic //----------------------------------------------------------------------------- @@ -131,7 +78,7 @@ FESlidingSurfaceBiphasic::FESlidingSurfaceBiphasic(FEModel* pfem) : FEBiphasicCo //! create material point data FEMaterialPoint* FESlidingSurfaceBiphasic::CreateMaterialPoint() { - return new FESlidingSurfaceBiphasic::Data; + return new FEBiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -180,7 +127,7 @@ void FESlidingSurfaceBiphasic::InitSlidingSurface() int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // Store current surface projection values as previous data.m_rsp = data.m_rs; data.m_pmep = data.m_pme; @@ -328,7 +275,7 @@ double FESlidingSurfaceBiphasic::GetContactArea() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors @@ -411,7 +358,7 @@ void FESlidingSurfaceBiphasic::GetVectorGap(int nface, vec3d& pg) pg = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_dg; } pg /= ni; @@ -425,7 +372,7 @@ void FESlidingSurfaceBiphasic::GetContactPressure(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_Ln; } pg /= ni; @@ -439,7 +386,7 @@ void FESlidingSurfaceBiphasic::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt += data.m_tr; } pt /= ni; @@ -453,7 +400,7 @@ void FESlidingSurfaceBiphasic::GetSlipTangent(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (!data.m_bstick) pt += data.m_s1; } pt /= ni; @@ -467,7 +414,7 @@ void FESlidingSurfaceBiphasic::GetMuEffective(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_mueff; } pg /= ni; @@ -481,7 +428,7 @@ void FESlidingSurfaceBiphasic::GetLocalFLS(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_fls; } pg /= ni; @@ -495,7 +442,7 @@ void FESlidingSurfaceBiphasic::GetNodalVectorGap(int nface, vec3d* pg) vec3d gi[FEElement::MAX_INTPOINTS]; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); gi[k] = data.m_dg; } el.project_to_nodes(gi, pg); @@ -517,7 +464,7 @@ void FESlidingSurfaceBiphasic::GetStickStatus(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (data.m_bstick) pg += 1.0; } pg /= ni; @@ -632,7 +579,7 @@ void FESlidingInterfaceBiphasic::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = pt.m_pme; if (pe != 0) { @@ -717,7 +664,7 @@ void FESlidingInterfaceBiphasic::CalcAutoPenalty(FESlidingSurfaceBiphasic& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -739,7 +686,7 @@ void FESlidingInterfaceBiphasic::CalcAutoPressurePenalty(FESlidingSurfaceBiphasi int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -882,7 +829,7 @@ void FESlidingInterfaceBiphasic::ProjectSurface(FESlidingSurfaceBiphasic& ss, FE for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point vec3d r = ss.Local2Global(el, j); @@ -1139,7 +1086,7 @@ vec3d FESlidingInterfaceBiphasic::SlipTangent(FESlidingSurfaceBiphasic& ss, cons FESurfaceElement& se = ss.Element(nel); // get integration point data - FESlidingSurfaceBiphasic::Data& data = static_cast(*se.GetMaterialPoint(nint)); + FEBiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(nint)); double g = data.m_gap; vec3d nu = data.m_nu; @@ -1206,7 +1153,7 @@ vec3d FESlidingInterfaceBiphasic::ContactTraction(FESlidingSurfaceBiphasic& ss, FESurfaceElement& se = ss.Element(nel); // get the integration point data - FESlidingSurfaceBiphasic::Data& data = static_cast(*se.GetMaterialPoint(n)); + FEBiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(n)); // penalty double eps = m_epsn*data.m_epsn*psf; @@ -1519,7 +1466,7 @@ void FESlidingInterfaceBiphasic::LoadVector(FEGlobalVector& R, const FETimeInfo& for (int j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact pressure and account for stick double pn; @@ -1706,7 +1653,7 @@ void FESlidingInterfaceBiphasic::StiffnessMatrix(FELinearSystem& LS, const FETim for (int j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact pressure and account for stick double pn; @@ -2308,7 +2255,7 @@ void FESlidingInterfaceBiphasic::UpdateContactPressures() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& sd = static_cast(*el.GetMaterialPoint(i)); // evaluate traction on primary surface double eps = m_epsn*sd.m_epsn*psf; if (sd.m_bstick) { @@ -2336,7 +2283,7 @@ void FESlidingInterfaceBiphasic::UpdateContactPressures() vec3d ti[MI]; for (int j=0; j(*pme->GetMaterialPoint(j)); + FEBiphasicContactPoint& md = static_cast(*pme->GetMaterialPoint(j)); // evaluate traction on secondary surface double eps = m_epsn*md.m_epsn*psf; @@ -2397,7 +2344,7 @@ bool FESlidingInterfaceBiphasic::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); if (ds.m_bstick) normL0 += ds.m_Lmt*ds.m_Lmt; else @@ -2409,7 +2356,7 @@ bool FESlidingInterfaceBiphasic::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); if (dm.m_bstick) normL0 += dm.m_Lmt*dm.m_Lmt; else @@ -2431,7 +2378,7 @@ bool FESlidingInterfaceBiphasic::Augment(int naug, const FETimeInfo& tp) if (m_bsmaug) m_ss.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface double eps = m_epsn*ds.m_epsn*psf; @@ -2493,7 +2440,7 @@ bool FESlidingInterfaceBiphasic::Augment(int naug, const FETimeInfo& tp) if (m_bsmaug) m_ms.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface double eps = m_epsn*dm.m_epsn*psf; diff --git a/FEBioMix/FESlidingInterfaceBiphasic.h b/FEBioMix/FESlidingInterfaceBiphasic.h index 27a4998e4..ab34abce0 100644 --- a/FEBioMix/FESlidingInterfaceBiphasic.h +++ b/FEBioMix/FESlidingInterfaceBiphasic.h @@ -33,32 +33,6 @@ SOFTWARE.*/ //----------------------------------------------------------------------------- class FEBIOMIX_API FESlidingSurfaceBiphasic : public FEBiphasicContactSurface { -public: - //! Integration point data - class Data : public FEBiphasicContactPoint - { - public: - Data(); - - void Init() override; - - void Serialize(DumpStream& ar) override; - - public: - 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 - }; - public: //! constructor FESlidingSurfaceBiphasic(FEModel* pfem); diff --git a/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp b/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp index 9bfc3a152..6cc3e99cb 100644 --- a/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp +++ b/FEBioMix/FESlidingInterfaceBiphasicMixed.cpp @@ -67,40 +67,6 @@ BEGIN_FECORE_CLASS(FESlidingInterfaceBiphasicMixed, FEContactInterface) ADD_PARAMETER(m_bshellbm , "shell_bottom_secondary"); END_FECORE_CLASS(); -//----------------------------------------------------------------------------- -FESlidingSurfaceBiphasicMixed::Data::Data() -{ - 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_nu = m_s1 = m_dg = vec3d(0,0,0); - m_rs = m_rsp = vec2d(0,0); - m_bstick = false; -} - -//----------------------------------------------------------------------------- -void FESlidingSurfaceBiphasicMixed::Data::Serialize(DumpStream& ar) -{ - FEBiphasicContactPoint::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; -} - //----------------------------------------------------------------------------- // FESlidingSurfaceBiphasic //----------------------------------------------------------------------------- @@ -114,7 +80,7 @@ FESlidingSurfaceBiphasicMixed::FESlidingSurfaceBiphasicMixed(FEModel* pfem) : FE //! create material point data FEMaterialPoint* FESlidingSurfaceBiphasicMixed::CreateMaterialPoint() { - return new FESlidingSurfaceBiphasicMixed::Data; + return new FEBiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -168,7 +134,7 @@ void FESlidingSurfaceBiphasicMixed::InitSlidingSurface() int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); // Store current surface projection values as previous data.m_rsp = data.m_rs; data.m_pmep = data.m_pme; @@ -316,7 +282,7 @@ double FESlidingSurfaceBiphasicMixed::GetContactArea() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); if (data.m_Ln > 0) { // get the base vectors @@ -403,7 +369,7 @@ void FESlidingSurfaceBiphasicMixed::GetVectorGap(int nface, vec3d& pg) pg = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_dg; } pg /= ni; @@ -417,7 +383,7 @@ void FESlidingSurfaceBiphasicMixed::GetContactPressure(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_Ln; } pg /= ni; @@ -431,7 +397,7 @@ void FESlidingSurfaceBiphasicMixed::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt += data.m_tr; } pt /= ni; @@ -445,7 +411,7 @@ void FESlidingSurfaceBiphasicMixed::GetSlipTangent(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (!data.m_bstick) pt += data.m_s1; } pt /= ni; @@ -459,7 +425,7 @@ void FESlidingSurfaceBiphasicMixed::GetMuEffective(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_mueff; } pg /= ni; @@ -473,7 +439,7 @@ void FESlidingSurfaceBiphasicMixed::GetLocalFLS(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_fls; } pg /= ni; @@ -487,7 +453,7 @@ void FESlidingSurfaceBiphasicMixed::GetNodalVectorGap(int nface, vec3d* pg) vec3d gi[FEElement::MAX_INTPOINTS]; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); gi[k] = data.m_dg; } el.project_to_nodes(gi, pg); @@ -513,7 +479,7 @@ void FESlidingSurfaceBiphasicMixed::GetStickStatus(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (data.m_bstick) pg += 1.0; } pg /= ni; @@ -660,7 +626,7 @@ void FESlidingInterfaceBiphasicMixed::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = pt.m_pme; if (pe != 0) { @@ -745,7 +711,7 @@ void FESlidingInterfaceBiphasicMixed::CalcAutoPenalty(FESlidingSurfaceBiphasicMi int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -767,7 +733,7 @@ void FESlidingInterfaceBiphasicMixed::CalcAutoPressurePenalty(FESlidingSurfaceBi int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -915,7 +881,7 @@ void FESlidingInterfaceBiphasicMixed::ProjectSurface(FESlidingSurfaceBiphasicMix for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point vec3d r = ss.Local2Global(el, j); @@ -1177,7 +1143,7 @@ vec3d FESlidingInterfaceBiphasicMixed::SlipTangent(FESlidingSurfaceBiphasicMixed FESurfaceElement& se = ss.Element(nel); // get integration point data - FESlidingSurfaceBiphasicMixed::Data& data = static_cast(*se.GetMaterialPoint(nint)); + FEBiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(nint)); double g = data.m_gap; vec3d nu = data.m_nu; @@ -1248,7 +1214,7 @@ vec3d FESlidingInterfaceBiphasicMixed::ContactTraction(FESlidingSurfaceBiphasicM FESurfaceElement& se = ss.Element(nel); // get the integration point data - FESlidingSurfaceBiphasicMixed::Data& data = static_cast(*se.GetMaterialPoint(n)); + FEBiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(n)); // penalty double eps = m_epsn*data.m_epsn*psf; @@ -1566,7 +1532,7 @@ void FESlidingInterfaceBiphasicMixed::LoadVector(FESlidingSurfaceBiphasicMixed& for (int j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact pressure and account for stick double pn; @@ -1779,7 +1745,7 @@ void FESlidingInterfaceBiphasicMixed::StiffnessMatrix(FESlidingSurfaceBiphasicMi for (int j=0; j(*se.GetMaterialPoint(j)); + FEBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact pressure and account for stick double pn; @@ -2402,7 +2368,7 @@ void FESlidingInterfaceBiphasicMixed::UpdateContactPressures() for (int i=0; i(*el.GetMaterialPoint(i)); + FEBiphasicContactPoint& sd = static_cast(*el.GetMaterialPoint(i)); // evaluate traction on primary surface double eps = m_epsn*sd.m_epsn*psf; if (sd.m_bstick) { @@ -2430,7 +2396,7 @@ void FESlidingInterfaceBiphasicMixed::UpdateContactPressures() vec3d ti[MI]; for (int j=0; j(*pme->GetMaterialPoint(j)); + FEBiphasicContactPoint& md = static_cast(*pme->GetMaterialPoint(j)); // evaluate traction on secondary surface double eps = m_epsn*md.m_epsn*psf; @@ -2491,7 +2457,7 @@ bool FESlidingInterfaceBiphasicMixed::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); if (ds.m_bstick) normL0 += ds.m_Lmt*ds.m_Lmt; else @@ -2503,7 +2469,7 @@ bool FESlidingInterfaceBiphasicMixed::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); if (dm.m_bstick) normL0 += dm.m_Lmt*dm.m_Lmt; else @@ -2525,7 +2491,7 @@ bool FESlidingInterfaceBiphasicMixed::Augment(int naug, const FETimeInfo& tp) if (m_bsmaug) m_ss.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface double eps = m_epsn*ds.m_epsn*psf; @@ -2587,7 +2553,7 @@ bool FESlidingInterfaceBiphasicMixed::Augment(int naug, const FETimeInfo& tp) if (m_bsmaug) m_ms.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface double eps = m_epsn*dm.m_epsn*psf; diff --git a/FEBioMix/FESlidingInterfaceBiphasicMixed.h b/FEBioMix/FESlidingInterfaceBiphasicMixed.h index 0d883431f..770dda249 100644 --- a/FEBioMix/FESlidingInterfaceBiphasicMixed.h +++ b/FEBioMix/FESlidingInterfaceBiphasicMixed.h @@ -33,30 +33,6 @@ SOFTWARE.*/ //----------------------------------------------------------------------------- class FEBIOMIX_API FESlidingSurfaceBiphasicMixed : public FEBiphasicContactSurface { -public: - //! Integration point data - class Data : public FEBiphasicContactPoint - { - public: - Data(); - - void Serialize(DumpStream& ar) override; - - public: - 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 - }; - public: //! constructor FESlidingSurfaceBiphasicMixed(FEModel* pfem); diff --git a/FEBioMix/FESlidingInterfaceMP.cpp b/FEBioMix/FESlidingInterfaceMP.cpp index 2fad95816..eeb466c77 100644 --- a/FEBioMix/FESlidingInterfaceMP.cpp +++ b/FEBioMix/FESlidingInterfaceMP.cpp @@ -83,65 +83,13 @@ BEGIN_FECORE_CLASS(FESlidingInterfaceMP, FEContactInterface) END_FECORE_CLASS(); //----------------------------------------------------------------------------- -FESlidingSurfaceMP::Data::Data() -{ - 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; -} - -//----------------------------------------------------------------------------- -void FESlidingSurfaceMP::Data::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; -} - -//----------------------------------------------------------------------------- -void FESlidingSurfaceMP::Data::Serialize(DumpStream& ar) +void FEMultiphasicContactPoint::Serialize(DumpStream& ar) { FEBiphasicContactPoint::Serialize(ar); - ar & m_gap; - ar & m_dg; - ar & m_nu; - ar & m_s1; - ar & m_rs; - ar & m_rsp; - ar & m_Lmd; - ar & m_Lmt; - ar & m_Lmp; ar & m_Lmc; - ar & m_epsn; - ar & m_epsp; ar & m_epsc; - ar & m_pg; ar & m_cg; - ar & m_p1; ar & m_c1; - ar & m_Ln; - ar & m_bstick; - ar & m_tr; - ar & m_mueff; - ar & m_fls; } //----------------------------------------------------------------------------- @@ -164,7 +112,7 @@ FESlidingSurfaceMP::~FESlidingSurfaceMP() //! create material point data FEMaterialPoint* FESlidingSurfaceMP::CreateMaterialPoint() { - return new FESlidingSurfaceMP::Data; + return new FEMultiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -267,7 +215,7 @@ bool FESlidingSurfaceMP::Init() int nint = el.GaussPoints(); if (nsol) { for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); data.m_Lmc.resize(nsol); data.m_epsc.resize(nsol); data.m_cg.resize(nsol); @@ -293,7 +241,7 @@ void FESlidingSurfaceMP::InitSlidingSurface() for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); data.m_rsp = data.m_rs; data.m_pmep = data.m_pme; } @@ -440,7 +388,7 @@ double FESlidingSurfaceMP::GetContactArea() for (int i=0; i(*el.GetMaterialPoint(i)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(i)); double s = (data.m_Ln > 0) ? 1 : 0; // get the base vectors @@ -528,7 +476,7 @@ void FESlidingSurfaceMP::GetVectorGap(int nface, vec3d& pg) pg = vec3d(0,0,0); for (int k=0; k(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_dg; } pg /= ni; @@ -542,7 +490,7 @@ void FESlidingSurfaceMP::GetContactPressure(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_Ln; } pg /= ni; @@ -556,7 +504,7 @@ void FESlidingSurfaceMP::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt += data.m_tr; } pt /= ni; @@ -570,7 +518,7 @@ void FESlidingSurfaceMP::GetSlipTangent(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k=0; k(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (!data.m_bstick) pt += data.m_s1; } pt /= ni; @@ -584,7 +532,7 @@ void FESlidingSurfaceMP::GetMuEffective(int nface, double& pg) pg = 0; for (int k=0; k(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_mueff; } pg /= ni; @@ -598,7 +546,7 @@ void FESlidingSurfaceMP::GetLocalFLS(int nface, double& pg) pg = 0; for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_fls; } pg /= ni; @@ -612,7 +560,7 @@ void FESlidingSurfaceMP::GetNodalVectorGap(int nface, vec3d* pg) vec3d gi[FEElement::MAX_INTPOINTS]; for (int k=0; k(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); gi[k] = data.m_dg; } el.project_to_nodes(gi, pg); @@ -634,7 +582,7 @@ void FESlidingSurfaceMP::GetStickStatus(int nface, double& pg) pg = 0; for (int k=0; k(*el.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); if (data.m_bstick) pg += 1.0; } pg /= ni; @@ -775,7 +723,7 @@ void FESlidingInterfaceMP::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FEMultiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = data.m_pme; if (pe != 0) @@ -864,7 +812,7 @@ void FESlidingInterfaceMP::CalcAutoPenalty(FESlidingSurfaceMP& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -886,7 +834,7 @@ void FESlidingInterfaceMP::CalcAutoPressurePenalty(FESlidingSurfaceMP& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -1024,7 +972,7 @@ void FESlidingInterfaceMP::CalcAutoConcentrationPenalty(FESlidingSurfaceMP& s, int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsc[isol] = eps; } } @@ -1190,7 +1138,7 @@ void FESlidingInterfaceMP::ProjectSurface(FESlidingSurfaceMP& ss, FESlidingSurfa for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point r = ss.Local2Global(el, j); @@ -1492,7 +1440,7 @@ vec3d FESlidingInterfaceMP::SlipTangent(FESlidingSurfaceMP& ss, const int nel, c FESurfaceElement& se = ss.Element(nel); // get integration point data - FESlidingSurfaceMP::Data& data = static_cast(*se.GetMaterialPoint(nint)); + FEMultiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(nint)); double g = data.m_gap; vec3d nu = data.m_nu; @@ -1563,7 +1511,7 @@ vec3d FESlidingInterfaceMP::ContactTraction(FESlidingSurfaceMP& ss, const int ne FESurfaceElement& se = ss.Element(nel); // get the integration point data - FESlidingSurfaceMP::Data& data = static_cast(*se.GetMaterialPoint(n)); + FEMultiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(n)); // penalty double eps = m_epsn*data.m_epsn*psf; @@ -1898,7 +1846,7 @@ void FESlidingInterfaceMP::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) for (int j=0; j(*se.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact pressure and account for stick double pn; @@ -2136,7 +2084,7 @@ void FESlidingInterfaceMP::StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& for (j=0; j(*se.GetMaterialPoint(j)); + FEMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // calculate contact traction and account for stick double pn; @@ -2794,7 +2742,7 @@ void FESlidingInterfaceMP::UpdateContactPressures() for (i=0; i(*el.GetMaterialPoint(i)); + FEMultiphasicContactPoint& sd = static_cast(*el.GetMaterialPoint(i)); // evaluate traction on primary surface double eps = m_epsn*sd.m_epsn*psf; if (sd.m_bstick) { @@ -2819,7 +2767,7 @@ void FESlidingInterfaceMP::UpdateContactPressures() vec3d ti[MI]; double pi[MI]; for (j=0; j(*pme->GetMaterialPoint(j)); + FEMultiphasicContactPoint& md = static_cast(*pme->GetMaterialPoint(j)); // evaluate traction on secondary surface double eps = m_epsn*md.m_epsn*psf; @@ -2884,7 +2832,7 @@ bool FESlidingInterfaceMP::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); if (ds.m_bstick) normL0 += ds.m_Lmt*ds.m_Lmt; else @@ -2896,7 +2844,7 @@ bool FESlidingInterfaceMP::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); if (dm.m_bstick) normL0 += dm.m_Lmt*dm.m_Lmt; else @@ -2917,7 +2865,7 @@ bool FESlidingInterfaceMP::Augment(int naug, const FETimeInfo& tp) if (m_bsmaug) m_ss.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface eps = m_epsn*ds.m_epsn*psf; @@ -2986,7 +2934,7 @@ bool FESlidingInterfaceMP::Augment(int naug, const FETimeInfo& tp) vec3d tn[FEElement::MAX_INTPOINTS]; if (m_bsmaug) m_ms.GetGPSurfaceTraction(i, tn); for (int j=0; j(*el.GetMaterialPoint(j)); + FEMultiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on master surface double eps = m_epsn*dm.m_epsn*psf; diff --git a/FEBioMix/FESlidingInterfaceMP.h b/FEBioMix/FESlidingInterfaceMP.h index ae50da89f..7f7ea477b 100644 --- a/FEBioMix/FESlidingInterfaceMP.h +++ b/FEBioMix/FESlidingInterfaceMP.h @@ -34,38 +34,25 @@ SOFTWARE.*/ #include //----------------------------------------------------------------------------- -class FEBIOMIX_API FESlidingSurfaceMP : public FEBiphasicContactSurface +class FEBIOMIX_API FEMultiphasicContactPoint : public FEBiphasicContactPoint { public: - //! integration point data - class Data : public FEBiphasicContactPoint + vector m_Lmc; //!< Lagrange multipliers for solute concentrations + vector m_epsc; //!< concentration penalty factors + vector m_cg; //!< concentration "gap" + vector m_c1; //!< solute concentration + + void Init() override { - public: - Data(); - - void Init() override; - - void Serialize(DumpStream& ar) override; - - public: - vec3d m_dg; //!< vector gap - double m_Lmd; //!< Lagrange multipliers for displacements - vec3d m_Lmt; //!< Lagrange multipliers for vector traction - double m_epsn; //!< displacement penalty factors - double m_epsp; //!< pressure penalty factors - 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 of integration point - vec2d m_rsp; //!< m_rs at the previous time step - bool m_bstick; //!< stick flag - vector m_Lmc; //!< Lagrange multipliers for solute concentrations - vector m_epsc; //!< concentration penalty factors - vector m_cg; //!< concentration "gap" - vector m_c1; //!< solute concentration - }; + FEBiphasicContactPoint::Init(); + } + + void Serialize(DumpStream& ar) override; +}; +//----------------------------------------------------------------------------- +class FEBIOMIX_API FESlidingSurfaceMP : public FEBiphasicContactSurface +{ public: //! constructor FESlidingSurfaceMP(FEModel* pfem); diff --git a/FEBioMix/FETiedBiphasicInterface.cpp b/FEBioMix/FETiedBiphasicInterface.cpp index f434d4734..4a321c9cf 100644 --- a/FEBioMix/FETiedBiphasicInterface.cpp +++ b/FEBioMix/FETiedBiphasicInterface.cpp @@ -56,22 +56,7 @@ BEGIN_FECORE_CLASS(FETiedBiphasicInterface, FEContactInterface) END_FECORE_CLASS(); //----------------------------------------------------------------------------- -FETiedBiphasicSurface::Data::Data() -{ - m_Gap = vec3d(0,0,0); - m_dg = vec3d(0,0,0); - m_nu = vec3d(0,0,0); - m_rs = vec2d(0,0); - m_Lmd = vec3d(0,0,0); - m_tr = vec3d(0,0,0); - m_Lmp = 0.0; - m_epsn = 1.0; - m_epsp = 1.0; - m_pg = 0.0; -} - -//----------------------------------------------------------------------------- -void FETiedBiphasicSurface::Data::Serialize(DumpStream& ar) +void FETiedBiphasicContactPoint::Serialize(DumpStream& ar) { FEBiphasicContactPoint::Serialize(ar); ar & m_Gap; @@ -97,7 +82,7 @@ FETiedBiphasicSurface::FETiedBiphasicSurface(FEModel* pfem) : FEBiphasicContactS //! create material point data FEMaterialPoint* FETiedBiphasicSurface::CreateMaterialPoint() { - return new FETiedBiphasicSurface::Data; + return new FETiedBiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -189,7 +174,7 @@ void FETiedBiphasicSurface::GetVectorGap(int nface, vec3d& pg) pg = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FETiedBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pg += data.m_dg; } pg /= ni; @@ -203,7 +188,7 @@ void FETiedBiphasicSurface::GetContactTraction(int nface, vec3d& pt) pt = vec3d(0,0,0); for (int k = 0; k < ni; ++k) { - Data& data = static_cast(*el.GetMaterialPoint(k)); + FETiedBiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(k)); pt += data.m_tr; } pt /= ni; @@ -293,7 +278,7 @@ void FETiedBiphasicInterface::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FETiedBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = pt.m_pme; if (pe != 0) { @@ -381,7 +366,7 @@ void FETiedBiphasicInterface::CalcAutoPenalty(FETiedBiphasicSurface& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -403,7 +388,7 @@ void FETiedBiphasicInterface::CalcAutoPressurePenalty(FETiedBiphasicSurface& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -498,7 +483,7 @@ void FETiedBiphasicInterface::InitialProjection(FETiedBiphasicSurface& ss, FETie // find the intersection point with the secondary surface pme = np.Project2(r, nu, rs); - FETiedBiphasicSurface::Data& pt = static_cast(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_pme = pme; pt.m_rs[0] = rs[0]; pt.m_rs[1] = rs[1]; @@ -547,7 +532,7 @@ void FETiedBiphasicInterface::ProjectSurface(FETiedBiphasicSurface& ss, FETiedBi for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point r = ss.Local2Global(el, j); @@ -568,6 +553,7 @@ void FETiedBiphasicInterface::ProjectSurface(FETiedBiphasicSurface& ss, FETiedBi // calculate the gap function vec3d g = q - r; pt.m_dg = g - pt.m_Gap; +// pt.m_gap = pt.m_dg.unit(); // calculate the pressure gap function bool mporo = ms.m_poro[pme->m_lid]; @@ -582,6 +568,7 @@ void FETiedBiphasicInterface::ProjectSurface(FETiedBiphasicSurface& ss, FETiedBi { // the node is not tied pt.m_dg = vec3d(0,0,0); +// pt.m_gap = 0; if (sporo) pt.m_pg = 0; } } @@ -659,7 +646,7 @@ void FETiedBiphasicInterface::LoadVector(FEGlobalVector& R, const FETimeInfo& tp // note that we are integrating over the current surface for (j=0; j(*se.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -868,7 +855,7 @@ void FETiedBiphasicInterface::StiffnessMatrix(FELinearSystem& LS, const FETimeIn // loop over all integration points for (j=0; j(*se.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -1226,7 +1213,7 @@ bool FETiedBiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); normL0 += ds.m_Lmd*ds.m_Lmd; } } @@ -1235,7 +1222,7 @@ bool FETiedBiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); normL0 += dm.m_Lmd*dm.m_Lmd; } } @@ -1252,7 +1239,7 @@ bool FETiedBiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface eps = m_epsn*ds.m_epsn; @@ -1280,7 +1267,7 @@ bool FETiedBiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedBiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface eps = m_epsn*dm.m_epsn; diff --git a/FEBioMix/FETiedBiphasicInterface.h b/FEBioMix/FETiedBiphasicInterface.h index 13d76c17e..ffbfdf806 100644 --- a/FEBioMix/FETiedBiphasicInterface.h +++ b/FEBioMix/FETiedBiphasicInterface.h @@ -27,32 +27,38 @@ SOFTWARE.*/ #pragma once -#include "FEBioMech/FEContactInterface.h" +#include #include "FEBiphasicContactSurface.h" //----------------------------------------------------------------------------- -class FEBIOMIX_API FETiedBiphasicSurface : public FEBiphasicContactSurface +class FEBIOMIX_API FETiedBiphasicContactPoint: public FEBiphasicContactPoint { public: - //! Integration point data - class Data : public FEBiphasicContactPoint + vec3d m_Gap; //!< initial gap in reference configuration + vec3d m_dg; //!< gap function at integration points + vec3d m_nu; //!< normal at integration points + vec2d m_rs; //!< natural coordinates of projection of integration point + vec3d m_Lmd; //!< lagrange multipliers for displacements + vec3d m_tr; //!< contact traction + double m_epsn; //!< penalty factors + double m_epsp; //!< pressure penalty factors + + void Init() override { - public: - Data(); - - void Serialize(DumpStream& ar) override; - - public: - vec3d m_Gap; //!< initial gap in reference configuration - vec3d m_dg; //!< gap function at integration points - vec3d m_nu; //!< normal at integration points - vec2d m_rs; //!< natural coordinates of projection of integration point - vec3d m_Lmd; //!< lagrange multipliers for displacements - vec3d m_tr; //!< contact traction - double m_epsn; //!< penalty factors - double m_epsp; //!< pressure penalty factors - }; + FEBiphasicContactPoint::Init(); + m_Gap = m_dg = m_nu = m_Lmd = m_tr = vec3d(0,0,0); + m_rs = vec2d(0,0); + m_epsn = 1.0; + m_epsp = 1.0; + } + void Serialize(DumpStream& ar) override; +}; + +//----------------------------------------------------------------------------- +class FEBIOMIX_API FETiedBiphasicSurface : public FEBiphasicContactSurface +{ +public: //! constructor FETiedBiphasicSurface(FEModel* pfem); diff --git a/FEBioMix/FETiedMultiphasicInterface.cpp b/FEBioMix/FETiedMultiphasicInterface.cpp index b3a67f8a6..670386579 100644 --- a/FEBioMix/FETiedMultiphasicInterface.cpp +++ b/FEBioMix/FETiedMultiphasicInterface.cpp @@ -61,29 +61,9 @@ BEGIN_FECORE_CLASS(FETiedMultiphasicInterface, FEContactInterface) END_FECORE_CLASS(); //----------------------------------------------------------------------------- -FETiedMultiphasicSurface::Data::Data() +void FETiedMultiphasicContactPoint::Serialize(DumpStream& ar) { - m_Gap = vec3d(0,0,0); - m_dg = vec3d(0,0,0); - m_nu = vec3d(0,0,0); - m_rs = vec2d(0,0); - m_Lmd = vec3d(0,0,0); - m_Lmp = 0.0; - m_epsn = 1.0; - m_epsp = 1.0; - m_pg = 0.0; -} - -void FETiedMultiphasicSurface::Data::Serialize(DumpStream& ar) -{ - FEBiphasicContactPoint::Serialize(ar); - ar & m_Gap; - ar & m_dg; - ar & m_nu; - ar & m_rs; - ar & m_Lmd; - ar & m_epsn; - ar & m_epsp; + FETiedBiphasicContactPoint::Serialize(ar); ar & m_Lmc; ar & m_epsc; ar & m_cg; @@ -103,7 +83,7 @@ FETiedMultiphasicSurface::FETiedMultiphasicSurface(FEModel* pfem) : FEBiphasicCo //! create material point data FEMaterialPoint* FETiedMultiphasicSurface::CreateMaterialPoint() { - return new FETiedMultiphasicSurface::Data; + return new FETiedMultiphasicContactPoint; } //----------------------------------------------------------------------------- @@ -229,7 +209,7 @@ bool FETiedMultiphasicSurface::Init() int nint = el.GaussPoints(); if (nsol) { for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& data = static_cast(*el.GetMaterialPoint(j)); data.m_Lmc.resize(nsol); data.m_epsc.resize(nsol); data.m_cg.resize(nsol); @@ -393,7 +373,7 @@ void FETiedMultiphasicInterface::BuildMatrixProfile(FEGlobalMatrix& K) int* sn = &se.m_node[0]; for (k=0; k(*se.GetMaterialPoint(k)); + FETiedMultiphasicContactPoint& data = static_cast(*se.GetMaterialPoint(k)); FESurfaceElement* pe = data.m_pme; if (pe != 0) { @@ -484,7 +464,7 @@ void FETiedMultiphasicInterface::CalcAutoPenalty(FETiedMultiphasicSurface& s) int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsn = eps; } } @@ -570,7 +550,7 @@ void FETiedMultiphasicInterface::CalcAutoPressurePenalty(FETiedMultiphasicSurfac int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsp = eps; } } @@ -647,7 +627,7 @@ void FETiedMultiphasicInterface::CalcAutoConcentrationPenalty(FETiedMultiphasicS int nint = el.GaussPoints(); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_epsc[isol] = eps; } } @@ -742,7 +722,7 @@ void FETiedMultiphasicInterface::InitialProjection(FETiedMultiphasicSurface& ss, // find the intersection point with the secondary surface pme = np.Project2(r, nu, rs); - FETiedMultiphasicSurface::Data& pt = static_cast(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); pt.m_pme = pme; pt.m_rs[0] = rs[0]; pt.m_rs[1] = rs[1]; @@ -801,7 +781,7 @@ void FETiedMultiphasicInterface::ProjectSurface(FETiedMultiphasicSurface& ss, FE for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*el.GetMaterialPoint(j)); // calculate the global position of the integration point r = ss.Local2Global(el, j); @@ -932,7 +912,7 @@ void FETiedMultiphasicInterface::LoadVector(FEGlobalVector& R, const FETimeInfo& // integration weights w[j] = se.GaussWeights()[j]; - FETiedMultiphasicSurface::Data& pt = static_cast(*se.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // contact traction double eps = m_epsn*pt.m_epsn; // penalty @@ -955,7 +935,7 @@ void FETiedMultiphasicInterface::LoadVector(FEGlobalVector& R, const FETimeInfo& // note that we are integrating over the current surface for (j=0; j(*se.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; if (pme) @@ -1159,7 +1139,7 @@ void FETiedMultiphasicInterface::StiffnessMatrix(FELinearSystem& LS, const FETim // integration weights w[j] = se.GaussWeights()[j]; - FETiedMultiphasicSurface::Data& pd = static_cast(*se.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pd = static_cast(*se.GetMaterialPoint(j)); // contact traction double eps = m_epsn*pd.m_epsn; // penalty @@ -1185,7 +1165,7 @@ void FETiedMultiphasicInterface::StiffnessMatrix(FELinearSystem& LS, const FETim // loop over all integration points for (j=0; j(*se.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& pt = static_cast(*se.GetMaterialPoint(j)); // get the secondary surface element FESurfaceElement* pme = pt.m_pme; @@ -1485,7 +1465,7 @@ bool FETiedMultiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); normL0 += ds.m_Lmd*ds.m_Lmd; } } @@ -1494,7 +1474,7 @@ bool FETiedMultiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j=0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); normL0 += dm.m_Lmd*dm.m_Lmd; } } @@ -1512,7 +1492,7 @@ bool FETiedMultiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ss.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& ds = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on primary surface eps = m_epsn*ds.m_epsn; @@ -1550,7 +1530,7 @@ bool FETiedMultiphasicInterface::Augment(int naug, const FETimeInfo& tp) FESurfaceElement& el = m_ms.Element(i); for (int j = 0; j(*el.GetMaterialPoint(j)); + FETiedMultiphasicContactPoint& dm = static_cast(*el.GetMaterialPoint(j)); // update Lagrange multipliers on secondary surface eps = m_epsn*dm.m_epsn; diff --git a/FEBioMix/FETiedMultiphasicInterface.h b/FEBioMix/FETiedMultiphasicInterface.h index 7397f3b53..4796b2875 100644 --- a/FEBioMix/FETiedMultiphasicInterface.h +++ b/FEBioMix/FETiedMultiphasicInterface.h @@ -27,35 +27,33 @@ SOFTWARE.*/ #pragma once -#include "FEBioMech/FEContactInterface.h" -#include "FEBiphasicContactSurface.h" +#include "FETiedBiphasicInterface.h" #include "FESolute.h" //----------------------------------------------------------------------------- -class FEBIOMIX_API FETiedMultiphasicSurface : public FEBiphasicContactSurface +class FEBIOMIX_API FETiedMultiphasicContactPoint: public FETiedBiphasicContactPoint { public: - //! Integration point data - class Data : public FEBiphasicContactPoint + vector m_Lmc; //!< Lagrange multipliers for solute concentrations + vector m_epsc; //!< concentration penatly factors + vector m_cg; //!< concentration "gap" + + void Init() override { - public: - Data(); - - void Serialize(DumpStream& ar) override; - - public: - vec3d m_Gap; //!< initial gap in reference configuration - vec3d m_dg; //!< gap function at integration points - vec3d m_nu; //!< normal at integration points - vec2d m_rs; //!< natural coordinates of projection of integration point - vec3d m_Lmd; //!< lagrange multipliers for displacements - double m_epsn; //!< penalty factors - double m_epsp; //!< pressure penalty factors - vector m_Lmc; //!< Lagrange multipliers for solute concentrations - vector m_epsc; //!< concentration penatly factors - vector m_cg; //!< concentration "gap" - }; + FEBiphasicContactPoint::Init(); + m_Gap = m_dg = m_nu = m_Lmd = m_tr = vec3d(0,0,0); + m_rs = vec2d(0,0); + m_epsn = 1.0; + m_epsp = 1.0; + } + void Serialize(DumpStream& ar) override; +}; + +//----------------------------------------------------------------------------- +class FEBIOMIX_API FETiedMultiphasicSurface : public FEBiphasicContactSurface +{ +public: //! constructor FETiedMultiphasicSurface(FEModel* pfem);