Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2013-05-03 15:47:52 +01:00
259 changed files with 13699 additions and 3213 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -0,0 +1,3 @@
laplacianFoam.C
EXE = $(FOAM_USER_APPBIN)/laplacianFoam-communicators

View File

@ -0,0 +1,3 @@
EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = -lfiniteVolume

View File

@ -0,0 +1,37 @@
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT
(
transportProperties.lookup("DT")
);

View File

@ -0,0 +1,817 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Application
laplacianFoam
Description
Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "simpleControl.H"
#include "globalIndex.H"
#include "lduPrimitiveMesh.H"
#include "processorGAMGInterface.H"
#include "GAMGInterfaceField.H"
#include "processorLduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void checkUpperTriangular
(
const label size,
const labelUList& l,
const labelUList& u
)
{
forAll(l, faceI)
{
if (u[faceI] < l[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Reversed face. Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
}
if (l[faceI] < 0 || u[faceI] < 0 || u[faceI] >= size)
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Illegal cell label. Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
}
}
for (label faceI=1; faceI < l.size(); faceI++)
{
if (l[faceI-1] > l[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Lower not in incremental cell order."
<< " Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI]
<< " previous l:" << l[faceI-1] << abort(FatalError);
}
else if (l[faceI-1] == l[faceI])
{
// Same cell.
if (u[faceI-1] > u[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Upper not in incremental cell order."
<< " Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI]
<< " previous u:" << u[faceI-1] << abort(FatalError);
}
}
}
}
void print(const string& msg, const lduMesh& mesh)
{
const lduAddressing& addressing = mesh.lduAddr();
const lduInterfacePtrsList interfaces = mesh.interfaces();
Pout<< "Mesh:" << msg.c_str() << nl
<< " cells:" << addressing.size() << nl
<< " faces:" << addressing.lowerAddr().size() << nl
<< " patches:" << interfaces.size() << nl;
const labelUList& l = addressing.lowerAddr();
const labelUList& startAddr = addressing.losortStartAddr();
const labelUList& addr = addressing.losortAddr();
forAll(addressing, cellI)
{
Pout<< " cell:" << cellI << nl;
label start = startAddr[cellI];
label end = startAddr[cellI+1];
for (label index = start; index < end; index++)
{
Pout<< " nbr:" << l[addr[index]] << nl;
}
}
Pout<< " Patches:" << nl;
forAll(interfaces, i)
{
if (interfaces.set(i))
{
if (isA<processorLduInterface>(interfaces[i]))
{
const processorLduInterface& pldui =
refCast<const processorLduInterface>(interfaces[i]);
Pout<< " " << i
<< " me:" << pldui.myProcNo()
<< " nbr:" << pldui.neighbProcNo()
<< " comm:" << pldui.comm()
<< " tag:" << pldui.tag()
<< nl;
}
{
Pout<< " " << i << " addressing:" << nl;
const labelUList& faceCells = interfaces[i].faceCells();
forAll(faceCells, i)
{
Pout<< "\t\t" << i << '\t' << faceCells[i] << nl;
}
}
}
}
}
template<class ProcPatch>
lduSchedule nonBlockingSchedule
(
const lduInterfacePtrsList& interfaces
)
{
lduSchedule schedule(2*interfaces.size());
label slotI = 0;
forAll(interfaces, i)
{
if (interfaces.set(i) && !isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = true;
slotI++;
schedule[slotI].patch = i;
schedule[slotI].init = false;
slotI++;
}
}
forAll(interfaces, i)
{
if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = true;
slotI++;
}
}
forAll(interfaces, i)
{
if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = false;
slotI++;
}
}
return schedule;
}
void sendReceive
(
const label comm,
const label tag,
const globalIndex& offsets,
const scalarField& field,
scalarField& allField
)
{
label nProcs = Pstream::nProcs(comm);
if (Pstream::master(comm))
{
allField.setSize(offsets.size());
// Assign master slot
SubList<scalar>
(
allField,
offsets.localSize(0),
offsets.offset(0)
).assign(field);
// Assign slave slots
for (label procI = 1; procI < nProcs; procI++)
{
SubList<scalar> procSlot
(
allField,
offsets.localSize(procI),
offsets.offset(procI)
);
Pout<< "Receiving allField from " << procI
<< " at offset:" << offsets.offset(procI)
<< " size:" << offsets.size()
<< endl;
IPstream::read
(
Pstream::nonBlocking,
procI,
reinterpret_cast<char*>(procSlot.begin()),
procSlot.byteSize(),
tag,
comm
);
}
}
else
{
OPstream::write
(
Pstream::nonBlocking,
0, // master
reinterpret_cast<const char*>(field.begin()),
field.byteSize(),
tag,
comm
);
}
}
void sendReceive
(
const label comm,
const label tag,
const globalIndex& offsets,
const FieldField<Field, scalar>& field,
FieldField<Field, scalar>& allField
)
{
PstreamBuffers pBufs(Pstream::nonBlocking, Pstream::msgType(), comm);
if (!Pstream::master(comm))
{
UOPstream toMaster(Pstream::masterNo(), pBufs);
Pout<< "To 0 sending " << field.size()
<< " fields." << endl;
forAll(field, intI)
{
toMaster << field[intI];
}
}
pBufs.finishedSends();
if (Pstream::master(comm))
{
allField.setSize(offsets.size());
forAll(allField, i)
{
allField.set(i, new scalarField(0));
}
// Insert master values
forAll(field, intI)
{
allField[intI] = field[intI];
}
// Receive and insert slave values
label nProcs = Pstream::nProcs(comm);
for (label procI = 1; procI < nProcs; procI++)
{
UIPstream fromSlave(procI, pBufs);
label nSlaveInts = offsets.localSize(procI);
Pout<< "From " << procI << " receiving "
<< nSlaveInts << " fields." << endl;
for (label intI = 0; intI < nSlaveInts; intI++)
{
label slotI = offsets.toGlobal(procI, intI);
Pout<< " int:" << intI << " goes into slot " << slotI
<< endl;
fromSlave >> allField[slotI];
}
}
}
}
void collect
(
const label comm,
const globalIndex& cellOffsets,
const globalIndex& faceOffsets,
const scalarField& diagonal,
const scalarField& upper,
const scalarField& lower,
scalarField& allDiagonal,
scalarField& allUpper,
scalarField& allLower
)
{
label nOutstanding = Pstream::nRequests();
int allDiagonalTag = Pstream::allocateTag("allDiagonal:" __FILE__);
int allUpperTag = Pstream::allocateTag("allUpper:" __FILE__);
int allLowerTag = Pstream::allocateTag("allLower:" __FILE__);
sendReceive
(
comm,
allDiagonalTag,
cellOffsets,
diagonal,
allDiagonal
);
sendReceive
(
comm,
allUpperTag,
faceOffsets,
upper,
allUpper
);
sendReceive
(
comm,
allLowerTag,
faceOffsets,
lower,
allLower
);
Pstream::waitRequests(nOutstanding);
Pstream::freeTag("allDiagonal:" __FILE__, allDiagonalTag);
Pstream::freeTag("allUpper:" __FILE__, allUpperTag);
Pstream::freeTag("allLower:" __FILE__, allLowerTag);
}
void setCommunicator(fvMesh& mesh, const label newComm)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// The current state is consistent with the mesh so check where the new
// communicator is and adjust accordingly.
forAll(pbm, patchI)
{
if (isA<processorPolyPatch>(pbm[patchI]))
{
processorPolyPatch& ppp = const_cast<processorPolyPatch&>
(
refCast
<
const processorPolyPatch
>(pbm[patchI])
);
label thisRank = UPstream::procNo
(
newComm,
ppp.comm(),
ppp.myProcNo()
);
label nbrRank = UPstream::procNo
(
newComm,
ppp.comm(),
ppp.neighbProcNo()
);
//ppp.comm() = newComm;
ppp.myProcNo() = thisRank;
ppp.neighbProcNo() = nbrRank;
}
}
mesh.polyMesh::comm() = newComm;
}
namespace Foam
{
typedef UPtrList<const GAMGInterfaceField> GAMGInterfaceFieldPtrsList;
}
// Gather matrices from processors procIDs[1..] on procIDs[0]
void gatherMatrices
(
const labelList& procIDs,
const PtrList<lduMesh>& procMeshes,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
PtrList<GAMGInterfaceFieldPtrsList>& otherInterfaces
)
{
const label meshComm = mat.mesh().comm();
//lduInterfacePtrsList interfaces(mesh.interfaces());
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Master.
otherMats.setSize(procIDs.size()-1);
otherBouCoeffs.setSize(procIDs.size()-1);
otherIntCoeffs.setSize(procIDs.size()-1);
otherInterfaces.setSize(procIDs.size()-1);
for (label i = 1; i < procIDs.size(); i++)
{
const lduMesh& procMesh = procMeshes[i];
const lduInterfacePtrsList procInterfaces = procMesh.interfaces();
IPstream fromSlave
(
Pstream::scheduled,
procIDs[i],
0, // bufSize
Pstream::msgType(),
meshComm
);
otherMats.set(i-1, new lduMatrix(procMesh, fromSlave));
// Receive number of/valid interfaces
boolList validTransforms(fromSlave);
List<int> validRanks(fromSlave);
// Size coefficients
otherBouCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherIntCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherInterfaces.set
(
i-1,
new GAMGInterfaceFieldPtrsList(validTransforms.size())
);
forAll(validTransforms, intI)
{
if (validTransforms[intI])
{
const processorGAMGInterface& interface =
refCast<const processorGAMGInterface>
(
procInterfaces[intI]
);
otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave));
otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave));
otherInterfaces[i-1].set
(
intI,
GAMGInterfaceField::New
(
interface, //procInterfaces[intI],
validTransforms[intI],
validRanks[intI]
).ptr()
);
}
}
}
}
else
{
// Send to master
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
Pstream::msgType(),
meshComm
);
// Count valid interfaces
boolList validTransforms(interfaceBouCoeffs.size(), false);
List<int> validRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
interfaces[intI]
);
validTransforms[intI] = interface.doTransform();
validRanks[intI] = interface.rank();
}
}
toMaster << mat << validTransforms << validRanks;
forAll(validTransforms, intI)
{
if (validTransforms[intI])
{
toMaster
<< interfaceBouCoeffs[intI]
<< interfaceIntCoeffs[intI];
}
}
}
}
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
simpleControl simple(mesh);
//const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating temperature distribution\n" << endl;
// Get a subset of processors
labelList subProcs(3);
subProcs[0] = 0;
subProcs[1] = 1;
subProcs[2] = 2;
const UPstream::communicator newComm
(
UPstream::worldComm,
subProcs,
true
);
Pout<< "procIDs world :" << UPstream::procID(UPstream::worldComm) << endl;
Pout<< "procIDs newComm:" << UPstream::procID(newComm) << endl;
//// On ALL processors: get the interfaces:
//lduInterfacePtrsList interfaces(mesh.interfaces());
//PtrList<lduMesh> procMeshes;
//
//if (Pstream::myProcNo(newComm) != -1)
//{
// print("InitialMesh", mesh);
//
// labelList procIDs(3);
// procIDs[0] = 0;
// procIDs[1] = 1;
// procIDs[2] = 2;
//
////XXXXXX
// // Collect meshes from procs 0,1 (in newComm) on 1.
// lduPrimitiveMesh::gather(mesh, procIDs, procMeshes);
//
// if (Pstream::myProcNo(newComm) == procIDs[0])
// {
// // Print a bit
// forAll(procMeshes, i)
// {
// const lduMesh& pMesh = procMeshes[i];
// print("procMesh" + Foam::name(i), pMesh);
//
// const lduAddressing& addr = pMesh.lduAddr();
// checkUpperTriangular
// (
// addr.size(),
// addr.lowerAddr(),
// addr.upperAddr()
// );
// }
//
//
// // Combine
// labelList cellOffsets;
// labelListList faceMap;
// labelListList boundaryMap;
// labelListListList boundaryFaceMap;
// //autoPtr<lduPrimitiveMesh> allMeshPtr = combineMeshes
// //(
// // newComm,
// // procIDs,
// // procMeshes,
// //
// // cellOffsets, // per mesh the starting cell
// // faceMap, // per mesh the starting face
// // boundaryMap, // per mesh,per interface the starting face
// // boundaryFaceMap
// //);
// //const lduPrimitiveMesh& allMesh = allMeshPtr();
// const lduPrimitiveMesh allMesh
// (
// newComm,
// procIDs,
// procMeshes,
//
// cellOffsets,
// faceMap,
// boundaryMap,
// boundaryFaceMap
// );
//
//
// print("ALLMESH", allMesh);
//
// forAll(procMeshes, procMeshI)
// {
// const lduMesh& pMesh = procMeshes[procMeshI];
// //const lduAddressing& pAddressing = pMesh.lduAddr();
//
// Pout<< "procMesh" << procMeshI << endl
// << " cells start at:" << cellOffsets[procMeshI] << endl
// << " faces to to:" << faceMap[procMeshI] << endl;
//
// lduInterfacePtrsList interfaces = pMesh.interfaces();
// forAll(interfaces, intI)
// {
// Pout<< " patch:" << intI
// << " becomes patch:" << boundaryMap[procMeshI][intI]
// << endl;
//
// Pout<< " patch:" << intI
// << " faces become faces:"
// << boundaryFaceMap[procMeshI][intI]
// << endl;
// }
// }
// }
//
//
// // Construct test data
// // ~~~~~~~~~~~~~~~~~~~
//
// GAMGInterfaceFieldPtrsList interfaces(interfaces.size());
// FieldField<Field, scalar> interfaceBouCoeffs(interfaces.size());
// FieldField<Field, scalar> interfaceIntCoeffs(interfaces.size());
//
// forAll(interfaces, intI)
// {
// if (interfaces.set(intI))
// {
// label size = interfaces[intI].size();
//
// interfaces.set
// (
// intI,
// GAMGInterfaceField::New
// (
// mesh.interfaces()[intI],
// interfaces[intI]
// )
// );
// interfaceBouCoeffs.set(intI, new scalarField(size, 111));
// interfaceIntCoeffs.set(intI, new scalarField(size, 222));
// }
// }
//
//
// PtrList<lduMatrix> otherMats;
// PtrList<FieldField<Field, scalar> > otherBouCoeffs;
// PtrList<FieldField<Field, scalar> > otherIntCoeffs;
// PtrList<GAMGInterfaceFieldPtrsList> otherInterfaces;
// gatherMatrices
// (
// procIDs,
// procMeshes,
//
// mat,
// interfaceBouCoeffs,
// interfaceIntCoeffs,
// interfaces,
//
// otherMats,
// otherBouCoeffs,
// otherIntCoeffs,
// otherInterfaces
// );
////XXXXXX
//}
{
Pout<< "World:" << UPstream::worldComm
<< " procID:" << 2
<< " subComm:" << newComm
<< " rank1:" << UPstream::procNo(newComm, UPstream::worldComm, 1)
<< " rank2:" << UPstream::procNo(newComm, UPstream::worldComm, 2)
<< endl;
}
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
while (simple.correctNonOrthogonal())
{
fvScalarMatrix Teqn
(
//fvm::ddt(T) - fvm::laplacian(DT, T)
fvm::laplacian(DT, T)
);
{
label oldWarn = UPstream::warnComm;
UPstream::warnComm = newComm;
label oldComm = mesh.comm();
setCommunicator(mesh, newComm);
Pout<< "** oldcomm:" << oldComm
<< " newComm:" << mesh.comm() << endl;
if (Pstream::myProcNo(mesh.comm()) != -1)
{
solve(Teqn);
}
setCommunicator(mesh, oldComm);
Pout<< "** reset mesh to:" << mesh.comm() << endl;
UPstream::warnComm = oldWarn;
}
}
#include "write.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Pout<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,46 @@
if (runTime.outputTime())
{
volVectorField gradT(fvc::grad(T));
volScalarField gradTx
(
IOobject
(
"gradTx",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::X)
);
volScalarField gradTy
(
IOobject
(
"gradTy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Y)
);
volScalarField gradTz
(
IOobject
(
"gradTz",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Z)
);
runTime.write();
}

View File

@ -0,0 +1,3 @@
Test-parallel-communicators.C
EXE = $(FOAM_USER_APPBIN)/Test-parallel-communicators

View File

@ -0,0 +1,203 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Application
Test-parallel-communicators
Description
Checks communication using user-defined communicators
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "IPstream.H"
#include "OPstream.H"
#include "vector.H"
#include "IOstreams.H"
#include "PstreamReduceOps.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar sumReduce
(
const label comm,
const scalar localValue
)
{
scalar sum;
if (Pstream::parRun())
{
if (UPstream::master(comm))
{
// Add master value and all slaves
sum = localValue;
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave(comm);
slave++
)
{
scalar slaveValue;
UIPstream::read
(
Pstream::blocking,
slave,
reinterpret_cast<char*>(&slaveValue),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
sum += slaveValue;
}
// Send back to slaves
for
(
int slave=UPstream::firstSlave();
slave<=UPstream::lastSlave(comm);
slave++
)
{
UOPstream::write
(
UPstream::blocking,
slave,
reinterpret_cast<const char*>(&sum),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
}
else
{
{
UOPstream::write
(
UPstream::blocking,
UPstream::masterNo(),
reinterpret_cast<const char*>(&localValue),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
{
UIPstream::read
(
UPstream::blocking,
UPstream::masterNo(),
reinterpret_cast<char*>(&sum),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
}
}
return sum;
}
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Allocate a communicator
label n = Pstream::nProcs(UPstream::worldComm);
DynamicList<label> bottom;
DynamicList<label> top;
for (label i = 0; i < n/2; i++)
{
bottom.append(i);
}
for (label i = n/2; i < n; i++)
{
top.append(i);
}
//Pout<< "bottom:" << bottom << endl;
Pout<< "top :" << top << endl;
scalar localValue = 111*UPstream::myProcNo(UPstream::worldComm);
Pout<< "localValue :" << localValue << endl;
label comm = Pstream::allocateCommunicator
(
UPstream::worldComm,
top
);
Pout<< "allocated comm :" << comm << endl;
Pout<< "comm myproc :" << Pstream::myProcNo(comm)
<< endl;
if (Pstream::myProcNo(comm) != -1)
{
//scalar sum = sumReduce(comm, localValue);
//scalar sum = localValue;
//reduce
//(
// UPstream::treeCommunication(comm),
// sum,
// sumOp<scalar>(),
// Pstream::msgType(),
// comm
//);
scalar sum = returnReduce
(
localValue,
sumOp<scalar>(),
Pstream::msgType(),
comm
);
Pout<< "sum :" << sum << endl;
}
Pstream::freeCommunicator(comm);
Pout<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -493,6 +493,7 @@ int main(int argc, char *argv[])
mesh.nFaces(), // start
patchI, // index
mesh.boundaryMesh(),// polyBoundaryMesh
Pstream::worldComm, // communicator
Pstream::myProcNo(),// myProcNo
nbrProcI // neighbProcNo
)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -448,6 +448,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI]
);
@ -475,6 +476,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI],
referPatch,

View File

@ -98,6 +98,7 @@ int main(int argc, char *argv[])
cci.write();
}
volScalarField V
(
IOobject
@ -110,10 +111,19 @@ int main(int argc, char *argv[])
false
),
mesh,
dimensionedScalar("V", mesh.V().dimensions(), 0.0),
calculatedFvPatchField<scalar>::typeName
);
mesh.V().setInstance(runTime.timeName());
mesh.V().write();
V.dimensionedInternalField() = mesh.V();
forAll(V.boundaryField(), patchI)
{
V.boundaryField()[patchI] =
V.boundaryField()[patchI].patch().magSf();
}
Info<< "Writing cellVolumes and patch faceAreas to " << V.name()
<< " in " << runTime.timeName() << endl;
V.write();
}
Info<< "\nEnd\n" << endl;

