Modified fluid surface force log output to use vector output

Author: gateshian

FEBioFluid/FEBioFluid.cpp

@@ -294,9 +294,6 @@ REGISTER_FECORE_CLASS(FELogFluidRateOfDefZZ    , "fdzz");
 REGISTER_FECORE_CLASS(FELogFluidRateOfDefXY    , "fdxy");
 REGISTER_FECORE_CLASS(FELogFluidRateOfDefYZ    , "fdyz");
 REGISTER_FECORE_CLASS(FELogFluidRateOfDefXZ    , "fdxz");
-REGISTER_FECORE_CLASS(FELogFluidSurfaceForceX  , "fsfx");
-REGISTER_FECORE_CLASS(FELogFluidSurfaceForceY  , "fsfy");
-REGISTER_FECORE_CLASS(FELogFluidSurfaceForceZ  , "fsfz");
 REGISTER_FECORE_CLASS(FELogFluidSurfaceForce   , "fsf" );
 
 //-----------------------------------------------------------------------------

FEBioFluid/FEBioFluidData.cpp

@@ -525,12 +525,12 @@ double FELogFluidRateOfDefXZ::value(FEElement& el)
 //=============================================================================
 
 //-----------------------------------------------------------------------------
-//! Fluid surface force along X
+//! Fluid surface force
 //-----------------------------------------------------------------------------
-double FELogFluidSurfaceForceX::value(FESurface& surface)
+vec3d FELogFluidSurfaceForce::typedValue(FESurface& surface)
 {
     FESurface* pcs = &surface;
-    if (pcs == 0) return false;
+    if (pcs == 0) return vec3d(0,0,0);
 
     int NF = pcs->Elements();
     vec3d fn(0,0,0);    // initialize
@@ -582,188 +582,5 @@ double FELogFluidSurfaceForceX::value(FESurface& surface)
         }
     }
 
-    return fn.x;
-}
-
-//-----------------------------------------------------------------------------
-//! Fluid surface force along Y
-//-----------------------------------------------------------------------------
-double FELogFluidSurfaceForceY::value(FESurface& surface)
-{
-    FESurface* pcs = &surface;
-    if (pcs == 0) return false;
-    
-    int NF = pcs->Elements();
-    vec3d fn(0,0,0);    // initialize
-    
-    // calculate the vectorial area of each surface element and to identify solid element associated with this surface element
-    m_area.resize(NF);
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        m_area[j] = pcs->SurfaceNormal(el,0,0)*pcs->FaceArea(el);
-    }
-    
-    // calculate net fluid force
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        
-        // get the element this surface element belongs to
-        FEElement* pe = el.m_elem[0].pe;
-        if (pe)
-        {
-            // get the material
-            FEMaterial* pm = GetFEModel()->GetMaterial(pe->GetMatID());
-            FEFluidMaterial* pfluid = pm->ExtractProperty<FEFluidMaterial>();
-            
-            if (!pfluid) {
-                pe = el.m_elem[1].pe;
-                if (pe) pfluid = GetFEModel()->GetMaterial(pe->GetMatID())->ExtractProperty<FEFluidMaterial>();
-            }
-            
-            // see if this is a fluid element
-            if (pfluid) {
-                // evaluate the average stress in this element
-                int nint = pe->GaussPoints();
-                mat3d s(mat3dd(0));
-                for (int n=0; n<nint; ++n)
-                {
-                    FEMaterialPoint& mp = *pe->GetMaterialPoint(n);
-                    FEFluidMaterialPoint& pt = *(mp.ExtractData<FEFluidMaterialPoint>());
-                    s += pt.m_sf;
-                }
-                s /= nint;
-                
-                // Evaluate contribution to net force on surface.
-                // Negate the fluid traction since we want the traction on the surface,
-                // which is the opposite of the traction on the fluid.
-                fn -= s*m_area[j];
-            }
-        }
-    }
-    
-    return fn.y;
-}
-
-//-----------------------------------------------------------------------------
-//! Fluid surface force along Z
-//-----------------------------------------------------------------------------
-double FELogFluidSurfaceForceZ::value(FESurface& surface)
-{
-    FESurface* pcs = &surface;
-    if (pcs == 0) return false;
-    
-    int NF = pcs->Elements();
-    vec3d fn(0,0,0);    // initialize
-    
-    // calculate the vectorial area of each surface element and to identify solid element associated with this surface element
-    m_area.resize(NF);
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        m_area[j] = pcs->SurfaceNormal(el,0,0)*pcs->FaceArea(el);
-    }
-    
-    // calculate net fluid force
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        
-        // get the element this surface element belongs to
-        FEElement* pe = el.m_elem[0].pe;
-        if (pe)
-        {
-            // get the material
-            FEMaterial* pm = GetFEModel()->GetMaterial(pe->GetMatID());
-            FEFluidMaterial* pfluid = pm->ExtractProperty<FEFluidMaterial>();
-            
-            if (!pfluid) {
-                pe = el.m_elem[1].pe;
-                if (pe) pfluid = GetFEModel()->GetMaterial(pe->GetMatID())->ExtractProperty<FEFluidMaterial>();
-            }
-            
-            // see if this is a fluid element
-            if (pfluid) {
-                // evaluate the average stress in this element
-                int nint = pe->GaussPoints();
-                mat3d s(mat3dd(0));
-                for (int n=0; n<nint; ++n)
-                {
-                    FEMaterialPoint& mp = *pe->GetMaterialPoint(n);
-                    FEFluidMaterialPoint& pt = *(mp.ExtractData<FEFluidMaterialPoint>());
-                    s += pt.m_sf;
-                }
-                s /= nint;
-                
-                // Evaluate contribution to net force on surface.
-                // Negate the fluid traction since we want the traction on the surface,
-                // which is the opposite of the traction on the fluid.
-                fn -= s*m_area[j];
-            }
-        }
-    }
-    
-    return fn.z;
-}
-
-//-----------------------------------------------------------------------------
-//! Fluid surface force magnitude
-//-----------------------------------------------------------------------------
-double FELogFluidSurfaceForce::value(FESurface& surface)
-{
-    FESurface* pcs = &surface;
-    if (pcs == 0) return false;
-    
-    int NF = pcs->Elements();
-    vec3d fn(0,0,0);    // initialize
-    
-    // calculate the vectorial area of each surface element and to identify solid element associated with this surface element
-    m_area.resize(NF);
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        m_area[j] = pcs->SurfaceNormal(el,0,0)*pcs->FaceArea(el);
-    }
-    
-    // calculate net fluid force
-    for (int j=0; j<NF; ++j)
-    {
-        FESurfaceElement& el = pcs->Element(j);
-        
-        // get the element this surface element belongs to
-        FEElement* pe = el.m_elem[0].pe;
-        if (pe)
-        {
-            // get the material
-            FEMaterial* pm = GetFEModel()->GetMaterial(pe->GetMatID());
-            FEFluidMaterial* pfluid = pm->ExtractProperty<FEFluidMaterial>();
-            
-            if (!pfluid) {
-                pe = el.m_elem[1].pe;
-                if (pe) pfluid = GetFEModel()->GetMaterial(pe->GetMatID())->ExtractProperty<FEFluidMaterial>();
-            }
-            
-            // see if this is a fluid element
-            if (pfluid) {
-                // evaluate the average stress in this element
-                int nint = pe->GaussPoints();
-                mat3d s(mat3dd(0));
-                for (int n=0; n<nint; ++n)
-                {
-                    FEMaterialPoint& mp = *pe->GetMaterialPoint(n);
-                    FEFluidMaterialPoint& pt = *(mp.ExtractData<FEFluidMaterialPoint>());
-                    s += pt.m_sf;
-                }
-                s /= nint;
-                
-                // Evaluate contribution to net force on surface.
-                // Negate the fluid traction since we want the traction on the surface,
-                // which is the opposite of the traction on the fluid.
-                fn -= s*m_area[j];
-            }
-        }
-    }
-    
-    return fn.norm();
+    return fn;
 }

