| 197 | } |
| 198 | |
| 199 | void TimeStepPF::solvePDConstraints() |
| 200 | { |
| 201 | Simulation *sim = Simulation::getCurrent(); |
| 202 | const unsigned int nFluids = sim->numberOfFluidModels(); |
| 203 | if (nFluids == 0) |
| 204 | return; |
| 205 | |
| 206 | // total number of active fluid particles |
| 207 | m_numActiveParticlesTotal = m_simulationData.getParticleOffset(nFluids - 1) + sim->getFluidModel(nFluids - 1)->numActiveParticles(); |
| 208 | |
| 209 | VectorXr x(3 * m_numActiveParticlesTotal); |
| 210 | VectorXr b(3 * m_numActiveParticlesTotal); |
| 211 | |
| 212 | for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) |
| 213 | { |
| 214 | FluidModel *model = sim->getFluidModel(fluidModelIndex); |
| 215 | const unsigned int offset = m_simulationData.getParticleOffset(fluidModelIndex); |
| 216 | |
| 217 | #pragma omp parallel for schedule(static) |
| 218 | for (int i = 0; i < (int)model->numActiveParticles(); i++) |
| 219 | { |
| 220 | ////////////////////////////////////////////////////////////////////////// |
| 221 | // initialize positions |
| 222 | ////////////////////////////////////////////////////////////////////////// |
| 223 | x.Vec3Block(offset + i) = m_simulationData.getS(fluidModelIndex, i); |
| 224 | |
| 225 | ////////////////////////////////////////////////////////////////////////// |
| 226 | // count number of fluid neighbors for relaxation |
| 227 | ////////////////////////////////////////////////////////////////////////// |
| 228 | unsigned int nNeighbors = 0; |
| 229 | for (unsigned int pid = 0; pid < nFluids; pid++) |
| 230 | nNeighbors += sim->numberOfNeighbors(fluidModelIndex, pid, i); |
| 231 | m_simulationData.setNumFluidNeighbors(fluidModelIndex, i, nNeighbors + 1u); |
| 232 | } |
| 233 | } |
| 234 | |
| 235 | ////////////////////////////////////////////////////////////////////////// |
| 236 | // Init linear system solver and preconditioner |
| 237 | ////////////////////////////////////////////////////////////////////////// |
| 238 | MatrixReplacement A(3 * m_numActiveParticlesTotal, matrixVecProd, (void*) this); |
| 239 | #ifdef PD_USE_DIAGONAL_PRECONDITIONER |
| 240 | preparePreconditioner(); |
| 241 | m_solver.preconditioner().init(m_numActiveParticlesTotal, diagonalMatrixElement, (void*)this); |
| 242 | #endif |
| 243 | m_solver.setMaxIterations(m_maxIterations); |
| 244 | m_solver.compute(A); |
| 245 | |
| 246 | for (m_iterations = 0u; m_iterations < m_maxIterations; m_iterations++) |
| 247 | { |
| 248 | ////////////////////////////////////////////////////////////////////////// |
| 249 | // Compute RHS |
| 250 | ////////////////////////////////////////////////////////////////////////// |
| 251 | matrixFreeRHS(x,b); |
| 252 | |
| 253 | ////////////////////////////////////////////////////////////////////////// |
| 254 | // Solve linear system |
| 255 | ////////////////////////////////////////////////////////////////////////// |
| 256 | #ifdef PD_USE_DIAGONAL_PRECONDITIONER |
nothing calls this directly
no test coverage detected