Removed linesearch timer and added timings to update timer. Overrode some rigid domain functions to avoid performance penalty. Made performance improvements to FEOctree and FENormalProjection. Added omp pragmas in certain places.

Author: SteveMaas1978

FEBioLib/FEBioModel.cpp

@@ -1785,7 +1785,6 @@ TimingInfo FEBioModel::GetTimingInfo()
 	ti.total_stiff        = GetTimer(TimerID::Timer_Stiffness       )->GetExclusiveTime(); total += ti.total_stiff;
 	ti.total_rhs          = GetTimer(TimerID::Timer_Residual        )->GetExclusiveTime(); total += ti.total_rhs;
 	ti.total_update       = GetTimer(TimerID::Timer_Update          )->GetExclusiveTime(); total += ti.total_update;
-	ti.total_line_search  = GetTimer(TimerID::Timer_LineSearch      )->GetExclusiveTime(); total += ti.total_line_search;
 	ti.total_qn           = GetTimer(TimerID::Timer_QNUpdate        )->GetExclusiveTime(); total += ti.total_qn;
 	ti.total_serialize    = GetTimer(TimerID::Timer_Serialize       )->GetExclusiveTime(); total += ti.total_serialize;
 	ti.total_callback     = GetTimer(TimerID::Timer_Callback        )->GetExclusiveTime(); total += ti.total_callback;
@@ -1858,7 +1857,6 @@ void FEBioModel::on_cb_solved()
 		Timer::time_str(ti.total_stiff      , sztime); feLog("\t   evaluating stiffness ......... : %s (%lg sec)\n\n", sztime, ti.total_stiff);
 		Timer::time_str(ti.total_rhs        , sztime); feLog("\t   evaluating residual .......... : %s (%lg sec)\n\n", sztime, ti.total_rhs);
 		Timer::time_str(ti.total_update     , sztime); feLog("\t   model update ................. : %s (%lg sec)\n\n", sztime, ti.total_update);
-		Timer::time_str(ti.total_line_search, sztime); feLog("\t   line search .................. : %s (%lg sec)\n\n", sztime, ti.total_line_search);
 		Timer::time_str(ti.total_qn         , sztime); feLog("\t   QN updates ................... : %s (%lg sec)\n\n", sztime, ti.total_qn);
 		Timer::time_str(linsol_time         , sztime); feLog("\t   time in linear solver ........ : %s (%lg sec)\n\n", sztime, linsol_time);
 		Timer::time_str(ti.total_time       , sztime); feLog("\tTotal elapsed time .............. : %s (%lg sec)\n\n", sztime, ti.total_time);

FEBioLib/febiolib_types.h

@@ -46,7 +46,6 @@ struct TimingInfo {
 	double total_rhs = 0;
 	double total_update = 0;
 	double total_qn = 0;
-	double total_line_search = 0;
 	double total_serialize = 0;
 	double total_callback = 0;
 	double total_other = 0;

FEBioMech/FERigidShellDomain.cpp

@@ -359,6 +359,22 @@ void FERigidShellDomain::MassMatrix(FELinearSystem& LS, double scale)
 	// nothing to do here
 }
 
+void FERigidShellDomain::Serialize(DumpStream& ar)
+{
+	// we bypass FEDomain::Serialize since we don't need to serialize anything for the rigid shells
+	FEMeshPartition::Serialize(ar);
+}
+
+void FERigidShellDomain::PreSolveUpdate(const FETimeInfo& timeInfo)
+{
+	// nothing to do here
+}
+
+void FERigidShellDomain::BuildMatrixProfile(FEGlobalMatrix& M)
+{
+	// nothing to do here
+}
+
 double FERigidShellDomain::detJ0(FEShellElement& el, int n)
 {
 	vector<vec3d> X(FEElement::MAX_NODES);

FEBioMech/FERigidShellDomain.h

@@ -98,6 +98,12 @@ class FERigidShellDomain : public FEShellDomain, public FEElasticDomain
 	//! get the material (overridden from FEDomain)
 	FEMaterial* GetMaterial() override;
 
+	void Serialize(DumpStream& ar) override;
+
+	void PreSolveUpdate(const FETimeInfo& timeInfo) override;
+
+	void BuildMatrixProfile(FEGlobalMatrix& M) override;
+
 public: // from FEShellDomain
 	// get a shell element
 	FEShellElement& Element(int i) override { return m_Elem[i]; }

FEBioMech/FERigidSolidDomain.cpp

@@ -64,6 +64,22 @@ void FERigidSolidDomain::Reset()
 	// nothing to reset here
 }
 