FEBioFluid/FEBioFluidData.h

@@ -324,50 +324,12 @@ class FELogFluidRateOfDefXZ : public FELogElemData
 //=============================================================================
 
 //-----------------------------------------------------------------------------
-//! Fluid surface force along X
-//!
-class FELogFluidSurfaceForceX : public FELogSurfaceData
+class FELogFluidSurfaceForce : public FELogSurfaceVec3dData
 {
 private:
     vector<vec3d>       m_area;
 
 public:
-    FELogFluidSurfaceForceX(FEModel* pfem) : FELogSurfaceData(pfem) {}
-    double value(FESurface& surface) override;
-};
-
-//! Fluid surface force along Y
-//!
-class FELogFluidSurfaceForceY : public FELogSurfaceData
-{
-private:
-    vector<vec3d>       m_area;
-    
-public:
-    FELogFluidSurfaceForceY(FEModel* pfem) : FELogSurfaceData(pfem) {}
-    double value(FESurface& surface) override;
-};
-
-//! Fluid surface force along Z
-//!
-class FELogFluidSurfaceForceZ : public FELogSurfaceData
-{
-private:
-    vector<vec3d>       m_area;
-    
-public:
-    FELogFluidSurfaceForceZ(FEModel* pfem) : FELogSurfaceData(pfem) {}
-    double value(FESurface& surface) override;
-};
-
-//! Fluid surface force magnitude
-//!
-class FELogFluidSurfaceForce : public FELogSurfaceData
-{
-private:
-    vector<vec3d>       m_area;
-    
-public:
-    FELogFluidSurfaceForce(FEModel* pfem) : FELogSurfaceData(pfem) {}
-    double value(FESurface& surface) override;
+    FELogFluidSurfaceForce(FEModel* fem) : FELogSurfaceVec3dData(fem) {}
+    vec3d typedValue(FESurface& surf) override;
 };