53{
54 constexpr unsigned D = DESCRIPTOR::d;
55
56 using namespace descriptors;
57 using namespace eler;
58 static_assert(PARTICLETYPE::template providesNested<MOBILITY, FLUIDVEL>(), "Field MOBILITY:FLUIDVEL has to be provided");
59
60 Vector<T,D> position = particle.template getField<GENERAL,POSITION>();
61 Vector<T,D> velocity = particle.template getField<MOBILITY,VELOCITY>();
62 T* positionArray = position.data();
63
64
65 const auto& cuboid = _blockGeometry.
getCuboid();
67 _blockLattice, this->_converter, cuboid);
68
69
70
71
72
73
74 T fluidVelArray[D] = {0.};
75
76
77 bool inside = true;
79 inside = cuboid.checkPoint(position);
80 }
81 if (inside){
82
83
84 blockInterpPhysVelF(fluidVelArray, positionArray);
85 particle.template setField<MOBILITY, FLUIDVEL>(fluidVelArray);
86 Vector<T,3> fvel = particle.template getField<MOBILITY, FLUIDVEL>();
87
88 Vector<T,9> rotmat = particle.template getField<NUMERICPROPERTIES,ROT_MATRIX>();
89
91Khathat = eler::matrixMultiply(eler::inverseTransformMatrix(rotmat),
92eler::matrixMultiply(particle.template getField<NUMERICPROPERTIES,TRANSLATIONAL_DYADIC>(), rotmat));
93Vector<T,3> Khathatvelo (0.);
94Khathatvelo = eler::matrixVectorMultiply(Khathat, fvel - velocity);
95
96
97
98
99
100
101Vector<T,3> F12, F13, F23, F21, F31, F32;
102Vector<T,3> force = Vector<T,3> (0.,0.,0.);
104Vector<T,9> invB12, invB13, invB23, invB21, invB31, invB32;
106invB12 = B12;
117
118Vector<T,9> L (0.0501, 0.0329, 0.0, 0.0182, 0.0173, 0.00, 0.00, 0.00, 0.0373);
119
120
121 Vector <T,3> posVector (0.);
122 posVector = positionArray;
123 Vector<T,3> x (1.,0.,0.);
124 Vector<T,3> y (0.,1.,0.);
125 Vector<T,3> z (0.,0.,1.);
126 Vector<T,3> gradx (0.);
127 Vector<T,3> grady (0.);
128 Vector<T,3> gradz (0.);
129 veloGradientInterpFD<T,DESCRIPTOR>(gradx,blockInterpPhysVelF, x,posVector, this->_converter.getPhysDeltaX()/100.);
130 veloGradientInterpFD<T,DESCRIPTOR>(grady,blockInterpPhysVelF, y,posVector, this->_converter.getPhysDeltaX()/100.);
131 veloGradientInterpFD<T,DESCRIPTOR>(gradz,blockInterpPhysVelF, z,posVector, this->_converter.getPhysDeltaX()/100.);
132
133
134
135
136
138 gradvel[0] = gradx[0];
139 gradvel[3] = gradx[1];
140 gradvel[6] = gradx[2];
141 gradvel[1] = grady[0];
142 gradvel[4] = grady[1];
143 gradvel[7] = grady[2];
144 gradvel[2] = gradz[0];
145 gradvel[5] = gradz[1];
146 gradvel[8] = gradz[2];
147
148using namespace eler;
149using namespace util;
150
151T lim = std::numeric_limits<T>::min();
152
153T Const =
M_PI*
M_PI*particle.template getField<PHYSPROPERTIES,RADIUS>()*particle.template getField<PHYSPROPERTIES,RADIUS>()*this->_converter.getPhysViscosity()*this->_converter.getPhysDensity()/sqrt(this->_converter.getPhysViscosity());
154
155if(std::fabs(gradvel[1]) < lim )
156F12 = Vector<T,3> (0.,0.,0.);
157else
158F12 = Const*gradvel[1]/sqrt(std::fabs(gradvel[1]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B12,matrixMultiply(L,matrixMultiply(invB12,Khathat)))),fvel - velocity);
159if(std::fabs(gradvel[2]) < lim )
160F13 = Vector<T,3> (0.,0.,0.);
161else
162F13 = Const*gradvel[2]/sqrt(std::fabs(gradvel[2]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B13,matrixMultiply(L,matrixMultiply(invB13,Khathat)))),fvel - velocity) ;
163if(std::fabs(gradvel[3]) < lim )
164F21 = Vector<T,3> (0.,0.,0.);
165else
166F21 = Const*gradvel[3]/sqrt(std::fabs(gradvel[3]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B21,matrixMultiply(L,matrixMultiply(invB21,Khathat)))),fvel - velocity) ;
167if(std::fabs(gradvel[5]) < lim )
168F23 = Vector<T,3> (0.,0.,0.);
169else
170F23 = Const*gradvel[5]/sqrt(std::fabs(gradvel[5]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B23,matrixMultiply(L,matrixMultiply(invB23,Khathat)))),fvel - velocity) ;
171if(std::fabs(gradvel[6]) < lim )
172F31 = Vector<T,3> (0.,0.,0.);
173else
174F31 = Const*gradvel[6]/sqrt(std::fabs(gradvel[6]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B31,matrixMultiply(L,matrixMultiply(invB31,Khathat)))),fvel - velocity) ;
175if(std::fabs(gradvel[7]) < lim )
176F32 = Vector<T,3> (0.,0.,0.);
177else
178F32 = Const*gradvel[7]/sqrt(std::fabs(gradvel[7]))*matrixVectorMultiply(matrixMultiply(Khathat,matrixMultiply(B32,matrixMultiply(L,matrixMultiply(invB32,Khathat)))),fvel - velocity) ;
179
180
181force = F12 + F13 + F21 + F23 + F31 + F32;
182
183 particle.template setField<FORCING,FORCE> ( particle.template getField<FORCING,FORCE>()+force);
184
185 }
186}
Cuboid< T, D > & getCuboid()
constexpr bool providesParallelization()
Provides group PARALLELIZATION.
std::conditional_t< DESCRIPTOR::d==2, BlockLatticeInterpPhysVelocity2D< T, DESCRIPTOR >, BlockLatticeInterpPhysVelocity3D< T, DESCRIPTOR > > BlockLatticeInterpPhysVelocity
Vector(T &&t, Ts &&... ts) -> Vector< std::remove_cvref_t< T >, 1+sizeof...(Ts)>