Merge branch 'develop' of https://github.com/febiosoftware/FEBio into…

… develop

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);
 };
 

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"          );

FEBioMix/FEBioMixPlot.cpp

@@ -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>
@@ -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);
 	});
 
@@ -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
 //=============================================================================

FEBioMix/FEBioMixPlot.h

@@ -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);
 };
 
@@ -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;
+};
+

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;
 }
 

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;