1404 {
1405 V uNew[DESCRIPTOR::d] { };
1406 const V epsilon = cell.template getField<descriptors::EPSILON>();
1407 const V nu = cell.template getField<descriptors::NU>();
1408 const V k = cell.template getField<descriptors::K>();
1409 const auto bodyF = cell.template getFieldPointer<descriptors::BODY_FORCE>();
1410
1411 const V uMag =
util::sqrt( util::normSqr<V,DESCRIPTOR::d>(u) );
1412 const V Fe = 0.;
1413 const V c_0 = 0.5*(1 + 0.5*epsilon*nu/k);
1415
1416 for (unsigned iD=0; iD < DESCRIPTOR::d; ++iD) {
1417 uNew[iD] = u[iD] * (c_0 +
util::sqrt(c_0*c_0 + c_1*uMag)) - V{0.5} * bodyF[iD] * epsilon;
1418 }
1419 STRESS().template compute<TYPE>(cell, rho, uNew, pi);
1421
1422 int iPi = 0;
1423 for (int iAlpha=0; iAlpha < DESCRIPTOR::d; ++iAlpha) {
1424 for (int iBeta=iAlpha; iBeta < DESCRIPTOR::d; ++iBeta) {
1425 forceTensor[iPi] = V{0.5} * rho * (bodyF[iAlpha]*uNew[iBeta] + uNew[iAlpha]*bodyF[iBeta]);
1426 ++iPi;
1427 }
1428 }
1429
1430 for (int iPi=0; iPi < util::TensorVal<DESCRIPTOR>::n; ++iPi) {
1431 pi[iPi] += forceTensor[iPi];
1432 }
1433 }
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
static constexpr int n
result stored in n