1882 {
1883 V uNew[DESCRIPTOR::d] { };
1884 const V epsilon = cell.template getField<descriptors::EPSILON>();
1885 const V nu = cell.template getField<descriptors::NU>();
1886 const V k = cell.template getField<descriptors::K>();
1887 const auto bodyF = cell.template getFieldPointer<descriptors::BODY_FORCE>();
1888
1890 const V Fe = 0.;
1891 const V c_0 = 0.5*(1 + 0.5*epsilon*nu/k);
1893
1894 for (unsigned iD=0; iD < DESCRIPTOR::d; ++iD) {
1895 uNew[iD] = u[iD] * (c_0 +
util::sqrt(c_0*c_0 + c_1*uMag)) - V{0.5} * bodyF[iD] * epsilon;
1896 }
1899
1900 int iPi = 0;
1901 for (int iAlpha=0; iAlpha < DESCRIPTOR::d; ++iAlpha) {
1902 for (int iBeta=iAlpha; iBeta < DESCRIPTOR::d; ++iBeta) {
1903 forceTensor[iPi] = V{0.5} * rho * (bodyF[iAlpha]*uNew[iBeta] + uNew[iAlpha]*bodyF[iBeta]);
1904 ++iPi;
1905 }
1906 }
1907
1908 for (int iPi=0; iPi < util::TensorVal<DESCRIPTOR>::n; ++iPi) {
1909 pi[iPi] += forceTensor[iPi];
1910 }
1911 }
auto normSqr(const ARRAY_LIKE &u) any_platform
Compute norm square of a d-dimensional vector.
void compute(CELL &cell, const RHO &rho, const U &u, PI &pi) any_platform
static constexpr int n
result stored in n