+void FERigidSolidDomain::Serialize(DumpStream& ar)
+{
+	// we bypass FEDomain::Serialize since we don't need to serialize anything for the rigid solids
+	FEMeshPartition::Serialize(ar);
+}
+
+void FERigidSolidDomain::PreSolveUpdate(const FETimeInfo& timeInfo)
+{
+	// nothing to do here
+}
+
+void FERigidSolidDomain::BuildMatrixProfile(FEGlobalMatrix& M)
+{
+	// nothing to do here
+}
+
 //-----------------------------------------------------------------------------
 //! Calculates the stiffness matrix for 3D rigid elements.
 //! Rigid elements don't generate stress, so there is nothing to do here

FEBioMech/FERigidSolidDomain.h

@@ -44,6 +44,13 @@ class FERigidSolidDomain : public FEElasticSolidDomain
 	//! reset data
 	void Reset() override;
 
+	//! serialization
+	void Serialize(DumpStream& ar) override;
+
+	void PreSolveUpdate(const FETimeInfo& timeInfo) override;
+
+	void BuildMatrixProfile(FEGlobalMatrix& M) override;
+
 public:
 
 	//! calculates the global stiffness matrix for this domain

FEBioMech/FERigidSolver.cpp

@@ -1514,6 +1514,7 @@ void FERigidSolverNew::UpdateRigidBodies(vector<double>& Ui, vector<double>& ui)
 	// Since the rigid nodes are repositioned we need to update the displacement DOFS
 	FEMesh& mesh = m_fem->GetMesh();
 	int N = mesh.Nodes();
