Line data Source code
1 : // Copyright (C) 2015 Technische Universitaet Muenchen
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder, or at
4 : // www5.in.tum.de/mamico
5 : #ifndef _MOLECULARDYNAMICS_COUPLING_CELLMAPPINGS_VELOCITYGRADIENTRELAXATIONMAPPING_H_
6 : #define _MOLECULARDYNAMICS_COUPLING_CELLMAPPINGS_VELOCITYGRADIENTRELAXATIONMAPPING_H_
7 :
8 : #include "coupling/CouplingMDDefinitions.h"
9 : #include "coupling/datastructures/CouplingCell.h"
10 : #include "coupling/indexing/IndexingService.h"
11 : #include "coupling/interface/MDSolverInterface.h"
12 : #include "coupling/interface/Molecule.h"
13 : #include "tarch/la/Matrix.h"
14 : #include <iostream>
15 :
16 : namespace coupling {
17 : namespace cellmappings {
18 : template <class LinkedCell, unsigned int dim> class VelocityGradientRelaxationMapping;
19 : template <class LinkedCell, unsigned int dim> class VelocityGradientRelaxationTopOnlyMapping;
20 : } // namespace cellmappings
21 : } // namespace coupling
22 :
23 : /** This class relaxes velocities of molecules towards an interpolated avergaged
24 : *velocity value. Only molecules within a three cell-wide boundary strip are
25 : *considered if they are inside the one cell-wide boundary strip that is spanned
26 : *by the outermost non-ghost coupling cell midpoints. Currently, a
27 : *second-order interpolation of the molecules in the outer boundary strip is
28 : *used (that is the first layer with three coupling cells needs to be
29 : *initialised with valid LB velocities).
30 : * @brief This class relaxes velocities of molecules towards an
31 : *interpolated avergaged velocity value.
32 : * @tparam LinkedCell cell type
33 : * @tparam dim Number of dimensions; it can be 1, 2 or 3
34 : * @author Philipp Neumann
35 : * @note ONLY SUPPORTS 3D AT THE MOMENT!
36 : * \todo Philipp please take a look on this class
37 : */
38 : template <class LinkedCell, unsigned int dim> class coupling::cellmappings::VelocityGradientRelaxationMapping {
39 : public:
40 : /** obtains relaxation factor and current velocity (the velocity towards which
41 : *the relaxation shall be done is later obtained from the
42 : *microscopicmomentum-buffers).
43 : * @param velocityRelaxationFactor
44 : * @param currentVelocity
45 : * @param cellIndex
46 : * @param mdSolverInterface
47 : * @param couplingCells
48 : */
49 0 : VelocityGradientRelaxationMapping(const double& velocityRelaxationFactor, const tarch::la::Vector<dim, double>& currentVelocity, const I02& cellIndex,
50 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface,
51 : const coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>* const couplingCells)
52 0 : : _mdSolverInterface(mdSolverInterface), _couplingCells(couplingCells), _velocityRelaxationFactor(velocityRelaxationFactor),
53 0 : _currentVelocity(currentVelocity), _innerLowerLeft(getInnerLowerLeftCorner()), _innerUpperRight(getInnerUpperRightCorner()),
54 0 : _outerLowerLeft(getOuterLowerLeftCorner()), _outerUpperRight(getOuterUpperRightCorner()), _cellIdx(cellIndex),
55 0 : _ignoreThisCell(ignoreThisCell(cellIndex)) {}
56 :
57 : /** Destructor */
58 0 : virtual ~VelocityGradientRelaxationMapping() {}
59 :
60 : /** empty function
61 : */
62 : void beginCellIteration() {}
63 :
64 : /** empty function
65 : */
66 : void endCellIteration() {}
67 :
68 : /** computes the current velocity directly from coupling cell and the new
69 : *velocity (microscopicMomentum) with second-order interpolation multiplies
70 : *the difference between the two with the velocity relaxation factor add it to
71 : *the velocity of the molecule and applies it to the molecule.
72 : * @param cell
73 : */
74 0 : void handleCell(LinkedCell& cell) {
75 : // if this coupling cell is not interesting, skip it immediately
76 0 : if (_ignoreThisCell) {
77 : return;
78 : }
79 :
80 : // std::cout << "Handle VelocityGradientRelaxation in cell " <<
81 : // I03{_cellIdx}.get() << "=" << _cellIdx <<
82 : // std::endl;
83 0 : coupling::interface::MoleculeIterator<LinkedCell, dim>* it = _mdSolverInterface->getMoleculeIterator(cell);
84 0 : it->begin();
85 0 : while (it->continueIteration()) {
86 0 : coupling::interface::Molecule<dim>& wrapper(it->get());
87 0 : tarch::la::Vector<dim, double> velocity = wrapper.getVelocity();
88 0 : const tarch::la::Vector<dim, double> position(wrapper.getPosition());
89 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
90 : std::cout << "VelocityGradientRelaxationMapping: Process cell " << _cellIdx << ", molecule " << position << std::endl;
91 : #endif
92 0 : if (relaxMolecule(position)) {
93 0 : tarch::la::Vector<dim, double> normalisedPosition(0.0);
94 0 : bool secondCake = false;
95 : // index of coupling cell values that are involved in interpolation
96 : // process
97 0 : I02 couplingCellIndex[10];
98 : // determine cell indices and the "correct cake" for interpolation
99 0 : createCouplingCellIndex4SecondOrderInterpolation(position, normalisedPosition, secondCake, couplingCellIndex);
100 :
101 0 : tarch::la::Vector<dim, double> newVelocity(0.0);
102 0 : tarch::la::Vector<dim, double> oldVelocity(0.0);
103 :
104 : // get newVelocity (microscopicMomentum) with second-order interpolation
105 0 : interpolateVelocitySecondOrder(couplingCellIndex, normalisedPosition, secondCake, newVelocity);
106 : // get current velocity directly from coupling cell
107 0 : oldVelocity = _couplingCells[_cellIdx.get()].getCurrentVelocity();
108 :
109 0 : velocity += _velocityRelaxationFactor * (newVelocity - oldVelocity);
110 0 : wrapper.setVelocity(velocity);
111 :
112 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
113 : for (unsigned int d = 0; d < dim; d++) {
114 : if (normalisedPosition[d] > 2.0) {
115 : std::cout << "ERROR cellmappings::VelocityGradientRelaxationMapping: "
116 : "normalisedPosition>2.0! "
117 : << "Local vector index=" << I03{_cellIdx} << " , global vector index=" << I01{_cellIdx} << ", currentLocalCellIndex=" << _cellIdx
118 : << ", Norm. position=" << normalisedPosition << ", Actual position=" << position << std::endl;
119 : for (unsigned int i = 0; i < 10; i++) {
120 : std::cout << couplingCellIndex[i] << " ";
121 : }
122 : std::cout << std::endl;
123 : exit(EXIT_FAILURE);
124 : }
125 : if (normalisedPosition[d] < 0.0) {
126 : std::cout << "ERROR cellmappings::VelocityGradientRelaxationMapping: "
127 : "normalisedPosition<0.0! "
128 : << I03{_cellIdx} << " ," << normalisedPosition << std::endl;
129 : for (unsigned int i = 0; i < 10; i++) {
130 : std::cout << couplingCellIndex[i] << " ";
131 : }
132 : std::cout << std::endl;
133 : exit(EXIT_FAILURE);
134 : }
135 : }
136 : if ((tarch::la::dot(newVelocity, newVelocity) > 1000.0) || (tarch::la::dot(oldVelocity, oldVelocity) > 1000.0)) {
137 : for (unsigned int i = 0; i < 10; i++) {
138 : std::cout << couplingCellIndex[i] << " ";
139 : }
140 : std::cout << std::endl;
141 : std::cout << "ERROR: velocity to high! " << oldVelocity << " " << newVelocity << " " << I03{_cellIdx} << std::endl;
142 : exit(EXIT_FAILURE);
143 : }
144 : #endif
145 : }
146 0 : it->next();
147 : }
148 0 : delete it;
149 : }
150 :
151 : protected:
152 : /**
153 : * @param position
154 : * @returns true, if the position 'position' is within the respective
155 : *boundary strip that is spanned by the two outermost coupling cell
156 : *midpoints
157 : */
158 0 : virtual bool relaxMolecule(const tarch::la::Vector<dim, double>& position) const {
159 : // check if molecule is in the respective boundary strip
160 0 : bool isInsideOuter = true;
161 0 : bool isOutsideInner = false;
162 0 : for (unsigned int d = 0; d < dim; d++) {
163 0 : isInsideOuter = isInsideOuter && (position[d] > _outerLowerLeft[d]) && (position[d] < _outerUpperRight[d]);
164 0 : isOutsideInner = isOutsideInner || (position[d] < _innerLowerLeft[d]) || (position[d] > _innerUpperRight[d]);
165 : }
166 0 : return isInsideOuter && isOutsideInner;
167 : }
168 :
169 : /**
170 : * @param idx
171 : * @returns true if all molecules in this coupling cell do not require
172 : *any velocity relaxation (-> used to speed up the code)
173 : */
174 0 : bool ignoreThisCell(const I02& idx) const {
175 0 : const tarch::la::Vector<dim, unsigned int> globalIndex{I01{idx}.get()};
176 0 : bool innerCell = true;
177 0 : for (unsigned int d = 0; d < dim; d++) {
178 0 : innerCell = innerCell && ((globalIndex[d] > 2) && (globalIndex[d] < I01::numberCellsInDomain[d] - 3));
179 : }
180 0 : return innerCell;
181 : }
182 :
183 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const _mdSolverInterface;
184 : const coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>* const _couplingCells;
185 : /** relaxation factor for velocity relaxation */
186 : const double _velocityRelaxationFactor;
187 : /** current velocity in this coupling cell */
188 : const tarch::la::Vector<dim, double> _currentVelocity;
189 :
190 : /** boundary points of the boundary strip under consideration */
191 : const tarch::la::Vector<dim, double> _innerLowerLeft;
192 : const tarch::la::Vector<dim, double> _innerUpperRight;
193 : const tarch::la::Vector<dim, double> _outerLowerLeft;
194 : const tarch::la::Vector<dim, double> _outerUpperRight;
195 :
196 : I02 _cellIdx;
197 : /** true, if all molecules in this cell do not require any velocity relaxation */
198 : bool _ignoreThisCell;
199 :
200 : private:
201 : /**
202 : * @returns the inner lower left corner of the cell
203 : */
204 0 : tarch::la::Vector<dim, double> getInnerLowerLeftCorner() const { return IDXS.getGlobalMDDomainOffset() + 1.5 * IDXS.getCouplingCellSize(); }
205 : /**
206 : * @returns the inner upper right corner of the cell
207 : */
208 0 : tarch::la::Vector<dim, double> getInnerUpperRightCorner() const {
209 0 : return IDXS.getGlobalMDDomainOffset() + IDXS.getGlobalMDDomainSize() - 1.5 * IDXS.getCouplingCellSize();
210 : }
211 : /**
212 : * @returns the outer lower left corner of the cell
213 : */
214 0 : tarch::la::Vector<dim, double> getOuterLowerLeftCorner() const { return IDXS.getGlobalMDDomainOffset() + 0.5 * IDXS.getCouplingCellSize(); }
215 : /**
216 : * @returns the outer upper right corner of the cell
217 : */
218 0 : tarch::la::Vector<dim, double> getOuterUpperRightCorner() const {
219 0 : return IDXS.getGlobalMDDomainOffset() + IDXS.getGlobalMDDomainSize() - 0.5 * IDXS.getCouplingCellSize();
220 : }
221 :
222 : /** creates a coupling cell index for the second order interpolation.
223 : * @param position
224 : * @param normalisedPosition
225 : * @param secondCake
226 : * @param res
227 : */
228 0 : void createCouplingCellIndex4SecondOrderInterpolation(const tarch::la::Vector<dim, double>& position, tarch::la::Vector<dim, double>& normalisedPosition,
229 : bool& secondCake, I02* res) const {
230 : // compute lower left cell index and normalised position;
231 : // the normalised position is chosen such that the origin (0,0...,0)
232 : // coincides with the lower left... coupling cell's midpoint
233 :
234 0 : tarch::la::Vector<dim, double> cellMidpoint = _cellIdx.getCellMidPoint();
235 0 : auto cellSize = IDXS.getCouplingCellSize();
236 0 : tarch::la::Vector<dim, unsigned int> globalCellIndex{I01{_cellIdx}.get()};
237 0 : tarch::la::Vector<dim, unsigned int> lowerLeftCellIndex{I03{_cellIdx}.get()};
238 :
239 : // determine cell mid point and normalised position in this loop
240 0 : for (unsigned int d = 0; d < dim; d++) {
241 : // determine mid point of current coupling cell
242 :
243 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
244 : if ((position[d] < cellMidpoint[d] - 0.5 * cellSize[d]) || (position[d] > cellMidpoint[d] + 0.5 * cellSize[d])) {
245 : std::cout << "ERROR "
246 : "VelocityGradientRelaxationMapping::"
247 : "createCouplingCellIndex4SecondOrderInterpolation: "
248 : "Input position "
249 : << position;
250 : std::cout << " out of range of cell " << globalCellIndex << "!" << std::endl;
251 : exit(EXIT_FAILURE);
252 : }
253 : #endif
254 :
255 : // put lowerLeftCellIndex such that it's right below the position
256 0 : if (position[d] < cellMidpoint[d]) {
257 0 : lowerLeftCellIndex[d]--;
258 0 : globalCellIndex[d]--;
259 : }
260 : // if this is on the left/ lower boundary, but still to far up, reduce
261 : // height
262 0 : if (lowerLeftCellIndex[d] == 1) {
263 0 : lowerLeftCellIndex[d]--;
264 0 : globalCellIndex[d]--;
265 0 : } else if (lowerLeftCellIndex[d] == I03::numberCellsInDomain[d] - 2) {
266 0 : lowerLeftCellIndex[d]--;
267 0 : globalCellIndex[d]--;
268 : }
269 :
270 : // determine normalised position
271 0 : normalisedPosition[d] = (position[d] - (IDXS.getGlobalMDDomainOffset()[d] - 0.5 * cellSize[d]) - globalCellIndex[d] * cellSize[d]) / cellSize[d];
272 : }
273 0 : unsigned int lowerLeftIndex = I02{I03{tarch::la::Vector<3, int>{lowerLeftCellIndex}}}.get();
274 :
275 : if (dim == 2) {
276 : std::cout << "Not implemented correctly yet!" << std::endl;
277 : exit(EXIT_FAILURE);
278 : } else if (dim == 3) {
279 0 : const unsigned int localNumberCells01 = I03::numberCellsInDomain[0] * I03::numberCellsInDomain[1];
280 0 : const unsigned int localNumber2Cells0 = 2 * I03::numberCellsInDomain[0];
281 0 : const unsigned int localNumber2Cells01Plus2Cells0 = 2 * localNumberCells01 + 2 * I03::numberCellsInDomain[0];
282 :
283 0 : tarch::la::Vector<dim, double> normal;
284 0 : normal[0] = 1.0;
285 0 : normal[1] = 1.0;
286 0 : normal[2] = 0.0;
287 0 : tarch::la::Vector<dim, double> relPos(normalisedPosition);
288 0 : relPos[0] -= 2.0;
289 :
290 : // THE CAKE IS A LIE
291 :
292 : // first cake:
293 0 : if (tarch::la::dot(relPos, normal) < 0.0) {
294 0 : secondCake = false;
295 0 : res[0] = I02{lowerLeftIndex};
296 0 : res[1] = I02{lowerLeftIndex + 1};
297 0 : res[2] = I02{lowerLeftIndex + 2};
298 0 : res[3] = I02{lowerLeftIndex + I03::numberCellsInDomain[0]};
299 0 : res[4] = I02{lowerLeftIndex + localNumber2Cells0};
300 0 : res[5] = I02{lowerLeftIndex + localNumberCells01};
301 0 : res[6] = I02{lowerLeftIndex + localNumberCells01 + I03::numberCellsInDomain[0] + 1};
302 0 : res[7] = I02{lowerLeftIndex + 2 * localNumberCells01};
303 0 : res[8] = I02{lowerLeftIndex + 2 * localNumberCells01 + 2};
304 0 : res[9] = I02{lowerLeftIndex + localNumber2Cells01Plus2Cells0};
305 : // second cake
306 : } else {
307 0 : secondCake = true;
308 0 : res[0] = I02{lowerLeftIndex + 2};
309 0 : res[1] = I02{lowerLeftIndex + I03::numberCellsInDomain[0] + 2};
310 0 : res[2] = I02{lowerLeftIndex + localNumber2Cells0};
311 0 : res[3] = I02{lowerLeftIndex + localNumber2Cells0 + 1};
312 0 : res[4] = I02{lowerLeftIndex + localNumber2Cells0 + 2};
313 0 : res[5] = I02{lowerLeftIndex + localNumberCells01 + I03::numberCellsInDomain[0] + 1};
314 0 : res[6] = I02{lowerLeftIndex + localNumberCells01 + 2 * I03::numberCellsInDomain[0] + 2};
315 0 : res[7] = I02{lowerLeftIndex + 2 * localNumberCells01 + 2};
316 0 : res[8] = I02{lowerLeftIndex + localNumber2Cells01Plus2Cells0};
317 0 : res[9] = I02{lowerLeftIndex + localNumber2Cells01Plus2Cells0 + 2};
318 : }
319 : } else {
320 : std::cout << "ERROR createCouplingCellIndex4Interpolation(position): "
321 : "only 2D/3D supported!"
322 : << std::endl;
323 : exit(EXIT_FAILURE);
324 : }
325 0 : }
326 :
327 : /** carries out second order interpolation of the velocity value
328 : *(microscopicmomentum-buffer) at position 'position'
329 : * @param cellIndex
330 : * @param normalisedPosition
331 : * @param secondCake
332 : * @param newVelocity
333 : */
334 0 : void interpolateVelocitySecondOrder(I02* cellIndex, const tarch::la::Vector<dim, double>& normalisedPosition, const bool& secondCake,
335 : tarch::la::Vector<dim, double>& newVelocity) const {
336 : if (dim == 2) {
337 : // TODO
338 : } else if (dim == 3) {
339 : // compute interpolation coefficients
340 0 : tarch::la::Vector<dim, double> coefficients[10];
341 0 : tarch::la::Vector<dim, double> velocity[10];
342 :
343 : // extract velocities from cells
344 0 : for (int i = 0; i < 10; i++) {
345 0 : velocity[i] = _couplingCells[cellIndex[i].get()].getMicroscopicMomentum();
346 : }
347 :
348 0 : if (secondCake) {
349 0 : const tarch::la::Vector<dim, double> v78 = velocity[7] + velocity[8];
350 0 : coefficients[0] =
351 0 : 2.0 * (velocity[0] + velocity[2] + velocity[9]) + 4.0 * (velocity[5] - velocity[1] - velocity[6] - velocity[3]) + 5.0 * velocity[4] - v78;
352 0 : coefficients[1] =
353 0 : 0.5 * (v78 - velocity[0]) + 2.0 * (velocity[1] - velocity[2] - velocity[5] + velocity[6]) + 4.0 * velocity[3] - 3.5 * velocity[4] - velocity[9];
354 0 : coefficients[2] =
355 0 : 2.0 * (velocity[3] - velocity[0] - velocity[5] + velocity[6]) + 4.0 * velocity[1] - 3.5 * velocity[4] + 0.5 * (v78 - velocity[2]) - velocity[9];
356 0 : coefficients[3] = 0.5 * (v78 - velocity[0] - velocity[2] - velocity[4]) + 2.0 * velocity[6] - 1.5 * velocity[9];
357 0 : coefficients[4] = 0.5 * (velocity[2] + velocity[4]) - velocity[3];
358 0 : coefficients[5] = 0.5 * (velocity[0] + velocity[4]) - velocity[1];
359 0 : coefficients[6] = 0.5 * (velocity[4] + velocity[9]) - velocity[6];
360 0 : coefficients[7] =
361 0 : 0.25 * (velocity[0] + velocity[2] - v78) - velocity[1] - velocity[3] + 1.5 * velocity[4] + velocity[5] - velocity[6] + 0.5 * velocity[9];
362 0 : coefficients[8] = 0.25 * (velocity[2] - velocity[4] - velocity[8] + velocity[9]);
363 0 : coefficients[9] = 0.25 * (velocity[0] - velocity[4] - velocity[7] + velocity[9]);
364 : } else {
365 0 : const tarch::la::Vector<dim, double> minus15V0 = -1.5 * velocity[0];
366 0 : coefficients[0] = velocity[0];
367 0 : coefficients[1] = minus15V0 + 2.0 * velocity[1] - 0.5 * velocity[2];
368 0 : coefficients[2] = minus15V0 + 2.0 * velocity[3] - 0.5 * velocity[4];
369 0 : coefficients[3] = minus15V0 + 2.0 * velocity[5] - 0.5 * velocity[7];
370 0 : coefficients[4] = 0.5 * (velocity[0] + velocity[2]) - velocity[1];
371 0 : coefficients[5] = 0.5 * (velocity[0] + velocity[4]) - velocity[3];
372 0 : coefficients[6] = 0.5 * (velocity[0] + velocity[7]) - velocity[5];
373 0 : coefficients[7] = 1.5 * velocity[0] + velocity[6] - velocity[1] - velocity[3] - velocity[5] + 0.5 * velocity[7] +
374 0 : 0.25 * (velocity[2] + velocity[4] - velocity[8] - velocity[9]);
375 0 : coefficients[8] = 0.25 * (velocity[0] - velocity[2] - velocity[7] + velocity[8]);
376 0 : coefficients[9] = 0.25 * (velocity[0] - velocity[4] - velocity[7] + velocity[9]);
377 : }
378 :
379 : // evaluate polynomial
380 0 : const tarch::la::Vector<dim, double> contribution0 =
381 : normalisedPosition[0] *
382 0 : (coefficients[1] + coefficients[4] * normalisedPosition[0] + coefficients[7] * normalisedPosition[1] + coefficients[8] * normalisedPosition[2]);
383 0 : const tarch::la::Vector<dim, double> contribution1 =
384 0 : normalisedPosition[1] * (coefficients[2] + coefficients[5] * normalisedPosition[1] + coefficients[9] * normalisedPosition[2]);
385 0 : const tarch::la::Vector<dim, double> contribution2 = normalisedPosition[2] * (coefficients[3] + coefficients[6] * normalisedPosition[2]);
386 :
387 0 : newVelocity = coefficients[0] + contribution0 + contribution1 + contribution2;
388 : }
389 0 : }
390 : };
391 :
392 : /** This is the same as the class
393 : *coupling::cellmappings::VelocityGradientRelaxationMapping, but relaxes only
394 : *those molecules which are located in the top boundary strip. derived from the
395 : *class coupling::cellmappings::VelocityGradientRelaxationMapping.
396 : * @brief This is the same as the class
397 : *coupling::cellmappings::VelocityGradientRelaxationMapping, but relaxes only
398 : *those molecules which are located in the top boundary strip.
399 : * @tparam LinkedCell cell type
400 : * @tparam dim Number of dimensions; it can be 1, 2 or 3
401 : * @author Philipp Neumann
402 : * @note ONLY SUPPORTS 3D AT THE MOMENT!
403 : * \todo Philipp please take a look on this class
404 : */
405 : template <class LinkedCell, unsigned int dim>
406 0 : class coupling::cellmappings::VelocityGradientRelaxationTopOnlyMapping : public coupling::cellmappings::VelocityGradientRelaxationMapping<LinkedCell, dim> {
407 : public:
408 : /** obtains relaxation factor and current velocity (the velocity towards which
409 : *the relaxation shall be done is later obtained from the
410 : *microscopicmomentum-buffers).
411 : * @param velocityRelaxationFactor
412 : * @param currentVelocity
413 : * @param cellIndex
414 : * @param mdSolverInterface
415 : * @param couplingCells
416 : * @note the method ignoreThisCell() of the class
417 : *coupling::cellmappings::VelocityGradientRelaxationMapping is replaceed.
418 : *Since this function is called in the constructor of the based class, we
419 : *cannot overwrite it; hence, we solve it by overwriting the respective
420 : *variable from within this constructor (called after base object is
421 : *constructed)
422 : */
423 0 : VelocityGradientRelaxationTopOnlyMapping(const double& velocityRelaxationFactor, const tarch::la::Vector<dim, double>& currentVelocity, const I02& cellIndex,
424 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface,
425 : const coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>* const couplingCells)
426 : : coupling::cellmappings::VelocityGradientRelaxationMapping<LinkedCell, dim>(velocityRelaxationFactor, currentVelocity, cellIndex, mdSolverInterface,
427 0 : couplingCells) {
428 :
429 : // the following snippet is basically a replacement of the method
430 : // ignoreThisCell(). Since this function is called in the constructor of the
431 : // based class, we cannot overwrite it; hence, we solve it by overwriting
432 : // the respective variable from within this constructor (called after base
433 : // object is constructed)
434 0 : coupling::cellmappings::VelocityGradientRelaxationMapping<LinkedCell, dim>::_ignoreThisCell =
435 0 : (I01{cellIndex}.get()[dim - 1] < (int)(I01::numberCellsInDomain[dim - 1]) - 3);
436 0 : }
437 :
438 : protected:
439 : /** returns true, if the position 'position' is within the respective boundary
440 : *strip that is spanned by the two outermost coupling cell midpoints
441 : * @param position
442 : * @return true, if the position 'position' is within the respective boundary
443 : *strip that is spanned by the two outermost coupling cell midpoints
444 : */
445 0 : virtual bool relaxMolecule(const tarch::la::Vector<dim, double>& position) const {
446 : // check if molecule is in the respective boundary strip (only upper part is
447 : // considered)
448 0 : return (position[dim - 1] > coupling::cellmappings::VelocityGradientRelaxationMapping<LinkedCell, dim>::_innerUpperRight[dim - 1]) &&
449 0 : (position[dim - 1] < coupling::cellmappings::VelocityGradientRelaxationMapping<LinkedCell, dim>::_outerUpperRight[dim - 1]);
450 : }
451 : };
452 :
453 : #endif // _MOLECULARDYNAMICS_COUPLING_CELLMAPPINGS_VELOCITYGRADIENTRELAXATIONMAPPING_H_
|
