calculateMeshToMesh0Weights.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "meshToMesh0.H"
27 #include "tetOverlapVolume.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 void Foam::meshToMesh0::calculateInverseDistanceWeights() const
32 {
33  if (debug)
34  {
36  << "Calculating inverse distance weighting factors" << endl;
37  }
38 
39  if (inverseDistanceWeightsPtr_)
40  {
42  << "weighting factors already calculated"
43  << exit(FatalError);
44  }
45 
46  //- Initialise overlap volume to zero
47  V_ = 0.0;
48 
49  inverseDistanceWeightsPtr_ = new scalarListList(toMesh_.nCells());
50  scalarListList& invDistCoeffs = *inverseDistanceWeightsPtr_;
51 
52  // get reference to source mesh data
53  const labelListList& cc = fromMesh_.cellCells();
54  const vectorField& centreFrom = fromMesh_.C();
55  const vectorField& centreTo = toMesh_.C();
56 
57  forAll(cellAddressing_, celli)
58  {
59  if (cellAddressing_[celli] != -1)
60  {
61  const vector& target = centreTo[celli];
62  scalar m = mag(target - centreFrom[cellAddressing_[celli]]);
63 
64  const labelList& neighbours = cc[cellAddressing_[celli]];
65 
66  // if the nearest cell is a boundary cell or there is a direct hit,
67  // pick up the value
68  label directCelli = -1;
69  if (m < directHitTol || neighbours.empty())
70  {
71  directCelli = celli;
72  }
73  else
74  {
75  forAll(neighbours, ni)
76  {
77  scalar nm = mag(target - centreFrom[neighbours[ni]]);
78  if (nm < directHitTol)
79  {
80  directCelli = neighbours[ni];
81  break;
82  }
83  }
84  }
85 
86 
87  if (directCelli != -1)
88  {
89  // Direct hit
90  invDistCoeffs[directCelli].setSize(1);
91  invDistCoeffs[directCelli][0] = 1.0;
92  V_ += fromMesh_.V()[cellAddressing_[directCelli]];
93  }
94  else
95  {
96  invDistCoeffs[celli].setSize(neighbours.size() + 1);
97 
98  // The first coefficient corresponds to the centre cell.
99  // The rest is ordered in the same way as the cellCells list.
100  scalar invDist = 1.0/m;
101  invDistCoeffs[celli][0] = invDist;
102  scalar sumInvDist = invDist;
103 
104  // now add the neighbours
105  forAll(neighbours, ni)
106  {
107  invDist = 1.0/mag(target - centreFrom[neighbours[ni]]);
108  invDistCoeffs[celli][ni + 1] = invDist;
109  sumInvDist += invDist;
110  }
111 
112  // divide by the total inverse-distance
113  forAll(invDistCoeffs[celli], i)
114  {
115  invDistCoeffs[celli][i] /= sumInvDist;
116  }
117 
118 
119  V_ +=
120  invDistCoeffs[celli][0]
121  *fromMesh_.V()[cellAddressing_[celli]];
122  for (label i = 1; i < invDistCoeffs[celli].size(); i++)
123  {
124  V_ +=
125  invDistCoeffs[celli][i]*fromMesh_.V()[neighbours[i-1]];
126  }
127  }
128  }
129  }
130 }
131 
132 
133 void Foam::meshToMesh0::calculateInverseVolumeWeights() const
134 {
135  if (debug)
136  {
138  << "Calculating inverse volume weighting factors" << endl;
139  }
140 
141  if (inverseVolumeWeightsPtr_)
142  {
144  << "weighting factors already calculated"
145  << exit(FatalError);
146  }
147 
148  //- Initialise overlap volume to zero
149  V_ = 0.0;
150 
151  inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
152  scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
153 
154  const labelListList& cellToCell = cellToCellAddressing();
155 
156  tetOverlapVolume overlapEngine;
157 
158  forAll(cellToCell, celli)
159  {
160  const labelList& overlapCells = cellToCell[celli];
161 
162  if (overlapCells.size() > 0)
163  {
164  invVolCoeffs[celli].setSize(overlapCells.size());
165 
166  forAll(overlapCells, j)
167  {
168  label cellFrom = overlapCells[j];
169  treeBoundBox bbFromMesh
170  (
171  pointField
172  (
173  fromMesh_.points(),
174  fromMesh_.cellPoints()[cellFrom]
175  )
176  );
177 
178  scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp
179  (
180  toMesh_,
181  celli,
182 
183  fromMesh_,
184  cellFrom,
185  bbFromMesh
186  );
187  invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
188 
189  V_ += v;
190  }
191  }
192  }
193 }
194 
195 
196 void Foam::meshToMesh0::calculateCellToCellAddressing() const
197 {
198  if (debug)
199  {
201  << "Calculating cell to cell addressing" << endl;
202  }
203 
204  if (cellToCellAddressingPtr_)
205  {
207  << "addressing already calculated"
208  << exit(FatalError);
209  }
210 
211  //- Initialise overlap volume to zero
212  V_ = 0.0;
213 
214  tetOverlapVolume overlapEngine;
215 
216  cellToCellAddressingPtr_ = new labelListList(toMesh_.nCells());
217  labelListList& cellToCell = *cellToCellAddressingPtr_;
218 
219 
220  forAll(cellToCell, iTo)
221  {
222  const labelList overLapCells =
223  overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo);
224  if (overLapCells.size() > 0)
225  {
226  cellToCell[iTo].setSize(overLapCells.size());
227  forAll(overLapCells, j)
228  {
229  cellToCell[iTo][j] = overLapCells[j];
230  V_ += fromMesh_.V()[overLapCells[j]];
231  }
232  }
233  }
234 }
235 
236 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
237 
238 const Foam::scalarListList& Foam::meshToMesh0::inverseDistanceWeights() const
239 {
240  if (!inverseDistanceWeightsPtr_)
241  {
242  calculateInverseDistanceWeights();
243  }
244 
245  return *inverseDistanceWeightsPtr_;
246 }
247 
248 
249 const Foam::scalarListList& Foam::meshToMesh0::inverseVolumeWeights() const
250 {
251  if (!inverseVolumeWeightsPtr_)
252  {
253  calculateInverseVolumeWeights();
254  }
255 
256  return *inverseVolumeWeightsPtr_;
257 }
258 
259 
260 const Foam::labelListList& Foam::meshToMesh0::cellToCellAddressing() const
261 {
262  if (!cellToCellAddressingPtr_)
263  {
264  calculateCellToCellAddressing();
265  }
266 
267  return *cellToCellAddressingPtr_;
268 }
269 
270 
271 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< scalarList > scalarListList
Definition: scalarList.H:51
List< label > labelList
A List of labels.
Definition: labelList.H:56
void setSize(const label)
Reset size of List.
Definition: List.C:281
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const volVectorField & C() const
Return cell centres as volVectorField.
const labelListList & cellCells() const
const labelListList & cellPoints() const
#define InfoInFunction
Report an information message using Foam::Info.