64 void Foam::cyclicFaPatch::calcTransforms()
66 const label sizeby2 = this->size()/2;
71 for (label edgei=0; edgei < sizeby2; ++edgei)
74 half1Ctrs[edgei] = this->
edgeCentres()[edgei + sizeby2];
80 const vectorField eN(edgeNormals()*magEdgeLengths());
82 scalar maxMatchError = 0;
85 for (label edgei = 0; edgei < sizeby2; ++edgei)
87 half0Normals[edgei] = eN[edgei];
88 half1Normals[edgei] = eN[edgei + sizeby2];
90 scalar magLe =
mag(half0Normals[edgei]);
91 scalar nbrMagLe =
mag(half1Normals[edgei]);
92 scalar avLe = 0.5*(magLe + nbrMagLe);
94 if (magLe < ROOTVSMALL && nbrMagLe < ROOTVSMALL)
99 half0Normals[edgei] =
point(1, 0, 0);
100 half1Normals[edgei] = half0Normals[edgei];
102 else if (
mag(magLe - nbrMagLe)/avLe > matchTol_)
111 half0Normals[edgei] /= magLe;
112 half1Normals[edgei] /= nbrMagLe;
117 if (maxMatchError > matchTol_)
119 label nbrEdgei = errorEdge + sizeby2;
120 scalar magLe =
mag(half0Normals[errorEdge]);
121 scalar nbrMagLe =
mag(half1Normals[errorEdge]);
122 scalar avLe = 0.5*(magLe + nbrMagLe);
125 <<
"edge " << errorEdge
126 <<
" area does not match neighbour " 127 << nbrEdgei <<
" by " 128 << 100*
mag(magLe - nbrMagLe)/avLe
129 <<
"% -- possible edge ordering problem." <<
endl 130 <<
"patch:" <<
name()
131 <<
" my area:" << magLe
132 <<
" neighbour area:" << nbrMagLe
133 <<
" matching tolerance:" << matchTol_
135 <<
"Mesh edge:" << start() + errorEdge
137 <<
"Neighbour edge:" << start() + nbrEdgei
139 <<
"Other errors also exist, only the largest is reported. " 140 <<
"Please rerun with cyclic debug flag set" 156 if (forwardT().size() > 1 || reverseT().size() > 1)
159 <<
"Transformation tensor is not constant for the cyclic " 160 <<
"patch. Please reconsider your setup and definition of " 161 <<
"cyclic boundaries." <<
endl;
173 const label sizeby2 = deltas.size()/2;
175 scalar maxMatchError = 0;
176 label errorEdge = -1;
178 for (label edgei = 0; edgei < sizeby2; ++edgei)
180 scalar avL = 0.5*(magL[edgei] + magL[edgei + sizeby2]);
184 mag(magL[edgei] - magL[edgei + sizeby2])/avL
193 mag(magL[edgei] - magL[edgei + sizeby2])/avL
199 scalar di = deltas[edgei];
200 scalar dni = deltas[edgei + sizeby2];
202 w[edgei] = dni/(di + dni);
203 w[edgei + sizeby2] = 1 - w[edgei];
207 if (maxMatchError > matchTol_)
209 scalar avL = 0.5*(magL[errorEdge] + magL[errorEdge + sizeby2]);
212 <<
"edge " << errorEdge <<
" and " << errorEdge + sizeby2
213 <<
" areas do not match by " 214 << 100*
mag(magL[errorEdge] - magL[errorEdge + sizeby2])/avL
215 <<
"% -- possible edge ordering problem." <<
nl 216 <<
"Cyclic area match tolerance = " 217 << matchTol_ <<
" patch: " <<
name()
226 const label sizeby2 = deltas.size()/2;
228 for (label edgei = 0; edgei < sizeby2; ++edgei)
230 scalar di = deltas[edgei];
231 scalar dni = deltas[edgei + sizeby2];
233 dc[edgei] = 1.0/(di + dni);
234 dc[edgei + sizeby2] = dc[edgei];
276 const label sizeby2 = patchD.size()/2;
279 auto& pdv = tpdv.ref();
284 for (label edgei = 0; edgei < sizeby2; ++edgei)
286 const vector& ddi = patchD[edgei];
287 const vector& dni = patchD[edgei + sizeby2];
289 pdv[edgei] = ddi - dni;
290 pdv[edgei + sizeby2] = -pdv[edgei];
295 for (label edgei = 0; edgei < sizeby2; ++edgei)
297 const vector& ddi = patchD[edgei];
298 const vector& dni = patchD[edgei + sizeby2];
300 pdv[edgei] = ddi -
transform(forwardT()[0], dni);
301 pdv[edgei + sizeby2] = -
transform(reverseT()[0], pdv[edgei]);
314 return patchInternalField(internalData);
325 patchInternalField(internalData, edgeFaces, tpfld.ref());
337 auto& pnf = tpnf.ref();
339 const label sizeby2 = this->size()/2;
341 for (label edgei=0; edgei < sizeby2; ++edgei)
343 pnf[edgei] = interfaceData[edgei + sizeby2];
344 pnf[edgei + sizeby2] = interfaceData[edgei];
357 return internalFieldTransfer(commsType, iF, this->
faceCells());
369 auto& pnf = tpnf.ref();
371 const label sizeby2 = this->size()/2;
373 for (label edgei=0; edgei < sizeby2; ++edgei)
375 pnf[edgei] = iF[edgeCells[edgei + sizeby2]];
376 pnf[edgei + sizeby2] = iF[edgeCells[edgei]];
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
virtual tmp< labelField > interfaceInternalField(const labelUList &internalData) const
Return the values of the given internal data adjacent to the interface as a field.
static const scalar matchTol_
Relative tolerance (for geometric matching). Is factor of.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
commsTypes
Communications types.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual tmp< labelField > internalFieldTransfer(const Pstream::commsTypes commsType, const labelUList &internalData) const
Return neighbour field.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual void movePoints(PstreamBuffers &, const pointField &)
Correct patches after moving points.
constexpr char nl
The newline '\n' character (0x0a)
void makeDeltaCoeffs(scalarField &) const
Make patch face - neighbour cell distances.
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
Smooth ATC in cells next to a set of patches supplied by type.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
virtual tmp< vectorField > delta() const =0
Return delta (P to N) vectors across coupled patch.
Macros for easy insertion into run-time selection tables.
UList< label > labelUList
A UList of labels.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
virtual void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ const Field< vector > edgeCentres(faMeshTools::flattenEdgeField(aMesh.edgeCentres(), true))
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void makeWeights(scalarField &) const
Make patch weighting factors.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
errorManip< error > abort(error &err)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
defineTypeNameAndDebug(combustionModel, 0)
Buffers for inter-processor communications streams (UOPstream, UIPstream).
cyclicFaPatch(const word &name, const dictionary &dict, const label index, const faBoundaryMesh &bm, const word &patchType)
Construct from dictionary.
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
vector point
Point is a vector.
coupledFaPatch is an abstract base class for patches that couple regions of the computational domain ...
Finite area boundary mesh.
Field< vector > vectorField
Specialisation of Field<T> for vector.
virtual tmp< labelField > transfer(const Pstream::commsTypes commsType, const labelUList &interfaceData) const
Transfer and return neighbour field.
A class for managing temporary objects.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
virtual void movePoints(PstreamBuffers &, const pointField &)
Correct patch after moving points.
virtual tmp< vectorField > delta() const
Return delta (P to N) vectors across coupled patch.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.