+#pragma omp parallel for schedule(dynamic, 64)
 	for (int i = 0; i<N; ++i)
 	{
 		FENode& node = mesh.Node(i);

FEBioMech/FESolidSolver2.cpp

@@ -663,37 +663,39 @@ void FESolidSolver2::UpdateIncrements(vector<double>& Ui, vector<double>& ui, bo
 	m_rigidSolver.UpdateIncrements(Ui, ui, emap);
 
 	// update flexible nodes
-	int n;
-	for (int i=0; i<mesh.Nodes(); ++i)
+	int NN = mesh.Nodes();
+#pragma omp parallel for 
+	for (int i=0; i<NN; ++i)
 	{
 		FENode& node = mesh.Node(i);
+		const std::vector<int>& id = node.m_ID;
+		int n;
 
 		// displacement dofs
 		// current position = initial + total at prev conv step + total increment so far + current increment
-		if ((n = node.m_ID[m_dofU[0]]) >= 0) Ui[n] += ui[n];
-		if ((n = node.m_ID[m_dofU[1]]) >= 0) Ui[n] += ui[n];
-		if ((n = node.m_ID[m_dofU[2]]) >= 0) Ui[n] += ui[n];
+		if ((n = id[m_dofU[0]]) >= 0) Ui[n] += ui[n];
+		if ((n = id[m_dofU[1]]) >= 0) Ui[n] += ui[n];
+		if ((n = id[m_dofU[2]]) >= 0) Ui[n] += ui[n];
 
 		// beam rotations
+		if ((id[m_dofQ[0]] >= 0) || (id[m_dofQ[1]] >= 0) || (id[m_dofQ[2]] >= 0))
 		{
 			vec3d ri, Ri;
-			if ((n = node.m_ID[m_dofQ[0]]) >= 0) { ri.x = ui[n]; Ri.x = Ui[n]; }
-			if ((n = node.m_ID[m_dofQ[1]]) >= 0) { ri.y = ui[n]; Ri.y = Ui[n]; }
-			if ((n = node.m_ID[m_dofQ[2]]) >= 0) { ri.z = ui[n]; Ri.z = Ui[n]; }
+			if ((n = id[m_dofQ[0]]) >= 0) { ri.x = ui[n]; Ri.x = Ui[n]; }
+			if ((n = id[m_dofQ[1]]) >= 0) { ri.y = ui[n]; Ri.y = Ui[n]; }
+			if ((n = id[m_dofQ[2]]) >= 0) { ri.z = ui[n]; Ri.z = Ui[n]; }
 			quatd qi(ri), Qi(Ri);
 			quatd Qn = qi * Qi;
 			vec3d rn = Qn.GetRotationVector();
-			if ((n = node.m_ID[m_dofQ[0]]) >= 0) { Ui[n] = rn.x; }
-			if ((n = node.m_ID[m_dofQ[1]]) >= 0) { Ui[n] = rn.y; }
-			if ((n = node.m_ID[m_dofQ[2]]) >= 0) { Ui[n] = rn.z; }
+			if ((n = id[m_dofQ[0]]) >= 0) { Ui[n] = rn.x; }
+			if ((n = id[m_dofQ[1]]) >= 0) { Ui[n] = rn.y; }
+			if ((n = id[m_dofQ[2]]) >= 0) { Ui[n] = rn.z; }
 		}
 
 		// shell dofs
-		{
-			if ((n = node.m_ID[m_dofSU[0]]) >= 0) Ui[n] += ui[n];
-			if ((n = node.m_ID[m_dofSU[1]]) >= 0) Ui[n] += ui[n];
-			if ((n = node.m_ID[m_dofSU[2]]) >= 0) Ui[n] += ui[n];
-		}
+		if ((n = id[m_dofSU[0]]) >= 0) Ui[n] += ui[n];
+		if ((n = id[m_dofSU[1]]) >= 0) Ui[n] += ui[n];
+		if ((n = id[m_dofSU[2]]) >= 0) Ui[n] += ui[n];
 	}
 
 	for (int i = 0; i < fem.NonlinearConstraints(); ++i)
@@ -727,23 +729,26 @@ void FESolidSolver2::UpdateIncrements(vector<double>& Ui, vector<double>& ui, bo
 //! Updates the current state of the model
 void FESolidSolver2::Update(vector<double>& ui)
 {
-    FEModel& fem = *GetFEModel();
-    FETimeInfo& tp = fem.GetTime();
-    tp.currentIteration = m_niter;
-    
-    // update EAS
-    UpdateEAS(ui);
-    UpdateIncrementsEAS(ui, true);
+	{
+		TRACK_TIME(Timer_Update)
+			FEModel& fem = *GetFEModel();
+		FETimeInfo& tp = fem.GetTime();
+		tp.currentIteration = m_niter;
 
-	// update kinematics
-	UpdateKinematics(ui);
+		// update EAS
+		UpdateEAS(ui);
+		UpdateIncrementsEAS(ui, true);
 
-	// update domains 
-	FEMesh& mesh = fem.GetMesh();
-	for (int i = 0; i < mesh.Domains(); ++i)
-	{
-		FEDomain& dom = mesh.Domain(i);
-		dom.IncrementalUpdate(ui, false);
+		// update kinematics
+		UpdateKinematics(ui);
+
+		// update domains 
+		FEMesh& mesh = fem.GetMesh();
+		for (int i = 0; i < mesh.Domains(); ++i)
+		{
+			FEDomain& dom = mesh.Domain(i);
+			dom.IncrementalUpdate(ui, false);
+		}
 	}
 
 	// update model state
@@ -864,7 +869,9 @@ void FESolidSolver2::PrepStep()
 
 		// store previous mesh state
 		// we need them for velocity and acceleration calculations
-		for (int i = 0; i < mesh.Nodes(); ++i)
+		int NN = mesh.Nodes();
+#pragma omp parallel for
+		for (int i = 0; i < NN; ++i)
 		{
 			FENode& ni = mesh.Node(i);
 			ni.m_rp = ni.m_rt;
@@ -1616,25 +1623,27 @@ void FESolidSolver2::ExternalForces(FEGlobalVector& RHS)
 
 	// set the nodal reaction forces
 	// TODO: Is this a good place to do this?
-	for (int i = 0; i<mesh.Nodes(); ++i)
+	int NN = mesh.Nodes();
+#pragma omp parallel for
+	for (int i = 0; i<NN; ++i)
 	{
 		FENode& node = mesh.Node(i);
-		node.set_load(m_dofU[0], 0);
-		node.set_load(m_dofU[1], 0);
-		node.set_load(m_dofU[2], 0);
 
 		int n;
 		if ((n = node.m_ID[m_dofU[0]]) >= 0) node.set_load(m_dofU[0], -m_Fr[n]);
-		if ((n = -node.m_ID[m_dofU[0]] - 2) >= 0) node.set_load(m_dofU[0], -m_Fr[n]);
+		else if ((n = -node.m_ID[m_dofU[0]] - 2) >= 0) node.set_load(m_dofU[0], -m_Fr[n]);
+		else node.set_load(m_dofU[0], 0);
 
 		if ((n = node.m_ID[m_dofU[1]]) >= 0) node.set_load(m_dofU[1], -m_Fr[n]);
-		if ((n = -node.m_ID[m_dofU[1]] - 2) >= 0) node.set_load(m_dofU[1], -m_Fr[n]);
+		else if ((n = -node.m_ID[m_dofU[1]] - 2) >= 0) node.set_load(m_dofU[1], -m_Fr[n]);
+		else node.set_load(m_dofU[1], 0);
 
 		if ((n = node.m_ID[m_dofU[2]]) >= 0) node.set_load(m_dofU[2], -m_Fr[n]);
-		if ((n = -node.m_ID[m_dofU[2]] - 2) >= 0) node.set_load(m_dofU[2], -m_Fr[n]);
+		else if ((n = -node.m_ID[m_dofU[2]] - 2) >= 0) node.set_load(m_dofU[2], -m_Fr[n]);
+		else node.set_load(m_dofU[2], 0);
 
 		// add nodal loads
-		double s = (m_arcLength>0 ? m_al_lam : 1.0);
+//		double s = (m_arcLength>0 ? m_al_lam : 1.0);
 //		if ((n = node.m_ID[m_dofU[0]]) >= 0) node.set_load(m_dofU[0], -m_Fn[n]*s);
 //		if ((n = node.m_ID[m_dofU[1]]) >= 0) node.set_load(m_dofU[1], -m_Fn[n]*s);
 //		if ((n = node.m_ID[m_dofU[2]]) >= 0) node.set_load(m_dofU[2], -m_Fn[n]*s);

FECore/FENewtonSolver.cpp

@@ -882,7 +882,6 @@ double FENewtonSolver::DoLineSearch()
 	m_ls = 1.0;
 	if (m_lineSearch && (m_lineSearch->m_LStol > 0.0))
 	{
-		TRACK_TIME(TimerID::Timer_LineSearch);
 		m_ls = m_lineSearch->DoLineSearch();
 	}
 	else

FECore/FENormalProjection.cpp

@@ -50,6 +50,8 @@ void FENormalProjection::Init()
 //! of the intersection point. It searches for the closest patch based on
 //! algebraic value of the gap function
 //!
+//!
+/*
 FESurfaceElement* FENormalProjection::Project(vec3d r, vec3d n, double rs[2])
 {
 	// let's find all the candidate surface elements
@@ -91,6 +93,44 @@ FESurfaceElement* FENormalProjection::Project(vec3d r, vec3d n, double rs[2])
 	// we did not find a surface
 	return 0;
 }
+*/
+
+FESurfaceElement* FENormalProjection::Project(vec3d r, vec3d n, double rs[2])
+{
+	bool found = false;
+
+	// now that we found candidate surface elements, lets see if we can find 
+	// those that intersect the ray, then pick the closest intersection
+	double rsl[2], gl, g = 0;
+	FESurfaceElement* pei = 0;
+	m_OT.VisitIntersectedLeaves(r, n, m_rad, [&](int j) {
+		// get the surface element
+		// project the node on the element
+		FESurfaceElement* pe = &m_surf.Element(j);
+		if (pe->isActive())
+		{
+			if (m_surf.Intersect(*pe, r, n, rsl, gl, m_tol)) {
+				if ((!found) && (gl > -m_rad)) {
+					found = true;
+					g = gl;
+					rs[0] = rsl[0];
+					rs[1] = rsl[1];
+					pei = pe;
+				}
+				else if ((gl < g) && (gl > -m_rad)) {
+					g = gl;
+					rs[0] = rsl[0];
+					rs[1] = rsl[1];
+					pei = pe;
+				}
+			}
+		}
+		});
+	if (found) return pei;
+
+	// we did not find a surface
+	return 0;
+}
 
 //-----------------------------------------------------------------------------
 //! This function finds the element which is intersected by the ray (r,n).