Merge branch 'volPointInterpolation'

This commit is contained in:
mattijs 2010-02-17 15:20:56 +00:00
commit 827d7d7a8f
55 changed files with 1299 additions and 2332 deletions

View File

@ -666,7 +666,7 @@ int main(int argc, char *argv[])
)
);
pointMesh procPMesh(procMesh, true);
pointMesh procPMesh(procMesh);
pointFieldDecomposer fieldDecomposer
(

View File

@ -109,7 +109,6 @@ Foam::domainDecomposition::domainDecomposition(const IOobject& io)
procNeighbourProcessors_(nProcs_),
procProcessorPatchSize_(nProcs_),
procProcessorPatchStartIndex_(nProcs_),
globallySharedPoints_(0),
cyclicParallel_(false)
{
if (decompositionDict_.found("distributed"))
@ -132,15 +131,6 @@ bool Foam::domainDecomposition::writeDecomposition()
{
Info<< "\nConstructing processor meshes" << endl;
// Make a lookup map for globally shared points
Map<label> sharedPointLookup(2*globallySharedPoints_.size());
forAll(globallySharedPoints_, pointi)
{
sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
}
// Mark point/faces/cells that are in zones.
// -1 : not in zone
// -2 : in multiple zones

View File

@ -103,9 +103,6 @@ class domainDecomposition
//- Start indices for inter-processor patches
labelListList procProcessorPatchStartIndex_;
//- List of globally shared point labels
labelList globallySharedPoints_;
//- Are there cyclic-parallel faces
bool cyclicParallel_;

View File

@ -643,79 +643,6 @@ void Foam::domainDecomposition::decomposeMesh()
// Reset the size of used points
procPointLabels.setSize(nUsedPoints);
}
// Gather data about globally shared points
// Memory management
{
labelList pointsUsage(nPoints(), 0);
// Globally shared points are the ones used by more than 2 processors
// Size the list approximately and gather the points
labelHashSet gSharedPoints
(
min(100, nPoints()/1000)
);
// Loop through all the processors and mark up points used by
// processor boundaries. When a point is used twice, it is a
// globally shared point
for (label procI = 0; procI < nProcs_; procI++)
{
// Get list of face labels
const labelList& curFaceLabels = procFaceAddressing_[procI];
// Get start of processor faces
const labelList& curProcessorPatchStarts =
procProcessorPatchStartIndex_[procI];
const labelList& curProcessorPatchSizes =
procProcessorPatchSize_[procI];
// Reset the lookup list
pointsUsage = 0;
forAll(curProcessorPatchStarts, patchi)
{
const label curStart = curProcessorPatchStarts[patchi];
const label curEnd = curStart + curProcessorPatchSizes[patchi];
for
(
label facei = curStart;
facei < curEnd;
facei++
)
{
// Mark the original face as used
// Remember to decrement the index by one (turning index)
//
const label curF = mag(curFaceLabels[facei]) - 1;
const face& f = fcs[curF];
forAll(f, pointi)
{
if (pointsUsage[f[pointi]] == 0)
{
// Point not previously used
pointsUsage[f[pointi]] = patchi + 1;
}
else if (pointsUsage[f[pointi]] != patchi + 1)
{
// Point used by some other patch = global point!
gSharedPoints.insert(f[pointi]);
}
}
}
}
}
// Grab the result from the hash list
globallySharedPoints_ = gSharedPoints.toc();
sort(globallySharedPoints_);
}
}
// ************************************************************************* //

View File

@ -26,7 +26,6 @@ License
#include "pointFieldDecomposer.H"
#include "processorPointPatchFields.H"
#include "globalPointPatchFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -40,12 +39,8 @@ Foam::pointFieldDecomposer::decomposeField
// Create and map the internal field values
Field<Type> internalField(field.internalField(), pointAddressing_);
// Create a list of pointers for the patchFields including one extra
// for the global patch
PtrList<pointPatchField<Type> > patchFields
(
boundaryAddressing_.size() + 1
);
// Create a list of pointers for the patchFields
PtrList<pointPatchField<Type> > patchFields(boundaryAddressing_.size());
// Create and map the patch field values
forAll(boundaryAddressing_, patchi)
@ -78,17 +73,6 @@ Foam::pointFieldDecomposer::decomposeField
}
}
// Add the global patch
patchFields.set
(
boundaryAddressing_.size(),
new globalPointPatchField<Type>
(
procMesh_.boundary().globalPatch(),
DimensionedField<Type, pointMesh>::null()
)
);
// Create the field for the processor
return tmp<GeometricField<Type, pointPatchField, pointMesh> >
(

View File

@ -37,7 +37,6 @@ SourceFiles
#ifndef lagrangianWriter_H
#define lagrangianWriter_H
#include "globalPointPatch.H"
#include "OFstream.H"
#include "Cloud.H"
#include "volFields.H"

View File

@ -643,7 +643,6 @@ DebugSwitches
perfectInterface 0;
pointIndexHitList 0;
pointPatchField 0;
pointPatchInterpolation 0;
pointScalarField 0;
pointScalarField::DimensionedInternalField 0;
pointSet 0;

View File

@ -453,7 +453,9 @@ $(constraintPointPatches)/processor/processorPointPatch.C
derivedPointPatches = $(pointPatches)/derived
$(derivedPointPatches)/coupled/coupledFacePointPatch.C
/*
$(derivedPointPatches)/global/globalPointPatch.C
*/
$(derivedPointPatches)/wall/wallPointPatch.C
pointBoundaryMesh = $(pointMesh)/pointBoundaryMesh
@ -508,7 +510,9 @@ $(constraintPointPatchFields)/processor/processorPointPatchFields.C
derivedPointPatchFields = $(pointPatchFields)/derived
$(derivedPointPatchFields)/slip/slipPointPatchFields.C
/*
$(derivedPointPatchFields)/global/globalPointPatchFields.C
*/
$(derivedPointPatchFields)/uniformFixedValue/uniformFixedValuePointPatchFields.C
$(derivedPointPatchFields)/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchFields.C

View File

@ -121,11 +121,19 @@ public:
) = 0;
//- Initialise swap of patch point values
virtual void initSwapAdd(Field<Type>&) const
virtual void initSwapAddSeparated
(
const Pstream::commsTypes,
Field<Type>&
) const
{}
//- Complete swap of patch point values and add to local values
virtual void swapAdd(Field<Type>&) const = 0;
virtual void swapAddSeparated
(
const Pstream::commsTypes,
Field<Type>&
) const = 0;
};

View File

@ -122,7 +122,11 @@ cyclicPointPatchField<Type>::cyclicPointPatchField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void cyclicPointPatchField<Type>::swapAdd(Field<Type>& pField) const
void cyclicPointPatchField<Type>::swapAddSeparated
(
const Pstream::commsTypes,
Field<Type>& pField
) const
{
Field<Type> pf(this->patchInternalField(pField));
@ -145,7 +149,7 @@ void cyclicPointPatchField<Type>::swapAdd(Field<Type>& pField) const
}
}
addToInternalField(pField, pf);
addToInternalField(pField, pf, cyclicPatch_.separatedPoints());
}

View File

@ -159,7 +159,11 @@ public:
{}
//- Complete swap of patch point values and add to local values
virtual void swapAdd(Field<Type>&) const;
virtual void swapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>&
) const;
};

View File

@ -96,16 +96,28 @@ processorPointPatchField<Type>::~processorPointPatchField()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void processorPointPatchField<Type>::initSwapAdd(Field<Type>& pField) const
void processorPointPatchField<Type>::initSwapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>& pField
)
const
{
if (Pstream::parRun())
{
// Get internal field into my point order
Field<Type> pf(this->patchInternalField(pField));
// Get internal field into correct order for opposite side
Field<Type> pf
(
this->patchInternalField
(
pField,
procPatch_.reverseMeshPoints()
)
);
OPstream::write
(
Pstream::blocking,
commsType,
procPatch_.neighbProcNo(),
reinterpret_cast<const char*>(pf.begin()),
pf.byteSize()
@ -115,7 +127,11 @@ void processorPointPatchField<Type>::initSwapAdd(Field<Type>& pField) const
template<class Type>
void processorPointPatchField<Type>::swapAdd(Field<Type>& pField) const
void processorPointPatchField<Type>::swapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>& pField
) const
{
if (Pstream::parRun())
{
@ -123,7 +139,7 @@ void processorPointPatchField<Type>::swapAdd(Field<Type>& pField) const
IPstream::read
(
Pstream::blocking,
commsType,
procPatch_.neighbProcNo(),
reinterpret_cast<char*>(pnf.begin()),
pnf.byteSize()
@ -140,22 +156,20 @@ void processorPointPatchField<Type>::swapAdd(Field<Type>& pField) const
}
else
{
const labelList& nonGlobalPatchPoints =
procPatch_.nonGlobalPatchPoints();
const labelListList& pointFaces = ppp.pointFaces();
forAll(nonGlobalPatchPoints, pfi)
forAll(pointFaces, pfi)
{
pnf[pfi] = transform
(
forwardT[pointFaces[nonGlobalPatchPoints[pfi]][0]],
forwardT[pointFaces[pfi][0]],
pnf[pfi]
);
}
}
}
addToInternalField(pField, pnf);
addToInternalField(pField, pnf, procPatch_.separatedPoints());
}
}

View File

@ -169,11 +169,19 @@ public:
)
{}
//- Initialise swap of patch point values
virtual void initSwapAdd(Field<Type>&) const;
//- Initialise swap of non-collocated patch point values
virtual void initSwapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>&
) const;
//- Complete swap of patch point values and add to local values
virtual void swapAdd(Field<Type>&) const;
virtual void swapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>&
) const;
};

View File

