71 for (
int i = 0; i < this->getTargetDim(); ++i) {
75 T epsilon = 1. - this->_blockLattice.get(input).template getField<descriptors::POROSITY>();
78 if ((epsilon > 1e-5)) {
79 T rho, u[DESCRIPTOR::d], u_s[DESCRIPTOR::d];
81 for (
int i = 0; i < DESCRIPTOR::d; i++) {
82 u_s[i] = this->_blockLattice.get(input).template getFieldComponent<descriptors::VELOCITY_SOLID>(i);
84 T paramA = this->_converter.getLatticeRelaxationTime() - 0.5;
86 T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
89 T omega = 1. / this->_converter.getLatticeRelaxationTime();
91 this->_blockLattice.get(input).computeRhoU(rho, u);
93 const T uSqr_s = util::normSqr<T,DESCRIPTOR::d>(u_s);
94 T uSqr = util::normSqr<T,DESCRIPTOR::d>(u);
95 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
99 - this->_blockLattice.get(input[0], input[1], input[2])[iPop])
101 * (this->_blockLattice.get(input[0], input[1], input[2])[iPop]
113 for (
int i = 0; i < this->getTargetDim(); ++i) {
114 output[i] -= descriptors::c<DESCRIPTOR>(iPop,i) * omega_s;
118 for (
int i = 0; i < this->getTargetDim(); ++i) {
119 output[i] = this->_converter.getPhysForce(output[i] * paramB);