Fixed calculation of eigenvectors of moment of inertia. Provided ability to calculate area-based MOI. Added tool to align eigenvectors of MOI with coordinate system axes.

Author: gateshian

FEBioStudio/MeasureCOMTool.cpp

@@ -44,47 +44,3 @@ bool CMeasureCOMTool::OnApply()
 	m_com = CalculateCOM(*mesh);
 	return true;
 }
-
-// This algorithm calculates the center of mass approximately,
-// by weighing the element centers by the element volumes. 
-// TODO: I should look into bringing the numerical integration tools
-//       from FEBio to FBS so this can be calculated accurately. 
-vec3d CalculateCOM(FSMesh& mesh)
-{
-	vec3d rc(0, 0, 0);
-	double Vtotal = 0.0;
-
-	int selectedElements = mesh.CountSelectedElements();
-
-	// loop over all elements
-	int NE = mesh.Elements();
-	for (int iel = 0; iel < NE; ++iel)
-	{
-		FSElement& el = mesh.Element(iel);
-		if ((selectedElements == 0) || el.IsSelected())
-		{
-			// Volume is calculated in global frame
-			double V0 = mesh.ElementVolume(el);
-
-			// center of mass is calculated in local frame
-			vec3d c_local = mesh.ElementCenter(el);
-			vec3d c_global = mesh.LocalToGlobal(c_local);
-
-			rc += c_global * V0;
-			Vtotal += V0;
-		}
-	}
-	if (Vtotal != 0.0)
-	rc /= Vtotal;
-
-	const double eps = 1e-12;
-	double L = rc.Length();
-	if (L > eps)
-	{
-		if (fabs(rc.x) < eps) rc.x = 0;
-		if (fabs(rc.y) < eps) rc.y = 0;
-		if (fabs(rc.z) < eps) rc.z = 0;
-	}
-
-	return rc;
-}

FEBioStudio/MeasureCOMTool.h

@@ -25,6 +25,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.*/
 
 #include "Tool.h"
+#include "MeasureTools.h"
 
 class CMainWindow;
 

FEBioStudio/MeasureMOITool.cpp

@@ -32,12 +32,14 @@ mat3d CalculateMOI(FSMesh& mesh);
 // constructor
 CMeasureMOITool::CMeasureMOITool(CMainWindow* wnd) : CBasicTool(wnd, "Moment of inertia", CBasicTool::HAS_APPLY_BUTTON)
 {
+    m_area = false;
 	m_moi.zero();
     m_evec.zero();
     m_eval = vec3d(0,0,0);
 	addMat3Property(&m_moi, "moment of inertia:")->setFlags(CProperty::Visible);
     addMat3Property(&m_evec, "eigenvectors of MOI:")->setFlags(CProperty::Visible);
     addVec3Property(&m_eval, "eigenvalues of MOI:")->setFlags(CProperty::Visible);
+    addBoolProperty(&m_area, "area");
 
 	SetInfo("Calculates the moment of inertia of a meshed object or element selection. The MOI is calculated with respect to the center of mass of the selection.");
 }
@@ -47,54 +49,16 @@ bool CMeasureMOITool::OnApply()
 {
 	FSMesh* mesh = GetActiveMesh();
 	if (mesh == nullptr) return false;
-	m_moi = CalculateMOI(*mesh);
+    m_moi = (m_area) ? CalculateAreaMOI(*mesh) : CalculateMOI(*mesh);
     double eval[3];
     vec3d evec[3];
     mat3ds mois = m_moi.sym();
     mois.eigen2(eval,evec);
-    m_evec = mat3d(evec[0].x, evec[1].x, evec[2].x,
-                   evec[0].y, evec[1].y, evec[2].y,
-                   evec[0].z, evec[1].z, evec[2].z);
+    for (int i=0; i<3; ++i) evec[i].Normalize();
+    m_evec = mat3d(evec[0].x, evec[0].y, evec[0].z,
+                   evec[1].x, evec[1].y, evec[1].z,
+                   evec[2].x, evec[2].y, evec[2].z);
     m_eval = vec3d(eval[0],eval[1],eval[2]);
+    m_evec = CleanUp(m_evec, 1e-12);
 	return true;
 }
-
-// This calculates the moment of inertia, but only approximately!
-// Requires a fine mesh for accurate results. 
-mat3d CalculateMOI(FSMesh& mesh)
-{
-	mat3dd I(1);        // identity tensor
-
-	mat3d moi; moi.zero();
-
-	// calculate the COM
-	vec3d com = CalculateCOM(mesh);
-
-	// loop over all elements
-	int NE = mesh.Elements();
-	for (int iel = 0; iel < NE; ++iel)
-	{
-		FSElement& el = mesh.Element(iel);
-
-		// ElementCenter is calculated in local frame
-		vec3d c_local = mesh.ElementCenter(el);
-		vec3d c_global = mesh.LocalToGlobal(c_local);
-
-		// Element volume is calculated in global frame
-		double V0 = mesh.ElementVolume(el);
-
-		vec3d r = c_global - com;
-		mat3d Iij = (r * r) * I - (r & r);
-		moi += Iij * V0;
-	}
-
-	const double eps = 1e-12;
-	double L = moi.norm();
-	if (L > eps)
-	{
-		for (int i = 0; i < 3; ++i)
-			for (int j = 0; j < 3; ++j)
-				if (fabs(moi[i][j]) < eps) moi[i][j] = 0.0;
-	}
-	return moi;
-}

