nonConformalCyclicLagrangianPatch.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) 2025-2026 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 
27 #include "LagrangianFields.H"
28 #include "meshSearch.H"
29 #include "RemoteData.H"
30 #include "tracking.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  template<class Type>
38  struct maxFirstOp
39  {
40  const Type& operator()(const Type& a, const Type& b) const
41  {
42  return a.first() > b.first() ? a : b;
43  }
44  };
45 }
46 
47 
48 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
49 
50 namespace Foam
51 {
53 
55  (
58  polyPatch
59  );
60 }
61 
62 
63 
64 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
65 
67 (
68  const polyPatch& poly,
69  const LagrangianBoundaryMesh& boundaryMesh
70 )
71 :
72  LagrangianPatch(poly, boundaryMesh),
73  nonConformalCyclicPoly_(refCast<const nonConformalCyclicPolyPatch>(poly)),
74  isNbrPatchMesh_(false),
75  nPositionalErrors_(0),
76  maxPositionalErrorSqr_(-vGreat),
77  maxPositionalErrorReceivePosition_(point::uniform(NaN))
78 {}
79 
80 
81 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
82 
84 {}
85 
86 
87 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
88 
91 {
92  return
93  boundaryMesh()
94  [
95  isNbrPatchMesh_
96  ? nonConformalCyclicPoly_.nbrPatchIndex()
97  : poly().index()
98  ].LagrangianPatch::mesh();
99 }
100 
101 
103 (
107 ) const
108 {
109  const LagrangianSubMesh& patchMesh = this->mesh();
110 
111  const meshSearch& searchEngine = meshSearch::New(mesh.poly());
112 
113  // Sub-set the geometry and topology of the elements
114  SubField<barycentric> patchCoordinates = patchMesh.sub(mesh.coordinates());
115  SubField<label> patchCelli = patchMesh.sub(mesh.celli());
116  SubField<label> patchFacei = patchMesh.sub(mesh.facei());
117  SubField<label> patchFaceTrii = patchMesh.sub(mesh.faceTrii());
118 
119  // Sub-set the receiving information
120  SubField<label> receivePatchFace =
121  patchMesh.sub(mesh.receivePatchFacePtr_());
122  SubField<vector> receivePosition =
123  patchMesh.sub(mesh.receivePositionPtr_());
124 
125  // Get a reference to the receiving original patch
126  const polyPatch& receivePp =
127  nonConformalCyclicPoly_.nbrPatch().origPatch();
128 
129  // Search for the elements on the receiving side
130  forAll(patchMesh, i)
131  {
132  patchCelli[i] =
133  mesh.poly().faceOwner()[receivePatchFace[i] + receivePp.start()];
134 
135  if
136  (
138  (
139  searchEngine,
140  receivePosition[i],
141  patchCoordinates[i],
142  patchCelli[i],
143  patchFacei[i],
144  patchFaceTrii[i],
145  fraction[i + patchMesh.start()]
146  )
147  )
148  {
149  // The search hit a boundary. Compute the positional error.
150  const scalar positionalErrorSqr =
151  magSqr
152  (
153  receivePosition[i]
155  (
156  mesh.poly(),
157  patchCoordinates[i],
158  patchCelli[i],
159  patchFacei[i],
160  patchFaceTrii[i],
161  fraction[i + patchMesh.start()]
162  )
163  );
164 
165  // If the positional error is significant, relative to the size of
166  // the receiving face, then register as a failure
167  if
168  (
169  sqr(positionalErrorSqr)
170  > sqr(small)*magSqr(receivePp.faceAreas()[receivePatchFace[i]])
171  )
172  {
173  nPositionalErrors_ ++;
174 
175  if (positionalErrorSqr > maxPositionalErrorSqr_)
176  {
177  maxPositionalErrorSqr_ = positionalErrorSqr;
178  maxPositionalErrorReceivePosition_ = receivePosition[i];
179  }
180  }
181  }
182  }
183 
184  // Set the elements to be in the adjacent cell
185  patchMesh.sub(mesh.states()) = LagrangianState::inCell;
186 
187  // The elements are now on the neighbour patch
188  isNbrPatchMesh_ = true;
189 
190  // Invalidate the now used receiving information
191  receivePatchFace = -1;
192  receivePosition = point::nan;
193 }
194 
195 
197 {
199 
200  // The elements are back on this patch
201  isNbrPatchMesh_ = false;
202 
203  // Check for positional errors and report if any are found
204  const label nPositionalErrors =
205  returnReduce(nPositionalErrors_, sumOp<label>());
206 
207  if (nPositionalErrors != 0)
208  {
209  const label nTransfers =
210  returnReduce(this->mesh().size(), sumOp<label>());
211 
212  const Tuple2<scalar, vector> maxPositionalErrorInfo =
214  (
216  (
217  maxPositionalErrorSqr_,
218  maxPositionalErrorReceivePosition_
219  ),
221  );
222 
223  maxPositionalErrorSqr_ = maxPositionalErrorInfo.first();
224  maxPositionalErrorReceivePosition_ = maxPositionalErrorInfo.second();
225 
227  << nPositionalErrors << "/" << nTransfers << " elements "
228  << "transferring to patch "
229  << nonConformalCyclicPoly_.nbrPatch().name()
230  << " were not accurately located. The largest positional error "
231  << "was " << sqrt(maxPositionalErrorSqr_)
232  << " from " << maxPositionalErrorReceivePosition_ << "."
233  << nl;
234  }
235 
236  nPositionalErrors_ = 0;
237  maxPositionalErrorSqr_ = -vGreat;
238  maxPositionalErrorReceivePosition_ = point::uniform(NaN);
239 }
240 
241 
242 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:449
Macros for easy insertion into run-time selection tables.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Boundary part of a Lagrangian mesh. Just a list of Lagrangian patches with some added convenience fun...
Class containing Lagrangian geometry and topology.
Base class for Lagrangian patches.
virtual void partition() const
Update for mesh changes.
Mesh that relates to a sub-section of a Lagrangian mesh. This is used to construct fields that relate...
word sub(const word &fieldName) const
Return the name of a field corresponding to this sub-mesh.
label start() const
Return start.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
Pre-declare related SubField type.
Definition: SubField.H:63
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:66
const Type2 & second() const
Return second.
Definition: Tuple2.H:131
const Type1 & first() const
Return first.
Definition: Tuple2.H:119
static const Form nan
Definition: VectorSpace.H:124
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
const polyMesh & poly() const
Return reference to polyMesh.
Definition: fvMesh.H:456
Mesh object that implements searches within the local cells and faces.
Definition: meshSearch.H:59
static const meshSearch & New(const polyMesh &mesh, const pointInCellShapes=pointInCellShapes::tets)
Lookup or construct from mesh and cell decomposition option.
Definition: meshSearch.C:61
virtual void evaluate(PstreamBuffers &, LagrangianMesh &, const LagrangianInternalScalarDynamicField &fraction) const
Evaluate changes in elements that have tracked to this patch.
nonConformalCyclicLagrangianPatch(const polyPatch &, const LagrangianBoundaryMesh &)
Construct from a patch and a boundary mesh.
virtual void partition() const
Update following partitioning of the mesh.
virtual const LagrangianSubMesh & mesh() const
Return the sub-mesh associated with this patch.
Non-conformal cyclic poly patch. As nonConformalCoupledPolyPatch, but the neighbouring patch is local...
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1321
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:71
const vectorField::subField faceAreas() const
Return face areas.
Definition: polyPatch.C:233
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:277
Foam::fvMesh mesh(Foam::IOobject(regionName, runTime.name(), runTime, Foam::IOobject::MUST_READ), false)
volScalarField & b
Definition: createFields.H:27
#define WarningInFunction
Report a warning using Foam::Warning.
point position(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label facei, const label faceTrii, const scalar stepFraction)
Return the position given the coordinates and tet topology.
Definition: trackingI.H:224
bool locate(const polyMesh &mesh, const point &position, barycentric &coordinates, label &celli, label &facei, label &faceTrii, const scalar stepFraction, const string &debugPrefix=NullObjectRef< string >())
Initialise the location at the given position. Returns whether or not a.
Definition: tracking.C:1592
const unitSet fraction
Namespace for OpenFOAM.
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
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:141
tmp< DimensionedField< typename outerProduct< Type, Type >::type, GeoMesh, Field >> sqr(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void sqrt(LagrangianPatchField< scalar > &f, const LagrangianPatchField< scalar > &f1)
tmp< DimensionedField< scalar, GeoMesh, Field > > magSqr(const DimensionedField< Type, GeoMesh, PrimitiveField > &df)
defineTypeNameAndDebug(atmosphericBoundaryLayer, 0)
static const char nl
Definition: Ostream.H:297
const Type & operator()(const Type &a, const Type &b) const