@ -1,169 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "globalPointPatchField.H"
#include "PstreamCombineReduceOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
template<class Type>
globalPointPatchField<Type>::globalPointPatchField
(
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF
)
:
coupledPointPatchField<Type>(p, iF),
globalPointPatch_(refCast<const globalPointPatch>(p))
{}
template<class Type>
globalPointPatchField<Type>::globalPointPatchField
(
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const dictionary& dict
)
:
coupledPointPatchField<Type>(p, iF, dict),
globalPointPatch_(refCast<const globalPointPatch>(p))
{
if (!isType<globalPointPatch>(p))
{
FatalIOErrorIn
(
"globalPointPatchField<Type>::globalPointPatchField\n"
"(\n"
" const pointPatch& p,\n"
" const Field<Type>& field,\n"
" const dictionary& dict\n"
")\n",
dict
) << "patch " << this->patch().index()
<< " not processorPoint type. "
<< "Patch type = " << p.type()
<< exit(FatalIOError);
}
}
template<class Type>
globalPointPatchField<Type>::globalPointPatchField
(
const globalPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper& mapper
)
:
coupledPointPatchField<Type>(ptf, p, iF, mapper),
globalPointPatch_(refCast<const globalPointPatch>(ptf.patch()))
{
if (!isType<globalPointPatch>(this->patch()))
{
FatalErrorIn
(
"globalPointPatchField<Type>::globalPointPatchField\n"
"(\n"
" const globalPointPatchField<Type>& ptf,\n"
" const pointPatch& p,\n"
" const DimensionedField<Type, pointMesh>& iF,\n"
" const pointPatchFieldMapper& mapper\n"
")\n"
) << "Field type does not correspond to patch type for patch "
<< this->patch().index() << "." << endl
<< "Field type: " << typeName << endl
<< "Patch type: " << this->patch().type()
<< exit(FatalError);
}
}
template<class Type>
globalPointPatchField<Type>::globalPointPatchField
(
const globalPointPatchField<Type>& ptf,
const DimensionedField<Type, pointMesh>& iF
)
:
coupledPointPatchField<Type>(ptf, iF),
globalPointPatch_(refCast<const globalPointPatch>(ptf.patch()))
{}
// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
globalPointPatchField<Type>::~globalPointPatchField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class Type>
void globalPointPatchField<Type>::swapAdd(Field<Type>& pField) const
{
// Create the global list and insert local values
if (globalPointPatch_.globalPointSize() > 0)
{
Field<Type> lpf = patchInternalField(pField);
const labelList& addr = globalPointPatch_.sharedPointAddr();
Field<Type> gpf
(
globalPointPatch_.globalPointSize(),
pTraits<Type>::zero
);
forAll(addr, i)
{
gpf[addr[i]] += lpf[i];
}
combineReduce(gpf, plusEqOp<Field<Type> >());
// Extract local data
forAll (addr, i)
{
lpf[i] = gpf[addr[i]];
}
setInInternalField(pField, lpf);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,166 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::globalPointPatchField
Description
Foam::globalPointPatchField
SourceFiles
globalPointPatchField.C
\*---------------------------------------------------------------------------*/
#ifndef globalPointPatchField_H
#define globalPointPatchField_H
#include "coupledPointPatchField.H"
#include "globalPointPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class globalPointPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class globalPointPatchField
:
public coupledPointPatchField<Type>
{
// Private data
//- Local reference to processorPoint patch
const globalPointPatch& globalPointPatch_;
public:
//- Runtime type information
TypeName("global");
// Constructors
//- Construct from patch and internal field
globalPointPatchField
(
const pointPatch&,
const DimensionedField<Type, pointMesh>&
);
//- Construct from patch, internal field and dictionary
globalPointPatchField
(
const pointPatch&,
const DimensionedField<Type, pointMesh>&,
const dictionary&
);
//- Construct by mapping given patchField<Type> onto a new patch
globalPointPatchField
(
const globalPointPatchField<Type>&,
const pointPatch&,
const DimensionedField<Type, pointMesh>&,
const pointPatchFieldMapper&
);
//- Construct and return a clone
virtual autoPtr<pointPatchField<Type> > clone() const
{
return autoPtr<pointPatchField<Type> >
(
new globalPointPatchField<Type>
(
*this
)
);
}
//- Construct as copy setting internal field reference
globalPointPatchField
(
const globalPointPatchField<Type>&,
const DimensionedField<Type, pointMesh>&
);
//- Construct and return a clone setting internal field reference
virtual autoPtr<pointPatchField<Type> > clone
(
const DimensionedField<Type, pointMesh>& iF
) const
{
return autoPtr<pointPatchField<Type> >
(
new globalPointPatchField
<Type>
(
*this,
iF
)
);
}
// Destructor
~globalPointPatchField();
// Member functions
// Evaluation functions
//- Evaluate the patch field
virtual void evaluate
(
const Pstream::commsTypes commsType=Pstream::blocking
)
{}
//- Complete swap of patch point values and add to local values
virtual void swapAdd(Field<Type>&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "globalPointPatchField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,44 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "globalPointPatchFields.H"
#include "pointPatchFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
makePointPatchFields(global);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,51 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#ifndef globalPointPatchFields_H
#define globalPointPatchFields_H
#include "globalPointPatchField.H"
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePointPatchFieldTypedefs(global);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -124,7 +124,8 @@ template<class Type>
template<class Type1>
tmp<Field<Type1> > pointPatchField<Type>::patchInternalField
(
const Field<Type1>& iF
const Field<Type1>& iF,
const labelList& meshPoints
) const
{
// Check size
@ -141,9 +142,6 @@ tmp<Field<Type1> > pointPatchField<Type>::patchInternalField
<< abort(FatalError);
}
// get addressing
const labelList& meshPoints = patch().meshPoints();
tmp<Field<Type1> > tvalues(new Field<Type1>(meshPoints.size()));
Field<Type1>& values = tvalues();
@ -156,6 +154,17 @@ tmp<Field<Type1> > pointPatchField<Type>::patchInternalField
}
template<class Type>
template<class Type1>
tmp<Field<Type1> > pointPatchField<Type>::patchInternalField
(
const Field<Type1>& iF
) const
{
return patchInternalField(iF, patch().meshPoints());
}
template<class Type>
template<class Type1>
void pointPatchField<Type>::addToInternalField
@ -201,6 +210,53 @@ void pointPatchField<Type>::addToInternalField
}
template<class Type>
template<class Type1>
void pointPatchField<Type>::addToInternalField
(
Field<Type1>& iF,
const Field<Type1>& pF,
const labelList& points
) const
{
// Check size
if (iF.size() != internalField().size())
{
FatalErrorIn
(
"void pointPatchField<Type>::"
"addToInternalField("
"Field<Type1>& iF, const Field<Type1>& iF, const labelList&) const"
) << "given internal field does not correspond to the mesh. "
<< "Field size: " << iF.size()
<< " mesh size: " << internalField().size()
<< abort(FatalError);
}
if (pF.size() != size())
{
FatalErrorIn
(
"void pointPatchField<Type>::"
"addToInternalField("
"Field<Type1>& iF, const Field<Type1>& iF, const labelList&) const"
) << "given patch field does not correspond to the mesh. "
<< "Field size: " << pF.size()
<< " mesh size: " << size()
<< abort(FatalError);
}
// Get the addressing
const labelList& mp = patch().meshPoints();
forAll(points, i)
{
label pointI = points[i];
iF[mp[pointI]] += pF[pointI];
}
}
template<class Type>
template<class Type1>
void pointPatchField<Type>::setInInternalField

View File

@ -291,6 +291,15 @@ public:
const Field<Type1>& iF
) const;
//- Return field created from selected internal field values
// given internal field reference
template<class Type1>
tmp<Field<Type1> > patchInternalField
(
const Field<Type1>& iF,
const labelList& meshPoints
) const;
//- Given the internal field and a patch field,
// add the patch field to the internal field
template<class Type1>
@ -300,6 +309,16 @@ public:
const Field<Type1>& pF
) const;
//- Given the internal field and a patch field,
// add selected elements of the patch field to the internal field
template<class Type1>
void addToInternalField
(
Field<Type1>& iF,
const Field<Type1>& pF,
const labelList& points
) const;
//- Given the internal field and a patch field,
// set the patch field in the internal field
template<class Type1>

View File

