Function computeMatrixL

Tools/MeshSkinning/main.cpp:640–683

void computeMatrixL()
{
	const int numParticles = (int) x0.size();

	#pragma omp parallel default(shared)
	{
		#pragma omp for schedule(static)  
		for (int i = 0; i < numParticles; i++)
		{
			const Vector3r &xi0 = x0[i];
			Matrix3r L_;
			L_.setZero();

			const size_t numNeighbors = initialNeighbors[i].size();

			//////////////////////////////////////////////////////////////////////////
			// Fluid
			//////////////////////////////////////////////////////////////////////////
			for (unsigned int j = 0; j < numNeighbors; j++)
			{
				const unsigned int neighborIndex = initialNeighbors[i][j];

				const Vector3r &xj0 = x0[neighborIndex];
				const Vector3r xj_xi_0 = xj0 - xi0;
				const Vector3r gradW = gradKernelFct(xj_xi_0);

				// minus because gradW(xij0) == -gradW(xji0)
				L_ -= restVolumes[neighborIndex] * gradW * xj_xi_0.transpose();
			}

			//// add 1 to z-component. otherwise we get a singular matrix in 2D
			//if (sim->is2DSimulation())
			//	L(2, 2) = 1.0;

			bool invertible = false;
			L_.computeInverseWithCheck(L[i], invertible, 1e-9);
			if (!invertible)
			{
				//MathFunctions::pseudoInverse(L, m_L[i]);
				L[i] = Matrix3r::Identity();
			}
		}
	}
}

#ifdef USE_AVX