FEBioStudio/MeasureMOITool.h

@@ -25,6 +25,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.*/
 
 #include "Tool.h"
+#include "MeasureTools.h"
 
 class CMainWindow;
 
@@ -42,4 +43,5 @@ class CMeasureMOITool : public CBasicTool
 	mat3d	m_moi;	// center of mass
     mat3d   m_evec; // eigenvectors of MOI
     vec3d   m_eval; // eigenvalues of MOI
+    bool    m_area; // boolean flag to report area MOI
 };

FEBioStudio/MeshPanel.cpp

@@ -274,6 +274,7 @@ void ModifierThread::stop()
 REGISTER_CLASS(FEAddNode              , CLASS_FEMODIFIER, "Add Node"       , EDIT_MESH);
 REGISTER_CLASS(FEAddTriangle          , CLASS_FEMODIFIER, "Add Triangle"   , EDIT_MESH);
 REGISTER_CLASS(FEAlignNodes           , CLASS_FEMODIFIER, "Align"          , EDIT_NODE);
+REGISTER_CLASS(FEAlignMeshMOI         , CLASS_FEMODIFIER, "Align MOI to CS", EDIT_MESH);
 REGISTER_CLASS(FEAutoPartition        , CLASS_FEMODIFIER, "Auto Partition" , EDIT_MESH);
 REGISTER_CLASS(FEBoundaryLayerMesher  , CLASS_FEMODIFIER, "Boundary Layer" , EDIT_FACE | EDIT_SAFE);
 REGISTER_CLASS(FEConvertMesh		  , CLASS_FEMODIFIER, "Convert Mesh"   , EDIT_MESH | EDIT_SAFE);

MeshTools/FEModifier.cpp

@@ -39,6 +39,8 @@ SOFTWARE.*/
 #include <MeshLib/Intersect.h>
 #include <GeomLib/GGroup.h>
 #include <MeshTools/LaplaceSolver.h>
+#include <FEBioStudio/MeasureTools.h>
+#include <numeric>
 
 FEModifier::FEModifier(const char* sz) { SetName(sz); }
 FEModifier::~FEModifier() {}
@@ -2048,3 +2050,56 @@ FSMesh* FEDetachElements::Apply(FSMesh* pm)
 
 	return newMesh;		
 }
+
+//=============================================================================
+// FEAlignMesh
+//-----------------------------------------------------------------------------
+
+FEAlignMeshMOI::FEAlignMeshMOI() : FEModifier("Align Mesh MOI")
+{
+    AddBoolParam(true, "Use area MOI");
+}
+
+FSMesh* FEAlignMeshMOI::Apply(FSMesh* pm)
+{
+    GObject* po = GObject::GetActiveObject();
+    if (po == nullptr) return nullptr;
+
+    bool use_area = GetBoolValue(0);
+    mat3d moi = (use_area) ? CalculateAreaMOI(*pm) : CalculateMOI(*pm);
+    double eval[3];
+    vec3d evec[3];
+    mat3ds mois = moi.sym();
+    mois.eigen2(eval,evec);
+    // sort eigenvalues in ascending order
+    std::vector<int> indices(3);
+    std::iota(indices.begin(), indices.end(),0);
+    std::vector<double> ev(eval,eval+3);
+    std::sort(indices.begin(), indices.end(), [&](int A, int B)->bool {return ev[A] < ev[B];});
+    // check handedness of eigenvectors and swap if needed
+    if (((evec[indices[0]] ^ evec[indices[1]])*evec[indices[2]]) < 0)
+        evec[indices[2]] = -evec[indices[2]];
+        
+        
+    mat3d fevec = mat3d(evec[indices[0]].x, evec[indices[0]].y, evec[indices[0]].z,
+                        evec[indices[1]].x, evec[indices[1]].y, evec[indices[1]].z,
+                        evec[indices[2]].x, evec[indices[2]].y, evec[indices[2]].z);
+    vec3d feval = vec3d(eval[indices[0]],eval[indices[1]],eval[indices[2]]);
+
+    // roundoff small numbers
+    const double eps = 1e-12;
+    fevec = CleanUp(fevec, eps);
+
+    // create a quaternion from the matrix of eigenvectors
+    quatd q(fevec);
+    q.MakeUnit();
+    
+    po->GetTransform().Rotate(q,vec3d(0,0,0));
+    vec3d com = use_area ? CalculateAreaCOM(*pm) : CalculateCOM(*pm);
+    po->GetTransform().Translate(-com);
+
+    FSMesh* newMesh = new FSMesh(*pm);
+    FSMeshBuilder meshBuilder(*newMesh);
+    
+    return newMesh;
+}

MeshTools/FEModifier.h

@@ -551,9 +551,19 @@ class FEInflateMesh: public FEModifier
 	FSMesh* Apply(FSMesh* pm);
 };
 
+//-----------------------------------------------------------------------------
 class FEDeleteElements : public FEModifier
 {
 public:
 	FEDeleteElements();
 	FSMesh* Apply(FSMesh* pm);
 };
+
+//-----------------------------------------------------------------------------
+class FEAlignMeshMOI: public FEModifier
+{
+public:
+    FEAlignMeshMOI();
+    FSMesh* Apply(FSMesh* pm);
+};
+