@ -39,6 +39,7 @@ SourceFiles
#include "labelField.H"
#include "typeInfo.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -39,6 +39,7 @@ SourceFiles
#include "lduInterface.H"
#include "primitiveFieldsFwd.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -27,7 +27,7 @@ License
#include "pointBoundaryMesh.H"
#include "polyBoundaryMesh.H"
#include "facePointPatch.H"
#include "globalPointPatch.H"
#include "pointMesh.H"
#include "PstreamBuffers.H"
#include "lduSchedule.H"
#include "globalMeshData.H"
@ -105,31 +105,6 @@ void Foam::pointBoundaryMesh::calcGeometry()
}
const Foam::globalPointPatch&
Foam::pointBoundaryMesh::globalPatch() const
{
const pointPatchList& patches = *this;
forAll (patches, patchI)
{
if (isType<globalPointPatch>(patches[patchI]))
{
return refCast<const globalPointPatch>(patches[patchI]);
}
}
FatalErrorIn
(
"const pointBoundaryMesh::"
"globalPointPatch& globalPatch() const"
) << "patch not found."
<< abort(FatalError);
// Dummy return
return refCast<const globalPointPatch>(patches[0]);
}
void Foam::pointBoundaryMesh::movePoints(const pointField& p)
{
PstreamBuffers pBufs(Pstream::defaultCommsType);

View File

@ -46,7 +46,6 @@ namespace Foam
// Forward declaration of classes
class pointMesh;
class polyBoundaryMesh;
class globalPointPatch;
/*---------------------------------------------------------------------------*\
Class pointBoundaryMesh Declaration
@ -98,9 +97,6 @@ public:
return mesh_;
}
//- Return reference to globalPointPatch
const globalPointPatch& globalPatch() const;
//- Correct polyBoundaryMesh after moving points
void movePoints(const pointField&);

View File

@ -26,7 +26,6 @@ License
#include "pointMesh.H"
#include "globalMeshData.H"
#include "globalPointPatch.H"
#include "pointMeshMapper.H"
#include "pointFields.H"
#include "MapGeometricFields.H"
@ -56,36 +55,12 @@ void Foam::pointMesh::mapFields(const mapPolyMesh& mpm)
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pointMesh::pointMesh
(
const polyMesh& pMesh,
bool alwaysConstructGlobalPatch
)
Foam::pointMesh::pointMesh(const polyMesh& pMesh)
:
MeshObject<polyMesh, pointMesh>(pMesh),
GeoMesh<polyMesh>(pMesh),
boundary_(*this, pMesh.boundaryMesh())
{
// Add the globalPointPatch if there are global points
if
(
alwaysConstructGlobalPatch
|| GeoMesh<polyMesh>::mesh_.globalData().nGlobalPoints()
)
{
boundary_.setSize(boundary_.size() + 1);
boundary_.set
(
boundary_.size() - 1,
new globalPointPatch
(
boundary_,
boundary_.size() - 1
)
);
}
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
}

View File

@ -79,11 +79,7 @@ public:
// Constructors
//- Construct from polyMesh
explicit pointMesh
(
const polyMesh& pMesh,
bool alwaysConstructGlobalPatch = false
);
explicit pointMesh(const polyMesh& pMesh);
// Member Functions

View File

@ -26,7 +26,8 @@ Class
Foam::genericPointPatch
Description
DirectMapped patch.
Substitute for unknown patches. Used for postprocessing when only
basic polyPatch info is needed.
SourceFiles
genericPointPatch.C

View File

@ -28,7 +28,6 @@ License
#include "pointBoundaryMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "pointMesh.H"
#include "globalPointPatch.H"
#include "edgeList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -51,104 +50,31 @@ addToRunTimeSelectionTable
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::cyclicPointPatch::initGeometry(PstreamBuffers&)
{
transformPairs_.setSize(0);
}
{}
void Foam::cyclicPointPatch::calcGeometry(PstreamBuffers&)
{
const edgeList& cp = cyclicPolyPatch_.coupledPoints();
const labelList& mp = cyclicPolyPatch_.meshPoints();
const pointField& points = cyclicPolyPatch_.points();
// If there are no global points create a 1->1 map
if (!boundaryMesh().mesh().globalData().nGlobalPoints())
DynamicList<label> separated;
forAll(cp, i)
{
nonGlobalPatchPoints_.setSize(mp.size());
forAll(nonGlobalPatchPoints_, i)
{
nonGlobalPatchPoints_[i] = i;
}
const edge& coupledSet = cp[i];
meshPoints_ = cyclicPolyPatch_.meshPoints();
transformPairs_ = cp;
// Assume all points are separated.
separated.append(coupledSet[0]);
separated.append(coupledSet[1]);
}
else
separatedPoints_.transfer(separated);
if (debug)
{
// Get reference to shared points
const labelList& sharedPoints =
boundaryMesh().globalPatch().meshPoints();
nonGlobalPatchPoints_.setSize(mp.size());
meshPoints_.setSize(mp.size());
labelList pointMap(mp.size(), -1);
label noFiltPoints = 0;
forAll (mp, pointI)
{
label curP = mp[pointI];
bool found = false;
forAll (sharedPoints, sharedI)
{
if (sharedPoints[sharedI] == curP)
{
found = true;
break;
}
}
if (!found)
{
pointMap[pointI] = noFiltPoints;
nonGlobalPatchPoints_[noFiltPoints] = pointI;
meshPoints_[noFiltPoints] = curP;
noFiltPoints++;
}
}
nonGlobalPatchPoints_.setSize(noFiltPoints);
meshPoints_.setSize(noFiltPoints);
transformPairs_.setSize(cp.size());
label noFiltPointPairs = 0;
forAll(cp, i)
{
if (pointMap[cp[i][0]] != -1 && pointMap[cp[i][1]] != -1)
{
transformPairs_[noFiltPointPairs][0] = pointMap[cp[i][0]];
transformPairs_[noFiltPointPairs][1] = pointMap[cp[i][1]];
noFiltPointPairs++;
}
else if (pointMap[cp[i][0]] == -1 && pointMap[cp[i][1]] != -1)
{
FatalErrorIn
(
"cyclicPointPatch::calcGeometry(PstreamBuffers&) const"
) << "Point " << cp[i][0] << "of point-pair " << i
<< " is a global point but the other point "
<< cp[i][1] << " is not"
<< exit(FatalError);
}
else if (pointMap[cp[i][0]] != -1 && pointMap[cp[i][1]] == -1)
{
FatalErrorIn
(
"cyclicPointPatch::calcGeometry(PstreamBuffers&) const"
) << "Point " << cp[i][1] << "of point-pair " << i
<< " is a global point but the other point "
<< cp[i][0] << " is not"
<< exit(FatalError);
}
}
transformPairs_.setSize(noFiltPointPairs);
Pout<< "cyclic:" << cyclicPolyPatch_.name()
<< " separated:" << separatedPoints_.size()
<< " out of points:" << mp.size() << endl;
}
}
@ -198,7 +124,13 @@ cyclicPointPatch::~cyclicPointPatch()
const edgeList& cyclicPointPatch::transformPairs() const
{
return transformPairs_;
return cyclicPolyPatch_.coupledPoints();
}
const labelList& cyclicPointPatch::separatedPoints() const
{
return separatedPoints_;
}

View File

@ -57,6 +57,8 @@ class cyclicPointPatch
//- Local reference cast into the cyclic patch
const cyclicPolyPatch& cyclicPolyPatch_;
//- List of local points that are not collocated
mutable labelList separatedPoints_;
// Private Member Functions
@ -69,10 +71,6 @@ class cyclicPointPatch
// Demand driven private data
//- The set of pairs of points that require transformation
// and/or mapping
edgeList transformPairs_;
//- Initialise the calculation of the patch geometry
virtual void initGeometry(PstreamBuffers&);
@ -147,6 +145,9 @@ public:
//- Return the set of pairs of points that require transformation
// and/or mapping
virtual const edgeList& transformPairs() const;
//- List of separated coupled points
virtual const labelList& separatedPoints() const;
};

View File

@ -28,7 +28,6 @@ License
#include "pointBoundaryMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "pointMesh.H"
#include "globalPointPatch.H"
#include "faceList.H"
#include "primitiveFacePatch.H"
#include "emptyPolyPatch.H"
@ -58,34 +57,22 @@ void Foam::processorPointPatch::initGeometry(PstreamBuffers& pBufs)
// Depending on whether the patch is a master or a slave, get the primitive
// patch points and filter away the points from the global patch.
if (isMaster())
// Create the reversed patch and pick up its points
// so that the order is correct
const polyPatch& pp = patch();
faceList masterFaces(pp.size());
forAll (pp, faceI)
{
meshPoints_ = procPolyPatch_.meshPoints();
}
else
{
// Slave side. Create the reversed patch and pick up its points
// so that the order is correct
const polyPatch& pp = patch();
faceList masterFaces(pp.size());
forAll (pp, faceI)
{
masterFaces[faceI] = pp[faceI].reverseFace();
}
meshPoints_ = primitiveFacePatch
(
masterFaces,
pp.points()
).meshPoints();
masterFaces[faceI] = pp[faceI].reverseFace();
}
if (Pstream::parRun())
{
initPatchPatchPoints(pBufs);
}
reverseMeshPoints_ = primitiveFacePatch
(
masterFaces,
pp.points()
).meshPoints();
}
@ -93,261 +80,46 @@ void Foam::processorPointPatch::calcGeometry(PstreamBuffers& pBufs)
{
if (Pstream::parRun())
{
calcPatchPatchPoints(pBufs);
}
const boolList& collocated = procPolyPatch_.collocated();
// If it is not runing parallel or there are no global points
// create a 1->1 map
if
(
!Pstream::parRun()
|| !boundaryMesh().mesh().globalData().nGlobalPoints()
)
{
nonGlobalPatchPoints_.setSize(meshPoints_.size());
forAll(nonGlobalPatchPoints_, i)
if (collocated.size() == 0)
{
nonGlobalPatchPoints_[i] = i;
separatedPoints_.setSize(0);
}
}
else
{
// Get reference to shared points
const labelList& sharedPoints =
boundaryMesh().globalPatch().meshPoints();
nonGlobalPatchPoints_.setSize(meshPoints_.size());
label noFiltPoints = 0;
forAll (meshPoints_, pointI)
else if (collocated.size() == 1)
{
label curP = meshPoints_[pointI];
bool found = false;
forAll (sharedPoints, sharedI)
// Uniformly
if (collocated[0])
{
if (sharedPoints[sharedI] == curP)
separatedPoints_.setSize(0);
}
else
{
separatedPoints_ = identity(size());
}
}
else
{
// Per face collocated or not.
const labelListList& pointFaces = procPolyPatch_.pointFaces();
DynamicList<label> separated;
forAll(pointFaces, pfi)
{
if (!collocated[pointFaces[pfi][0]])
{
found = true;
break;
separated.append(pfi);
}
}
if (!found)
{
nonGlobalPatchPoints_[noFiltPoints] = pointI;
meshPoints_[noFiltPoints] = curP;
noFiltPoints++;
}
}
nonGlobalPatchPoints_.setSize(noFiltPoints);
meshPoints_.setSize(noFiltPoints);
}
}
void processorPointPatch::initPatchPatchPoints(PstreamBuffers& pBufs)
{
if (debug)
{
Info<< "processorPointPatch::initPatchPatchPoints(PstreamBuffers&) : "
<< "constructing patch-patch points"
<< endl;
}
const polyBoundaryMesh& bm = boundaryMesh().mesh()().boundaryMesh();
// Get the mesh points for this patch corresponding to the faces
const labelList& ppmp = meshPoints();
// Create a HashSet of the point labels for this patch
Map<label> patchPointSet(2*ppmp.size());
forAll (ppmp, ppi)
{
patchPointSet.insert(ppmp[ppi], ppi);
}
// Create the lists of patch-patch points
labelListList patchPatchPoints(bm.size());
// Create the lists of patch-patch point normals
List<List<vector> > patchPatchPointNormals(bm.size());
// Loop over all patches looking for other patches that share points
forAll(bm, patchi)
{
if
(
patchi != index() // Ignore self-self
&& !isA<emptyPolyPatch>(bm[patchi]) // Ignore empty
&& !bm[patchi].coupled() // Ignore other couples
)
{
// Get the meshPoints for the other patch
const labelList& meshPoints = bm[patchi].meshPoints();
// Get the normals for the other patch
const vectorField& normals = bm[patchi].pointNormals();
label pppi = 0;
forAll(meshPoints, pointi)
{
label ppp = meshPoints[pointi];
// Check to see if the point of the other patch is shared with
// this patch
Map<label>::iterator iter = patchPointSet.find(ppp);
if (iter != patchPointSet.end())
{
// If it is shared initialise the patchPatchPoints for this
// patch
if (!patchPatchPoints[patchi].size())
{
patchPatchPoints[patchi].setSize(ppmp.size());
patchPatchPointNormals[patchi].setSize(ppmp.size());
}
// and add the entry
patchPatchPoints[patchi][pppi] = iter();
patchPatchPointNormals[patchi][pppi] = normals[pointi];
pppi++;
}
}
// Resise the list of shared points and normals for the patch
// being considerd
patchPatchPoints[patchi].setSize(pppi);
patchPatchPointNormals[patchi].setSize(pppi);
separatedPoints_.transfer(separated);
}
}
// Send the patchPatchPoints to the neighbouring processor
UOPstream toNeighbProc(neighbProcNo(), pBufs);
toNeighbProc
<< ppmp.size() // number of points for checking
<< patchPatchPoints
<< patchPatchPointNormals;
if (debug)
{
Info<< "processorPointPatch::initPatchPatchPoints() : "
<< "constructed patch-patch points"
<< endl;
}
}
void Foam::processorPointPatch::calcPatchPatchPoints(PstreamBuffers& pBufs)
{
// Get the patchPatchPoints from the neighbouring processor
UIPstream fromNeighbProc(neighbProcNo(), pBufs);
label nbrNPoints(readLabel(fromNeighbProc));
labelListList patchPatchPoints(fromNeighbProc);
List<List<vector> > patchPatchPointNormals(fromNeighbProc);
pointBoundaryMesh& pbm = const_cast<pointBoundaryMesh&>(boundaryMesh());
const labelList& ppmp = meshPoints();
// Simple check for the very rare situation when not the same number
// of points on both sides. This can happen with decomposed cyclics.
// If on one side the cyclic shares a point with proc faces coming from
// internal faces it will have a different number of points from
// the situation where the cyclic and the 'normal' proc faces are fully
// separate.
if (nbrNPoints != ppmp.size())
{
WarningIn("processorPointPatch::calcPatchPatchPoints(PstreamBuffers&)")
<< "Processor patch " << name()
<< " has " << ppmp.size() << " points; coupled patch has "
<< nbrNPoints << " points." << endl
<< " (usually due to decomposed cyclics)."
<< " This might give problems" << endl
<< " when using point fields (interpolation, mesh motion)."
<< endl;
}
// Loop over the patches looking for other patches that share points
forAll(patchPatchPoints, patchi)
{
const labelList& patchPoints = patchPatchPoints[patchi];
const List<vector>& patchPointNormals = patchPatchPointNormals[patchi];
// If there are potentially shared points for the patch being considered
if (patchPoints.size())
{
// Get the current meshPoints list for the patch
facePointPatch& fpp = refCast<facePointPatch>(pbm[patchi]);
const labelList& fmp = fpp.meshPoints();
labelList& mp = fpp.meshPoints_;
const vectorField& fnormals = fpp.pointNormals();
vectorField& normals = fpp.pointNormals_;
// Create a HashSet of the point labels for the patch
Map<label> patchPointSet(2*fmp.size());
forAll (fmp, ppi)
{
patchPointSet.insert(fmp[ppi], ppi);
}
label nPoints = mp.size();
label lpi = 0;
bool resized = false;
// For each potentially shared point...
forAll(patchPoints, ppi)
{
// Check if it is not already in the patch,
// i.e. not part of a face of the patch
if (!patchPointSet.found(ppmp[patchPoints[ppi]]))
{
// If it isn't already in the patch check if the local
// meshPoints is already set and if not initialise the
// meshPoints_ and pointNormals_
if (!resized)
{
if (!mp.size() && fmp.size())
{
mp = fmp;
normals = fnormals;
nPoints = mp.size();
}
mp.setSize(nPoints + patchPoints.size());
loneMeshPoints_.setSize(patchPoints.size());
normals.setSize(nPoints + patchPoints.size());
resized = true;
}
// Add the new point to the patch
mp[nPoints] = ppmp[patchPoints[ppi]];
loneMeshPoints_[lpi++] = ppmp[patchPoints[ppi]];
normals[nPoints++] = patchPointNormals[ppi];
}
}
// If the lists have been resized points have been added.
// Shrink the lists to the current size.
if (resized)
{
mp.setSize(nPoints);
loneMeshPoints_.setSize(lpi);
normals.setSize(nPoints);
}
}
Pout<< "processor:" << name()
<< " separated:" << separatedPoints_.size()
<< " out of points:" << size() << endl;
}
}
@ -393,6 +165,20 @@ processorPointPatch::~processorPointPatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const labelList& processorPointPatch::reverseMeshPoints() const
{
return reverseMeshPoints_;
}
const labelList& processorPointPatch::separatedPoints() const
{
return separatedPoints_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -64,6 +64,9 @@ class processorPointPatch
const processorPolyPatch& procPolyPatch_;
mutable labelList reverseMeshPoints_;
mutable labelList separatedPoints_;
// Private Member Functions
@ -73,14 +76,6 @@ class processorPointPatch
//- Calculate the patch geometry
virtual void calcGeometry(PstreamBuffers&);
//- Initialise the points on this patch which are should also be
// on a neighbouring patch but are not part of faces of that patch
void initPatchPatchPoints(PstreamBuffers&);
//- Calculate the points on this patch which are should also be
// on a neighbouring patch but are not part of faces of that patch
void calcPatchPatchPoints(PstreamBuffers&);
//- Initialise the patches for moving points
virtual void initMovePoints(PstreamBuffers&, const pointField&);
@ -165,6 +160,13 @@ public:
{
return procPolyPatch_;
}
//- Return mesh points in the correct order for the receiving side
const labelList& reverseMeshPoints() const;
//- List of separated coupled points
virtual const labelList& separatedPoints() const;
};

View File

@ -57,25 +57,6 @@ coupledFacePointPatch::~coupledFacePointPatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const labelList& coupledFacePointPatch::nonGlobalPatchPoints() const
{
return nonGlobalPatchPoints_;
}
const labelList& coupledFacePointPatch::loneMeshPoints() const
{
return loneMeshPoints_;
}
const vectorField& coupledFacePointPatch::pointNormals() const
{
notImplemented("coupledFacePointPatch::pointNormals() const");
return Field<vector>::null();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -73,16 +73,6 @@ class coupledFacePointPatch
protected:
// Demand driven private data
//- The set of labels of the pointPatch points which are
// non-global, i.e. present only in this coupledPointPatch.
// MUST be initialised by calcGeometry()!
labelList nonGlobalPatchPoints_;
labelList loneMeshPoints_;
// Construction of demand-driven data
//- Calculate mesh points
@ -120,20 +110,8 @@ public:
return true;
}
// Access functions for demand driven data
//- Return the set of labels of the pointPatch points which are
// non-global, i.e. present in this coupledFacePointPatch
virtual const labelList& nonGlobalPatchPoints() const;
//- Return the set of labels of the pointPatch points which are
// lone, i.e. present in this coupledFacePointPatch but not
// associated with any faces
virtual const labelList& loneMeshPoints() const;
//- Return point unit normals. Not implemented.
virtual const vectorField& pointNormals() const;
//- List of separated coupled points
virtual const labelList& separatedPoints() const = 0;
};

View File

@ -1,58 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "globalPointPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::globalPointPatch, 0);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::globalPointPatch::globalPointPatch
(
const pointBoundaryMesh& bm,
const label index
)
:
pointPatch(bm),
coupledPointPatch(bm),
index_(index)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::globalPointPatch::~globalPointPatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -1,209 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::globalPointPatch
Description
Foam::globalPointPatch
SourceFiles
globalPointPatch.C
\*---------------------------------------------------------------------------*/
#ifndef globalPointPatch_H
#define globalPointPatch_H
#include "pointPatch.H"
#include "coupledPointPatch.H"
#include "globalMeshData.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class globalPointPatch Declaration
\*---------------------------------------------------------------------------*/
class globalPointPatch
:
public pointPatch,
public coupledPointPatch
{
// Private data
// //- Reference to the basic globalMeshData
// const globalMeshData& globalMeshData_;
//- Index in the boundary mesh
label index_;
// Protected Member Functions
//- Initialise the calculation of the patch geometry
virtual void initGeometry(PstreamBuffers&)
{}
//- Calculate the patch geometry
virtual void calcGeometry(PstreamBuffers&)
{}
//- Initialise the patches for moving points
virtual void initMovePoints(PstreamBuffers&, const pointField&)
{}
//- Correct patches after moving points
virtual void movePoints(PstreamBuffers&, const pointField&)
{}
//- Initialise the update of the patch topology
virtual void initUpdateMesh(PstreamBuffers&)
{}
//- Update of the patch topology
virtual void updateMesh(PstreamBuffers&)
{}
// Private Member Functions
//- Disallow default construct as copy
globalPointPatch
(
const globalPointPatch&
);
//- Disallow default assignment
void operator=(const globalPointPatch&);
public:
//- Runtime type information
TypeName("global");
// Constructors
//- Construct from components
globalPointPatch
(
const pointBoundaryMesh&,
const label index
);
// Destructor
virtual ~globalPointPatch();
// Member functions
//- Return name
virtual const word& name() const
{
// There can only be a single patch of this type - therefore
// its name is hard-coded.
return type();
}
//- Return size
virtual label size() const
{
return meshPoints().size();
}
//- Return true if running parallel
virtual bool coupled() const
{
if (Pstream::parRun())
{
return true;
}
else
{
return false;
}
}
//- Return number of faces
virtual label nFaces() const
{
return 0;
}
//- Return the index of this patch in the pointBoundaryMesh
virtual label index() const
{
return index_;
}
//- Return mesh points
virtual const labelList& meshPoints() const
{
return boundaryMesh().mesh().globalData().sharedPointLabels();
}
//- Return local points. Not implemented
virtual const pointField& localPoints() const
{
notImplemented("globalPointPatch::localPoints() const");
return pointField::null();
}
//- Return point normals. Not implemented
virtual const vectorField& pointNormals() const
{
notImplemented("globalPointPatch::pointNormals() const");
return vectorField::null();
}
//- Return total number of shared points
virtual label globalPointSize() const
{
return boundaryMesh().mesh().globalData().nGlobalPoints();
}
//- Return addressing into the global points list
const labelList& sharedPointAddr() const
{
return boundaryMesh().mesh().globalData().sharedPointAddr();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -52,11 +52,7 @@ addToRunTimeSelectionTable
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void facePointPatch::initGeometry(PstreamBuffers&)
{
meshPoints_.setSize(0);
localPoints_.setSize(0);
pointNormals_.setSize(0);
}
{}
void facePointPatch::calcGeometry(PstreamBuffers&)
@ -94,60 +90,6 @@ facePointPatch::facePointPatch
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const labelList& facePointPatch::meshPoints() const
{
if (meshPoints_.size())
{
return meshPoints_;
}
else
{
return polyPatch_.meshPoints();
}
}
const pointField& facePointPatch::localPoints() const
{
if (meshPoints_.size())
{
if (localPoints_.size() != meshPoints_.size())
{
const labelList& meshPts = meshPoints();
localPoints_.setSize(meshPts.size());
const pointField& points = polyPatch_.points();
forAll (meshPts, pointi)
{
localPoints_[pointi] = points[meshPts[pointi]];
}
}
return localPoints_;
}
else
{
return polyPatch_.localPoints();
}
}
const vectorField& facePointPatch::pointNormals() const
{
if (pointNormals_.size())
{
return pointNormals_;
}
else
{
return polyPatch_.pointNormals();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -30,7 +30,6 @@ Description
SourceFiles
facePointPatch.C
facePointPatchM.C
newPointPatch.C
\*---------------------------------------------------------------------------*/
@ -65,13 +64,6 @@ protected:
//- Reference to the underlying polyPatch
const polyPatch& polyPatch_;
// Optional data used if the pointPatch has points not associated
// with faces, i.e. not accessible via polyPatch
mutable labelList meshPoints_;
mutable pointField localPoints_;
mutable vectorField pointNormals_;
// Protected Member Functions
@ -182,13 +174,24 @@ public:
}
//- Return mesh points
virtual const labelList& meshPoints() const;
virtual const labelList& meshPoints() const
{
return polyPatch_.meshPoints();
}
//- Return pointField of points in patch
virtual const pointField& localPoints() const;
virtual const pointField& localPoints() const
{
return polyPatch_.localPoints();
}
//- Return point unit normals
virtual const vectorField& pointNormals() const;
virtual const vectorField& pointNormals() const
{
return polyPatch_.pointNormals();
}
};

View File

@ -319,6 +319,17 @@ class globalMeshData
void calcGlobalEdgeAllSlaves() const;
//- Synchronise pointwise data
template<class Type, class CombineOp>
void syncPointData
(
List<Type>& pointData,
const labelListList& slaves,
const mapDistribute& slavesMap,
const CombineOp& cop
) const;
//- Disallow default bitwise copy construct
globalMeshData(const globalMeshData&);
@ -505,6 +516,13 @@ public:
// distributed by below map.
const labelListList& globalPointSlaves() const;
const mapDistribute& globalPointSlavesMap() const;
//- Helper to synchronise mesh data
template<class Type, class CombineOp>
void syncPointData
(
List<Type>& pointData,
const CombineOp& cop
) const;
// Coupled edge to coupled edges.
@ -539,6 +557,13 @@ public:
const globalIndex& globalPointAllNumbering()const;
const labelListList& globalPointAllSlaves() const;
const mapDistribute& globalPointAllSlavesMap() const;
//- Helper to synchronise mesh data
template<class Type, class CombineOp>
void syncPointAllData
(
List<Type>& pointData,
const CombineOp& cop
) const;
// Coupled edge to all coupled edges (same numbering as
// collocated)
@ -600,6 +625,13 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "globalMeshDataTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View File

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "globalMeshData.H"
#include "polyMesh.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class CombineOp>
void Foam::globalMeshData::syncPointData
(
List<Type>& pointData,
const labelListList& slaves,
const mapDistribute& slavesMap,
const CombineOp& cop
) const
{
if (pointData.size() != mesh_.nPoints())
{
FatalErrorIn("globalMeshData::syncPointData(..)")
<< "Number of elements in data:" << pointData.size()
<< " differs from number of points in mesh:" << mesh_.nPoints()
<< abort(FatalError);
}
const indirectPrimitivePatch& cpp = coupledPatch();
const labelList& meshPoints = cpp.meshPoints();
// Copy mesh (point)data to coupled patch (point)data
Field<Type> cppFld(slavesMap.constructSize());
forAll(meshPoints, patchPointI)
{
cppFld[patchPointI] = pointData[meshPoints[patchPointI]];
}
// Pull slave data onto master
slavesMap.distribute(cppFld);
// Combine master data with slave data
forAll(slaves, patchPointI)
{
const labelList& slavePoints = slaves[patchPointI];
// Combine master with slave data
forAll(slavePoints, i)
{
cop(cppFld[patchPointI], cppFld[slavePoints[i]]);
}
// Copy result back to slave slots
forAll(slavePoints, i)
{
cppFld[slavePoints[i]] = cppFld[patchPointI];
}
}
// Push master data back to slaves
slavesMap.reverseDistribute(meshPoints.size(), cppFld);
// Update mesh (point)data from coupled patch (point)data
forAll(meshPoints, patchPointI)
{
pointData[meshPoints[patchPointI]] = cppFld[patchPointI];
}
}
template<class Type, class CombineOp>
void Foam::globalMeshData::syncPointData
(
List<Type>& pointData,
const CombineOp& cop
) const
{
const labelListList& slaves = globalPointSlaves();
const mapDistribute& map = globalPointSlavesMap();
syncPointData
(
pointData,
slaves,
map,
cop
);
}
template<class Type, class CombineOp>
void Foam::globalMeshData::syncPointAllData
(
List<Type>& pointData,
const CombineOp& cop
) const
{
syncPointData
(
pointData,
globalPointAllSlaves(),
globalPointAllSlavesMap(),
cop
);
}
// ************************************************************************* //

View File

@ -35,141 +35,50 @@ defineTypeNameAndDebug(Foam::globalPoints, 0);
const Foam::label Foam::globalPoints::fromCollocated = labelMax/2;
const Foam::scalar Foam::globalPoints::mergeDist = ROOTVSMALL;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Routines to handle global indices
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool Foam::globalPoints::noTransform(const tensor& tt, const scalar mergeDist)
{
return
(mag(tt.xx()-1) < mergeDist)
&& (mag(tt.yy()-1) < mergeDist)
&& (mag(tt.zz()-1) < mergeDist)
&& (mag(tt.xy()) < mergeDist)
&& (mag(tt.xz()) < mergeDist)
&& (mag(tt.yx()) < mergeDist)
&& (mag(tt.yz()) < mergeDist)
&& (mag(tt.zx()) < mergeDist)
&& (mag(tt.zy()) < mergeDist);
}
// Calculates per face whether couple is collocated.
Foam::PackedBoolList Foam::globalPoints::collocatedFaces
(
const coupledPolyPatch& pp,
const scalar mergeDist
)
{
// Initialise to false
PackedBoolList collocated(pp.size());
const vectorField& separation = pp.separation();
const tensorField& forwardT = pp.forwardT();
if (forwardT.size() == 0)
{
// Parallel.
if (separation.size() == 0)
{
collocated = 1u;
}
else if (separation.size() == 1)
{
// Fully separate. Do not synchronise.
}
else
{
// Per face separation.
forAll(pp, faceI)
{
if (mag(separation[faceI]) < mergeDist)
{
collocated[faceI] = 1u;
}
}
}
}
else if (forwardT.size() == 1)
{
// Fully transformed.
}
else
{
// Per face transformation.
forAll(pp, faceI)
{
if (noTransform(forwardT[faceI], mergeDist))
{
collocated[faceI] = 1u;
}
}
}
return collocated;
}
Foam::PackedBoolList Foam::globalPoints::collocatedPoints
(
const coupledPolyPatch& pp,
const scalar mergeDist
const coupledPolyPatch& pp
)
{
// Initialise to false
PackedBoolList collocated(pp.nPoints());
PackedBoolList isCollocated(pp.nPoints());
const vectorField& separation = pp.separation();
const tensorField& forwardT = pp.forwardT();
const boolList& collocated = pp.collocated();
if (forwardT.size() == 0)
if (collocated.size() == 0)
{
// Parallel.
if (separation.size() == 0)
{
collocated = 1u;
}
else if (separation.size() == 1)
{
// Fully separate.
}
else
{
// Per face separation.
for (label pointI = 0; pointI < pp.nPoints(); pointI++)
{
label faceI = pp.pointFaces()[pointI][0];
if (mag(separation[faceI]) < mergeDist)
{
collocated[pointI] = 1u;
}
}
}
isCollocated = 1;
}
else if (forwardT.size() == 1)
else if (collocated.size() == 1)
{
// Fully transformed.
// Uniform.
if (collocated[0])
{
isCollocated = 1;
}
}
else
{
// Per face transformation.
for (label pointI = 0; pointI < pp.nPoints(); pointI++)
{
label faceI = pp.pointFaces()[pointI][0];
// Per face collocated or not.
const labelListList& pointFaces = pp.pointFaces();
if (noTransform(forwardT[faceI], mergeDist))
forAll(pointFaces, pfi)
{
if (collocated[pointFaces[pfi][0]])
{
collocated[pointI] = 1u;
isCollocated[pfi] = 1;
}
}
}
return collocated;
return isCollocated;
}
Foam::label Foam::globalPoints::toGlobal
(
@ -467,8 +376,7 @@ void Foam::globalPoints::initOwnPoints
(
collocatedPoints
(
refCast<const coupledPolyPatch>(pp),
mergeDist
refCast<const coupledPolyPatch>(pp)
)
);
@ -563,8 +471,7 @@ void Foam::globalPoints::sendPatchPoints
(
collocatedPoints
(
procPatch,
mergeDist
procPatch
)
);
@ -663,8 +570,7 @@ void Foam::globalPoints::receivePatchPoints
(
collocatedPoints
(
procPatch,
mergeDist
procPatch
)
);
@ -726,8 +632,7 @@ void Foam::globalPoints::receivePatchPoints
(
collocatedPoints
(
cycPatch,
mergeDist
cycPatch
)
);
@ -1233,8 +1138,7 @@ void Foam::globalPoints::receiveSharedPoints
(
collocatedPoints
(
cycPatch,
mergeDist
cycPatch
)
);

View File

@ -127,9 +127,6 @@ class globalPoints
// collocated coupled points.
static const label fromCollocated;
//- Distance to check whether points/faces are collocated.
static const scalar mergeDist;
// Private data
@ -162,22 +159,9 @@ class globalPoints
// Private Member Functions
//- Is identity transform?
static bool noTransform(const tensor&, const scalar mergeDist);
//- Return per face collocated status
static PackedBoolList collocatedFaces
(
const coupledPolyPatch&,
const scalar mergeDist
);
//- Return per point collocated status
static PackedBoolList collocatedPoints
(
const coupledPolyPatch&,
const scalar mergeDist
);
static PackedBoolList collocatedPoints(const coupledPolyPatch&);
// Wrappers around global point numbering to add collocated bit

View File

@ -294,13 +294,13 @@ Foam::polyMesh::polyMesh(const IOobject& io)
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
// Warn if global empty mesh (constructs globalData!)
if (globalData().nTotalPoints() == 0)
// Warn if global empty mesh
if (returnReduce(nPoints(), sumOp<label>()) == 0)
{
WarningIn("polyMesh(const IOobject&)")
<< "no points in mesh" << endl;
}
if (globalData().nTotalCells() == 0)
if (returnReduce(nCells(), sumOp<label>()) == 0)
{
WarningIn("polyMesh(const IOobject&)")
<< "no cells in mesh" << endl;
@ -743,8 +743,12 @@ void Foam::polyMesh::resetPrimitives
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
// Warn if global empty mesh (constructs globalData!)
if (globalData().nTotalPoints() == 0 || globalData().nTotalCells() == 0)
// Warn if global empty mesh
if
(
(returnReduce(nPoints(), sumOp<label>()) == 0)
|| (returnReduce(nCells(), sumOp<label>()) == 0)
)
{
FatalErrorIn
(

View File

@ -285,6 +285,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
separation_.setSize(0);
forwardT_ = I;
reverseT_ = I;
collocated_.setSize(0);
}
else
{
@ -299,10 +300,14 @@ void Foam::coupledPolyPatch::calcTransformTensors
{
// Rotation, no separation
// Assume per-face differening transformation, correct later
separation_.setSize(0);
forwardT_.setSize(Cf.size());
reverseT_.setSize(Cf.size());
collocated_.setSize(Cf.size());
collocated_ = false;
forAll (forwardT_, facei)
{
@ -321,6 +326,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
{
forwardT_.setSize(1);
reverseT_.setSize(1);
collocated_.setSize(1);
if (debug)
{
@ -332,11 +338,15 @@ void Foam::coupledPolyPatch::calcTransformTensors
}
else
{
// No rotation, possible separation
forwardT_.setSize(0);
reverseT_.setSize(0);
separation_ = (nf&(Cr - Cf))*nf;
collocated_.setSize(separation_.size());
// Three situations:
// - separation is zero. No separation.
// - separation is same. Single separation vector.
@ -344,15 +354,23 @@ void Foam::coupledPolyPatch::calcTransformTensors
// Check for different separation per face
bool sameSeparation = true;
bool doneWarning = false;
forAll(separation_, facei)
{
scalar smallSqr = sqr(smallDist[facei]);
collocated_[facei] = (magSqr(separation_[facei]) < smallSqr);
// Check if separation differing w.r.t. face 0.
if (magSqr(separation_[facei] - separation_[0]) > smallSqr)
{
if (debug)
sameSeparation = false;
if (!doneWarning && debug)
{
doneWarning = true;
Pout<< " separation " << separation_[facei]
<< " at " << facei
<< " differs from separation[0] " << separation_[0]
@ -360,15 +378,13 @@ void Foam::coupledPolyPatch::calcTransformTensors
<< smallDist[facei]
<< ". Assuming non-uniform separation." << endl;
}
sameSeparation = false;
break;
}
}
if (sameSeparation)
{
// Check for zero separation (at 0 so everywhere)
if (magSqr(separation_[0]) < sqr(smallDist[0]))
if (collocated_[0])
{
if (debug)
{
@ -378,6 +394,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
}
separation_.setSize(0);
collocated_ = boolList(1, true);
}
else
{
@ -389,6 +406,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
}
separation_.setSize(1);
collocated_ = boolList(1, false);
}
}
}

View File

@ -64,6 +64,9 @@ class coupledPolyPatch
//- Neighbour-cell transformation tensor
mutable tensorField reverseT_;
//- Are faces collocated. Either size 0,1 or length of patch.
mutable boolList collocated_;
public:
// Static data members
@ -261,7 +264,6 @@ public:
return separation_;
}
//- Are the cyclic planes parallel
bool parallel() const
{
@ -280,6 +282,12 @@ public:
return reverseT_;
}
//- Are faces collocated. Either size 0,1 or length of patch
const boolList& collocated() const
{
return collocated_;
}
//- Initialize ordering for primitivePatch. Does not
// refer to *this (except for name() and type() etc.)

View File

@ -26,8 +26,8 @@ Class
Foam::genericPolyPatch
Description
Determines a mapping between patch face centres and mesh cell centres and
processors they're on.
Substitute for unknown patches. Used for postprocessing when only
basic polyPatch info is needed.
Note
Storage is not optimal. It stores all face centres and cells on all

View File

@ -74,6 +74,17 @@ public:
//- The various text representations for a switch value.
// These also correspond to the entries in names.
# undef FALSE
# undef TRUE
# undef OFF
# undef ON
# undef NO
# undef YES
# undef NO_1
# undef YES_1
# undef NONE
# undef PLACEHOLDER
# undef INVALID
enum switchType
{
FALSE = 0, TRUE = 1,

View File

@ -82,6 +82,8 @@ EqOp(orEq, x = (x || y))
EqOp(eqMinus, x = -y)
EqOp(nopEq, (void)x)
#undef EqOp

View File

@ -27,7 +27,6 @@ License
#include "motionSmoother.H"
#include "meshTools.H"
#include "processorPointPatchFields.H"
#include "globalPointPatchFields.H"
#include "pointConstraint.H"
#include "syncTools.H"

View File

@ -187,7 +187,7 @@ $(interpolation)/interpolationCellPointWallModified/cellPointWeightWallModified/
$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
volPointInterpolation = interpolation/volPointInterpolation
$(volPointInterpolation)/pointPatchInterpolation/pointPatchInterpolation.C
/* $(volPointInterpolation)/pointPatchInterpolation/pointPatchInterpolation.C */
$(volPointInterpolation)/volPointInterpolation.C
surfaceInterpolation = interpolation/surfaceInterpolation

View File

@ -1,208 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "pointPatchInterpolation.H"
#include "volFields.H"
#include "pointFields.H"
#include "emptyFvPatch.H"
#include "valuePointPatchField.H"
#include "coupledPointPatchField.H"
#include "coupledFacePointPatch.H"
#include "transform.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
void pointPatchInterpolation::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf,
bool overrideFixedValue
) const
{
if (debug)
{
Info<< "pointPatchInterpolation::interpolate("
<< "const GeometricField<Type, fvPatchField, volMesh>&, "
<< "GeometricField<Type, pointPatchField, pointMesh>&) : "
<< "interpolating field from cells to points"
<< endl;
}
// Interpolate patch values: over-ride the internal values for the points
// on the patch with the interpolated point values from the faces of the
// patch
const fvBoundaryMesh& bm = fvMesh_.boundary();
const pointBoundaryMesh& pbm = pointMesh::New(fvMesh_).boundary();
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
pointPatchField<Type>& ppf = pf.boundaryField()[patchi];
// Only map the values corresponding to the points associated with
// faces, not "lone" points due to decomposition
ppf.setInInternalField
(
pf.internalField(),
patchInterpolators_[patchi]
.faceToPointInterpolate(vf.boundaryField()[patchi])(),
bm[patchi].patch().meshPoints()
);
if
(
overrideFixedValue
&& isA<valuePointPatchField<Type> >(ppf)
)
{
refCast<valuePointPatchField<Type> >(ppf) = ppf;
}
}
else if (bm[patchi].coupled())
{
// Initialise the "lone" points on the coupled patch to zero,
// these values are obtained from the couple-transfer
const labelList& loneMeshPoints =
refCast<const coupledFacePointPatch>(pbm[patchi])
.loneMeshPoints();
forAll(loneMeshPoints, i)
{
pf[loneMeshPoints[i]] = pTraits<Type>::zero;
}
}
}
// Correct patch-patch boundary points by interpolation "around" corners
const labelListList& PointFaces = fvMesh_.pointFaces();
forAll(patchPatchPoints_, pointi)
{
const label curPoint = patchPatchPoints_[pointi];
const labelList& curFaces = PointFaces[curPoint];
label fI = 0;
// Reset the boundary value before accumulation
pf[curPoint] = pTraits<Type>::zero;
// Go through all the faces
forAll(curFaces, facei)
{
if (!fvMesh_.isInternalFace(curFaces[facei]))
{
label patchi =
fvMesh_.boundaryMesh().whichPatch(curFaces[facei]);
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
label faceInPatchi =
bm[patchi].patch().whichFace(curFaces[facei]);
pf[curPoint] +=
patchPatchPointWeights_[pointi][fI]
*vf.boundaryField()[patchi][faceInPatchi];
fI++;
}
}
}
}
// Update coupled boundaries
forAll(pf.boundaryField(), patchi)
{
if (pf.boundaryField()[patchi].coupled())
{
refCast<coupledPointPatchField<Type> >(pf.boundaryField()[patchi])
.initSwapAdd(pf.internalField());
}
}
forAll(pf.boundaryField(), patchi)
{
if (pf.boundaryField()[patchi].coupled())
{
refCast<coupledPointPatchField<Type> >(pf.boundaryField()[patchi])
.swapAdd(pf.internalField());
}
}
// Override constrained pointPatchField types with the constraint value.
// This relys on only constrained pointPatchField implementing the evaluate
// function
pf.correctBoundaryConditions();
// Apply multiple constraints on edge/corner points
applyCornerConstraints(pf);
if (debug)
{
Info<< "pointPatchInterpolation::interpolate("
<< "const GeometricField<Type, fvPatchField, volMesh>&, "
<< "GeometricField<Type, pointPatchField, pointMesh>&) : "
<< "finished interpolating field from cells to points"
<< endl;
}
}
template<class Type>
void pointPatchInterpolation::applyCornerConstraints
(
GeometricField<Type, pointPatchField, pointMesh>& pf
) const
{
forAll(patchPatchPointConstraintPoints_, pointi)
{
pf[patchPatchPointConstraintPoints_[pointi]] = transform
(
patchPatchPointConstraintTensors_[pointi],
pf[patchPatchPointConstraintPoints_[pointi]]
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,337 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "pointPatchInterpolation.H"
#include "fvMesh.H"
#include "volFields.H"
#include "pointFields.H"
#include "emptyFvPatch.H"
#include "demandDrivenData.H"
#include "coupledPointPatchFields.H"
#include "pointConstraint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(pointPatchInterpolation, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void pointPatchInterpolation::makePatchPatchAddressing()
{
if (debug)
{
Info<< "pointPatchInterpolation::makePatchPatchAddressing() : "
<< "constructing boundary addressing"
<< endl;
}
const fvBoundaryMesh& bm = fvMesh_.boundary();
const pointBoundaryMesh& pbm = pointMesh::New(fvMesh_).boundary();
// first count the total number of patch-patch points
label nPatchPatchPoints = 0;
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
nPatchPatchPoints += bm[patchi].patch().boundaryPoints().size();
}
}
// Go through all patches and mark up the external edge points
Map<label> patchPatchPointSet(2*nPatchPatchPoints);
patchPatchPoints_.setSize(nPatchPatchPoints);
List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
label pppi = 0;
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
const labelList& bp = bm[patchi].patch().boundaryPoints();
const labelList& meshPoints = bm[patchi].patch().meshPoints();
forAll(bp, pointi)
{
label ppp = meshPoints[bp[pointi]];
Map<label>::iterator iter = patchPatchPointSet.find(ppp);
if (iter == patchPatchPointSet.end())
{
patchPatchPointSet.insert(ppp, pppi);
patchPatchPoints_[pppi] = ppp;
pbm[patchi].applyConstraint
(
bp[pointi],
patchPatchPointConstraints[pppi]
);
pppi++;
}
else
{
pbm[patchi].applyConstraint
(
bp[pointi],
patchPatchPointConstraints[iter()]
);
}
}
}
}
nPatchPatchPoints = pppi;
patchPatchPoints_.setSize(nPatchPatchPoints);
patchPatchPointConstraints.setSize(nPatchPatchPoints);
patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
label nConstraints = 0;
forAll(patchPatchPointConstraints, i)
{
if (patchPatchPointConstraints[i].first() != 0)
{
patchPatchPointConstraintPoints_[nConstraints] =
patchPatchPoints_[i];
patchPatchPointConstraintTensors_[nConstraints] =
patchPatchPointConstraints[i].constraintTransformation();
nConstraints++;
}
}
patchPatchPointConstraintPoints_.setSize(nConstraints);
patchPatchPointConstraintTensors_.setSize(nConstraints);
patchInterpolators_.clear();
patchInterpolators_.setSize(bm.size());
forAll(bm, patchi)
{
patchInterpolators_.set
(
patchi,
new primitivePatchInterpolation(bm[patchi].patch())
);
}
if (debug)
{
Info<< "pointPatchInterpolation::makePatchPatchAddressing() : "
<< "finished constructing boundary addressing"
<< endl;
}
}
void pointPatchInterpolation::makePatchPatchWeights()
{
if (debug)
{
Info<< "pointPatchInterpolation::makePatchPatchWeights() : "
<< "constructing boundary weighting factors"
<< endl;
}
patchPatchPointWeights_.clear();
patchPatchPointWeights_.setSize(patchPatchPoints_.size());
const labelListList& pf = fvMesh_.pointFaces();
const volVectorField& centres = fvMesh_.C();
const fvBoundaryMesh& bm = fvMesh_.boundary();
pointScalarField sumWeights
(
IOobject
(
"sumWeights",
fvMesh_.polyMesh::instance(),
fvMesh_
),
pointMesh::New(fvMesh_),
dimensionedScalar("zero", dimless, 0)
);
forAll(patchPatchPoints_, pointi)
{
const label curPoint = patchPatchPoints_[pointi];
const labelList& curFaces = pf[curPoint];
patchPatchPointWeights_[pointi].setSize(curFaces.size());
scalarList& pw = patchPatchPointWeights_[pointi];
label nFacesAroundPoint = 0;
const vector& pointLoc = fvMesh_.points()[curPoint];
forAll(curFaces, facei)
{
if (!fvMesh_.isInternalFace(curFaces[facei]))
{
label patchi =
fvMesh_.boundaryMesh().whichPatch(curFaces[facei]);
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
vector d =
pointLoc
- centres.boundaryField()[patchi]
[bm[patchi].patch().whichFace(curFaces[facei])];
pw[nFacesAroundPoint] = 1.0/(mag(d)+VSMALL);
nFacesAroundPoint++;
}
}
}
// Reset the sizes of the local weights
pw.setSize(nFacesAroundPoint);
// Collect the sum of weights for parallel correction
sumWeights[curPoint] += sum(pw);
}
// Do parallel correction of weights
// Update coupled boundaries
forAll(sumWeights.boundaryField(), patchi)
{
if (sumWeights.boundaryField()[patchi].coupled())
{
refCast<coupledPointPatchScalarField>
(sumWeights.boundaryField()[patchi]).initSwapAdd
(
sumWeights.internalField()
);
}
}
forAll(sumWeights.boundaryField(), patchi)
{
if (sumWeights.boundaryField()[patchi].coupled())
{
refCast<coupledPointPatchScalarField>
(sumWeights.boundaryField()[patchi]).swapAdd
(
sumWeights.internalField()
);
}
}
// Re-scale the weights for the current point
forAll(patchPatchPoints_, pointi)
{
scalarList& pw = patchPatchPointWeights_[pointi];
scalar sumw = sumWeights[patchPatchPoints_[pointi]];
forAll(pw, facei)
{
pw[facei] /= sumw;
}
}
if (debug)
{
Info<< "pointPatchInterpolation::makePatchPatchWeights() : "
<< "finished constructing boundary weighting factors"
<< endl;
}
}
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
pointPatchInterpolation::pointPatchInterpolation(const fvMesh& vm)
:
fvMesh_(vm)
{
updateMesh();
}
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
pointPatchInterpolation::~pointPatchInterpolation()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void pointPatchInterpolation::updateMesh()
{
makePatchPatchAddressing();
makePatchPatchWeights();
}
bool pointPatchInterpolation::movePoints()
{
forAll(patchInterpolators_, patchi)
{
patchInterpolators_[patchi].movePoints();
}
makePatchPatchWeights();
return true;
}
// Specialisaion of applyCornerConstraints for scalars because
// no constraint need be applied
template<>
void pointPatchInterpolation::applyCornerConstraints<scalar>
(
GeometricField<scalar, pointPatchField, pointMesh>& pf
) const
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,169 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::pointPatchInterpolation
Description
Foam::pointPatchInterpolation
SourceFiles
pointPatchInterpolation.C
\*---------------------------------------------------------------------------*/
#ifndef pointPatchInterpolation_H
#define pointPatchInterpolation_H
#include "primitivePatchInterpolation.H"
#include "PtrList.H"
#include "volFieldsFwd.H"
#include "pointFieldsFwd.H"
#include "scalarList.H"
#include "className.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class fvMesh;
class pointMesh;
/*---------------------------------------------------------------------------*\
Class pointPatchInterpolation Declaration
\*---------------------------------------------------------------------------*/
class pointPatchInterpolation
{
// Private data
const fvMesh& fvMesh_;
//- Primitive patch interpolators
PtrList<primitivePatchInterpolation> patchInterpolators_;
//- List of patch-patch edge points that require special treatement
labelList patchPatchPoints_;
//- Weights for patch-patch boundary points
scalarListList patchPatchPointWeights_;
labelList patchPatchPointConstraintPoints_;
tensorField patchPatchPointConstraintTensors_;
// Private member functions
//- Construct addressing for patch-patch boundary points
void makePatchPatchAddressing();
//- Construct weights for patch-patch boundary points
void makePatchPatchWeights();
//- Disallow default bitwise copy construct
pointPatchInterpolation(const pointPatchInterpolation&);
//- Disallow default bitwise assignment
void operator=(const pointPatchInterpolation&);
public:
// Declare name of the class and its debug switch
ClassName("pointPatchInterpolation");
// Constructors
//- Constructor given fvMesh and pointMesh.
pointPatchInterpolation(const fvMesh&);
// Destructor
~pointPatchInterpolation();
// Member functions
// Access
const fvMesh& mesh() const
{
return fvMesh_;
}
// Edit
//- Update mesh topology using the morph engine
void updateMesh();
//- Correct weighting factors for moving mesh.
bool movePoints();
// Interpolation functions
template<class Type>
void interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, pointPatchField, pointMesh>&,
bool overrideFixedValue
) const;
template<class Type>
void applyCornerConstraints
(
GeometricField<Type, pointPatchField, pointMesh>& pf
) const;
};
template<>
void pointPatchInterpolation::applyCornerConstraints<scalar>
(
GeometricField<scalar, pointPatchField, pointMesh>& pf
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "pointPatchInterpolate.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -27,7 +27,11 @@ License
#include "volPointInterpolation.H"
#include "volFields.H"
#include "pointFields.H"
#include "globalPointPatch.H"
#include "emptyFvPatch.H"
#include "mapDistribute.H"
#include "coupledPointPatchField.H"
#include "valuePointPatchField.H"
#include "transform.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -36,6 +40,48 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
void volPointInterpolation::addSeparated
(
GeometricField<Type, pointPatchField, pointMesh>& pf
) const
{
if (debug)
{
Pout<< "volPointInterpolation::addSeparated" << endl;
}
forAll(pf.boundaryField(), patchI)
{
if (pf.boundaryField()[patchI].coupled())
{
refCast<coupledPointPatchField<Type> >
(pf.boundaryField()[patchI]).initSwapAddSeparated
(
Pstream::blocking, //Pstream::nonBlocking,
pf.internalField()
);
}
}
// Block for any outstanding requests
//Pstream::waitRequests();
forAll(pf.boundaryField(), patchI)
{
if (pf.boundaryField()[patchI].coupled())
{
refCast<coupledPointPatchField<Type> >
(pf.boundaryField()[patchI]).swapAddSeparated
(
Pstream::blocking, //Pstream::nonBlocking,
pf.internalField()
);
}
}
}
template<class Type>
void volPointInterpolation::interpolateInternalField
(
@ -45,7 +91,7 @@ void volPointInterpolation::interpolateInternalField
{
if (debug)
{
Info<< "volPointInterpolation::interpolateInternalField("
Pout<< "volPointInterpolation::interpolateInternalField("
<< "const GeometricField<Type, fvPatchField, volMesh>&, "
<< "GeometricField<Type, pointPatchField, pointMesh>&) : "
<< "interpolating field from cells to points"
@ -57,19 +103,166 @@ void volPointInterpolation::interpolateInternalField
// Multiply volField by weighting factor matrix to create pointField
forAll(pointCells, pointi)
{
const scalarList& pw = pointWeights_[pointi];
const labelList& ppc = pointCells[pointi];
pf[pointi] = pTraits<Type>::zero;
forAll(ppc, pointCelli)
if (!isPatchPoint_[pointi])
{
pf[pointi] += pw[pointCelli]*vf[ppc[pointCelli]];
const scalarList& pw = pointWeights_[pointi];
const labelList& ppc = pointCells[pointi];
pf[pointi] = pTraits<Type>::zero;
forAll(ppc, pointCelli)
{
pf[pointi] += pw[pointCelli]*vf[ppc[pointCelli]];
}
}
}
}
template<class Type>
tmp<Field<Type> > volPointInterpolation::flatBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
const fvMesh& mesh = vf.mesh();
const fvBoundaryMesh& bm = mesh.boundary();
tmp<Field<Type> > tboundaryVals
(
new Field<Type>(mesh.nFaces()-mesh.nInternalFaces())
);
Field<Type>& boundaryVals = tboundaryVals();
forAll(vf.boundaryField(), patchI)
{
label bFaceI = bm[patchI].patch().start() - mesh.nInternalFaces();
if (!isA<emptyFvPatch>(bm[patchI]) && !bm[patchI].coupled())
{
SubList<Type>
(
boundaryVals,
vf.boundaryField()[patchI].size(),
bFaceI
).assign(vf.boundaryField()[patchI]);
}
else
{
const polyPatch& pp = bm[patchI].patch();
forAll(pp, i)
{
boundaryVals[bFaceI++] = pTraits<Type>::zero;
}
}
}
return tboundaryVals;
}
template<class Type>
void volPointInterpolation::interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf,
const bool overrideFixedValue
) const
{
const primitivePatch& boundary = boundaryPtr_();
Field<Type>& pfi = pf.internalField();
// Get face data in flat list
tmp<Field<Type> > tboundaryVals(flatBoundaryField(vf));
const Field<Type>& boundaryVals = tboundaryVals();
// Do points on 'normal' patches from the surrounding patch faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(boundary.meshPoints(), i)
{
label pointI = boundary.meshPoints()[i];
if (isPatchPoint_[pointI])
{
const labelList& pFaces = boundary.pointFaces()[i];
const scalarList& pWeights = boundaryPointWeights_[i];
Type& val = pfi[pointI];
val = pTraits<Type>::zero;
forAll(pFaces, j)
{
if (boundaryIsPatchFace_[pFaces[j]])
{
val += pWeights[j]*boundaryVals[pFaces[j]];
}
}
}
}
// Sum collocated contributions
mesh().globalData().syncPointData(pfi, plusEqOp<Type>());
// And add separated contributions
addSeparated(pf);
// Push master data to slaves. It is possible (not sure how often) for
// a coupled point to have its master on a different patch so
// to make sure just push master data to slaves. Reuse the syncPointData
// structure.
mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
if (overrideFixedValue)
{
forAll(pf.boundaryField(), patchI)
{
pointPatchField<Type>& ppf = pf.boundaryField()[patchI];
if (isA<valuePointPatchField<Type> >(ppf))
{
refCast<valuePointPatchField<Type> >(ppf) =
ppf.patchInternalField();
}
}
}
// Override constrained pointPatchField types with the constraint value.
// This relys on only constrained pointPatchField implementing the evaluate
// function
pf.correctBoundaryConditions();
// Sync any dangling points
vf.mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
// Apply multiple constraints on edge/corner points
applyCornerConstraints(pf);
}
template<class Type>
void volPointInterpolation::applyCornerConstraints
(
GeometricField<Type, pointPatchField, pointMesh>& pf
) const
{
forAll(patchPatchPointConstraintPoints_, pointi)
{
pf[patchPatchPointConstraintPoints_[pointi]] = transform
(
patchPatchPointConstraintTensors_[pointi],
pf[patchPatchPointConstraintPoints_[pointi]]
);
}
}
template<class Type>
void volPointInterpolation::interpolate
(
@ -79,7 +272,7 @@ void volPointInterpolation::interpolate
{
if (debug)
{
Info<< "volPointInterpolation::interpolate("
Pout<< "volPointInterpolation::interpolate("
<< "const GeometricField<Type, fvPatchField, volMesh>&, "
<< "GeometricField<Type, pointPatchField, pointMesh>&) : "
<< "interpolating field from cells to points"
@ -89,7 +282,7 @@ void volPointInterpolation::interpolate
interpolateInternalField(vf, pf);
// Interpolate to the patches preserving fixed value BCs
boundaryInterpolator_.interpolate(vf, pf, false);
interpolateBoundaryField(vf, pf, false);
}
@ -101,23 +294,7 @@ volPointInterpolation::interpolate
const wordList& patchFieldTypes
) const
{
wordList types(patchFieldTypes);
const pointMesh& pMesh = pointMesh::New(vf.mesh());
// If the last patch of the pointBoundaryMesh is the global patch
// it must be added to the list of patchField types
if
(
isType<globalPointPatch>
(
pMesh.boundary()[pMesh.boundary().size() - 1]
)
)
{
types.setSize(types.size() + 1);
types[types.size()-1] = pMesh.boundary()[types.size()-1].type();
}
const pointMesh& pm = pointMesh::New(vf.mesh());
// Construct tmp<pointField>
tmp<GeometricField<Type, pointPatchField, pointMesh> > tpf
@ -128,18 +305,18 @@ volPointInterpolation::interpolate
(
"volPointInterpolate(" + vf.name() + ')',
vf.instance(),
pMesh.thisDb()
pm.thisDb()
),
pMesh,
pm,
vf.dimensions(),
types
patchFieldTypes
)
);
interpolateInternalField(vf, tpf());
// Interpolate to the patches overriding fixed value BCs
boundaryInterpolator_.interpolate(vf, tpf(), true);
interpolateBoundaryField(vf, tpf(), true);
return tpf;
}
@ -186,7 +363,7 @@ volPointInterpolation::interpolate
);
interpolateInternalField(vf, tpf());
boundaryInterpolator_.interpolate(vf, tpf(), false);
interpolateBoundaryField(vf, tpf(), false);
return tpf;
}

View File

@ -44,15 +44,114 @@ defineTypeNameAndDebug(volPointInterpolation, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void volPointInterpolation::makeWeights()
void volPointInterpolation::calcBoundaryAddressing()
{
if (debug)
{
Info<< "volPointInterpolation::makeWeights() : "
<< "constructing weighting factors"
Pout<< "volPointInterpolation::calcBoundaryAddressing() : "
<< "constructing boundary addressing"
<< endl;
}
boundaryPtr_.reset
(
new primitivePatch
(
SubList<face>
(
mesh().faces(),
mesh().nFaces()-mesh().nInternalFaces(),
mesh().nInternalFaces()
),
mesh().points()
)
);
const primitivePatch& boundary = boundaryPtr_();
boundaryIsPatchFace_.setSize(boundary.size());
boundaryIsPatchFace_ = false;
isPatchPoint_.setSize(mesh().nPoints());
isPatchPoint_ = false;
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (!isA<emptyPolyPatch>(pp) && !pp.coupled())
{
label bFaceI = pp.start()-mesh().nInternalFaces();
forAll(pp, i)
{
boundaryIsPatchFace_[bFaceI] = true;
const face& f = boundary[bFaceI++];
forAll(f, fp)
{
isPatchPoint_[f[fp]] = true;
}
}
}
}
// Make sure point status is synchronised so even processor that holds
// no face of a certain patch still can have boundary points marked.
if (debug)
{
boolList oldData(isPatchPoint_);
mesh().globalData().syncPointData(isPatchPoint_, orEqOp<bool>());
forAll(isPatchPoint_, pointI)
{
if (isPatchPoint_[pointI] != oldData[pointI])
{
Pout<< "volPointInterpolation::calcBoundaryAddressing():"
<< " added dangling mesh point:" << pointI
<< " at:" << mesh().points()[pointI]
<< endl;
}
}
label nPatchFace = 0;
forAll(boundaryIsPatchFace_, i)
{
if (boundaryIsPatchFace_[i])
{
nPatchFace++;
}
}
label nPatchPoint = 0;
forAll(isPatchPoint_, i)
{
if (isPatchPoint_[i])
{
nPatchPoint++;
}
}
Pout<< "boundary:" << nl
<< " faces :" << boundary.size() << nl
<< " of which on proper patch:" << nPatchFace << nl
<< " points:" << boundary.nPoints() << nl
<< " of which on proper patch:" << nPatchPoint << endl;
}
}
void volPointInterpolation::makeInternalWeights(scalarField& sumWeights)
{
if (debug)
{
Pout<< "volPointInterpolation::makeInternalWeights() : "
<< "constructing weighting factors for internal and non-coupled"
<< " points." << endl;
}
const pointField& points = mesh().points();
const labelListList& pointCells = mesh().pointCells();
const vectorField& cellCentres = mesh().cellCentres();
@ -61,11 +160,91 @@ void volPointInterpolation::makeWeights()
pointWeights_.clear();
pointWeights_.setSize(points.size());
forAll(pointWeights_, pointi)
// Calculate inverse distances between cell centres and points
// and store in weighting factor array
forAll(points, pointi)
{
pointWeights_[pointi].setSize(pointCells[pointi].size());
if (!isPatchPoint_[pointi])
{
const labelList& pcp = pointCells[pointi];
scalarList& pw = pointWeights_[pointi];
pw.setSize(pcp.size());
forAll(pcp, pointCelli)
{
pw[pointCelli] =
1.0/mag(points[pointi] - cellCentres[pcp[pointCelli]]);
sumWeights[pointi] += pw[pointCelli];
}
}
}
}
void volPointInterpolation::makeBoundaryWeights(scalarField& sumWeights)
{
if (debug)
{
Pout<< "volPointInterpolation::makeBoundaryWeights() : "
<< "constructing weighting factors for boundary points." << endl;
}
const pointField& points = mesh().points();
const pointField& faceCentres = mesh().faceCentres();
const primitivePatch& boundary = boundaryPtr_();
boundaryPointWeights_.clear();
boundaryPointWeights_.setSize(boundary.meshPoints().size());
forAll(boundary.meshPoints(), i)
{
label pointI = boundary.meshPoints()[i];
if (isPatchPoint_[pointI])
{
const labelList& pFaces = boundary.pointFaces()[i];
scalarList& pw = boundaryPointWeights_[i];
pw.setSize(pFaces.size());
sumWeights[pointI] = 0.0;
forAll(pFaces, i)
{
if (boundaryIsPatchFace_[pFaces[i]])
{
label faceI = mesh().nInternalFaces() + pFaces[i];
pw[i] = 1.0/mag(points[pointI] - faceCentres[faceI]);
sumWeights[pointI] += pw[i];
}
else
{
pw[i] = 0.0;
}
}
}
}
}
void volPointInterpolation::makeWeights()
{
if (debug)
{
Pout<< "volPointInterpolation::makeWeights() : "
<< "constructing weighting factors"
<< endl;
}
// Update addressing over all boundary faces
calcBoundaryAddressing();
// Running sum of weights
pointScalarField sumWeights
(
IOobject
@ -78,60 +257,344 @@ void volPointInterpolation::makeWeights()
dimensionedScalar("zero", dimless, 0)
);
// Calculate inverse distances between cell centres and points
// and store in weighting factor array
forAll(points, pointi)
// Create internal weights; add to sumWeights
makeInternalWeights(sumWeights);
// Create boundary weights; add to sumWeights
makeBoundaryWeights(sumWeights);
//forAll(boundary.meshPoints(), i)
//{
// label pointI = boundary.meshPoints()[i];
//
// if (isPatchPoint_[pointI])
// {
// Pout<< "Calculated Weight at boundary point:" << i
// << " at:" << mesh().points()[pointI]
// << " sumWeight:" << sumWeights[pointI]
// << " from:" << boundaryPointWeights_[i]
// << endl;
// }
//}
// Sum collocated contributions
mesh().globalData().syncPointData(sumWeights, plusEqOp<scalar>());
// And add separated contributions
addSeparated(sumWeights);
// Push master data to slaves. It is possible (not sure how often) for
// a coupled point to have its master on a different patch so
// to make sure just push master data to slaves. Reuse the syncPointData
// structure.
mesh().globalData().syncPointData(sumWeights, nopEqOp<scalar>());
// Normalise internal weights
forAll(pointWeights_, pointI)
{
scalarList& pw = pointWeights_[pointi];
const labelList& pcp = pointCells[pointi];
forAll(pcp, pointCelli)
scalarList& pw = pointWeights_[pointI];
// Note:pw only sized for !isPatchPoint
forAll(pw, i)
{
pw[pointCelli] =
1.0/mag(points[pointi] - cellCentres[pcp[pointCelli]]);
sumWeights[pointi] += pw[pointCelli];
pw[i] /= sumWeights[pointI];
}
}
forAll(sumWeights.boundaryField(), patchi)
// Normalise boundary weights
const primitivePatch& boundary = boundaryPtr_();
forAll(boundary.meshPoints(), i)
{
if (sumWeights.boundaryField()[patchi].coupled())
label pointI = boundary.meshPoints()[i];
scalarList& pw = boundaryPointWeights_[i];
// Note:pw only sized for isPatchPoint
forAll(pw, i)
{
refCast<coupledPointPatchScalarField>
(sumWeights.boundaryField()[patchi]).initSwapAdd
(
sumWeights.internalField()
);
pw[i] /= sumWeights[pointI];
}
}
forAll(sumWeights.boundaryField(), patchi)
if (debug)
{
if (sumWeights.boundaryField()[patchi].coupled())
{
refCast<coupledPointPatchScalarField>
(sumWeights.boundaryField()[patchi]).swapAdd
(
sumWeights.internalField()
);
}
Pout<< "volPointInterpolation::makeWeights() : "
<< "finished constructing weighting factors"
<< endl;
}
}
void volPointInterpolation::makePatchPatchAddressing()
{
if (debug)
{
Pout<< "volPointInterpolation::makePatchPatchAddressing() : "
<< "constructing boundary addressing"
<< endl;
}
forAll(points, pointi)
{
scalarList& pw = pointWeights_[pointi];
const fvBoundaryMesh& bm = mesh().boundary();
const pointBoundaryMesh& pbm = pointMesh::New(mesh()).boundary();
forAll(pw, pointCelli)
// first count the total number of patch-patch points
label nPatchPatchPoints = 0;
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
pw[pointCelli] /= sumWeights[pointi];
nPatchPatchPoints += bm[patchi].patch().boundaryPoints().size();
if (debug)
{
Pout<< "On patch:" << bm[patchi].patch()
<< " nBoundaryPoints:"
<< bm[patchi].patch().boundaryPoints().size() << endl;
}
}
}
if (debug)
{
Info<< "volPointInterpolation::makeWeights() : "
<< "finished constructing weighting factors"
Pout<< "Found nPatchPatchPoints:" << nPatchPatchPoints << endl;
}
// Go through all patches and mark up the external edge points
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// From meshpoint to index in patchPatchPointConstraints.
Map<label> patchPatchPointSet(2*nPatchPatchPoints);
// Constraints (initialised to unconstrained)
List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
// From constraint index to mesh point
labelList patchPatchPoints(nPatchPatchPoints);
label pppi = 0;
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
const labelList& bp = bm[patchi].patch().boundaryPoints();
const labelList& meshPoints = bm[patchi].patch().meshPoints();
forAll(bp, pointi)
{
label ppp = meshPoints[bp[pointi]];
Map<label>::iterator iter = patchPatchPointSet.find(ppp);
label constraintI = -1;
if (iter == patchPatchPointSet.end())
{
patchPatchPointSet.insert(ppp, pppi);
patchPatchPoints[pppi] = ppp;
constraintI = pppi++;
}
else
{
constraintI = iter();
}
// Apply to patch constraints
pbm[patchi].applyConstraint
(
bp[pointi],
patchPatchPointConstraints[constraintI]
);
}
}
}
if (debug)
{
Pout<< "Have (local) constrained points:"
<< nPatchPatchPoints << endl;
}
// Extend set with constraints across coupled points
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
const globalMeshData& gd = mesh().globalData();
const labelListList& globalPointSlaves = gd.globalPointSlaves();
const mapDistribute& globalPointSlavesMap = gd.globalPointSlavesMap();
const Map<label>& cpPointMap = gd.coupledPatch().meshPointMap();
const labelList& cpMeshPoints = gd.coupledPatch().meshPoints();
// Constraints on coupled points
List<pointConstraint> constraints
(
globalPointSlavesMap.constructSize()
);
// Copy from patchPatch constraints into coupledConstraints.
forAll(bm, patchi)
{
if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
{
const labelList& bp = bm[patchi].patch().boundaryPoints();
const labelList& meshPoints = bm[patchi].patch().meshPoints();
forAll(bp, pointi)
{
label ppp = meshPoints[bp[pointi]];
Map<label>::const_iterator fnd = cpPointMap.find(ppp);
if (fnd != cpPointMap.end())
{
// Can just copy (instead of apply) constraint
// will already be consistent across multiple patches.
constraints[fnd()] = patchPatchPointConstraints
[
patchPatchPointSet[ppp]
];
}
}
}
}
// Exchange data
globalPointSlavesMap.distribute(constraints);
// Combine master with slave constraints
forAll(globalPointSlaves, pointI)
{
const labelList& slaves = globalPointSlaves[pointI];
// Combine master constraint with slave constraints
forAll(slaves, i)
{
constraints[pointI].combine(constraints[slaves[i]]);
}
// Duplicate master constraint into slave slots
forAll(slaves, i)
{
constraints[slaves[i]] = constraints[pointI];
}
}
// Send back
globalPointSlavesMap.reverseDistribute
(
cpMeshPoints.size(),
constraints
);
// Add back into patchPatch constraints
forAll(constraints, coupledPointI)
{
if (constraints[coupledPointI].first() != 0)
{
label meshPointI = cpMeshPoints[coupledPointI];
Map<label>::iterator iter = patchPatchPointSet.find(meshPointI);
label constraintI = -1;
if (iter == patchPatchPointSet.end())
{
//Pout<< "on meshpoint:" << meshPointI
// << " coupled:" << coupledPointI
// << " at:" << mesh().points()[meshPointI]
// << " have new constraint:"
// << constraints[coupledPointI]
// << endl;
// Allocate new constraint
if (patchPatchPoints.size() <= pppi)
{
patchPatchPoints.setSize(pppi+100);
}
patchPatchPointSet.insert(meshPointI, pppi);
patchPatchPoints[pppi] = meshPointI;
constraintI = pppi++;
}
else
{
//Pout<< "on meshpoint:" << meshPointI
// << " coupled:" << coupledPointI
// << " at:" << mesh().points()[meshPointI]
// << " have possibly extended constraint:"
// << constraints[coupledPointI]
// << endl;
constraintI = iter();
}
// Combine (new or existing) constraint with one
// on coupled.
patchPatchPointConstraints[constraintI].combine
(
constraints[coupledPointI]
);
}
}
}
nPatchPatchPoints = pppi;
patchPatchPoints.setSize(nPatchPatchPoints);
patchPatchPointConstraints.setSize(nPatchPatchPoints);
if (debug)
{
Pout<< "Have (global) constrained points:"
<< nPatchPatchPoints << endl;
}
// Copy out all non-trivial constraints
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
label nConstraints = 0;
forAll(patchPatchPointConstraints, i)
{
if (patchPatchPointConstraints[i].first() != 0)
{
patchPatchPointConstraintPoints_[nConstraints] =
patchPatchPoints[i];
patchPatchPointConstraintTensors_[nConstraints] =
patchPatchPointConstraints[i].constraintTransformation();
nConstraints++;
}
}
if (debug)
{
Pout<< "Have non-trivial constrained points:"
<< nConstraints << endl;
}
patchPatchPointConstraintPoints_.setSize(nConstraints);
patchPatchPointConstraintTensors_.setSize(nConstraints);
if (debug)
{
Pout<< "volPointInterpolation::makePatchPatchAddressing() : "
<< "finished constructing boundary addressing"
<< endl;
}
}
@ -141,8 +604,7 @@ void volPointInterpolation::makeWeights()
volPointInterpolation::volPointInterpolation(const fvMesh& vm)
:
MeshObject<fvMesh, volPointInterpolation>(vm),
boundaryInterpolator_(vm)
MeshObject<fvMesh, volPointInterpolation>(vm)
{
updateMesh();
}
@ -159,19 +621,27 @@ volPointInterpolation::~volPointInterpolation()
void volPointInterpolation::updateMesh()
{
makeWeights();
boundaryInterpolator_.updateMesh();
}
bool volPointInterpolation::movePoints()
{
makeWeights();
boundaryInterpolator_.movePoints();
return true;
}
// Specialisaion of applyCornerConstraints for scalars because
// no constraint need be applied
template<>
void volPointInterpolation::applyCornerConstraints<scalar>
(
GeometricField<scalar, pointPatchField, pointMesh>& pf
) const
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -38,7 +38,9 @@ SourceFiles
#define volPointInterpolation_H
#include "MeshObject.H"
#include "pointPatchInterpolation.H"
#include "scalarList.H"
#include "volFields.H"
#include "pointFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -58,18 +60,80 @@ class volPointInterpolation
{
// Private data
//- Boundary interpolation engine.
pointPatchInterpolation boundaryInterpolator_;
//- Interpolation scheme weighting factor array.
scalarListList pointWeights_;
// Boundary handling
//- Boundary addressing
autoPtr<primitivePatch> boundaryPtr_;
//- Per boundary face whether is on non-coupled patch
boolList boundaryIsPatchFace_;
//- Per mesh(!) point whether is on non-coupled patch (on any
// processor)
boolList isPatchPoint_;
//- Per boundary point the weights per pointFaces.
scalarListList boundaryPointWeights_;
// Patch-patch constraints
//- Mesh points on which to apply special constraints
labelList patchPatchPointConstraintPoints_;
//- Special constraints
tensorField patchPatchPointConstraintTensors_;
// Private member functions
//- Construct point weighting factors
//- Construct addressing over all boundary faces
void calcBoundaryAddressing();
//- Make weights for internal and coupled-only boundarypoints
void makeInternalWeights(scalarField& sumWeights);
//- Make weights for points on uncoupled patches
void makeBoundaryWeights(scalarField& sumWeights);
//- Construct all point weighting factors
void makeWeights();
//- Make patch-patch constraints
void makePatchPatchAddressing();
//- Get boundary field in same order as boundary faces. Field is
// zero on all coupled and empty patches
template<class Type>
tmp<Field<Type> > flatBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const;
template<class Type>
void interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf,
const bool overrideFixedValue
) const;
template<class Type>
void applyCornerConstraints
(
GeometricField<Type, pointPatchField, pointMesh>& pf
) const;
//- Add separated contributions
template<class Type>
void addSeparated
(
GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Disallow default bitwise copy construct
volPointInterpolation(const volPointInterpolation&);
@ -96,14 +160,6 @@ public:
// Member functions
// Access
const fvMesh& mesh() const
{
return boundaryInterpolator_.mesh();
}
// Edit
//- Update mesh topology using the morph engine
@ -169,6 +225,13 @@ public:
};
template<>
void volPointInterpolation::applyCornerConstraints<scalar>
(
GeometricField<scalar, pointPatchField, pointMesh>& pf
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -44,7 +44,6 @@ License
#include "pointFields.H"
#include "slipPointPatchFields.H"
#include "fixedValuePointPatchFields.H"
#include "globalPointPatchFields.H"
#include "calculatedPointPatchFields.H"
#include "processorPointPatch.H"
#include "globalIndex.H"
@ -1409,11 +1408,7 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
forAll(pointPatches, patchI)
{
if (isA<globalPointPatch>(pointPatches[patchI]))
{
patchFieldTypes[patchI] = globalPointPatchVectorField::typeName;
}
else if (isA<processorPointPatch>(pointPatches[patchI]))
if (isA<processorPointPatch>(pointPatches[patchI]))
{
patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
}