Fixed problem with exporting fiber components to log file.

febiosoftware · Nov 27, 2024 · 9c2483a · 9c2483a
1 parent 84d85f8
commit 9c2483a
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 67 deletions.
diff --git a/FEBioMech/FEBioMechData.cpp b/FEBioMech/FEBioMechData.cpp
@@ -1672,69 +1672,39 @@ double FELogElemFiberStretch::value(FEElement& el)
 }
 
 //-----------------------------------------------------------------------------
-double FELogElemFiberVectorX::value(FEElement& el)
+double FELogElemFiberVector_::value(FEElement& el)
 {
 	int matID = el.GetMatID();
 	FEMaterial* mat = GetFEModel()->GetMaterial(matID);
 
-	int n = el.GaussPoints();
-	double l = 0.0;
-	for (int j = 0; j<n; ++j)
-	{
-		FEMaterialPoint& mp = *el.GetMaterialPoint(j);
-		FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
-		mat3d Q = mat->GetLocalCS(mp);
+	FEElasticMaterial* pme = mat->ExtractProperty<FEElasticMaterial>();
+	if (pme == nullptr) return 0.0;
 
-		vec3d ri = Q.col(0);
-		vec3d r = pt.m_F*ri;
-
-		l += r.x;
-	}
-	l /= (double)n;
-	return l;
-}
-
-//-----------------------------------------------------------------------------
-double FELogElemFiberVectorY::value(FEElement& el)
-{
-	int matID = el.GetMatID();
-	FEMaterial* mat = GetFEModel()->GetMaterial(matID);
+	FEVec3dValuator* vec = dynamic_cast<FEVec3dValuator*>(pme->GetProperty("fiber"));
+	if (vec == nullptr) return 0.0;
 
 	int n = el.GaussPoints();
 	double l = 0.0;
 	for (int j = 0; j<n; ++j)
 	{
 		FEMaterialPoint& mp = *el.GetMaterialPoint(j);
-		FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
-		mat3d Q = mat->GetLocalCS(mp);
-
-		vec3d ri = Q.col(0);
-		vec3d r = pt.m_F*ri;
-
-		l += r.y;
-	}
-	l /= (double)n;
-	return l;
-}
-
-//-----------------------------------------------------------------------------
-double FELogElemFiberVectorZ::value(FEElement& el)
-{
-	int matID = el.GetMatID();
-	FEMaterial* mat = GetFEModel()->GetMaterial(matID);
-
-	int n = el.GaussPoints();
-	double l = 0.0;
-	for (int j = 0; j<n; ++j)
-	{
-		FEMaterialPoint& mp = *el.GetMaterialPoint(j);
-		FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
-		mat3d Q = mat->GetLocalCS(mp);
-
-		vec3d ri = Q.col(0);
-		vec3d r = pt.m_F*ri;
+
+		const FEElasticMaterialPoint* pt = mp.ExtractData<const FEElasticMaterialPoint>();
+		if (pt)
+		{
+			mat3d Q = pme->GetLocalCS(mp);
+			mat3d F = pt->m_F;
+			vec3d a0 = vec->unitVector(mp);
+			vec3d ar = Q * a0;
+			vec3d a = F * ar; a.unit();
 
-		l += r.z;
+			switch (m_comp)
+			{
+			case 0: l += a.x; break;
+			case 1: l += a.y; break;
+			case 2: l += a.z; break;
+			}
+		}
 	}
 	l /= (double)n;
 	return l;

diff --git a/FEBioMech/FEBioMechData.h b/FEBioMech/FEBioMechData.h
@@ -971,27 +971,20 @@ class FELogElemFiberStretch : public FELogElemData
 };
 
 //-----------------------------------------------------------------------------
-class FELogElemFiberVectorX : public FELogElemData
+class FELogElemFiberVector_ : public FELogElemData
 {
 public:
-	FELogElemFiberVectorX(FEModel* pfem) : FELogElemData(pfem){}
+	FELogElemFiberVector_(FEModel* pfem, int comp) : FELogElemData(pfem), m_comp(comp) {}
 	double value(FEElement& el);
-};
 
-//-----------------------------------------------------------------------------
-class FELogElemFiberVectorY : public FELogElemData
-{
-public:
-	FELogElemFiberVectorY(FEModel* pfem) : FELogElemData(pfem){}
-	double value(FEElement& el);
+private:
+	int m_comp;
 };
 
-//-----------------------------------------------------------------------------
-class FELogElemFiberVectorZ : public FELogElemData
+template <int N> class FELogElemFiberVector_N : public FELogElemFiberVector_
 {
 public:
-	FELogElemFiberVectorZ(FEModel* pfem) : FELogElemData(pfem){}
-	double value(FEElement& el);
+	FELogElemFiberVector_N(FEModel* fem) : FELogElemFiberVector_(fem, N) {}
 };
 
 //-----------------------------------------------------------------------------

diff --git a/FEBioMech/FEBioMechModule.cpp b/FEBioMech/FEBioMechModule.cpp
@@ -1033,9 +1033,9 @@ void FEBioMech::InitModule()
 	REGISTER_FECORE_CLASS(FELogElemStrainEnergyDensity, "sed");
 	REGISTER_FECORE_CLASS(FELogElemDevStrainEnergyDensity, "devsed");
 	REGISTER_FECORE_CLASS(FELogElemFiberStretch, "fiber_stretch");
-	REGISTER_FECORE_CLASS(FELogElemFiberVectorX, "fiber_x");
-	REGISTER_FECORE_CLASS(FELogElemFiberVectorY, "fiber_y");
-	REGISTER_FECORE_CLASS(FELogElemFiberVectorZ, "fiber_z");
+	REGISTER_FECORE_CLASS_T(FELogElemFiberVector_N, 0, "fiber_x");
+	REGISTER_FECORE_CLASS_T(FELogElemFiberVector_N, 1, "fiber_y");
+	REGISTER_FECORE_CLASS_T(FELogElemFiberVector_N, 2, "fiber_z");
 	REGISTER_FECORE_CLASS(FELogDamage, "D");
 	REGISTER_FECORE_CLASS_T(FELogDamage_n, 0, "damage_1");
 	REGISTER_FECORE_CLASS_T(FELogDamage_n, 1, "damage_2");