View File

@ -0,0 +1,3 @@
surfaceHookUp.C
EXE = $(FOAM_APPBIN)/surfaceHookUp

View File

@ -0,0 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools \
-lfileFormats \
-ltriSurface

View File

@ -0,0 +1,525 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Application
surfaceHookUp
Description
Find close open edges and stitches the surface along them
Usage
- surfaceHookUp hookDistance [OPTION]
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "triSurfaceMesh.H"
#include "indexedOctree.H"
#include "treeBoundBox.H"
#include "PackedBoolList.H"
#include "unitConversion.H"
#include "searchableSurfaces.H"
using namespace Foam;
// Split faceI along edgeI at position newPointI
void greenRefine
(
const triSurface& surf,
const label faceI,
const label edgeI,
const label newPointI,
DynamicList<labelledTri>& newFaces
)
{
const labelledTri& f = surf.localFaces()[faceI];
const edge& e = surf.edges()[edgeI];
// Find index of edge in face.
label fp0 = findIndex(f, e[0]);
label fp1 = f.fcIndex(fp0);
label fp2 = f.fcIndex(fp1);
if (f[fp1] == e[1])
{
// Edge oriented like face
newFaces.append
(
labelledTri
(
f[fp0],
newPointI,
f[fp2],
f.region()
)
);
newFaces.append
(
labelledTri
(
newPointI,
f[fp1],
f[fp2],
f.region()
)
);
}
else
{
newFaces.append
(
labelledTri
(
f[fp2],
newPointI,
f[fp1],
f.region()
)
);
newFaces.append
(
labelledTri
(
newPointI,
f[fp0],
f[fp1],
f.region()
)
);
}
}
//scalar checkEdgeAngle
//(
// const triSurface& surf,
// const label edgeIndex,
// const label pointIndex,
// const scalar& angle
//)
//{
// const edge& e = surf.edges()[edgeIndex];
// vector eVec = e.vec(surf.localPoints());
// eVec /= mag(eVec) + SMALL;
// const labelList& pEdges = surf.pointEdges()[pointIndex];
//
// forAll(pEdges, eI)
// {
// const edge& nearE = surf.edges()[pEdges[eI]];
// vector nearEVec = nearE.vec(surf.localPoints());
// nearEVec /= mag(nearEVec) + SMALL;
// const scalar dot = eVec & nearEVec;
// const scalar minCos = degToRad(angle);
// if (mag(dot) > minCos)
// {
// return false;
// }
// }
// return true;
//}
void createBoundaryEdgeTrees
(
const PtrList<triSurfaceMesh>& surfs,
PtrList<indexedOctree<treeDataEdge> >& bEdgeTrees,
labelListList& treeBoundaryEdges
)
{
forAll(surfs, surfI)
{
const triSurface& surf = surfs[surfI];
// Boundary edges
treeBoundaryEdges[surfI] =
labelList
(
identity(surf.nEdges() - surf.nInternalEdges())
+ surf.nInternalEdges()
);
Random rndGen(17301893);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are less face/edge aligned items.
treeBoundBox bb
(
treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bEdgeTrees.set
(
surfI,
new indexedOctree<treeDataEdge>
(
treeDataEdge
(
false, // cachebb
surf.edges(), // edges
surf.localPoints(), // points
treeBoundaryEdges[surfI] // selected edges
),
bb, // bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"hook surfaces to other surfaces by moving and retriangulating their"
"boundary edges to match other surface boundary edges"
);
argList::noParallel();
argList::validArgs.append("hookTolerance");
# include "addDictOption.H"
# include "setRootCase.H"
# include "createTime.H"
const word dictName("surfaceHookUpDict");
# include "setSystemRunTimeDictionaryIO.H"
Info<< "Reading " << dictName << nl << endl;
const IOdictionary dict(dictIO);
const scalar dist(args.argRead<scalar>(1));
const scalar matchTolerance(SMALL);
Info<< "Hooking distance = " << dist << endl;
searchableSurfaces surfs
(
IOobject
(
"surfacesToHook",
runTime.constant(),
"triSurface",
runTime
),
dict
);
Info<< nl << "Reading surfaces: " << endl;
forAll(surfs, surfI)
{
Info<< incrIndent;
Info<< nl << indent << "Surface = " << surfs.names()[surfI] << endl;
const wordList& regions = surfs[surfI].regions();
forAll(regions, surfRegionI)
{
Info<< incrIndent;
Info<< indent << "Regions = " << regions[surfRegionI] << endl;
Info<< decrIndent;
}
Info<< decrIndent;
}
PtrList<indexedOctree<treeDataEdge> > bEdgeTrees(surfs.size());
labelListList treeBoundaryEdges(surfs.size());
List<DynamicList<labelledTri> > newFaces(surfs.size());
List<DynamicList<point> > newPoints(surfs.size());
List<PackedBoolList> visitedFace(surfs.size());
PtrList<triSurfaceMesh> newSurfaces(surfs.size());
forAll(surfs, surfI)
{
const triSurfaceMesh& surf =
refCast<const triSurfaceMesh>(surfs[surfI]);
newSurfaces.set
(
surfI,
new triSurfaceMesh
(
IOobject
(
"hookedSurface_" + surfs.names()[surfI],
runTime.constant(),
"triSurface",
runTime
),
surf
)
);
}
label nChanged = 0;
label nIters = 0;
do
{
Info<< nl << "Iteration = " << nIters++ << endl;
nChanged = 0;
createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);
forAll(newSurfaces, surfI)
{
const triSurface& newSurf = newSurfaces[surfI];
newFaces[surfI] = newSurf.localFaces();
newPoints[surfI] = newSurf.localPoints();
visitedFace[surfI] = PackedBoolList(newSurf.size(), false);
}
forAll(newSurfaces, surfI)
{
const triSurface& surf = newSurfaces[surfI];
List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
labelList bPointsHitTree(surf.boundaryPoints().size(), -1);
const labelListList& pointEdges = surf.pointEdges();
forAll(bPointsTobEdges, bPointI)
{
pointIndexHit& nearestHit = bPointsTobEdges[bPointI];
const label pointI = surf.boundaryPoints()[bPointI];
const point& samplePt = surf.localPoints()[pointI];
const labelList& pEdges = pointEdges[pointI];
// Add edges connected to the edge to the shapeMask
DynamicList<label> shapeMask;
shapeMask.append(pEdges);
forAll(bEdgeTrees, treeI)
{
const indexedOctree<treeDataEdge>& bEdgeTree =
bEdgeTrees[treeI];
pointIndexHit currentHit =
bEdgeTree.findNearest
(
samplePt,
sqr(dist),
treeDataEdge::findNearestOpSubset
(
bEdgeTree,
shapeMask
)
);
if
(
currentHit.hit()
&&
(
!nearestHit.hit()
||
(
magSqr(currentHit.hitPoint() - samplePt)
< magSqr(nearestHit.hitPoint() - samplePt)
)
)
)
{
nearestHit = currentHit;
bPointsHitTree[bPointI] = treeI;
}
}
scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);
if (nearestHit.hit())
{
// bool rejectEdge =
// checkEdgeAngle
// (
// surf,
// nearestHit.index(),
// pointI,
// 30
// );
if (dist2 > Foam::sqr(dist))
{
nearestHit.setMiss();
}
}
}
forAll(bPointsTobEdges, bPointI)
{
const pointIndexHit& eHit = bPointsTobEdges[bPointI];
if (eHit.hit())
{
const label hitSurfI = bPointsHitTree[bPointI];
const triSurface& hitSurf = newSurfaces[hitSurfI];
const label eIndex =
treeBoundaryEdges[hitSurfI][eHit.index()];
const edge& e = hitSurf.edges()[eIndex];
const label pointI = surf.boundaryPoints()[bPointI];
const labelList& eFaces = hitSurf.edgeFaces()[eIndex];
if (eFaces.size() != 1)
{
WarningIn(args.executable())
<< "Edge is attached to " << eFaces.size()
<< " faces." << endl;
continue;
}
const label faceI = eFaces[0];
if (visitedFace[hitSurfI][faceI])
{
continue;
}
DynamicList<labelledTri> newFacesFromSplit(2);
const point& pt = surf.localPoints()[pointI];
if
(
(
magSqr(pt - hitSurf.localPoints()[e.start()])
< matchTolerance
)
|| (
magSqr(pt - hitSurf.localPoints()[e.end()])
< matchTolerance
)
)
{
continue;
}
nChanged++;
// Keep the points in the same place and move the edge
newPoints[hitSurfI].append(newPoints[surfI][pointI]);
// Move the points to the edges
//newPoints[pointI] = eHit.hitPoint();
//newPoints.append(eHit.hitPoint());
visitedFace[hitSurfI][faceI] = true;
// Split the other face.
greenRefine
(
hitSurf,
faceI,
eIndex,
newPoints[hitSurfI].size() - 1,
newFacesFromSplit
);
forAll(newFacesFromSplit, newFaceI)
{
if (newFaceI == 0)
{
newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
}
else
{
newFaces[hitSurfI].append
(
newFacesFromSplit[newFaceI]
);
}
}
}
}
}
Info<< " Number of edges split = " << nChanged << endl;
forAll(newSurfaces, surfI)
{
newSurfaces.set
(
surfI,
new triSurfaceMesh
(
IOobject
(
"hookedSurface_" + surfs.names()[surfI],
runTime.constant(),
"triSurface",
runTime
),
triSurface
(
newFaces[surfI],
newSurfaces[surfI].patches(),
pointField(newPoints[surfI])
)
)
);
}
} while (nChanged > 0);
Info<< endl;
forAll(newSurfaces, surfI)
{
const triSurfaceMesh& newSurf = newSurfaces[surfI];
Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
<< endl;
newSurf.searchableSurface::write();
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,22 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object surfaceHookUpDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
surface1.stl {type triSurfaceMesh;}
surface2.stl {type triSurfaceMesh;}
// ************************************************************************* //

View File

@ -31,6 +31,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "triSurface.H"
#include "triSurfaceSearch.H"
#include "argList.H"
#include "OFstream.H"
#include "IFstream.H"
@ -40,6 +41,7 @@ Description
#include "indexedOctree.H"
#include "treeDataTriSurface.H"
#include "Random.H"
#include "volumeType.H"
using namespace Foam;
@ -242,26 +244,16 @@ int main(int argc, char *argv[])
// Read surface to select on
triSurface selectSurf(surfName);
// bb of surface
treeBoundBox bb(selectSurf.localPoints());
// Random number generator
Random rndGen(354543);
// search engine
indexedOctree<treeDataTriSurface> selectTree
triSurfaceSearch searchSelectSurf
(
treeDataTriSurface
(
selectSurf,
indexedOctree<treeDataTriSurface>::perturbTol()
),
bb.extend(rndGen, 1e-4), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
selectSurf,
indexedOctree<treeDataTriSurface>::perturbTol(),
8
);
const indexedOctree<treeDataTriSurface>& selectTree =
searchSelectSurf.tree();
// Check if face (centre) is in outside or inside.
forAll(facesToSubset, faceI)
{
@ -269,14 +261,13 @@ int main(int argc, char *argv[])
{
const point fc(surf1[faceI].centre(surf1.points()));
indexedOctree<treeDataTriSurface>::volumeType t =
selectTree.getVolumeType(fc);
volumeType t = selectTree.getVolumeType(fc);
if
(
outside
? (t == indexedOctree<treeDataTriSurface>::OUTSIDE)
: (t == indexedOctree<treeDataTriSurface>::INSIDE)
? (t == volumeType::OUTSIDE)
: (t == volumeType::INSIDE)
)
{
facesToSubset[faceI] = true;