deltaBoundary.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2007-2019 PCOpt/NTUA
9  Copyright (C) 2013-2019 FOSS GP
10  Copyright (C) 2019 OpenCFD Ltd.
11 -------------------------------------------------------------------------------
12 License
13  This file is part of OpenFOAM.
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "deltaBoundary.H"
31 #include "fvMesh.H"
32 #include "surfaceFields.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
40 
41 tensor deltaBoundary::tensorCrossVector(const tensor& T, const vector& v)
42 {
43  // The correct approach when T is not a diagonal tensor
44  tensor res(Zero);
45  vector vec1(T.xx(), T.yx(), T.zx());
46  vector res1(vec1 ^ v);
47  res.xx() = res1.x(); res.yx() = res1.y(); res.zx() = res1.z();
48 
49  vector vec2(T.xy(), T.yy(), T.zy());
50  vector res2(vec2 ^ v);
51  res.xy() = res2.x(); res.yy() = res2.y(); res.zy() = res2.z();
52 
53  vector vec3(T.xz(), T.yz(), T.zz());
54  vector res3(vec3 ^ v);
55  res.xz() = res3.x(); res.yz() = res3.y(); res.zz() = res3.z();
56 
57  return res;
58 }
59 
60 
61 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
62 
63 deltaBoundary::deltaBoundary(const fvMesh& mesh)
64 :
65  mesh_(mesh)
66 {}
67 
68 
69 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
70 
72 (
73  const pointField& p,
74  const pointField& p_d
75 )
76 {
77  vector fCtrs_d(Zero);
78  vector fAreas_d(Zero);
79  vector unitVector_d(Zero);
80 
81  // Container field to return results
82  vectorField deltaVecs(3, Zero);
83 
84  label nPoints = p.size();
85 
86  // If the face is a triangle, do a direct calculation for efficiency
87  // and to avoid round-off error-related problems
88  if (nPoints == 3)
89  {
90  //fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
91  vector fAreas = 0.5*((p[1] - p[0])^(p[2] - p[0]));
92 
93  fCtrs_d = (1.0/3.0)*(p_d[0] + p_d[1] + p_d[2]);
94  fAreas_d =
95  0.5*((p_d[1] - p_d[0])^(p[2] - p[0]))
96  + 0.5*((p[1] - p[0])^(p_d[2] - p_d[0]));
97  scalar ds = mag(fAreas);
98  unitVector_d = fAreas_d/ds - (fAreas*(fAreas&fAreas_d))/ds/ds/ds;
99 
100  deltaVecs[0] = fCtrs_d;
101  deltaVecs[1] = fAreas_d;
102  deltaVecs[2] = unitVector_d;
103  }
104  else
105  {
106  vector sumN(Zero);
107  vector sumN_d(Zero);
108  scalar sumA(0);
109  scalar sumA_d(0);
110  vector sumAc(Zero);
111  vector sumAc_d(Zero);
112 
113  point fCentre = p[0];
114  point fCentre_d = p_d[0];
115  for (label pi = 1; pi < nPoints; pi++)
116  {
117  fCentre += p[pi];
118  fCentre_d += p_d[pi];
119  }
120 
121  fCentre /= nPoints;
122  fCentre_d /= nPoints;
123 
124  for (label pi = 0; pi < nPoints; pi++)
125  {
126  const point& nextPoint = p[(pi + 1) % nPoints];
127  const point& nextPoint_d = p_d[(pi + 1) % nPoints];
128 
129  vector c = p[pi] + nextPoint + fCentre;
130  vector c_d = p_d[pi] + nextPoint_d + fCentre_d;
131 
132  vector n = (nextPoint - p[pi])^(fCentre - p[pi]);
133  vector n_d =
134  ((nextPoint_d - p_d[pi])^(fCentre - p[pi]))
135  + ((nextPoint - p[pi])^(fCentre_d - p_d[pi]));
136 
137  scalar a = mag(n);
138  if (a < ROOTVSMALL)
139  {
140  // This shouldn't happen in general.
141  // Manually zero contribution from zero area face for now
143  << "Zero area face sub triangle found " << endl
144  << p[pi] << " " << nextPoint << " " << fCentre << endl
145  << "Neglecting contributions of this element " << endl;
146  }
147  else
148  {
149  scalar a_d = (n&n_d)/mag(n);
150 
151  sumN += n;
152  sumN_d += n_d;
153 
154  sumA += a;
155  sumA_d += a_d;
156 
157  sumAc += a*c;
158  sumAc_d += a_d*c + a*c_d;
159  }
160  }
161 
162  // fCtrs[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
163  vector fAreas = 0.5*sumN;
164  fCtrs_d = (1.0/3.0)*(sumAc_d*sumA - sumAc*sumA_d)/sumA/sumA;
165  fAreas_d = 0.5*sumN_d;
166  scalar ds = mag(fAreas);
167  unitVector_d = fAreas_d/ds - (fAreas*(fAreas&fAreas_d))/ds/ds/ds;
168 
169  deltaVecs[0] = fCtrs_d;
170  deltaVecs[1] = fAreas_d;
171  deltaVecs[2] = unitVector_d;
172  }
173 
174  return deltaVecs;
175 }
176 
177 
179 (
180  const pointField& p,
181  const tensorField& p_d
182 )
183 {
184  label nPoints = p.size();
185  tensor fCtrs_d(Zero);
186  tensor fAreas_d(Zero);
187  tensor unitVector_d(Zero);
188 
189  // Container field to return results
190  tensorField deltaVecs(3, Zero);
191 
192  // If the face is a triangle, do a direct calculation for efficiency
193  // and to avoid round-off error-related problems
194  if (nPoints == 3)
195  {
196  vector fAreas = 0.5*((p[1] - p[0])^(p[2] - p[0]));
197 
198  fCtrs_d = (1.0/3.0)*(p_d[0] + p_d[1] + p_d[2]);
199  fAreas_d =
200  0.5*tensorCrossVector(p_d[1] - p_d[0], p[2] - p[0])
201  //minus sign since it is vector ^ tensor
202  - 0.5*tensorCrossVector(p_d[2] - p_d[0], p[1] - p[0]);
203  scalar ds = mag(fAreas);
204  unitVector_d = fAreas_d/ds - (fAreas*(fAreas&fAreas_d))/ds/ds/ds;
205 
206  deltaVecs[0] = fCtrs_d;
207  deltaVecs[1] = fAreas_d;
208  deltaVecs[2] = unitVector_d;
209  }
210  else
211  {
212  vector sumN(Zero);
213  tensor sumN_d(Zero);
214  scalar sumA(0);
215  vector sumA_d(Zero);
216  vector sumAc(Zero);
217  tensor sumAc_d(Zero);
218 
219  point fCentre = p[0];
220  tensor fCentre_d = p_d[0];
221  for (label pi = 1; pi < nPoints; pi++)
222  {
223  fCentre += p[pi];
224  fCentre_d += p_d[pi];
225  }
226 
227  fCentre /= nPoints;
228  fCentre_d /= nPoints;
229 
230  for (label pi = 0; pi < nPoints; pi++)
231  {
232  const point& nextPoint = p[(pi + 1) % nPoints];
233  const tensor& nextPoint_d = p_d[(pi + 1) % nPoints];
234 
235  vector c = p[pi] + nextPoint + fCentre;
236  tensor c_d = p_d[pi] + nextPoint_d + fCentre_d;
237 
238  vector n = (nextPoint - p[pi])^(fCentre - p[pi]);
239  tensor n_d =
240  tensorCrossVector(nextPoint_d - p_d[pi], fCentre - p[pi])
241  //minus sign since it is vector ^ tensor
242  - tensorCrossVector(fCentre_d - p_d[pi], nextPoint - p[pi]);
243 
244  scalar a = mag(n);
245  if (a < ROOTVSMALL)
246  {
247  // This shouldn't happen in general.
248  // Manually zero contribution from zero area face for now
250  << "Zero area face sub triangle found " << nl
251  << p[pi] << " " << nextPoint << " " << fCentre << nl
252  << "Neglecting contributions of this element " << endl;
253  }
254  else
255  {
256  vector a_d = (n & n_d)/a;
257 
258  sumN += n;
259  sumN_d += n_d;
260 
261  sumA += a;
262  sumA_d += a_d;
263 
264  sumAc += a*c;
265  // c*a_d since we need to get the correct outer product
266  sumAc_d += (c*a_d) + a*c_d;
267  }
268  }
269 
270  vector fAreas = 0.5*sumN;
271  fCtrs_d = (1.0/3.0)*(sumAc_d/sumA - (sumAc*sumA_d)/sqr(sumA));
272  fAreas_d = 0.5*sumN_d;
273  scalar ds = mag(fAreas);
274  unitVector_d = fAreas_d/ds - (fAreas*(fAreas&fAreas_d))/ds/ds/ds;
275 
276  deltaVecs[0] = fCtrs_d;
277  deltaVecs[1] = fAreas_d;
278  deltaVecs[2] = unitVector_d;
279  }
280 
281  return deltaVecs;
282 }
283 
284 
286 {
288  const labelList& pointCellsI(pointCells[pointI]);
289  const pointField& points(mesh_.points());
290  auto tC_d = tmp<tensorField>::New(pointCellsI.size(), Zero);
291  auto& C_d = tC_d.ref();
292 
293  const labelList& pointFaces(mesh_.pointFaces()[pointI]);
294  tensorField Cf_d(pointFaces.size(), Zero);
295  tensorField Sf_d(pointFaces.size(), Zero);
296 
297  forAll(pointFaces, pfI)
298  {
299  const label pointFaceI = pointFaces[pfI];
300  const face& faceI = mesh_.faces()[pointFaceI];
301  tensorField p_d(faceI.size(), Zero);
302  forAll(faceI, pI)
303  {
304  if (faceI[pI] == pointI)
305  {
306  p_d[pI] = tensor::I;
307  break;
308  }
309  }
310 
311  pointField facePoints(faceI.points(points));
312 
313  // Compute changes in the face
314  tensorField dFace(makeFaceCentresAndAreas_d(facePoints, p_d));
315  Cf_d[pfI] = dFace[0];
316  Sf_d[pfI] = dFace[1];
317  }
318 
319  // Face variations have now been computed. Now, compute cell contributions
320  forAll(pointCellsI, pcI)
321  {
322  const label pointCellI = pointCellsI[pcI];
323  const cell& cellI(mesh_.cells()[pointCellI]);
324  vectorField fAreas(cellI.size(), Zero);
325  vectorField fCtrs(cellI.size(), Zero);
326  tensorField fAreas_d(cellI.size(), Zero);
327  tensorField fCtrs_d(cellI.size(), Zero);
328  forAll(cellI, fI)
329  {
330  const label globalFaceI = cellI[fI];
331 
332  // Assign values to faceAreas and faceCtrs
333  if (globalFaceI < mesh_.nInternalFaces())
334  {
335  fAreas[fI] = mesh_.Sf()[globalFaceI];
336  fCtrs[fI] = mesh_.Cf()[globalFaceI];
337  }
338  else
339  {
340  const label whichPatch =
341  mesh_.boundaryMesh().whichPatch(globalFaceI);
342  const fvPatch& patch = mesh_.boundary()[whichPatch];
343  const label patchStart = patch.patch().start();
344  const label localFace = globalFaceI - patchStart;
345  fAreas[fI] = patch.Sf()[localFace];
346  fCtrs[fI] = patch.Cf()[localFace];
347  }
348 
349  // Assign values to differentiated face areas and centres
350  forAll(pointFaces, pfI)
351  {
352  if (pointFaces[pfI] == globalFaceI)
353  {
354  fAreas_d[fI] = Sf_d[pfI];
355  fCtrs_d[fI] = Cf_d[pfI];
356  }
357  }
358  }
359  C_d[pcI] = makeCellCentres_d(fAreas, fCtrs, fAreas_d, fCtrs_d);
360  }
361 
362  return tC_d;
363 }
364 
365 
366 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
367 
368 } // End namespace Foam
369 
370 // ************************************************************************* //
Foam::Tensor< scalar >
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::Tensor::I
static const Tensor I
Definition: Tensor.H:85
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::primitiveMesh::pointFaces
const labelListList & pointFaces() const
Definition: primitiveMeshPointFaces.C:34
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:138
Foam::deltaBoundary::mesh_
const fvMesh & mesh_
Reference to the mesh.
Definition: deltaBoundary.H:65
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
surfaceFields.H
Foam::surfaceFields.
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::deltaBoundary::makeCellCentres_d
pT makeCellCentres_d(const vectorField &fAreas, const vectorField &fCtrs, const Field< pT > &fAreas_d, const Field< pT > &fCtrs_d)
Definition: deltaBoundaryTemplates.C:41
Foam::deltaBoundary::cellCenters_d
tmp< tensorField > cellCenters_d(const label pointI)
Definition: deltaBoundary.C:285
Foam::fvPatch
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:65
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
Foam::polyBoundaryMesh::whichPatch
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Definition: polyBoundaryMesh.C:812
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::deltaBoundary::makeFaceCentresAndAreas_d
vectorField makeFaceCentresAndAreas_d(const pointField &p, const pointField &p_d)
Definition: deltaBoundary.C:72
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
fvMesh.H
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
deltaBoundary.H
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:685
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1094
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
Foam::Vector< scalar >
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::List< labelList >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::face::points
pointField points(const UList< point > &pts) const
Return the points corresponding to this face.
Definition: faceI.H:87
Foam::primitiveMesh::pointCells
const labelListList & pointCells() const
Definition: primitiveMeshPointCells.C:110
Foam::tmp::New
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::pointCells
Smooth ATC in cells having a point to a set of patches supplied by type.
Definition: pointCells.H:56
Foam::fvMesh::Cf
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
Definition: fvMeshGeometry.C:352
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
Foam::fvMesh::Sf
const surfaceVectorField & Sf() const
Return cell face area vectors.
Definition: fvMeshGeometry.C:319